anchor.diff:  
     logical :: debug = .false.
   
  +  real(dp) :: tmin, emin
  +  type(wall_id) :: imin
  +
     public :: escaped
     interface escaped
        module procedure escaped_photon
        module procedure escaped_cell
     end interface escaped
   
  +  public :: next_cell
  +  interface next_cell
  +     module procedure next_cell_int
  +     module procedure next_cell_wall_id
  +  end interface next_cell
  +
     public :: setup_grid_geometry
     type(grid_geometry_desc),public,target :: geo
  -
  -  ! The following list gives the ID of the walls, edges, and vertices
  -  ! xmin             = 1
  -  ! xmax             = 2
  -  ! ymin             = 3
  -  ! ymax             = 4
  -  ! zmin             = 5
  -  ! zmax             = 6
  -  ! xmin, ymin       = 7
  -  ! xmax, ymin       = 8
  -  ! xmin, ymax       = 9
  -  ! xmax, ymax       = 10
  -  ! xmin, zmin       = 11
  -  ! xmax, zmin       = 12
  -  ! xmin, zmax       = 13
  -  ! xmax, zmax       = 14
  -  ! ymin, zmin       = 15
  -  ! ymax, zmin       = 16
  -  ! ymin, zmax       = 17
  -  ! ymax, zmax       = 18
  -  ! xmin, ymin, zmin = 19
  -  ! xmax, ymin, zmin = 20
  -  ! xmin, ymax, zmin = 21
  -  ! xmax, ymax, zmin = 22
  -  ! xmin, ymin, zmax = 23
  -  ! xmax, ymin, zmax = 24
  -  ! xmin, ymax, zmax = 25
  -  ! xmax, ymax, zmax = 26
  -
  -  !    1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
  -  ! 1  x  x  7  9 11 13  x  x  x  x  x  x  x  x 19 21 23 25  x  x  x  x  x  x  x  x
  -  ! 2  x  x  8 10 12 14  x  x  x  x  x  x  x  x 20 22 24 26  x  x  x  x  x  x  x  x
  -  ! 3  7  8  x  x 15 17  x  x  x  x 19 20 23 24  x  x  x  x  x  x  x  x  x  x  x  x
  -  ! 4  9 10  x  x 16 18  x  x  x  x 21 22 25 26  x  x  x  x  x  x  x  x  x  x  x  x
  -  ! 5 11 12 15 16  x  x 19 20 21 22  x  x  x  x  x  x  x  x  x  x  x  x  x  x  x  x
  -  ! 6 13 14 17 18  x  x 23 24 25 26  x  x  x  x  x  x  x  x  x  x  x  x  x  x  x  x
  -
  -  integer,parameter :: combine_wall(0:26, 6) = (/ 1,  0,  0,  7,  9, 11, 13,  0,  0,  0,  0,  0,  0,  0,  0, 19, 21, 23, 25,  0,  0,  0,  0,  0,  0,  0,  0, &
  -       &                                          2,  0,  0,  8, 10, 12, 14,  0,  0,  0,  0,  0,  0,  0,  0, 20, 22, 24, 26,  0,  0,  0,  0,  0,  0,  0,  0, &
  -       &                                          3,  7,  8,  0,  0, 15, 17,  0,  0,  0,  0, 19, 20, 23, 24,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0, &
  -       &                                          4,  9, 10,  0,  0, 16, 18,  0,  0,  0,  0, 21, 22, 25, 26,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0, &
  -       &                                          5, 11, 12, 15, 16,  0,  0, 19, 20, 21, 22,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0, &
  -       &                                          6, 13, 14, 17, 18,  0,  0, 23, 24, 25, 26,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0/)
   
   contains
   
  
positive.code: 
      logical,intent(in) :: radial
  
      integer,intent(out) :: id_min
      ! ID of next wall
  
      real(dp),intent(out)  :: tmin
      ! tmin to nearest wall
  
      real(dp) :: tx, ty, tz
      real(dp) :: tx_s, ty_s, tz_s
  
      logical :: pos_vx, pos_vy, pos_vz
  
      ! Can store this in photon, because it only changes at each interaction
      pos_vx = p%v%x > 0._dp ! whether photon is moving in the +ve x direction
      pos_vy = p%v%y > 0._dp ! whether photon is moving in the +ve y direction
      pos_vz = p%v%z > 0._dp ! whether photon is moving in the +ve z direction
  
      ! Store inv_v to go faster
      if(pos_vx) then
         tx = ( geo%w1(p%icell%i1+1) - p%r%x ) / p%v%x
      else if(p%v%x < 0._dp) then
         tx = ( geo%w1(p%icell%i1)   - p%r%x ) / p%v%x
      else
         tx = huge(1._dp)
      end if
  
      if(pos_vy) then
         ty = ( geo%w2(p%icell%i2+1) - p%r%y ) / p%v%y
      else if(p%v%y < 0._dp) then
         ty = ( geo%w2(p%icell%i2)   - p%r%y ) / p%v%y
      else
         ty = huge(1._dp)
      end if
  
      if(pos_vz) then
         tz = ( geo%w3(p%icell%i3+1) - p%r%z ) / p%v%z
      else if(p%v%z < 0._dp) then
         tz = ( geo%w3(p%icell%i3)   - p%r%z ) / p%v%z
      else
         tz = huge(1._dp)
      end if
  
      ! Following is potential slowdown, in fact, could just test tmin after
      if(min(tx,ty,tz) .lt. 0._dp) call error("find_wall","negative t")
  
      ! Find the closest of the three walls. The following effectively
      ! finds the minimum of three values. A lot of code for such a
      ! small thing, but this runs much much faster than using a built
      ! in min function or any kind of loop. Each iteraction comprises only
      ! three if statements and two pointer assignements.
  
      ! Find out what the spacing between floating-point values is
      tx_s = spacing(tx) * 2._dp
      ty_s = spacing(ty) * 2._dp
      tz_s = spacing(tz) * 2._dp
  
      if(tx.lt.tz - tz_s) then
         if(tx.lt.ty - ty_s) then
            if(pos_vx) then
               id_min = 2
            else
               id_min = 1
            end if
            tmin = tx
         else if(tx.gt.ty + ty_s) then
            if(pos_vy) then
               id_min = 4
            else
               id_min = 3
            end if
            tmin = ty
         else ! tx == ty
            id_min = 7
            if(pos_vx) id_min = id_min + 1
            if(pos_vy) id_min = id_min + 2
            tmin = tx
         end if
      else if(tx.gt.tz + tz_s) then
         if(tz.lt.ty - ty_s) then
            if(pos_vz) then
               id_min = 6
            else
               id_min = 5
            end if
            tmin = tz
         else if(tz.gt.ty + ty_s) then
            if(pos_vy) then
               id_min = 4
            else
               id_min = 3
            end if
            tmin = ty
         else ! ty == tz
            id_min = 15
            if(pos_vy) id_min = id_min + 1
            if(pos_vz) id_min = id_min + 2
            tmin = ty
         end if
      else ! tx == tz
         if(tx.lt.ty - ty_s) then
            id_min = 11
            if(pos_vx) id_min = id_min + 1
            if(pos_vz) id_min = id_min + 2
            tmin = tx
         else if(tx.gt.ty + ty_s) then
            if(pos_vy) then
               id_min = 4
            else
               id_min = 3
            end if
            tmin = ty
         else
            id_min = 19
            if(pos_vx) id_min = id_min + 1
            if(pos_vy) id_min = id_min + 2
            if(pos_vz) id_min = id_min + 4
            tmin = tx
         end if
      end if
  
    end subroutine find_wall
  
  end module grid_geometry_specific
  
negative.code:     escaped_cell = .false.
    end function escaped_cell
  
    type(grid_cell) function next_cell_int(cell, direction, intersection)
      implicit none
      type(grid_cell),intent(in) :: cell
      integer,intent(in) :: direction
  