onyx.f90

计算光子晶体能带的程序

read(infile,*) izmax
if (izmax<1) call err(2)
write(logfile,'(A1,A12,I10)') '#','izmax=',izmax
write(outfile,'(A1,A12,I10)') '#','izmax=',izmax
write(fields,'(A1,A12,I10)') '#','izmax=',izmax
write(trans,'(A1,A12,I10)') '#','izmax=',izmax
write(ref,'(A1,A12,I10)') '#','izmax=',izmax

read(infile,*) Q(1)
if (Q(1)<=0.0) call err(2)
write(logfile,'(A1,A12,F15.5)') '#','Q(1)=',Q(1)
write(outfile,'(A1,A12,F15.5)') '#','Q(1)=',Q(1)
write(fields,'(A1,A12,F15.5)') '#','Q(1)=',Q(1)
write(trans,'(A1,A12,F15.5)') '#','Q(1)=',Q(1)
write(ref,'(A1,A12,F15.5)') '#','Q(1)=',Q(1)

read(infile,*) Q(2)
if (Q(2)<=0.0) call err(2)
write(logfile,'(A1,A12,F15.5)') '#','Q(2)=',Q(2)
write(outfile,'(A1,A12,F15.5)') '#','Q(2)=',Q(2)
write(fields,'(A1,A12,F15.5)') '#','Q(2)=',Q(2)
write(trans,'(A1,A12,F15.5)') '#','Q(2)=',Q(2)
write(ref,'(A1,A12,F15.5)') '#','Q(2)=',Q(2)

read(infile,*) Q(3)
if (Q(3)<=0.0) call err(2)
write(logfile,'(A1,A12,F15.5)') '#','Q(3)=',Q(3)
write(outfile,'(A1,A12,F15.5)') '#','Q(3)=',Q(3)
write(fields,'(A1,A12,F15.5)') '#','Q(3)=',Q(3)
write(trans,'(A1,A12,F15.5)') '#','Q(3)=',Q(3)
write(ref,'(A1,A12,F15.5)') '#','Q(3)=',Q(3)

read(infile,*) dt
if (dt<=0.0) call err(2)

write(logfile,'(A1,A12,F15.5)') '#','dt=',dt
write(outfile,'(A1,A12,F15.5)') '#','dt=',dt
write(fields,'(A1,A12,F15.5)') '#','dt=',dt
write(trans,'(A1,A12,F15.5)') '#','dt=',dt
write(ref,'(A1,A12,F15.5)') '#','dt=',dt

read(infile,*) akx
if (akx<0.0.or.akx>pi) call err(2)
write(logfile,'(A1,A12,F10.5)') '#','akx=',akx
write(outfile,'(A1,A12,F10.5)') '#','akx=',akx
write(fields,'(A1,A12,F10.5)') '#','akx=',akx
write(trans,'(A1,A12,F10.5)') '#','akx=',akx
write(ref,'(A1,A12,F10.5)') '#','akx=',akx

read(infile,*) aky
if (aky<0.0.or.aky>pi) call err(2)
write(logfile,'(A1,A12,F10.5)') '#','aky=',aky
write(outfile,'(A1,A12,F10.5)') '#','aky=',aky
write(fields,'(A1,A12,F10.5)') '#','aky=',aky
write(trans,'(A1,A12,F10.5)') '#','aky=',aky
write(ref,'(A1,A12,F10.5)') '#','aky=',aky

read(infile,*) akz
if (akz<0.0.or.akz>pi) call err(2)
write(logfile,'(A1,A12,F10.5)') '#','akz=',akz
write(outfile,'(A1,A12,F10.5)') '#','akz=',akz
write(fields,'(A1,A12,F10.5)') '#','akz=',akz
write(trans,'(A1,A12,F10.5)') '#','akz=',akz
write(ref,'(A1,A12,F10.5)') '#','akz=',akz

read(infile,*) ikmax
if (ikmax<0) call err(2)
write(logfile,'(A1,A12,I10)') '#','ikmax=',ikmax
write(outfile,'(A1,A12,I10)') '#','ikmax=',ikmax
write(fields,'(A1,A12,I10)') '#','ikmax=',ikmax
write(trans,'(A1,A12,I10)') '#','ikmax=',ikmax
write(ref,'(A1,A12,I10)') '#','ikmax=',ikmax

read(infile,*) n_block
if (n_block<1) call err(2)
write(logfile,'(A1,A12,I10)') '#','n_block=',n_block
write(outfile,'(A1,A12,I10)') '#','n_block=',n_block
write(fields,'(A1,A12,I10)') '#','n_block=',n_block
write(trans,'(A1,A12,I10)') '#','n_block=',n_block
write(ref,'(A1,A12,I10)') '#','n_block=',n_block

read(infile,*) block_size
if (block_size<1) call err(2)
write(logfile,'(A1,A12,I10)') '#','block_size=',block_size
write(outfile,'(A1,A12,I10)') '#','block_size=',block_size
write(fields,'(A1,A12,I10)') '#','block_size=',block_size
write(trans,'(A1,A12,I10)') '#','block_size=',block_size
write(ref,'(A1,A12,I10)') '#','block_size=',block_size

itmax=n_block*block_size
write(logfile,'(A1,A12,I10)') '#','itmax=',itmax
write(outfile,'(A1,A12,I10)') '#','itmax=',itmax
write(fields,'(A1,A12,I10)') '#','itmax=',itmax
write(trans,'(A1,A12,I10)') '#','itmax=',itmax
write(ref,'(A1,A12,I10)') '#','itmax=',itmax

read(infile,*) n_pts_store
if (n_pts_store<0) call err(2)
write(logfile,'(A1,A12,I10)') '#','n_pts_store=',n_pts_store
write(outfile,'(A1,A12,I10)') '#','n_pts_store=',n_pts_store
write(fields,'(A1,A12,I10)') '#','n_pts_store=',n_pts_store
write(trans,'(A1,A12,I10)') '#','n_pts_store=',n_pts_store
write(ref,'(A1,A12,I10)') '#','n_pts_store=',n_pts_store

read(infile,*) n_segment
if (n_segment<1) call err(2)
write(logfile,'(A1,A12,I10)') '#','n_segment=',n_segment
write(outfile,'(A1,A12,I10)') '#','n_segment=',n_segment
write(fields,'(A1,A12,I10)') '#','n_segment=',n_segment
write(trans,'(A1,A12,I10)') '#','n_segment=',n_segment
write(ref,'(A1,A12,I10)') '#','n_segment=',n_segment

read(infile,*) overlap
write(logfile,'(A1,A12,L5)') '#','overlap=',overlap
write(outfile,'(A1,A12,L5)') '#','overlap=',overlap
write(fields,'(A1,A12,L5)') '#','overlap=',overlap
write(trans,'(A1,A12,L5)') '#','overlap=',overlap
write(ref,'(A1,A12,L5)') '#','overlap=',overlap

if (overlap) then
  fft_size=itmax/(n_segment+1)
else
  fft_size=itmax/(2*n_segment)
endif

if (fft_size<1) call err(2)

nbits=bit_size(fft_size)
set_bits=0
i=0
do
   if (btest(fft_size,i)) set_bits=set_bits+1
   if (i==nbits.or.set_bits>1) exit
   i=i+1
enddo
if (set_bits>1) call err(2)

read(infile,*) damping
write(logfile,'(A1,A12,F10.5)') '#','damping=',damping
write(outfile,'(A1,A12,F10.5)') '#','damping=',damping
write(fields,'(A1,A12,F10.5)') '#','damping=',damping
write(trans,'(A1,A12,F10.5)') '#','damping=',damping
write(ref,'(A1,A12,F10.5)') '#','damping=',damping
if (damping<0.0) call err(2)

read(infile,*) nw
write(logfile,'(A1,A12,I10)') '#','nw=',nw
write(outfile,'(A1,A12,I10)') '#','nw=',nw
write(fields,'(A1,A12,I10)') '#','nw=',nw
write(trans,'(A1,A12,I10)') '#','nw=',nw
write(ref,'(A1,A12,I10)') '#','nw=',nw
if (nw<1) call err(2)

return
end subroutine setparam

! -----------------------------------------------------------------------------
! Subroutine initfields
! Set up the initial values of the E & H fields
! Should be able to read fields in from a file 
! or set them to some appropriate analytic form
! -----------------------------------------------------------------------------
subroutine initfields(g1,g2,g3,e,h,eps_hat,mu_hat,ix_cur,iy_cur,iz_cur,i_pol)
use interface3
use interface5
use matrices
use files
use parameters
use physconsts
implicit none

real,intent(in) :: g1(3),g2(3),g3(3)
complex,pointer :: e(:,:,:,:)
complex,pointer :: h(:,:,:,:)
complex,pointer :: eps_hat(:,:,:,:,:),mu_hat(:,:,:,:,:)
integer,intent(in) :: ix_cur,iy_cur,iz_cur,i_pol

integer :: ix,iy,iz,jx,jy,jz,i
real :: G(3),k(3),v(3),k_plus_G(3)
complex :: vec(3),h0(3),expo
complex,pointer :: div(:,:,:)
real :: dtcq2

write(logfile,*) 'Initialize fields in initfields'

dtcq2=(dt*c0/Q0)**2
e=(0.0,0.0)
h=(0.0,0.0)

! Set up initial fields

! Initial fields following Chan, Yu and Ho's Order N paper
! Can be used for Band structure

v=(/1.0,0.0,0.0/)

! Loops over jx,jy,jz loop over the required reciprocal lattice vectors

do jz=-5,5
do jy=-5,5
do jx=-5,5

   G=2.0*pi*(g1*(jx)/real(ixmax)+g2*(jy)/real(iymax)+g3*(jz)/real(izmax))
   k=(akx/real(ixmax))*g1+(aky/real(iymax))*g2+(akz/real(izmax))*g3
   k_plus_G=k+G
   vec(1)=exp(ci*(k_plus_G(1)))-1.0
   vec(2)=exp(ci*(k_plus_G(2)))-1.0
   vec(3)=exp(ci*(k_plus_G(3)))-1.0
   h0(1)=v(2)*vec(3)-v(3)*vec(2)
   h0(2)=v(3)*vec(1)-v(1)*vec(3)
   h0(3)=v(1)*vec(2)-v(2)*vec(1)

   do iz=1,izmax
   do iy=1,iymax
   do ix=1,ixmax

      expo=exp(ci*(k_plus_G(1)*ix+k_plus_G(2)*iy+k_plus_G(3)*iz))
      do i=1,3
        h(i,ix,iy,iz)=h(i,ix,iy,iz)+h0(i)*expo
      enddo

   enddo
   enddo
   enddo

enddo
enddo
enddo

! Gaussian field profile

!!$do ix=1,ixmax+1
!!$do iy=1,iymax+1
!!$do iz=1,izmax+1
!!$   e(1,ix,iy,iz)=(1.0,0.0)*exp(-((iz-iz_cur)**2)/50.0) 
!!$enddo
!!$enddo
!!$enddo

! Delta funtion
! Needed for DOS, can be used for band structure

!!$if (i_pol<4) then
!!$  e(i_pol,ix_cur,iy_cur,iz_cur)=-dt*(0.0,1.0)*Q(i_pol)
!!$else
!!$ h(i_pol-3,ix_cur,iy_cur,iz_cur)=-dt*(0.0,1.0)*Q(i_pol-3)
!!$endif

! Initialise the boundaries

!call bc_xmax_metal(e,h)
!call bc_ymax_metal(e,h)
!call bc_zmax_metal(e,h)

!call bc_xmin_metal(e,h)
!call bc_ymin_metal(e,h)
!call bc_zmin_metal(e,h)

call bc_xmax_bloch(e,h)
call bc_ymax_bloch(e,h)
call bc_zmax_bloch(e,h)

call bc_xmin_bloch(e,h)
call bc_ymin_bloch(e,h)
call bc_zmin_bloch(e,h)

! Test that the fields are divergence free here?

allocate(div(ixmax,iymax,izmax))
call calc_div_B(h,mu_hat,div)
write(logfile,*) 'Sum div B',sum(abs(div))
!if (maxval(abs(div))>10.0e+3*emach) call err(6)
call calc_div_D(e,eps_hat,div)
write(logfile,*) 'sum div D',sum(abs(div))
!if (maxval(abs(div))>10.0e+3*emach) call err(5)
deallocate(div)

return
end subroutine initfields

! -----------------------------------------------------------------------------
! Subroutine initfields_dos
! Set up the initial values of the E & H fields
! Should be able to read fields in from a file 
! or set them to some appropriate analytic form
! Specifically tailored for the DOS calculation
! -----------------------------------------------------------------------------
subroutine initfields_dos(g1,g2,g3,e,h,eps_hat,mu_hat, &
&                         ix_cur,iy_cur,iz_cur,i_pol)
use interface3
use interface5
use matrices
use files
use parameters
use physconsts
implicit none

real,intent(in) :: g1(3),g2(3),g3(3)
complex,pointer :: e(:,:,:,:)
complex,pointer :: h(:,:,:,:)
complex,pointer :: eps_hat(:,:,:,:,:),mu_hat(:,:,:,:,:)
integer,intent(in) :: ix_cur,iy_cur,iz_cur,i_pol

complex,pointer :: div(:,:,:)
real :: dtcq2

write(logfile,*) 'Initialize fields in initfields_dos'

dtcq2=(dt*c0/Q0)**2
e=(0.0,0.0)
h=(0.0,0.0)

! Delta funtion
! Needed for DOS, can be used for band structure

if (i_pol<4) then
  e(i_pol,ix_cur,iy_cur,iz_cur)=-dt*(0.0,1.0)*Q(i_pol)
else
 h(i_pol-3,ix_cur,iy_cur,iz_cur)=-dt*(0.0,1.0)*Q(i_pol-3)
endif

! Initialise the boundaries

!call bc_xmax_metal(e,h)
!call bc_ymax_metal(e,h)
!call bc_zmax_metal(e,h)

!call bc_xmin_metal(e,h)
!call bc_ymin_metal(e,h)
!call bc_zmin_metal(e,h)

call bc_xmax_bloch(e,h)
call bc_ymax_bloch(e,h)
call bc_zmax_bloch(e,h)

call bc_xmin_bloch(e,h)
call bc_ymin_bloch(e,h)
call bc_zmin_bloch(e,h)

! Test that the fields are divergence free here?

allocate(div(ixmax,iymax,izmax))
call calc_div_B(h,mu_hat,div)
write(logfile,*) 'Sum div B',sum(abs(div))
!if (maxval(abs(div))>10.0e+3*emach) call err(6)
call calc_div_D(e,eps_hat,div)
write(logfile,*) 'sum div D',sum(abs(div))
!if (maxval(abs(div))>10.0e+3*emach) call err(5)
deallocate(div)

return
end subroutine initfields_dos

! -----------------------------------------------------------------------------
! Subroutine initfields_trans
! Set up the initial values of the E & H fields
! Should be able to read fields in from a file 
! or set them to some appropriate analytic form
! Cannot use delta function to start the trans/ref calculation as need to
! know k parallel for incident beam. But can use a delta function in the 
! time domain (ie set the fields at a single layer of points)
! -----------------------------------------------------------------------------
subroutine initfields_trans(g1,g2,g3,e,h,eps_hat,mu_hat, &
&                         ix_cur,iy_cur,iz_cur,i_pol)
use interface3
use interface5
use matrices
use files
use parameters
use physconsts
implicit none

real,intent(in) :: g1(3),g2(3),g3(3)
complex,pointer :: e(:,:,:,:)
complex,pointer :: h(:,:,:,:)
complex,pointer :: eps_hat(:,:,:,:,:),mu_hat(:,:,:,:,:)
integer,intent(in) :: ix_cur,iy_cur,iz_cur,i_pol

integer :: ix,iy,iz
complex,pointer :: div(:,:,:)
real :: dtcq2

write(logfile,*) 'Initialize fields in initfields_trans'

dtcq2=(dt*c0/Q0)**2

e=(0.0,0.0)
h=(0.0,0.0)

! Set up initial fields
! Need to fix incident kx, ky and then have a pulse with
! all frequencies - ie. delta function at t=0
! Gaussian field profile
! Note if the pulse is too wide the normalization at higher
! frequences will fail as 0/0

! Note only one incident beam is considered in this version

do ix=1,ixmax
do iy=1,iymax
do iz=1,izmax
   if (i_pol==1) then
   e(1,ix,iy,iz)=(1.0,0.0)*exp(-((iz-iz_cur)**2)/1.0) &
&                         *exp(ci*(akx*ix/real(ixmax)+aky*iy/real(iymax)))
   h(2,ix,iy,iz)=(dt*c0/(Q0))*e(1,ix,iy,iz)
   else
   e(2,ix,iy,iz)=(1.0,0.0)*exp(-((iz-iz_cur)**2)/1.0) &
&                         *exp(ci*(akx*ix/real(ixmax)+aky*iy/real(iymax)))
   h(1,ix,iy,iz)=-(dt*c0/(Q0))*e(2,ix,iy,iz)
   endif
enddo
enddo
enddo

! Initialise the boundaries

call bc_zmax_bloch(e,h)
call bc_zmin_bloch(e,h)

call bc_xmax_bloch(e,h)
call bc_ymax_bloch(e,h)

call bc_xmin_bloch(e,h)
call bc_ymin_bloch(e,h)

! Test that the fields are divergence free here?

allocate(div(ixmax,iymax,izmax))
call calc_div_B(h,mu_hat,div)
write(logfile,*) 'Sum div B',sum(abs(div))
!if (maxval(abs(div))>10.0e+3*emach) call err(6)
call calc_div_D(e,eps_hat,div)
write(logfile,*) 'sum div D',sum(abs(div))
!if (maxval(abs(div))>10.0e+3*emach) call err(5)
deallocate(div)

return
end subroutine initfields_trans

! -----------------------------------------------------------------------------
! Subroutine defcell
! Define the eps, mu, sigma and sigma_m for our unit cell
! Either analytically or from a file
! Defines an empty cell
! -----------------------------------------------------------------------------
!!$subroutine defcell(eps,mu,sigma,sigma_m,u1,u2,u3)
!!$use files
!!$use parameters
!!$use physconsts
!!$implicit none
!!$
!!$complex,pointer :: eps(:,:,:),mu(:,:,:)
!!$real,pointer :: sigma(:,:,:),sigma_m(:,:,:)
!!$real :: u1(3),u2(3),u3(3)
!!$
!!$write(logfile,*) 'Define unit cell'
!!$
!!$bclayer1=40
!!$bclayer2=40
!!$epsref=(1.0,0.0)
!!$
!!$eps=epsref
!!$mu=(1.0,0.0)
!!$sigma=(0.0,0.0)
!!$sigma_m=(0.0,0.0)
!!$
!!$eps(1,1,200)=10.0
!!$eps(1,1,201)=10.0
!!$eps(2,1,200)=10.0
!!$eps(2,1,201)=10.0
!!$
!!$return
!!$end subroutine defcell

! -----------------------------------------------------------------------------
! Subroutine defcell
! Define the eps, mu, sigma and sigma_m for our unit cell
! Either analytically or from a file