onyx.f90

计算光子晶体能带的程序

subroutine postproc_dos(fft_data,spectrum)
use files
use physconsts
use parameters
use interface6
implicit none

complex,pointer :: fft_data(:,:)
real,pointer :: spectrum(:)

integer :: i
complex :: ave

write(logfile,*) 'Postprocessing in postproc_dos'

! Eliminate zero frequency part.

ave=sum(fft_data(1,:))/itmax
do i=1,itmax
  fft_data(1,i)=(fft_data(1,i)-ave)*exp(-(i-1)*damping)
enddo

! Call FFT routine to calculate the power spectrum
! For the DOS

call FFT(fft_data,1)
do i=1,fft_size

  spectrum(i)=spectrum(i)-aimag(fft_data(1,i))/pi

enddo

return
end subroutine postproc_dos

! -----------------------------------------------------------------------------
! Subroutine finalproc_dos
! Final post-process step to write out the complete density of states
! -----------------------------------------------------------------------------
subroutine finalproc_dos(spectrum)
use files
use physconsts
use parameters
implicit none

real,pointer :: spectrum(:)

integer :: i

write(logfile,*) 'Final procoessing in finalproc_dos'

do i=1,fft_size

! Write out frequencies in dimensionless units
 write(outfile,*) (i/(dt*2.0*fft_size))*(Q(1)*ixmax)/c0,spectrum(i)

! Write out frequencies in electron volts
!  write(outfile,*) (27.2*2.0*pi*i/(dt*2*fft_size)),spectrum(i)

enddo

write(logfile,*) 'Sum of DOS', sum(spectrum)

return
end subroutine finalproc_dos

! -----------------------------------------------------------------------------
! Subroutine postproc_trans
! Post-process the stored time dependent fields to extract the trans/ref at
! that k-point
! -----------------------------------------------------------------------------
subroutine postproc_trans(fft_data,eps_inv,mu_inv)
use files
use physconsts
use parameters
use interface6
implicit none

complex,pointer :: fft_data(:,:)
complex,pointer :: eps_inv(:,:,:,:,:),mu_inv(:,:,:,:,:)

complex,pointer :: rvec(:,:),lvec(:,:),kz(:)
complex :: sum_in,sum_ref,sum_trans,kx,ky,sum_err
real :: in,out_ref,out_trans,out_err
integer :: i,j,k,isign,iw,ikx,iky,ipol,ix,iy,jpol
complex :: ave
real :: w,wp
logical :: ifail

write(logfile,*) 'Postprocessing in postproc_trans'

allocate(rvec(4*ixmax*iymax,4))
allocate(lvec(4*ixmax*iymax,4))
allocate(kz(4*ixmax*iymax))

! Remove the zero frequency part
do j=1,n_pts_store
ave=sum(fft_data(j,:))/itmax
do i=1,itmax
  fft_data(j,i)=(fft_data(j,i)-ave)*exp(-(i-1)*damping)
enddo
enddo

! Call FFT routine to calculate the frequency spectrum

isign=+1
call fft(fft_data,isign)

! Loop over frequencies and calculate the transmission and reflection
do iw=2,nw,1

w=(2.0*pi*(iw-1)/real(itmax))

call defpw(iw,rvec,lvec,kz,eps_inv,mu_inv)

in=0.0
out_err=0.0
ifail=.true.

! Project out incident beam
! Here we sum over all possible incident beams but this is not necessary
do ikx=1,ixmax
do iky=1,iymax
do ipol=1,2

  kx=(akx+2.0*pi*(ikx-1))/real(ixmax)
  ky=(aky+2.0*pi*(iky-1))/real(iymax)

! j labels the incident mode
  j=3+(ipol-1)
! k labels the PML reflected mode
  k=1+(ipol-1)

  sum_in=0.0
  sum_err=0.0

! Project out incident beams
  if (abs(imag(kz(j+4*(ikx-1+ixmax*(iky-1)))))<10000.0*emach) then
    do ix=1,ixmax
    do iy=1,iymax
    do jpol=1,4
      i=jpol+4*(ix-1+ixmax*(iy-1))
      sum_in = sum_in + fft_data(i,iw) &
&     *lvec(4*(ikx-1+ixmax*(iky-1))+jpol,j)*exp(ci*(kx*(ix)+ky*(iy)))
    enddo
    enddo
    enddo
  endif

! Project out beams reflected by the PML
  if (abs(imag(kz(k+4*(ikx-1+ixmax*(iky-1)))))<10000.0*emach) then
    do ix=1,ixmax
    do iy=1,iymax
    do jpol=1,4
      i=jpol+4*(ix-1+ixmax*(iy-1))
      sum_err = sum_err + fft_data(i+4*ixmax*iymax,iw)&
&     *lvec(4*(ikx-1+ixmax*(iky-1))+jpol,k)*exp(ci*(kx*(ix)+ky*(iy)))
    enddo
    enddo
    enddo
  endif

  in=in+abs(sum_in*conjg(sum_in))
  out_err=out_err+abs(sum_err*conjg(sum_err))

! Cause error if incident current is zero
  if (in>0.0) ifail=.false.

enddo
enddo
enddo

! Error trap - no current carrrying incident modes
if (ifail) call err(10)

out_ref=0.0
out_trans=0.0

do ikx=1,ixmax
do iky=1,iymax
do ipol=1,2

  kx=akx*Q(1)+2.0*pi*(ikx-1)/real(ixmax)
  ky=aky*Q(2)+2.0*pi*(iky-1)/real(iymax)

! j labels the reflected mode
  j=1+(ipol-1)
! k labels the transmitted mode
  k=3+(ipol-1)

  sum_ref=0.0
  sum_trans=0.0

! Project out reflected beams
  if (abs(imag(kz(j+4*(ikx-1+ixmax*(iky-1)))))<10000.0*emach) then
    do ix=1,ixmax
    do iy=1,iymax
    do jpol=1,4
      i=jpol+4*(ix-1+ixmax*(iy-1))
      sum_ref = sum_ref + fft_data(i,iw) &
&     *lvec(4*(ikx-1+ixmax*(iky-1))+jpol,j)*exp(ci*(kx*(ix)+ky*(iy)))
    enddo
    enddo
    enddo
  endif

! Project out transmitted beams
  if (abs(imag(kz(k+4*(ikx-1+ixmax*(iky-1)))))<10000.0*emach) then
    do ix=1,ixmax
    do iy=1,iymax
    do jpol=1,4
      i=jpol+4*(ix-1+ixmax*(iy-1))
      sum_trans = sum_trans + fft_data(i+4*ixmax*iymax,iw)&
&     *lvec(4*(ikx-1+ixmax*(iky-1))+jpol,k)*exp(ci*(kx*(ix)+ky*(iy)))
    enddo
    enddo
    enddo
  endif
  out_ref=out_ref+abs(sum_ref*conjg(sum_ref))
  out_trans=out_trans+abs(sum_trans*conjg(sum_trans))
enddo
enddo
enddo

! Write transmission and reflection intensities

! Define the dimensionless frequency
wp=(iw/(dt*2.0*fft_size))*(Q(1)*ixmax)/c0

write(ref,*) wp,abs(out_ref)/abs(in)
write(trans,*) wp,abs(out_trans)/abs(in)
write(outfile,*) wp,abs(out_err)/abs(in)

! Check the current conservation
write(logfile,*) 'Total Current at freq',wp,abs(out_trans)/abs(in) & 
&                                         + abs(out_ref)/abs(in)

enddo

deallocate(rvec,lvec,kz)

return
end subroutine postproc_trans

! -----------------------------------------------------------------------------
! Subroutine defpw
! Initialize the plane wave basis set
! For the transmission / reflectrion calculation
! -----------------------------------------------------------------------------
subroutine defpw(iw,rvec,lvec,kz,eps_inv,mu_inv)
use files
use parameters
use physconsts
use matrices
implicit none

integer,intent(in) :: iw
complex,pointer :: rvec(:,:),lvec(:,:),kz(:)
complex,pointer :: eps_inv(:,:,:,:,:),mu_inv(:,:,:,:,:)

integer :: i,j,ipol,ikx,iky,imode,jpol
real :: w,dtcq2,kx,ky
real :: wpp
complex :: e(3),h(3),kp(3),km(3),wp,wm,u(3),mag,temp1,temp2,const
logical :: fail

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

! Fix unit normal 
u(1)=(0.0,0.0)
u(2)=(0.0,0.0)
u(3)=(1.0,0.0)

w=2.0*pi*(iw-1)/real(itmax)

! Loop over possible incident beams

do ikx=1,ixmax
do iky=1,iymax

! Use dispersion relation to fix kz for given w and kx,ky
kx=(akx*Q(1)+2.0*pi*(ikx-1))/real(ixmax)
ky=(aky*Q(2)+2.0*pi*(iky-1))/real(iymax)

temp1=1.0+((cos(w)-1.0)*epsref)/dtcq2+(1.0-cos(kx))+(1.0-cos(ky))
temp2=-ci*log(temp1+sqrt(temp1**2-1.0))

if (abs(abs(exp(ci*temp2*Q(3)))-1.0)<100.0*emach) then
  if (sin(real(temp2*Q(3)))<0.0) temp2=-temp2
else
  if (imag(temp2)<0.0) temp2=-temp2
endif

do ipol=1,4

  imode=ipol+4*(ikx-1+ixmax*(iky-1))

  if (ipol==1.or.ipol==2) then

! Waves from +infinity first
    kz(imode)=-temp2
  else

! Wave from -infinity next
    kz(imode)=temp2
  endif

! Write out the free space band structure as a test
!  wpp=((iw-1)/(dt*2.0*fft_size))*(Q(1))/c0
!  if (abs(imag(kz(imode)))<1.0e-11) then
!    write(outfile,*) real(kz(imode)),wpp
!  else
!    write(outfile,*) imag(kz(imode)),wpp
!  endif

! Define unit k vector and w
  kp(1)=(exp(+ci*kx)-1.0)/(ci)
  kp(2)=(exp(+ci*ky)-1.0)/(ci)
  kp(3)=(exp(+ci*kz(imode))-1.0)/(ci)

  km(1)=(1.0-exp(-ci*kx))/(ci)
  km(2)=(1.0-exp(-ci*ky))/(ci)
  km(3)=(1.0-exp(-ci*kz(imode)))/(ci)

  wp=(exp(-ci*w)-1.0)/(-ci)
  wm=(1.0-exp(+ci*w))/(-ci)

  if (ipol==1.or.ipol==3) then

! s polarization
    if (abs(kx)<emach.and.abs(ky)<emach) then

! Normal incidence is a special case
      e(1)=(1.0,0.0)
      e(2)=(0.0,0.0)
      e(3)=(0.0,0.0)
    else
      e=c_cross_prod(km,u)
    endif

! Normalize
  mag=sqrt(sum(e(:)*conjg(e(:))))
  e(:)=e(:)/mag
  h(1)=(dtcq2*(kp(2)*e(3)-kp(3)*e(2))/wm) * mu_inv(1,1,1,1,bclayer1+1)
  h(2)=(dtcq2*(kp(3)*e(1)-kp(1)*e(3))/wm) * mu_inv(2,2,1,1,bclayer1+1)
  h(3)=(dtcq2*(kp(1)*e(2)-kp(2)*e(1))/wm) * mu_inv(3,3,1,1,bclayer1+1)
  else

! p polarization
    if (abs(kx)<emach.and.abs(ky)<emach) then

! Normal incidence is a special case
      h(1)=(1.0,0.0)
      h(2)=(0.0,0.0)
      h(3)=(0.0,0.0)
    else
      h=c_cross_prod(kp,u)
    endif
! Normalize
  mag=sqrt(sum(h(:)*conjg(h(:))))
  h(:)=h(:)/mag
  e(1)=-((km(2)*h(3)-km(3)*h(2))/(wp)) * eps_inv(1,1,1,1,bclayer1+1)
  e(2)=-((km(3)*h(1)-km(1)*h(3))/(wp)) * eps_inv(2,2,1,1,bclayer1+1)
  e(3)=-((km(1)*h(2)-km(2)*h(1))/(wp)) * eps_inv(3,3,1,1,bclayer1+1)

  endif

! Use e and h to define the right eigenvectors

! Normalize to the current
  mag=(exp(ci*kz(imode))*e(1)*conjg(h(2))            &
&     -exp(ci*kz(imode))*e(2)*conjg(h(1))            &
&     +conjg(exp(ci*kz(imode))*e(1))*exp(ci*w)*h(2)  &
&     -conjg(exp(ci*kz(imode))*e(2))*exp(ci*w)*h(1))

  if (abs(real(mag))<1.0e-8) then
    mag=1.0
  else
    mag=abs(mag)
  endif

  rvec(imode,1)=e(1)/sqrt(abs(mag))
  rvec(imode,2)=e(2)/sqrt(abs(mag))
  rvec(imode,3)=h(1)/sqrt(abs(mag))
  rvec(imode,4)=h(2)/sqrt(abs(mag))

  enddo

! Now set up the left vectors - inverse of right set
! Use analytic form for left-vectors given by the transfer matrix

  const=-(wm/wp)/(epsref*dtcq2)
  imode=4*(ikx-1+ixmax*(iky-1))

  do ipol=1,4
    if (ipol==1) jpol=2
    if (ipol==2) jpol=1
    if (ipol==3) jpol=4
    if (ipol==4) jpol=3
    lvec(imode+1,ipol)=      -rvec(jpol+imode,2)
    lvec(imode+2,ipol)=      +rvec(jpol+imode,1)
    lvec(imode+3,ipol)=+const*rvec(jpol+imode,4)
    lvec(imode+4,ipol)=-const*rvec(jpol+imode,3)
  enddo

! Normalise the left-vectors

  do i=1,4
    mag=0.0
    do jpol=1,4
      mag=mag+lvec(imode+jpol,i)*rvec(imode+i,jpol)
    enddo
    do jpol=1,4
      lvec(imode+jpol,i)=lvec(imode+jpol,i)/abs(mag)
    enddo
  enddo

enddo
enddo

! Check left and right vectors are indeed orthonormal
do imode=1,ixmax*iymax
  do i=1,4
  do j=1,4
    mag=0.0
    do jpol=1,4
      mag=mag+lvec(4*(imode-1)+jpol,j)*rvec(4*(imode-1)+i,jpol)
    enddo
    if ((i/=j).and.(abs(mag))>(1.0E+6*emach) &
&   .or.(i==j).and.(abs(mag)-1.0)>(1.0E+6*emach)) then
      write(logfile,*) iw,i,j,abs(mag)
      call err(9)
    endif
  enddo
  enddo
enddo

return
end subroutine defpw

! -----------------------------------------------------------------------------
! Subroutine bc_xmin_bloch
! Set the boundary conditions at ix=1
! -----------------------------------------------------------------------------
subroutine bc_xmin_bloch(e,h)
use physconsts
use parameters
implicit none

complex,pointer :: e(:,:,:,:)
complex,pointer :: h(:,:,:,:)

integer :: iy,iz,i
complex :: fxm1

! Set the Bloch phase shift

fxm1=exp(-ci*akx)

! Loop over the surface at ix=1

do iy=1,iymax
   do iz=1,izmax
      do i=1,3
         e(i,0,iy,iz)=fxm1*e(i,ixmax,iy,iz)
         h(i,0,iy,iz)=fxm1*h(i,ixmax,iy,iz)
      enddo
   enddo
enddo

return
end subroutine bc_xmin_bloch

! -----------------------------------------------------------------------------
! Subroutine bc_