1st upload

ImageCreator.jl

@@ -0,0 +1,35 @@
+const c=299792458.
+
+function IC(l,p,h,X,Y,Z,tilt,azimut,phase,amplitude,order)
+  const Rx=vcat([-1 0 0],[0 1 0],[0 0 1])
+  const Ry=vcat([1 0 0],[0 -1 0],[0 0 1])
+  const Rz=vcat([1 0 0],[0 1 0],[0 0 -1])
+  const i0=amplitude*[sin(tilt).*cos(azimut), sin(tilt).*sin(azimut), cos(tilt)]
+  POS=zeros((order+1)^3*8,7)
+
+  count=1
+  for i=0:order
+    perc=round(i/order*1000)/10
+    println("$perc %")
+    for j=0:order
+      for k=0:order
+        n = int(i+j+k);
+          ip = (-1)^n*Rx^abs(i)*Ry^abs(j)*Rz^abs(k)*i0;
+          x = 2*i*l+(-1)^abs(i)*X;
+          y =	2*j*p+(-1)^abs(j)*Y;
+          z = 2*k*h+(-1)^abs(k)*Z;
+          POS[count,:]=[ip' x  y  z n ]
+          POS[count+1,:]=[(-Rz*ip)' x  y  -z n+1 ]
+          POS[count+2,:]=[(-Ry*ip)' x  -y  z n+1 ]
+          POS[count+3,:]=[(-Ry*-Rz*ip)' x  -y  -z n+2 ]
+          POS[count+4,:]=[(-Rx*ip)' -x  y  z n+1 ]
+          POS[count+5,:]=[(-Rx*-Rz*ip)' -x  y  -z n+2 ]
+          POS[count+6,:]=[(-Rx*-Ry*ip)' -x  -y  z n+2 ]
+          POS[count+7,:]=[(-Rx*-Ry*-Rz*ip)' -x  -y  -z n+3 ]
+	        count+=8
+      end
+    end
+  end
+  perm=sortperm(POS[:,7])
+  return POS[perm,:]
+end

License.md

@@ -0,0 +1,7 @@
+This model represents the work done during my Ph.D. The set of programs is provided under a license or a nondisclosure agreement, and it may be used only in accordance with terms of those agreements. By using this software, you agree to credit the authors: Emmanuel Amador, Philippe Besnier and Christophe Lemoine from the Institut d'Électronique et de Télécommunications de Rennes IETR and the main article describing this model:
+
+E. Amador, C. Lemoine, P. Besnier, and A. Laisné, “Reverberation chamber modeling based on image theory: Investigation in the pulse regime,” Electromagnetic Compatibility, IEEE Transactions on, vol. 52, no. 4, pp. 778 – 789, 2010. IEEExplore, Bibtex
+
+You can use the source-code and make modifications to suit your needs. You cannot distribute this software or an altered version of this software without authors' permissions.
+
+Comments and suggestions are greatly welcomed.

README.md

@@ -0,0 +1,28 @@
+Reverb_Julia
+============
+
+Overview
+--------
+
+This a Julia code that computes the electric field in an electromagnetic reverberation chamber in the time domain or in the frequency domain.
+This model is based on image theory.
+This model represents the work done during my Ph.D. The set of programs is provided under a license or a nondisclosure agreement, and it may be used only in accordance with terms of those agreements. By using this software, you agree to credit the authors: Emmanuel Amador, Philippe Besnier and Christophe Lemoine from the Institut d'Électronique et de Télécommunications de Rennes IETR and the main article describing this model:
+
+E. Amador, C. Lemoine, P. Besnier, and A. Laisné, “Reverberation chamber modeling based on image theory: Investigation in the pulse regime,” *Electromagnetic Compatibility, IEEE Transactions on, vol. 52, no. 4*, pp. 778 – 789, 2010.
+
+http://dx.doi.org/10.1109/TEMC.2010.2049576
+
+You can use the source-code and make modifications to suit your needs. You cannot distribute this software or an altered version of this software without authors' permissions.
+
+Comments and suggestions are greatly welcomed.
+
+visit http://projets.ietr.fr/imagetheorymodel/ for more details.
+
+Requirements
+------------
+
+* Julia 0.3
+* PyPlot
+* HFS5
+* A lot of memory...
+

TDRC.jl

@@ -0,0 +1,157 @@
+include("ImageCreator.jl")
+
+
+const c = 299792458.
+const mu0 = 4*pi*1e-7
+const eps0 = 1/(mu0*c^2)
+const epsr=1.
+
+#cavity
+const L=8.7
+const l=3.7
+const h=2.9
+const losses= 0.998
+const Lt=.3e-6
+const dmax=Lt*c
+const order=int(round(dmax/minimum([L;l;h]))+1)
+const freq=.5e9
+const w=2*pi*freq  # omega
+const k=w/c     # wave number
+
+using HDF5,JLD
+
+##dipole
+#dipole position
+const x0=1
+const y0=2
+const z0=1
+const tilt=pi/2-acos(sqrt(2/3));
+const azimut=pi/4
+const phi=pi/2
+#dipole moment
+#total time averaged radiated power P= 1 W dipole moment => |p|=sqrt(12πcP/µOω⁴)
+const Pow=1
+const amplitude=sqrt(12*pi*c*Pow/(mu0*(2*pi*freq)^4))
+
+
+#Images
+POS=IC(L,l,h,x0,y0,z0,tilt,azimut,phi,amplitude,order)
+const dist=sqrt(POS[:,4].^2+POS[:,5].^2+POS[:,6].^2)
+const perm=sortperm(dist)
+const U=find(dist[perm].>dmax)
+POS=POS[perm[1:U[1]-1],:]
+#@save "POS.jld" POS
+#@load "POS.jld" POS
+
+const numberofimages=length(POS[:,1])
+
+tic()
+
+
+#observation points
+const dx=0.05
+const dy=0.05
+const dz=0.05
+
+const x=linspace(0,L,int(L/dx))
+const y=linspace(0,l,int(l/dy))
+const z=2
+
+const nx=length(x)
+const ny=length(y)
+
+const N=int(Lt*freq)
+const t=linspace(Lt/N,Lt,N)
+
+println("Computing the radiation...")
+Z=zeros(nx,ny,N)
+Et=zeros(nx,ny,N)
+Bt=zeros(nx,ny,N)
+for i=1:nx
+  perc=round(i/nx*1000)/10
+  println("$perc %")
+  for j=1:ny
+    r=[x[i],y[j],z] #position of the receiver
+    E=zeros(Complex128,3,N)
+    B=zeros(Complex128,3,N)
+    for m=1:numberofimages
+      ord=POS[m,7]#order of the dipole
+      p=vec(POS[m,1:3])*losses^ord #image dipole moment
+      R=vec(POS[m,4:6]) #image dipole position
+      rprime=r-R  # r'=r-R
+      magrprime=sqrt(sum(rprime.^2)) # |r-R|
+      krp=k*magrprime  # k*|r-R|
+      rprime_cross_p = cross(rprime, p) # (r-R) x p
+      rp_c_p_c_rp = cross(rprime_cross_p, rprime) # ((r-R) x p) x (r-R)
+      ta=int(magrprime/c/Lt*N)
+      if ta == 0
+        ta=1
+      end
+      expfac=exp(1im*(-w*t+krp-phi))
+      #expfac[1:ta]=0
+      Ep=1/(4*pi*eps0*epsr)*(w^2/(c^2*magrprime^3)*rp_c_p_c_rp+(1/magrprime^3-w*im/(c*magrprime^2))*(3*rprime*dot(rprime,p)/magrprime^2-p))
+      Bp=1/(4*pi*eps0*epsr)*1/(magrprime*c^3)*(w^2*rprime_cross_p)/magrprime*(1-c/(im*w*magrprime))
+      u=1
+      for n=ta:N
+        E[:,n]+=expfac[u]*Ep
+        B[:,n]+=expfac[u]*Bp
+        u+=1
+      end
+    end
+    Z[i,j,:]=sqrt(sum(real(E).^2,1))./sqrt(sum(real(B).^2,1))*mu0 #real(E).^2#0.5*numpy.cross(E.T,conjugate(B.T)) #impedance
+    Et[i,j,:]=real(E[3,:]) #vertical E-field
+    Bt[i,j,:]=real(B[2,:]) #y B-field
+  end
+end
+
+save("Cartos.jld", "Et", Et, "Bt", Bt, "Z", Z, "t", t, "N", N)
+
+#Some figures
+using PyPlot
+pygui(false)
+for u=1:N
+  ts=round(t[u]/1e-9*1000)/1000
+
+  figure(num=1,figsize=(10,4)) #total E-field
+  title("\$E_z\$ (V/m), t=$ts ns")
+  pcolor(x,y,Et[:,:,u]',cmap="jet")
+  clim(-20,20)
+  axis("scaled")
+  xlim(0,L)
+  ylim(0,l)
+  grid()
+  colorbar(shrink=1,orientation="horizontal")
+  xlabel("\$x\$/m")
+  ylabel("\$y\$/m")
+  savefig("Ez_$u.png",bbox="tight")
+  clf()
+
+  figure(num=2,figsize=(10,4))  #total B-field
+  title("\$B_y\$ (A/m), t=$ts ns")
+  pcolor(x,y,Bt[:,:,u]',cmap="jet")
+  clim(-7e-8,7e-8)
+  axis("scaled")
+  xlim(0,L)
+  ylim(0,l)
+  grid()
+  colorbar(shrink=1,orientation="horizontal")
+  xlabel("\$x\$/m")
+  ylabel("\$y\$/m")
+  savefig("By_$u.png",bbox="tight")
+  clf()
+
+  figure(num=3,figsize=(10,4))  #Z
+  title("\$Z\$ (\$ \\Omega\$), t=$ts ns")
+  pcolor(x,y,Z[:,:,u]',cmap="jet")
+  clim(0,1000)
+  axis("scaled")
+  xlim(0,L)
+  ylim(0,l)
+  grid()
+  colorbar(shrink=1,orientation="horizontal")
+  xlabel("\$x\$/m")
+  ylabel("\$y\$/m")
+  savefig("Z_$u.png",bbox="tight")
+  clf()
+end
+

TDRC_stat.jl

@@ -0,0 +1,108 @@
+include("ImageCreator.jl")
+
+
+const c = 299792458.
+const mu0 = 4*pi*1e-7
+const eps0 = 1/(mu0*c^2)
+const epsr=1.
+
+#cavity
+const L=8.7
+const l=3.7
+const h=2.9
+const losses= 0.998
+
+#simulation parameters
+const Lt=1e-6
+const freq=.5e9
+
+const dmax=Lt*c
+const order=int(round(dmax/minimum([L;l;h]))+1)
+const w=2*pi*freq  # omega
+const k=w/c     # wave number
+const N=int(Lt*freq)
+const t=linspace(Lt/N,Lt,N)
+
+
+using HDF5,JLD
+
+##dipole
+#dipole position
+const x0=1
+const y0=2
+const z0=1
+const tilt=pi/2-acos(sqrt(2/3));
+const azimut=pi/4
+const phi=pi/2
+#dipole moment
+#total time averaged radiated power P= 1 W dipole moment => |p|=sqrt(12πcP/µOω⁴)
+const Pow=1
+const amplitude=sqrt(12*pi*c*Pow/(mu0*(2*pi*freq)^4))
+
+
+#Images
+println("Computing images...")
+POS=IC(L,l,h,x0,y0,z0,tilt,azimut,phi,amplitude,order)
+const dist=sqrt(POS[:,4].^2+POS[:,5].^2+POS[:,6].^2)
+const perm=sortperm(dist)
+const U=find(dist[perm].>dmax)
+POS=POS[perm[1:U[1]-1],:]
+#@save "POS.jld" POS
+#@load "POS.jld" POS
+
+const numberofimages=length(POS[:,1])
+
+tic()
+
+#observation points
+M=150
+x=rand(M)*(L-c/freq)+2c/freq
+y=rand(M)*(l-c/freq)+2c/freq
+z=rand(M)*(h-c/freq)+2c/freq
+
+
+
+println("Computing the fields...")
+E=zeros(Complex128,M,N,3)
+B=zeros(Complex128,M,N,3)
+Z=zeros(M,N)
+for i=1:M
+  r=[x[i],y[i],z[i]] #position of the receiver
+  Et=zeros(Complex128,3,N)
+  Bt=zeros(Complex128,3,N)
+  for m=1:numberofimages
+    ord=POS[m,7]#order of the dipole
+    p=vec(POS[m,1:3])*losses^ord #image dipole moment
+    R=vec(POS[m,4:6]) #image dipole position
+    rprime=r-R  # r'=r-R
+    magrprime=sqrt(sum(rprime.^2)) # |r-R|
+    krp=k*magrprime  # k*|r-R|
+    rprime_cross_p = cross(rprime, p) # (r-R) x p
+    rp_c_p_c_rp = cross(rprime_cross_p, rprime) # ((r-R) x p) x (r-R)
+    ta=int(magrprime/c/Lt*N)
+    if ta == 0
+      ta=1
+    end
+    expfac=exp(1im*(-w*t+krp-phi))
+    Ep=1/(4*pi*eps0*epsr)*(w^2/(c^2*magrprime^3)*rp_c_p_c_rp+(1/magrprime^3-w*im/(c*magrprime^2))*(3*rprime*dot(rprime,p)/magrprime^2-p))
+    Bp=1/(4*pi*eps0*epsr)*1/(magrprime*c^3)*(w^2*rprime_cross_p)/magrprime*(1-c/(im*w*magrprime))
+    u=1
+    for n=ta:N
+      Et[:,n]+=expfac[u]*Ep
+      Bt[:,n]+=expfac[u]*Bp
+      u+=1
+    end
+  end
+  E[i,:,:]=Et' # E-field
+  B[i,:,:]=Bt' # B-field
+  Z[i,:]=sqrt(sum(real(Et).^2,1))./sqrt(sum(real(Bt).^2,1))*mu0 #real(E).^2#0.5*numpy.cross(E.T,conjugate(B.T)) #impedance
+  perc=round(i/M*1000)/10
+  println("$perc %")
+end
+
+Er=real(E)
+Ei=imag(E)
+Br=real(B)
+Bi=imag(B)
+
+save("Stats150.jld", "Er", Er, "Ei",Ei, "Br", Br, "Bi", Bi,"Z",Z,"t",t,"N",N)