update_grids Subroutine

    subroutine update_grids(grid, pos, dir, d_sdf, packet, mua)
    !! record fluence using path length estimators. Uses voxel grid

        use vector_class
        use photonMod
        use gridMod
        use iarray,     only: phasor, jmean, absorb
        use constants , only : sp
        
        !> grid stores voxel grid information (voxel walls and etc)
        type(cart_grid), intent(IN)    :: grid
        !> dir is the current direction (0,0,1) is up
        type(vector),    intent(IN)    :: dir
        !> d_sdf is the distance to travel in voxel grid
        real(kind=wp),   intent(IN)    :: d_sdf
        !> absoprtion coefficent
        real(kind=wp), optional, intent(IN) :: mua
        !> pos is current position with origin in centre of medium (0,0,0)
        type(vector),    intent(INOUT) :: pos
        !> packet stores the photon related variables
        type(photon),    intent(INOUT) :: packet
        
        complex(kind=sp) :: phasec
        type(vector)  :: old_pos
        logical       :: ldir(3)
        integer       :: celli, cellj, cellk
        real(kind=wp) :: dcell, delta=1e-8_wp, d, mua_real

        if(present(mua))then
            mua_real = mua
        else
            mua_real = 1._wp
        end if

        !convert to different coordinate system. Origin is at lower left corner of fluence grid
        old_pos = vector(pos%x+grid%xmax, pos%y+grid%ymax, pos%z+grid%zmax)
        call update_voxels(grid, old_pos, celli, cellj, cellk)
        packet%xcell = celli
        packet%ycell = cellj
        packet%zcell = cellk

        d = 0._wp
        !if packet outside grid return
        if(celli == -1 .or. cellj == -1 .or. cellk == -1)then
            packet%tflag = .true.
            pos = vector(old_pos%x-grid%xmax, old_pos%y-grid%ymax, old_pos%z-grid%zmax)
            return
        end if
        !move photon through grid updating path length estimators
        do
            ldir = (/.FALSE., .FALSE., .FALSE./)

            dcell = wall_dist(grid, celli, cellj, cellk, old_pos, dir, ldir)
            if(d + dcell > d_sdf)then
                dcell = d_sdf - d
                d = d_sdf
! needs to be atomic so dont write to same array address with more than 1 thread at a time
                packet%phase = packet%phase + dcell
!$omp atomic
                    jmean(celli,cellj,cellk) = jmean(celli,cellj,cellk) + real(dcell, kind=sp)
                call update_pos(grid, old_pos, celli, cellj, cellk, dcell, .false., dir, ldir, delta)
                exit
            else
                d = d + dcell
                packet%phase = packet%phase + dcell
!$omp atomic
                    jmean(celli,cellj,cellk) = jmean(celli,cellj,cellk) + real(dcell, kind=sp)
                call update_pos(grid, old_pos, celli, cellj, cellk, dcell, .true., dir, ldir, delta)
            end if
            if(celli == -1 .or. cellj == -1 .or. cellk == -1)then
                packet%tflag = .true.
                exit
            end if
        end do
        pos = vector(old_pos%x-grid%xmax, old_pos%y-grid%ymax, old_pos%z-grid%zmax)
        packet%xcell = celli
        packet%ycell = cellj
        packet%zcell = cellk

    end subroutine update_grids