! This file is part of TammFTN. ! Copyright (C) 2013 Adam Hirst ! ! TammFTN is free software: you can redistribute it and/or modify ! it under the terms of the GNU General Public License as published by ! the Free Software Foundation, either version 3 of the License, or ! (at your option) any later version. ! ! TammFTN is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License ! along with TammFTN. If not, see . module photonics implicit none integer, parameter :: dp = selected_real_kind(15, 307) real(dp), parameter :: c = 299792458_dp real(dp), parameter :: hbar = 1.05457172647E-34_dp real(dp), parameter :: q = 1.60217656535E-19_dp ! elementary charge real(dp), parameter :: e0 = 8.854187817E-12_dp real(dp), parameter :: mu0 = 16*atan(1.0_dp)*1E-7 ! 4 pi 10^-7 real(dp), parameter :: pi = 4*atan(1.0_dp) contains !! subroutine readin(Es, func, materl, pol, res, struct, thetas) implicit none complex(dp), pointer, intent(out) :: Es(:) character(len=5), intent(out) :: func character(len=10), pointer, intent(out) :: materl(:) character(len=2), intent(out) :: pol integer, intent(out) :: res real(dp), pointer, intent(out) :: struct(:), thetas(:) integer :: a, b, i, ioerr, j, gap, n, tokens real(dp), pointer :: Etemp(:) character(len=30) :: filename character(len=10000) :: line character(len=10) :: resstr call getarg(1, filename) open(unit=10, file=trim(filename), status='old', iostat=ioerr) if (ioerr /= 0) then stop 'Error reading input file.' end if call getarg(2, func) ! function to perform, 'RT', 'EHS', 'SOLVE' or 'PHASE' if (index('RT EHS SOLVE PHASE', func) == 0) then stop 'Invalid function specified.' end if call getarg(3, pol) ! polarisation, either 'TE' or 'TM' if (index('TE TM', pol) == 0) then stop 'Invalid polarisation specified.' end if call getarg(4, resstr) ! number of points for graph-data generation, as string read(resstr,*) res ! convert to integer via internal read do i = 1, 4 ! first 4 lines have variable length, materials, thicknesses, energies, angles read(10,'(a)') line ! '(a)' reads entire line ! EXPAND n(A B) INTO A B A B ... if (i .le. 2) then ! only do for materials and thicknesses do a = scan(line, '(') b = scan(line, ')') if (a > b) then stop 'Formatting error, check parentheses.' else if ((a == 0) .or. (b == 0)) then exit ! no more () terms, done end if gap = scan(line(1:a), ' ', .true.) ! search backwards from ( for blank read(line(gap+1:a-1),*) n ! convert to integer using internal read line = line(1:gap) // repeat(line(a+1:b-1) // ' ', n) // line(b+2:) ! b+2 is start of first token after ')' end do end if ! REMOVE DUPLICATED SPACES do gap = index(line(1:len_trim(line)),' ') if (gap == 0) then exit ! no more duplicated spaces, done else line(gap:) = line(gap+1:) end if end do ! COUNT TOKENS IN LINE tokens = 0 gap = 1 do if (len_trim(line(gap:)) == 0) then exit ! break when no tokens in line end if gap = gap + scan(line(gap+1:), ' ') ! find next blank (end of current token) tokens = tokens + 1 end do ! LOAD INTO ARRAYS if (i == 1) then ! materials allocate(materl(tokens)) read(line,*) materl ! internal read else if (i == 2) then ! layer thicknesses, in nanometres if (tokens /= size(materl) - 1) then stop 'Mismatch between materials and thicknesses.' end if allocate(struct(tokens)) read(line,*) struct ! internal read ! converting thicknessses to positions deferred until later else if (i == 3) then ! energies, in eV if (mod(tokens,2) /= 0) then ! n complex numbers => 2n real numbers stop 'Mismatch between number of real and imaginary energy components.' end if allocate(Etemp(tokens)) ! use Etemp as temporary store for Ereal and Eimag read(line,*) Etemp ! internal read allocate(Es(tokens/2)) do j = 1, size(Es) Es(j) = cmplx(Etemp(2*j-1),Etemp(2*j),dp) ! cast real and imaginary components of energy into Es end do else if (i == 4) then ! angles, in degrees if ((size(Es) == 2) .and. (size(thetas) == 2)) then stop '2D contour plot of energy and angle not implemented at this time.' end if allocate(thetas(tokens)) read(line,*) thetas thetas = pi*thetas/180 end if end do close(10) end subroutine readin ! subroutine refindex(Ein, materl, n) implicit none complex(dp), intent(in) :: Ein character(len=*), intent(in) :: materl(:) complex(dp), pointer, intent(out) :: n(:) complex(dp) :: E integer :: i allocate(n(size(materl))) E = cmplx(real(Ein),0,dp) ! only real(E) needed for refractive index do i = 1, size(materl) ! STANDARD DIELECTRICS if (materl(i) == 'air') then n(i) = (1,0) else if (materl(i) == 'glass') then n(i) = (1.5,0) else if (materl(i) == 'prism') then n(i) = (2,0) else if (materl(i) == 'SiO2') then ! Silicon Dioxide n(i) = (1.47,0) else if (materl(i) == 'TiO2') then ! Titanium Dioxide n(i) = (2.37,0) ! MODEL METAL else if (materl(i) == 'modmet') then n(i) = sqrt(1 - 1/(E*(E + (0,0.01)))) ! STANDARD METALS else if (materl(i) == 'Au') then ! Gold n(i) = sqrt(1 - 73.1025/(E*(E + (0,0.0184)))) ! OTHER DRUDE-FORM MATERIALS else if (materl(i) == 'TiN') then !n(i) = sqrt(1 - 33.64/(E*(E + (0,1.1)))) n(i) = sqrt( 1 - 33.64/(E*(E + (0,0.3))) + 33.64/(15.21 - E*(E + (0,0.3))) ) ! 1.1 ! "SPP study of optical dielectric function of TiNx" else if (materl(i) == 'ITO') then n(i) = 1.94936 * sqrt(1 - 1.26113/(E * (E + (0,0.111)))) ! loss 0.111 ! Rhodes et al, full loss term else if (materl(i) == 'ITO_NL') then ! NL = No Losses n(i) = 1.94936 * sqrt(1 - 1.26113/(E * (E + (0,0)))) ! SPECIAL MATERIALS else if (materl(i) == 'TiN_F') then n(i) = sqrt(1.95 - 40.2083/(E**2 + E*(0,0.746)) + (2.16736)/(12.0409 - E**2 & - E*(0,0.91261)) + (229.85)/(33.4084 - E**2 - E*(0,5.85514))) ! TiN Fitted, Applied Surface Science 175-176, p270 else write(*,*) 'Unsupported material: ', materl(i), '.' stop end if end do end subroutine refindex ! subroutine wavevect(n, theta, k) implicit none complex(dp), intent(in) :: n(:) real(dp), intent(in) :: theta complex(dp), pointer, intent(out) :: k(:) integer :: i allocate(k(0:size(n))) ! k'z is at entry 0 k(0) = n(1) * sin(theta) ! fixed at left of structure do i = 1, size(n) k(i) = sqrt(n(i)**2 - k(0)**2) end do end subroutine wavevect ! subroutine interpos(struct) implicit none real(dp), intent(inout) :: struct(:) integer :: i if (size(struct) .gt. 1) then do i = 2, size(struct) struct(i) = struct(i) + struct(i-1) ! convert thicknesses into interface positions end do end if end subroutine interpos ! subroutine trnsmtrx(E, materl, pol, struct, theta, k, n, trnsfr) implicit none complex(dp), intent(in) :: E character(len=*), intent(in) :: materl(:), pol real(dp), intent(in) :: struct(:), theta complex(dp), pointer, intent(out) :: k(:), n(:), trnsfr(:,:,:) complex(dp) :: M(2,2) integer :: i call refindex(E, materl, n) call wavevect(n, theta, k) allocate(trnsfr(size(n),2,2)) trnsfr(1,:,:) = reshape((/1,0,0,1/),(/2,2/)) ! identity matrix, since transferring into layer 1 does not alter coefficients do i = 2, size(n) if (pol == 'TE') then M(1,1) = (k(i) + k(i-1)) * exp(E * cmplx(0,pi*struct(i-1)/620,dp) * (k(i) - k(i-1))) M(1,2) = (k(i) - k(i-1)) * exp(E * cmplx(0,pi*struct(i-1)/620,dp) * (k(i) + k(i-1))) M(2,1) = (k(i) - k(i-1)) * exp(E * cmplx(0,-pi*struct(i-1)/620,dp) * (k(i) + k(i-1))) M(2,2) = (k(i) + k(i-1)) * exp(E * cmplx(0,-pi*struct(i-1)/620,dp) * (k(i) - k(i-1))) ! interface matrix characterising coefficient transfer from (i-1) 'into' (i) ! struct(i-1) since (i-1)th interface precedes (i)th layer trnsfr(i,:,:) = M / (2 * k(i)) else if (pol == 'TM') then M(1,1) = (k(i-1)*n(i)**2 + k(i)*n(i-1)**2) * exp(E * cmplx(0,pi*struct(i-1)/620,dp) & * (k(i) - k(i-1))) M(1,2) = (k(i-1)*n(i)**2 - k(i)*n(i-1)**2) * exp(E * cmplx(0,pi*struct(i-1)/620,dp) & * (k(i-1) + k(i))) M(2,1) = (k(i-1)*n(i)**2 - k(i)*n(i-1)**2) * exp(E * cmplx(0,-pi*struct(i-1)/620,dp) & * (k(i-1) + k(i))) M(2,2) = (k(i-1)*n(i)**2 + k(i)*n(i-1)**2) * exp(E * cmplx(0,pi*struct(i-1)/620,dp) & * (k(i-1) - k(i))) trnsfr(i,:,:) = M / (2 * k(i-1) * n(i)**2) end if trnsfr(i,:,:) = matmul(trnsfr(i,:,:),trnsfr(i-1,:,:)) ! multiply interface matrix with previous transfer matrix end do end subroutine trnsmtrx ! subroutine fieldint(A, B, E, k, kz, n, pol, x, EH) implicit none complex(dp), intent(in) :: A, B, E, k, kz, n character(len=2), intent(in) :: pol real(dp), intent(in) :: x real(dp), intent(out) :: EH(2,2) ! E in row 1, H in row 2 if (pol == 'TE') then ! (:,1) is perpendicular to planar interfaces (x) ! (:,2) is parallel to planar interfaces (y or z) ! E_x, E_y EH(1,1) = 0_dp EH(1,2) = abs(A*exp(k*E*cmplx(0,-pi*x/620,dp)) + B*exp(k*E*cmplx(0,pi*x/620,dp)))**2 ! H_x, H_z EH(2,1) = abs(kz*(A*exp(k*E*cmplx(0,-pi*x/620,dp)) + B*exp(k*E*cmplx(0,pi*x/620,dp))) & / (c*mu0))**2 EH(2,2) = abs(k*(B*exp(k*E*cmplx(0,pi*x/620,dp)) - A*exp(k*E*cmplx(0,-pi*x/620,dp))) & /(c*mu0))**2 else if (pol == 'TM') then ! E_x, E_z EH(1,1) = abs((kz/k)*(A*exp(k*E*cmplx(0,-pi*x/620,dp)) + B*exp(k*E*cmplx(0,pi*x/620,dp))))**2 EH(1,2) = abs(A*exp(k*E*cmplx(0,-pi*x/620,dp)) + B*exp(k*E*cmplx(0,pi*x/620,dp)))**2 ! H_x, H_y EH(2,1) = 0_dp EH(2,2) = abs((A*exp(k*E*cmplx(0,-pi*x/620,dp)) - B*exp(k*E*cmplx(0,pi*x/620,dp))) & * (n**2 * c * e0)/k)**2 end if end subroutine fieldint ! subroutine poweravg(A, B, E, k, n, pol, x, S) implicit none complex(dp), intent(in) :: A, B, E, k, n character(len=*), intent(in) :: pol real(dp), intent(in) :: x real(dp), intent(out) :: S if (pol == 'TE') then S = real((A*exp(k*E*cmplx(0,-pi*x/620,dp)) + B*exp(k*E*cmplx(0,pi*x/620,dp))) & * conjg(k*(B*exp(k*E*cmplx(0,pi*x/620,dp)) - A*exp(k*E*cmplx(0,-pi*x/620,dp)))) / (2*c*mu0)) else if (pol == 'TM') then S = real((A*exp(k*E*cmplx(0,-pi*x/620,dp)) + B*exp(k*E*cmplx(0,pi*x/620,dp))) * (c*e0/2) & * conjg((n**2 / k) * (A*exp(k*E*cmplx(0,-pi*x/620,dp)) - B*exp(k*E*cmplx(0,pi*x/620,dp))))) end if end subroutine poweravg ! subroutine reftrans(E, materl, pol, struct, theta, Rc, Tc, phi) implicit none complex(dp), intent(in) :: E character(len=*), intent(in) :: materl(:), pol real(dp), intent(in) :: struct(:), theta real(dp), intent(out) :: phi, Rc, Tc complex(dp), pointer :: k(:), n(:), trnsfr(:,:,:) complex(dp) :: M(2,2), r, t real(dp) :: SI, ST call trnsmtrx(E, materl, pol, struct, theta, k, n, trnsfr) M = trnsfr(size(n),:,:) r = -M(1,2)/M(1,1) t = M(2,2) - M(1,2)*M(2,1)/M(1,1) Rc = abs(r)**2 call poweravg(cmplx(0,0,dp), t, E, k(size(n)), n(size(n)), pol, struct(1), ST) call poweravg(cmplx(0,0,dp), cmplx(1,0,dp), E, k(1), n(1), pol, struct(size(struct)), SI) Tc = ST/SI phi = 180*atan2(imag(r),real(r))/pi ! ENSURE -180 RENDERS AS 180 if (phi == -180) then phi = 180 end if end subroutine reftrans ! subroutine rtenergy(Ecs, materl, pol, res, struct, theta, output) implicit none complex(dp), intent(in) :: Ecs(2) character(len=*), intent(in) :: materl(:), pol integer, intent(in) :: res real(dp), intent(in) :: struct(:), theta real(dp), pointer, intent(out) :: output(:,:) complex(dp) :: E, Es(2) integer :: i write(*,*) 'Plotting R, T and phi against E...' do i = 1,2 Es(i) = cmplx(real(Ecs(i)),0,dp) ! varying E, so disregard imaginary components end do allocate(output(res,4)) ! E, R, T, phi do i = 1, res E = Es(1) + (Es(2) - Es(1))*(i-1)/res output(i,1) = real(E) call reftrans(E, materl, pol, struct, theta, output(i,2), output(i,3), output(i,4)) ! outputs are (Rc, Tc, phi) end do end subroutine rtenergy ! subroutine rttheta(E, materl, pol, res, struct, thetas, output) implicit none complex(dp), intent(in) :: E character(len=*), intent(in) :: materl(:), pol integer, intent(in) :: res real(dp), intent(in) :: struct(:), thetas(2) real(dp), pointer, intent(out) :: output(:,:) integer :: i real(dp) :: theta write(*,*) 'Plotting R, T and phi against theta...' allocate(output(res,4)) ! theta, R, T, phi do i = 1, res theta = thetas(1) + (thetas(2) - thetas(1))*(i-1)/res output(i,1) = 180*theta/pi call reftrans(E, materl, pol, struct, theta, output(i,2), output(i,3), output(i,4)) ! R, T, phi end do end subroutine rttheta ! subroutine ehsplot(E, func, materl, pol, res, struct, theta, output) implicit none complex(dp), intent(in) :: E character(len=*), intent(in) :: func, materl(:), pol integer, intent(in) :: res real(dp), intent(in) :: struct(:), theta real(dp), pointer, intent(out) :: output(:,:) complex(dp), pointer :: AB(:,:), k(:), n(:), trnsfr(:,:,:) real(dp) :: EH(2,2), S, x integer :: i, j complex(dp) :: M(2,2), r write(*,*) 'Plotting E, H and S through structure...' allocate(output(res,7)) ! x, n, E(1er), E(||), H(1er), H(||), S call trnsmtrx(E, materl, pol, struct, theta, k, n, trnsfr) allocate(AB(size(n),2)) ! SET LHS COEFFICIENTS TO (r,1) OR (1,0) FOR INCIDENT/SURFACE RESPECTIVELY if (func == 'SOLVE') then AB(1,:) = (/cmplx(1,0,dp), cmplx(0,0,dp)/) ! reflected = r, 0 (confined) else M = trnsfr(size(n),:,:) r = -M(1,2)/M(1,1) AB(1,:) = (/r, cmplx(1,0,dp)/) ! reflected, incident end if ! GENERATE ALL COEFFICIENTS do i = 2, size(n) ! no need to recalculate AB(1) AB(i,:) = matmul(trnsfr(i,:,:),AB(1,:)) end do ! CALCULATE E, H AND S FOR ALL x do i = 1, res x = struct(1) + (struct(size(struct)) - struct(1))*(i-1)/res output(i,1) = x do j = 1, size(struct) if (x < struct(j)) then ! .lt. interface i => in material i call fieldint(AB(j,1), AB(j,2), E, k(j), k(0), n(j), pol, x, EH) call poweravg(AB(j,1), AB(j,2), E, k(j), n(j), pol, x, S) output(i,2) = real(n(j)) exit end if end do output(i,3:4) = EH(1,:) ! store E output(i,5:6) = EH(2,:) ! store H output(i,7) = S ! store S end do end subroutine ehsplot ! subroutine surfsolv(E, materl, pol, res, struct, theta, output) implicit none complex(dp), intent(in) :: E character(len=*), intent(in) :: materl(:), pol integer, intent(in) :: res real(dp), intent(in) :: struct(:), theta real(dp), pointer, intent(out) :: output(:,:) integer :: best, i, ioerr complex(dp) :: Es(5) character(len=10) :: filename complex(dp), pointer :: k(:), n(:), trnsfr(:,:,:) real(dp) :: M(5), step write(*,*) 'Solving for defect/surface state...' call getarg(1, filename) open(unit=20, file= trim(filename) // '_log', status='replace', iostat=ioerr) if (ioerr /= 0) then stop 'Error writing to log file.' end if ! DETERMINE BEST SOLUTION OF E step = 0.01 ! eV Es(1) = E do ! GENERATE "CROSS" OF E "POINTS" Es(:) = (/Es(1), Es(1)+step*(0,1), Es(1)+step, Es(1)-step*(0,1), Es(1)-step/) ! centre, up (+i), right (+1), down (-i), left (-1) ! CALCULATE M(1,1) FOR EACH E IN CROSS do i = 1, 5 call trnsmtrx(Es(i), materl, pol, struct, theta, k, n, trnsfr) M(i) = abs(trnsfr(size(n),1,1)) if (i == 1) then write(*,*) real(E)*real(k(0)) end if end do best = minloc(M,1) write(20,'(4ES14.6)') real(Es(best)), imag(Es(best)), step, M(best) ! HALVE STEP SIZE IF BEST E IS THE "CENTRE", ELSE RECENTRE AROUND BEST E if (best == 1) then step = step/2 else Es(1) = Es(best) end if ! SMALL STEP SIZE IMPLIES SUFFICIENT CONVERGENCE if (step < epsilon(step)) then exit end if end do write(*,*) 'Solution found!' write(*,*) 'E = ', real(Es(1)), '-', -imag(Es(1)), 'i' write(*,*) 'Characteristic lifetime is ',-hbar/(q*imag(Es(1))),' seconds.' close(20) ! GENERATE EHS PLOT DATA AT OPTIMUM E call ehsplot(Es(1), 'SOLVE', materl, pol, res, struct, theta, output) end subroutine surfsolv ! subroutine phasmtch(Es, materl, pol, res, struct, theta, output) implicit none complex(dp), intent(in) :: Es(:) character(len=*), intent(in) :: materl(:), pol integer, intent(in) :: res real(dp), intent(in) :: struct(:), theta real(dp), pointer, intent(out) :: output(:,:) integer :: layer character(len=10), pointer :: matlef(:), matrig(:) complex(dp), pointer :: n(:) real(dp), pointer :: strlef(:), strrig(:), outlef(:,:), outrig(:,:) real(dp) :: thetin !write(*,*) 'Generating phase plots for interface matching condition...' allocate(output(res,4)) ! E, phiL, phiR, phiTOT write(*,*) 'First three materials: ',trim(materl(1)),', ',trim(materl(2)),', ',trim(materl(3)),'.' write(*,*) 'In which layer should the virtual interface be constructed (at the far left)?' read(*,*) layer ! read layer # in which to perform phase-matching, as string ! ALLOCATE LEFT AND RIGHT materl(:) AND struct(:) ARRAYS allocate(matlef(layer)) allocate(matrig(size(materl) - layer + 1)) allocate(strlef(layer - 1)) allocate(strrig(size(materl) - layer )) ! CONSTRUCT materl(:) ARRAYS matlef(:) = materl(layer:1:-1) ! third entry indicates "stride", negative reverses matrig(1) = materl(layer) matrig(2:size(matrig)) = materl(layer:size(materl)) ! CONSTRUCT struct(:) ARRAYS strlef(1) = 0 strlef(2:size(strlef)) = struct(layer-1:2:-1) strrig(1) = 0 strrig(2:size(strrig)) = struct(layer:size(struct)) ! CONVERT THICKNESSES TO POSITIONS call interpos(strlef) call interpos(strrig) ! DETERMINE INTERNAL ANGLE WITHIN CONSIDERED INTERFACE call refindex(cmplx(1,0,dp), materl(1:layer:layer-1), n) thetin = asin(real(n(1))*sin(theta)/real(n(2))) write(*,*) 'Internal angle is ',180*thetin/pi,' degrees.' ! PERFORM REFLECTION CALCULATION FOR LEFT AND RIGHT HALF-STRUCTURES call rtenergy(Es, matlef, pol, res, strlef, thetin, outlef) write(*,*) 'Left phase complete!' call rtenergy(Es, matrig, pol, res, strrig, thetin, outrig) write(*,*) 'Right phase complete!' ! EXTRACT PHASES AND CONSTRUCT FINAL OUTPUT ARRAY output(:,1) = outlef(:,1) ! populate Energy column output(:,2) = outlef(:,4) ! phiL output(:,3) = outrig(:,4) ! phiR output(:,4) = outlef(:,4) + outrig(:,4) ! phiL + phiR = phiTOT end subroutine phasmtch ! subroutine writeout(output) implicit none real(dp), intent(in) :: output(:,:) character(len=20) :: filename integer :: i, ioerr, j call getarg(1, filename) write(*,*) 'Writing data to file...' open(unit=30, file=trim(filename) // '_out', status='replace', iostat=ioerr) if (ioerr /= 0) then stop 'Error writing to output file.' end if do i = 1, size(output(:,1)) write(30,'(7ES13.5)') (output(i,j), j = 1, size(output(1,:))) ! format statement with 7 repeats, scientific, 14 width, 6 decimal ! implied do loop end do close(30) write(*,*) 'Data written!' end subroutine writeout ! subroutine flowctrl implicit none complex(dp), pointer :: Es(:) character(len=5) :: func character(len=10), pointer :: materl(:) character(len=2) :: pol integer :: res real(dp), pointer :: output(:,:), struct(:), thetas(:) call readin(Es, func, materl, pol, res, struct, thetas) if (func .ne. 'PHASE') then call interpos(struct) ! convert thicknesses to positions, except for phase-matching end if if ((func == 'RT') .and. (size(thetas) == 1)) then call rtenergy(Es, materl, pol, res, struct, thetas(1), output) else if ((func == 'RT') .and. (size(thetas) == 2)) then call rttheta(Es(1), materl, pol, res, struct, thetas, output) else if (func == 'EHS') then call ehsplot(Es(1), func, materl, pol, res, struct, thetas(1), output) else if (func == 'SOLVE') then call surfsolv(Es(1), materl, pol, res, struct, thetas(1), output) else if (func == 'PHASE') then call phasmtch(Es, materl, pol, res, struct, thetas(1), output) else stop 'Parameters incompatible with desired function, or function unavailable.' end if call writeout(output) end subroutine flowctrl !! end module photonics ! ! program main use photonics implicit none call flowctrl end program main