-----------------------------------------------------------------------------
! ONYX - 2
! This is the second release version of our Order N photonics code
! Modified for transmission and reflection calculations 
! (c) Andrew Ward, Imperial College, London. 1997, 1998, 1999.
! All comments, bugs and inquiries to a.j.ward@ic.ac.uk
!
! $Revision: 2.1 $
! $Date: 1999/07/30 14:44:39 $
! $Source: /home/andreww/f90/orderN/PML/onyx.f90,v $
! $Author: andreww $
! -----------------------------------------------------------------------------
module interface1
interface

subroutine init_store_pts_band(store_pts)
implicit none
integer,pointer :: store_pts(:,:)
end subroutine init_store_pts_band

subroutine init_store_pts_trans(store_pts)
implicit none
integer,pointer :: store_pts(:,:)
end subroutine init_store_pts_trans

subroutine init_store_pts_dos(store_pts,ix_cur,iy_cur,iz_cur,i_pol)
implicit none
integer,pointer :: store_pts(:,:)
integer,intent(in) :: ix_cur,iy_cur,iz_cur,i_pol
end subroutine init_store_pts_dos

subroutine initfields(g1,g2,g3,e,h,eps_hat,mu_hat,ix,iy,iz,i_pol)
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,iy,iz,i_pol
end subroutine initfields

subroutine initfields_trans(g1,g2,g3,e,h,eps_hat,mu_hat,ix,iy,iz,i_pol)
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,iy,iz,i_pol
end subroutine initfields_trans

subroutine initfields_dos(g1,g2,g3,e,h,eps_hat,mu_hat,ix,iy,iz,i_pol)
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,iy,iz,i_pol
end subroutine initfields_dos

subroutine defcell(eps,mu,sigma,sigma_m,u1,u2,u3)
implicit none
complex,pointer :: eps(:,:,:),mu(:,:,:)
real,pointer :: sigma(:,:,:),sigma_m(:,:,:)
real :: u1(3),u2(3),u3(3)
end subroutine defcell

subroutine renorm(eps,mu,sigma,sigma_m,eps_inv,mu_inv,eps_hat,mu_hat,g, &
&                 omega)
implicit none
complex,pointer :: eps(:,:,:),mu(:,:,:)
complex,pointer :: eps_inv(:,:,:,:,:)
complex,pointer :: mu_inv(:,:,:,:,:)
complex,pointer :: eps_hat(:,:,:,:,:)
complex,pointer :: mu_hat(:,:,:,:,:)
real,pointer :: sigma(:,:,:),sigma_m(:,:,:)
real,intent(in) :: g(3,3),omega
end subroutine renorm

subroutine driver(e,h,eps_inv,mu_inv,eps_hat,mu_hat,sigma,sigma_m, &
&                 fft_data,store_pts)
implicit none
complex,pointer :: e(:,:,:,:)
complex,pointer :: h(:,:,:,:)
complex,pointer :: eps_inv(:,:,:,:,:)
complex,pointer :: mu_inv(:,:,:,:,:)
complex,pointer :: eps_hat(:,:,:,:,:)
complex,pointer :: mu_hat(:,:,:,:,:)
complex,pointer :: fft_data(:,:)
real,pointer :: sigma(:,:,:),sigma_m(:,:,:)
integer,pointer :: store_pts(:,:)
end subroutine driver

subroutine postproc_band(fft_data,spectrum)
implicit none
complex,pointer :: fft_data(:,:)
real,pointer :: spectrum(:)
end subroutine postproc_band

subroutine postproc_trans(fft_data,eps_inv,mu_inv)
implicit none
complex,pointer :: fft_data(:,:)
complex,pointer :: eps_inv(:,:,:,:,:),mu_inv(:,:,:,:,:)
end subroutine postproc_trans

subroutine postproc_dos(fft_data,spectrum)
implicit none
complex,pointer :: fft_data(:,:)
real,pointer :: spectrum(:)
end subroutine postproc_dos

subroutine finalproc_dos(spectrum)
implicit none
real,pointer :: spectrum(:)
end subroutine finalproc_dos

end interface
end module interface1

! -----------------------------------------------------------------------------
module interface2
interface

subroutine int(e_cur,h_cur,e_prev,h_prev,eps_inv,mu_inv,sigma,sigma_m,   &
&              n_tsteps,iblock,fft_data,store_pts,intcurle_1,intcurle_2, &
&              intcurlh_1,intcurlh_2,wz_1,wz_2)
implicit none
integer,intent(in) :: n_tsteps,iblock
complex,pointer :: e_cur(:,:,:,:)
complex,pointer :: e_prev(:,:,:,:)
complex,pointer :: h_cur(:,:,:,:)
complex,pointer :: h_prev(:,:,:,:)
complex,pointer :: eps_inv(:,:,:,:,:)
complex,pointer :: mu_inv(:,:,:,:,:)
complex,pointer :: fft_data(:,:)
real,pointer :: sigma(:,:,:),sigma_m(:,:,:)
integer,intent(in) :: store_pts(:,:)
complex,pointer :: intcurle_1(:,:,:)
complex,pointer :: intcurlh_1(:,:,:)
complex,pointer :: intcurle_2(:,:,:)
complex,pointer :: intcurlh_2(:,:,:)
real,intent(in) :: wz_1(:),wz_2(:)
end subroutine int

subroutine set_wz(wz_1,wz_2)
implicit none
real,intent(out) :: wz_1(:)
real,intent(out) :: wz_2(:)
end subroutine set_wz

end interface
end module interface2

! -----------------------------------------------------------------------------
module interface3
interface

subroutine calc_div_B(h,mu_hat,div_B)
implicit none
complex,pointer :: h(:,:,:,:)
complex,pointer :: mu_hat(:,:,:,:,:)
complex,pointer :: div_B(:,:,:)
end subroutine calc_div_B

subroutine calc_div_D(e,eps_hat,div_D)
implicit none
complex,pointer :: e(:,:,:,:)
complex,pointer :: eps_hat(:,:,:,:,:)
complex,pointer :: div_D(:,:,:)
end subroutine calc_div_D

subroutine calc_energy_density(E_cur,E_prev,H_prev,rho,eps_hat,mu_hat)
implicit none
complex,pointer :: E_cur(:,:,:,:)
complex,pointer :: E_prev(:,:,:,:)
complex,pointer :: H_prev(:,:,:,:)
complex,pointer :: eps_hat(:,:,:,:,:)
complex,pointer :: mu_hat(:,:,:,:,:)
real,pointer :: rho(:,:,:)
end subroutine calc_energy_density

subroutine calc_current(J,e_cur,h_cur,h_prev,ix,iy,iz)
implicit none
complex,intent(out) :: J(3)
complex,pointer :: e_cur(:,:,:,:),h_cur(:,:,:,:)
complex,pointer :: h_prev(:,:,:,:)
integer,intent(in) :: ix,iy,iz
end subroutine calc_current

end interface
end module interface3

! -----------------------------------------------------------------------------
module interface4
interface

subroutine power_spec(f,spectrum)
implicit none
complex,pointer :: f(:,:)
real,pointer :: spectrum(:)
end subroutine power_spec


end interface
end module interface4

! -----------------------------------------------------------------------------
module interface5
interface

subroutine bc_xmin_bloch(e,h)
implicit none
complex,pointer :: e(:,:,:,:)
complex,pointer :: h(:,:,:,:)
end subroutine bc_xmin_bloch

subroutine bc_xmax_bloch(e,h)
implicit none
complex,pointer :: e(:,:,:,:)
complex,pointer :: h(:,:,:,:)
end subroutine bc_xmax_bloch

subroutine bc_ymin_bloch(e,h)
implicit none
complex,pointer :: e(:,:,:,:)
complex,pointer :: h(:,:,:,:)
end subroutine bc_ymin_bloch

subroutine bc_ymax_bloch(e,h)
implicit none
complex,pointer :: e(:,:,:,:)
complex,pointer :: h(:,:,:,:)
end subroutine bc_ymax_bloch

subroutine bc_zmin_bloch(e,h)
implicit none
complex,pointer :: e(:,:,:,:)
complex,pointer :: h(:,:,:,:)
end subroutine bc_zmin_bloch

subroutine bc_zmax_bloch(e,h)
implicit none
complex,pointer :: e(:,:,:,:)
complex,pointer :: h(:,:,:,:)
end subroutine bc_zmax_bloch

subroutine bc_xmin_metal(e,h)
implicit none
complex,pointer :: e(:,:,:,:)
complex,pointer :: h(:,:,:,:)
end subroutine bc_xmin_metal 

subroutine bc_xmax_metal(e,h)
implicit none
complex,pointer :: e(:,:,:,:)
complex,pointer :: h(:,:,:,:)
end subroutine bc_xmax_metal 

subroutine bc_ymin_metal(e,h)
implicit none
complex,pointer :: e(:,:,:,:)
complex,pointer :: h(:,:,:,:)
end subroutine bc_ymin_metal 

subroutine bc_ymax_metal(e,h)
implicit none
complex,pointer :: e(:,:,:,:)
complex,pointer :: h(:,:,:,:)
end subroutine bc_ymax_metal 

subroutine bc_zmin_metal(e,h)
implicit none
complex,pointer :: e(:,:,:,:)
complex,pointer :: h(:,:,:,:)
end subroutine bc_zmin_metal 

subroutine bc_zmax_metal(e,h)
implicit none
complex,pointer :: e(:,:,:,:)
complex,pointer :: h(:,:,:,:)
end subroutine bc_zmax_metal 

subroutine PML_e(e_cur,e_prev,h_prev,intcurlh_1,intcurlh_2,wz_1,wz_2,eps_inv)
implicit none
complex,pointer :: e_cur(:,:,:,:),e_prev(:,:,:,:),h_prev(:,:,:,:)
complex,pointer :: intcurlh_1(:,:,:),intcurlh_2(:,:,:)
complex,pointer :: eps_inv(:,:,:,:,:)
real,intent(in) :: wz_1(:),wz_2(:)
end subroutine PML_e

subroutine PML_h(e_cur,h_cur,h_prev,intcurle_1,intcurle_2,wz_1,wz_2,mu_inv)
implicit none
complex,pointer :: e_cur(:,:,:,:),h_cur(:,:,:,:),h_prev(:,:,:,:)
complex,pointer :: intcurle_1(:,:,:),intcurle_2(:,:,:)
complex,pointer :: mu_inv(:,:,:,:,:)
real,intent(in) :: wz_1(:),wz_2(:)
end subroutine PML_h

end interface
end module interface5

! -----------------------------------------------------------------------------
module interface6
interface

subroutine defpw(iw,rvec,lvec,kz,eps_inv,mu_inv)
implicit none
integer,intent(in) :: iw
complex,pointer :: rvec(:,:),lvec(:,:),kz(:)
complex,pointer :: eps_inv(:,:,:,:,:),mu_inv(:,:,:,:,:)
end subroutine defpw

subroutine FFT(f,isign)
implicit none
integer,intent(in) :: isign
complex,pointer :: f(:,:)
end subroutine FFT

subroutine FT(f,isign)
implicit none
integer,intent(in) :: isign
complex,pointer :: f(:,:)
end subroutine FT

end interface
end module interface6

! -----------------------------------------------------------------------------
module matrices
interface

function matinv3(A,emach,fail)
implicit none
complex,dimension(3,3) :: matinv3
complex,intent(in) :: A(3,3)
real,intent(in) :: emach
logical,intent(out) :: fail
end function matinv3

function cross_prod(a,b)
implicit none
real,dimension(3) :: cross_prod
real,intent(in) :: a(3),b(3)
end function cross_prod

function c_cross_prod(a,b)
implicit none
complex,dimension(3) :: c_cross_prod
complex,intent(in) :: a(3),b(3)
end function c_cross_prod

function matinv(A,emach,fail)
implicit none
complex,pointer :: A(:,:)
complex,dimension(ubound(A,1),ubound(A,2)) :: matinv
real,intent(in) :: emach
logical,intent(out) :: fail
end function matinv

end interface
end module matrices

! -----------------------------------------------------------------------------
module matrices2
interface

subroutine lu_decomp(A,D,indx,emach,fail)
implicit none
complex,pointer :: A(:,:)
real,intent(in) :: emach
logical,intent(out) :: fail
integer,intent(out) :: D
integer,pointer :: indx(:)
end subroutine lu_decomp

subroutine lu_bksub(A,indx,B)
implicit none
complex,pointer :: A(:,:)
complex,pointer :: B(:)
integer,pointer :: indx(:)
end subroutine lu_bksub

end interface
end module matrices2

! -----------------------------------------------------------------------------
module physconsts
! Use Atomic units
! c0=Speed of light
! pi=3.1415...
! emach=Machine accuracy
real,parameter :: pi=3.14159265358979323846
real,parameter :: c0=137.0360
real,parameter :: eps0=1.0/(4.0*pi)
real,parameter :: mu0=1.0/(eps0*c0**2)
real,parameter :: abohr=5.291772e-11
real :: x
real,parameter :: emach=epsilon(x)
complex,parameter :: ci=(0.0,1.0)
end module physconsts

! -----------------------------------------------------------------------------
module files
! Set the unit numbers for the files we want to access
integer,parameter :: logfile=10,infile=11,outfile=12,fields=13
integer,parameter :: ref=14,trans=15
end module files

! -----------------------------------------------------------------------------
module parameters
! Define the run-time parameter list. These will be read in from a file
! ixmax, iymax, izmax define the number of mesh points in each of the 
! principle directions
! itmax and dt are the number and size of timesteps respectively
! akx, aky, akz are the components of the k-vector expressed on the reciprocal
! lattice. They are normalised by the mesh spacing in the relevent direction
! ie. akx = (Q1*ixmax)*kx
! ikmax = No. of k-points
! Q(3) is the matrix holding the mesh spacing in each direction
! n_block : Break down itmax into `n_block' blocks each of `block_size' time 
! steps. Fields can be writen out, div's calcuated etc at the end of each block
! itmax=n_block*block_size
! n_pts_store : The number of real space points at which we store the fields
! For the FFT the itmax timesteps can be split up into segments which can be 
! either overlapping or not. For details see Numerical Receipes pg 425ff.
! If overlap=.true. then itmax=fft_size*(n_segment+1)
! If overlap=.false. then itmax=fft_size*(2*n_segment)
! In both cases fft_size is the size the segment which is fourier transformed
! so must be a power of 2. 
! nw = Number of frequencies to be analyized / written out
! bclayer1, bclayer2 = thickness of the absorbing layer at either end of system
! epsref is the dielectric constant of the host material
! Notice that bclayer1, bclayer2 and epsref will be over-ridden in defcell
integer :: ixmax,iymax,izmax,itmax,ikmax,n_block,block_size,n_pts_store
real :: Q0,Q(3),dt,akx,aky,akz
integer :: n_segment,fft_size,nw
logical :: overlap
real :: damping
integer :: bclayer1=2
integer :: bclayer2=2
complex :: epsref=(1.0,0.0)
end module parameters

! -----------------------------------------------------------------------------

program onyx
use interface1
use parameters
use physconsts
use files
implicit none

! Define necessary variables

integer :: ios,ik,ix_cur,iy_cur,iz_cur,i_pol
integer,pointer :: store_pts(:,:)
complex,pointer :: e(:,:,:,:)
complex,pointer :: h(:,:,:,:)
complex,pointer :: eps(:,:,:),mu(:,:,:)
complex,pointer :: eps_inv(:,:,:,:,:),mu_inv(:,:,:,:,:)
complex,pointer :: eps_hat(:,:,:,:,:),mu_hat(:,:,:,:,:)
complex,pointer :: fft_data(:,:)
real,pointer :: sigma(:,:,:),sigma_m(:,:,:)
real,pointer :: spectrum(:)
real :: u1(3),u2(3),u3(3)
real :: g1(3),g2(3),g3(3)
real :: g(3,3),omega

! Open the neccessary files

open(unit=logfile,file='log.dat',status='replace',action='write')
open(unit=outfile,file='out.dat',status='replace',action='write')
open(unit=fields,file='fields.dat',status='replace',action='write')
open(unit=ref,file='ref.dat',status='replace',action='write')
open(unit=trans,file='trans.dat',status='replace',action='write')
open(unit=infile,file='infile.dat',status='old',action='read',iostat=ios)
if (ios.ne.0) call err(1)

! Set the run time parameters

call setparam()