! This file is part of TammFTN.
! Copyright (C) 2013 Adam Hirst <[email protected]>
!
! 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 <http://www.gnu.org/licenses/>.
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
Comments
0 B
|👍
/👎
0 B
|👍
/👎