Dear dr. Steven and meep users,
I made a 2d-tm planewave near to far transformation code (Fortran)
But I found there is something wrong with the code.
First I use meep to calculate the near field code.(a poin source)
Attached is the code.
Then I combined the data to three files (ez.txt, hx.txt, hy.txt). The order
that I stored the data is upper surface-->right surface-->lower surface-->left
surface, because the direction of the integral of the surface current is
clockwise.
After the transformation, I think the magnitude of far field electric field
(y-axis) vs. angle (x-axis) is supposed to be a straight line. However, I found
it is not.
Attached is the ntf Fortran code.
Could you help me find the error of my code? I tested the code many times but I
still can not figure out what the error is.
Could you spend sometime on my code please?
I think the probable error is that the parameters are matched with meep.
I really appreciate your help.
Lei Zhu
program main
integer ::wl=40 !--resolution
in meep (40 grids per wavelength)
double precision,parameter :: eps0=1
double precision,parameter :: xmu0=1
double precision ::pi=3.1415926
complex(kind=8),dimension(:),allocatable ::Ezf,Hxf,Hyf !--far fields
real startangle, endangle,degreestep
complex(kind=8) ::fx,fy,fz,fmx,fmy,fmz !--electric
and magnetic current moment on the virtual surface
double precision ::fai,dl,cfai,sfai,kk,rr,Ez0
complex(kind=8),dimension(:),allocatable ::Ez,Hx,Hy !--electric
and magnetic field on the virtual surface
complex(kind=8) ::jj=(0,1)
integer ::Imin=-100,Imax=100,Jmin=-100,Jmax=100,Count !--virtual
surface
number=(Imax+Jmax-Imin-Jmin)*2
startangle=0
endangle=360
degreestep=1
open(20,file='ez.txt') !--file of
near field data which got from meep
open(21,file='hx.txt') !--file of
near field data which got from meep
open(22,file='hy.txt') !--file of
near field data which got from meep
allocate(Ez(1:number),Hx(1:number),Hy(1:number))
allocate(Ezf(0:361),Hxf(0:361),Hyf(0:361))
read(20,*) Ez
read(21,*) Hx
read(22,*) Hy
rr=100*wl !--the
distance of the far field
print *,Imin,Imax,Jmin,Jmax
close(20)
kk=2.0*pi/real(wl) !--the phase
change per space step (per grid)
dl=1.0
print *," angle/pi Ez"
open(unit=1,file='degree.dat',status='replace')
open(unit=2,file='farfield.dat',status='replace')
open(unit=3,file='ang.dat',status='replace')
Ez0=0.0
do N=startangle,endangle,degreestep
fai=pi*real(N)/180.0
fx=0;fy=0;fz=0;fmx=0;fmy=0;fmz=0
cfai=cos(fai)
sfai=sin(fai)
K=1
!upper surface
fz=fz+(-Hy(K)*exp(jj*kk*(real(Imin)*cfai+real(Jmax)*sfai))*dl)*0.5
fmy=fmy+(-Ez(K)*exp(jj*kk*(real(Imin)*cfai+real(Jmax)*sfai))*dl)*0.5
fz=fz+(-Hx(K)*exp(jj*kk*(real(Imin)*cfai+real(Jmax)*sfai))*dl)*0.5
fmx=fmx+(-Ez(K)*exp(jj*kk*(real(Imin)*cfai+real(Jmax)*sfai))*dl)*0.5
K=K+1
do I=Imin+1,Imax-1
fz=fz+(-Hx(K)*exp(jj*kk*(real(I)*cfai+real(Jmax)*sfai))*dl)
fmx=fmx+(-Ez(K)*exp(jj*kk*(real(I)*cfai+real(Jmax)*sfai))*dl)
K=K+1
end do
! right surface
fz=fz+(-Hx(K)*exp(jj*kk*(real(Imax)*cfai+real(Jmax)*sfai))*dl)*0.5
fmx=fmx+(-Ez(K)*exp(jj*kk*(real(Imax)*cfai+real(Jmax)*sfai))*dl)*0.5
fz=fz+(Hy(K)*exp(jj*kk*(real(Imax)*cfai+real(Jmax)*sfai))*dl)*0.5
fmy=fmy+(Ez(K)*exp(jj*kk*(real(Imax)*cfai+real(Jmax)*sfai))*dl)*0.5
K=K+1
do J=Jmax-1,Jmin+1,-1
fz=fz+(Hy(K)*exp(jj*kk*(real(Imax)*cfai+real(J)*sfai))*dl)
fmy=fmy+(Ez(K)*exp(jj*kk*(real(Imax)*cfai+real(J)*sfai))*dl)
K=K+1
end do
!lower surface
fz=fz+(Hy(K)*exp(jj*kk*(real(Imax)*cfai+real(Jmin)*sfai))*dl)*0.5
fmy=fmy+(Ez(K)*exp(jj*kk*(real(Imax)*cfai+real(Jmin)*sfai))*dl)*0.5
fz=fz+(Hx(K)*exp(jj*kk*(real(Jmax)*cfai+real(Jmin)*sfai))*dl)*0.5
fmx=fmx+(Ez(K)*exp(jj*kk*(real(Imax)*cfai+real(Jmin)*sfai))*dl)*0.5
K=K+1
do I=Imax-1,Imin+1,-1
fz=fz+(Hx(K)*exp(jj*kk*(real(I)*cfai+real(Jmin)*sfai))*dl)
fmx=fmx+(Ez(K)*exp(jj*kk*(real(I)*cfai+real(Jmin)*sfai))*dl)
K=K+1
end do
!left surface
fz=fz+(Hx(K)*exp(jj*kk*(real(Imin)*cfai+real(Jmin)*sfai))*dl)*0.5
fmx=fmx+(Ez(K)*exp(jj*kk*(real(Imin)*cfai+real(Jmin)*sfai))*dl)*0.5
fz=fz+(-Hy(K)*exp(jj*kk*(real(Imin)*cfai+real(Jmin)*sfai))*dl)*0.5
fmy=fmy+(-Ez(K)*exp(jj*kk*(real(Imin)*cfai+real(Jmin)*sfai))*dl)*0.5
K=K+1
do J=Jmin+1,Jmax-1
fz=fz+(-Hy(K)*exp(jj*kk*(real(Imin)*cfai+real(J)*sfai))*dl)
fmy=fmy+(-Ez(K)*exp(jj*kk*(real(Imin)*cfai+real(J)*sfai))*dl)
K=K+1
end do
!-------------calculate the far field Ez---------------
Ezf(N)=0.5*jj*kk*exp(-jj*kk*rr)*(-sqrt(xmu0/eps0)*fz-fmx*sfai+fmy*cfai)/sqrt(2.0*jj*pi*kk*rr)
Ez0=abs(Ezf(N))+Ez0
write(*,"(f16.4,e20.8)") real(N),abs(Ezf(N))
write(1,*) fai
write(2,*) abs(Ezf(N))
write(3,*) atan(imag(Ezf(N))/real(Ezf(N)))
end do
close(21)
close(22)
close(23)
close(1)
close(2)
close(3)
end
point.ctl
Description: point.ctl
_______________________________________________ meep-discuss mailing list [email protected] http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
