TammFTN

! 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