My old Master's code. I had planned to rewrite this, and re-reading it has given me a little more motivation to actually do so. Oh, how much I've improved in merely a year...

photonics.f90

! This file is part of TammFTN.
! Copyright (C) 2013 Adam Hirst <adam@aphirst.karoo.co.uk>
!
! 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)
    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)
    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)
    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)
    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)
    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)
    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)
    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)
    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
    ! atan2(r%im, r%re)
    ! ENSURE -180 RENDERS AS 180
    if (phi == -180.0) then
      phi = 180.0
    end if
  end subroutine reftrans

  subroutine rtenergy(Ecs, materl, pol, res, struct, theta, output)
    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)
    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)
    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)

    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)
    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)
    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
    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