c program 计算三维有限厚度TE-like透过率
***************************
**** YOU need alter these parameter according your promlem
parameter(number=334)
**** 晶体大小
parameter (nx1=-8,nx2=263,ny1=-138,ny2=138)
****** 积分位置 nny1,nny2 the input field integration position
****** nny3,nny4 one branch output field integration position
****** and nny5, nny6 the another
parameter (nny1=-27,nny2=+27,ny3=51,ny4=104,ny5=39,ny6=112)
***** 入射场设置 we should change the ny0 position
********** width1 is the input light source half width
parameter (nnx1=nx1-3,ny0=0,width1=30,nk2=4)
***************************************
****************************************
parameter (nk1=1,nnx2=nx2+3,width2=width1*1.5)
parameter (iimin=nx1-10,iimax=nx2+6)
parameter (jjmin=ny1-6,jjmax=ny2+6)
parameter (imin=iimin-12,imax=iimax+12)
parameter (jmin=jjmin-12,jmax=jjmax+12)
parameter (kkmin=0,kmin=0,kkmax=nk2+5,kmax=kkmax+12)
real ez(imin:imax,jmin:jmax,kmin:kmax-1)
real ex(imin:imax-1,jmin:jmax,kmin+1:kmax)
real ey(imin:imax,jmin:jmax-1,kmin+1:kmax)
real exy(imin:imax-1,jmin:jmax,kmin+1:kmax)
real exz(imin:imax-1,jmin:jmax,kmin+1:kmax)
real ezx(imin:imax,jmin:jmax,kmin:kmax-1)
real ezy(imin:imax,jmin:jmax,kmin:kmax-1)
real eyx(imin:imax,jmin:jmax-1,kmin+1:kmax)
real eyz(imin:imax,jmin:jmax-1,kmin+1:kmax)
real hx(imin:imax,jmin:jmax-1,kmin:kmax-1)
real hy(imin:imax-1,jmin:jmax,kmin:kmax-1)
real hz(imin:imax-1,jmin:jmax-1,kmin+1:kmax)
real hxy(imin:imax,jmin:jmax-1,kmin:kmax-1)
real hxz(imin:imax,jmin:jmax-1,kmin:kmax-1)
real hyx(imin:imax-1,jmin:jmax,kmin:kmax-1)
real hyz(imin:imax-1,jmin:jmax,kmin:kmax-1)
real hzx(imin:imax-1,jmin:jmax-1,kmin+1:kmax)
real hzy(imin:imax-1,jmin:jmax-1,kmin+1:kmax)
real epsx1(imin:imax,jmin:jmax,kmin:kmax)
real epsx2(imin:imax,jmin:jmax,kmin:kmax)
real epsx3(imin:imax,jmin:jmax,kmin:kmax)
real epsx4(imin:imax,jmin:jmax,kmin:kmax)
real epsx5(imin:imax,jmin:jmax,kmin:kmax)
real epsx6(imin:imax,jmin:jmax,kmin:kmax)
real epsx7(imin:imax,jmin:jmax,kmin:kmax)
real epsx8(imin:imax,jmin:jmax,kmin:kmax)
real epsxex(imin:imax,jmin:jmax,kmin:kmax)
real epsxey(imin:imax,jmin:jmax,kmin:kmax)
real epsxez(imin:imax,jmin:jmax,kmin:kmax)
real fax(imin:imax),famx(imin:imax)
real fay(jmin:jmax),famy(jmin:jmax)
real faz(kmin:kmax),famz(kmin:kmax)
complex hi(jmin:jmax,nk1:nk2,0:2000)
complex hi1(jmin:jmax,nk1:nk2,0:2000)
complex hi2(jmin:jmax,nk1:nk2,0:2000)
complex ho(jmin:jmax,nk1:nk2,0:2000)
complex heo(jmin:jmax,nk1:nk2,0:2000)
complex xubu
real thi(0:2000),tho(0:2000),theo(0:2000)
real thi1(0:2000),thi2(0:2000)
real input(number,3)
real ezinc(imin:imax,jmin:jmax,kmin:kmax)
real eyinc(imin:imax,jmin:jmax,kmin:kmax)
real fs,fs1,fs2,fs3,fs4,fsmin,fsmax
real eps0,xmu0,c0,abc,xx,yy,zz,e1,e0
real tt,ttt,pinlv,freq,ffmax,ffmin,f0,ff
real t0,omiga1,omiga2,omiga3,maxfax,maxfay,maxfaz
integer it,i,j,k,itmax,n0,nx1,nx2,ny1,ny2
integer nfreq,mfreq,hh
open(0,file='/userdata1/chengby/fdtd/han/test02/
$in.dat',status='old')
open(1,file='/userdata1/chengby/fdtd/han/test02/
$TE-d=42-r=12-h=22-bend-1.dat',status='unknown')
open(2,file='/userdata1/chengby/fdtd/han/test02/
$TE-d=42-r=12-h=22-bend-2.dat',status='unknown')
open(3,file='/userdata1/chengby/fdtd/han/test02/
$TE-d=42-r=12-h=22-bend-3.dat',status='unknown')
open(4,file='/userdata1/chengby/fdtd/han/test01/
$TE-d=42-r=12-h=22-bend-4.dat',status='unknown')
open(5,file='/userdata1/chengby/fdtd/han/test02/
$TE-d=42-r=12-h=22-bend-5.dat',status='unknown')
open(6,file='/userdata1/chengby/fdtd/han/test02/
$epsx-fenbu.dat',status='unknown')
open(7,file='/userdata1/chengby/fdtd/han/test02/
$TE-d=42-r=12-h=22-bend-7.dat',status='unknown')
open(8,file='/userdata1/chengby/fdtd/han/test02/
$TE-d=42-r=12-h=22-bend-8.dat',status='unknown')
open(9,file='/userdata1/chengby/fdtd/han/test02/
$TE-d=42-r=12-h=22-bend-9.dat',status='unknown')
open(10,file='/userdata1/chengby/fdtd/han/test02/
$TE-d=42-r=12-h=22-bend-10.dat',status='unknown')
open(11,file='/userdata1/chengby/fdtd/han/test02/
$TE-d=42-r=12-h=22-bend-11.dat',status='unknown')
open(12,file='/userdata1/chengby/fdtd/han/test02/
$TE-d=42-r=12-h=22-bend-12.dat',status='unknown')
open(13,file='/userdata1/chengby/fdtd/han/test02/
$TE-d=42-r=12-h=22-bend-13.dat',status='unknown')
open(14,file='/userdata1/chengby/fdtd/han/test02/
$TE-d=42-r=12-h=22-bend-14.dat',status='unknown')
open(15,file='/userdata1/chengby/fdtd/han/test02/
$TE-d=42-r=12-h=22-bend-15.dat',status='unknown')
open(16,file='/userdata1/chengby/fdtd/han/test02/
$TE-d=42-r=12-h=22-bend-16.dat',status='unknown')
eps0=8.854e-12
pi=3.1415926
xmu0=pi*4.0e-7
c0=1.0/sqrt(eps0*xmu0)
xubu=(0.0,1.0)
time=200.0
***********************
**************************
******************************
**** YOU need alter these parameter according your promlem
period=420.0
height=224.0
**************** 输入频率范围
ffmax=3.0e14
ffmin=1.0e14
*************** 请输入中心频率
f0=2.0e14
**************** 输出频率单位 THz unit
freq=1.0e12
*************** total time step
itmax=60000
*************** e0 为背景介电场数
**************** e1 为介质介电常数
e0=3.4**2
e1=1.0
**************************************
!!!!!!!!!!! 以下可以不作改动
**************************************
xx=period/150.0*1.0e-8
yy=period/150.0*1.0e-8
zz=period/150.0*1.0e-8
d=15.0
r=radius/period*d
h=height/period*d/2.
ratio=15.0*1.0/period !!!每个晶格常数中有15个网格
hh=int(h-0.5)
abc=1.0/sqrt(1.0/(xx*xx)+1.0/(yy*yy)+1.0/(zz*zz))
ttt=abc/c0
tt=0.99*ttt
t0=1.3/(f0*tt)
omiga1=2*pi/(t0*tt)
omiga2=2*pi*2/(t0*tt)
omiga3=2*pi*3/(t0*tt)
ff=1.0/(tt*itmax)
mfreq=int(ffmax/ff)
nfreq=int(ffmin/ff)
pinlv=ff/freq
****************
********************
********************
********************
********************
n0=2
maxfax=12*(n0+1)*eps0*c0/(2*(imax-iimax)*xx)
maxfay=12*(n0+1)*eps0*c0/(2*(jmax-jjmax)*yy)
maxfaz=12*(n0+1)*eps0*c0/(2*(kmax-kkmax)*zz)
do i=imin,imax
if(i.ge.iimin.and.i.le.iimax)then
fax(i)=0.0
famx(i)=0.0
end if
if(i.ge.iimax)then
famx(i)=maxfax*((i-iimax+0.5)/real(imax-iimax))**n0*xmu0/eps0
end if
if(i.gt.iimax)then
fax(i)=maxfax*((i-iimax)/real(imax-iimax))**n0
end if
if(i.lt.iimin)then
fax(i)=maxfax*((iimin-i)/real(iimin-imin))**n0
famx(i)=maxfax*((iimin-i-0.5)/real(iimin-imin))**n0*xmu0/eps0
end if
end do
do j=jmin,jmax
if(j.ge.jjmin.and.j.le.jjmax)then
fay(j)=0.0
famy(j)=0.0
end if
if(j.ge.jjmax)then
famy(j)=maxfay*((j-jjmax+0.5)/real(jmax-jjmax))**n0*xmu0/eps0
end if
if(j.gt.jjmax)then
fay(j)=maxfay*((j-jjmax)/real(jmax-jjmax))**n0
end if
if(j.lt.jjmin)then
fay(j)=maxfay*((jjmin-j)/real(jjmin-jmin))**n0
famy(j)=maxfay*((jjmin-j-0.5)/real(jjmin-jmin))**n0*xmu0/eps0
end if
end do
do k=kmin,kmax
if(k.ge.kkmin.and.k.le.kkmax)then
faz(k)=0.0
famz(k)=0.0
end if
if(k.ge.kkmax)then
famz(k)=maxfaz*((k-kkmax+0.5)/real(kmax-kkmax))**n0*xmu0/eps0
end if
if(k.gt.kkmax)then
faz(k)=maxfaz*((k-kkmax)/real(kmax-kkmax))**n0
end if
end do
******************
*****************
do m=1,number
read(0,*) input(m,1),input(m,2),input(m,3)
end do
!!!!!!!!!!! 场区介电参数设置 !!!!!!!!!!!!!!!
do i=imin,imax
do j=jmin,jmax
do k=kmin,kmax
epsx1(i,j,k)=1.0
epsx2(i,j,k)=1.0
epsx3(i,j,k)=1.0
epsx4(i,j,k)=1.0
epsx5(i,j,k)=1.0
epsx6(i,j,k)=1.0
epsx7(i,j,k)=1.0
epsx8(i,j,k)=1.0
if(i.ge.nx1.and.i.le.nx2)then
if(j.ge.ny1.and.j.le.ny2)then
if(k.le.hh)then
epsx1(i,j,k)=e0
epsx2(i,j,k)=e0
epsx3(i,j