tauint2 Subroutine

    subroutine tauint2(grid, packet, sdfs_array)
    !! optical depth integration subroutine
    !! Moves photons to interaction location
    !! Calculated is any reflection or refraction happens whilst moving
    !
        use gridMod,      only : cart_grid
        use photonMod,    only : photon
        use random,       only : ran2
        use sdfs,         only : sdf, calcNormal
        use surfaces,     only : reflect_refract
        use vector_class, only : vector

        type(cart_grid),   intent(in)    :: grid
        type(photon),      intent(inout) :: packet
        type(sdf),   intent(in)    :: sdfs_array(:)

        real(kind=wp) :: tau, d_sdf, t_sdf, taurun, ds(size(sdfs_array)), dstmp(size(sdfs_array))
        real(kind=wp) :: eps, dtot, old(size(sdfs_array)), new(size(sdfs_array)), n1, n2, Ri
        integer       :: i, oldlayer, new_layer
        type(vector)  :: pos, dir, oldpos, N
        logical       :: rflag

        !setup temp variables
        pos = packet%pos
        oldpos = pos
        dir = vector(packet%nxp, packet%nyp, packet%nzp)


        !round off distance
        eps = 1e-8_wp
        !get random tau
        tau = -log(ran2())
        taurun = 0.
        dtot = 0.
        do
            !setup sdf distance and current layer
            ds = 0.
            do i = 1, size(ds)
                ds(i) = abs(sdfs_array(i)%evaluate(pos))
            end do
            packet%cnts = packet%cnts + size(ds)
            d_sdf = minval(ds)

            if(d_sdf < eps)then
                packet%tflag=.true.
                exit
            end if

            do while(d_sdf > eps)
                t_sdf = d_sdf * sdfs_array(packet%layer)%getkappa()
                if(taurun + t_sdf <= tau)then
                    !move full distance to sdf surface
                    taurun = taurun + t_sdf
                    oldpos = pos
                    !comment out for phase screen
                    call update_grids(grid, oldpos, dir, d_sdf, packet, sdfs_array(packet%layer)%getmua())
                    pos = pos + d_sdf * dir
                    dtot = dtot + d_sdf
                else
                    !run out of tau so move remaining tau and exit
                    d_sdf = (tau - taurun) / sdfs_array(packet%layer)%getkappa()
                    dtot = dtot + d_sdf
                    taurun = tau
                    oldpos = pos
                    pos = pos + d_sdf * dir
                    !comment out for phase screen
                    call update_grids(grid, oldpos, dir, d_sdf, packet, sdfs_array(packet%layer)%getmua())
                    exit
                end if
                ! get distance to nearest sdf
                ds = 0._wp
                do i = 1, size(ds)
                    ds(i) = sdfs_array(i)%evaluate(pos)
                end do
                d_sdf = minval(abs(ds),dim=1)
                packet%cnts = packet%cnts + size(ds)

                !check if outside all sdfs
                if(minval(ds) >= 0._wp)then
                    packet%tflag = .true.
                    exit
                end if
            end do

            !exit early if conditions met
            if(taurun >= tau .or. packet%tflag)then
                exit
            end if
            
            ds = 0._wp
            do i = 1, size(ds)
                ds(i) = sdfs_array(i)%evaluate(pos)
            end do
            packet%cnts = packet%cnts + size(ds)
            
            dstmp = ds
            ds = abs(ds)
            
            !step a bit into next sdf to get n2
            d_sdf = minval(ds) + 2._wp*eps
            oldpos = pos
            pos = pos + d_sdf*dir
            ds = 0._wp
            do i = 1, size(ds)
                ds(i) = sdfs_array(i)%evaluate(pos)
            end do
            packet%cnts = packet%cnts + size(ds)
            
            new = 0._wp
            old = 0._wp
            do i = 1, size(ds)
                if(dstmp(i) < 0.)then
                    old(i)=-1._wp
                    exit
                end if
            end do
            do i = 1, size(ds)
                if(ds(i) < 0.)then
                    new(i)=-1._wp     
                    exit
                end if
            end do

            !check for fresnel reflection
            n1 = sdfs_array(packet%layer)%getn()
            new_layer = minloc(new, dim=1)

            n2 = sdfs_array(new_layer)%getn()
            !carry out refelction/refraction
            if (n1 /= n2)then
                !get correct sdf normal
                if(ds(packet%layer) < 0._wp .and. ds(new_layer) < 0._wp)then
                    oldlayer = minloc(abs([ds(packet%layer), ds(new_layer)]), dim=1)
                elseif(dstmp(packet%layer) < 0._wp .and. dstmp(new_layer) < 0._wp)then
                    oldlayer=maxloc([dstmp(packet%layer), dstmp(new_layer)],dim=1)
                elseif(ds(packet%layer) > 0._wp .and. ds(new_layer) < 0._wp)then
                    oldlayer = packet%layer
                elseif(ds(packet%layer) > 0._wp .and. ds(new_layer) > 0._wp)then
                    packet%tflag = .true.
                    exit
                else
                    error stop "This should not be reached!"
                end if
                if(oldlayer == 1)then
                    oldlayer = packet%layer
                else
                    oldlayer = new_layer
                end if
                N = calcNormal(pos, sdfs_array(oldlayer))

                rflag = .false.
                call reflect_refract(dir, N, n1, n2, rflag, Ri)
                packet%weight = packet%weight * Ri
                tau = -log(ran2())
                taurun = 0._wp
                if(.not.rflag)then
                    packet%layer = new_layer
                else
                    !step back inside original sdf
                    pos = oldpos
                    !reflect so incrment bounce counter
                    packet%bounces = packet%bounces + 1
                    if(packet%bounces > 1000)then
                        packet%tflag=.true.
                        exit
                    end if
                end if
            else
                packet%layer = new_layer
            end if
            if(packet%tflag)exit
        end do

        packet%pos = pos
        packet%nxp = dir%x
        packet%nyp = dir%y
        packet%nzp = dir%z

        packet%phi = atan2(dir%y, dir%x)
        packet%sinp = sin(packet%phi)
        packet%cosp = cos(packet%phi)

        packet%cost = dir%z
        packet%sint = sqrt(1._wp - packet%cost**2)

        ! packet%step = dtot
        if(abs(packet%pos%x) > grid%xmax)then
            packet%tflag = .true.
        end if
        if(abs(packet%pos%y) > grid%ymax)then
            packet%tflag = .true.
        end if
        if(abs(packet%pos%z) > grid%zmax)then
            packet%tflag = .true.
        end if
    end subroutine tauint2