anchor.diff:  
           real(dp) :: dx(2), xi, yi, zi, xip1, yip1, zip1, xj(1), yj(1), zj(1)
           integer(ip) :: i, i0, j0, np, j
  +        real(force_double) :: n1, n2
           character(len=charlen) :: burnt
   
           if(present(burn_type)) then
  
positive.code:                 zip1 = z(i+1)
  
                  ! Interpolate np+1 points along the line segment between points i and (i+1)
                  np = ceiling(max((xip1-xi)/dx(1), (yip1-yi)/dx(2))) + 1
                  do j = 0, np
                      xj = xi + (j*(xip1-xi))/np
                      yj = yi + (j*(yip1-yi))/np
                      zj = zi + (j*(zip1-zi))/np
                      call burn_xyz_into_grid(xj, yj, zj, grid, lower_left, upper_right, burnt)
                  end do
              else
  
negative.code:             print*, 'FAIL'
          end if
  
          !print*, 'grid(1:6, 4): ', grid(1:6, 4)
          !print*, 'expected_grid(1:6, 4): ', expected_grid(1:6, 4) 
          !print*, 'grid(6, 5:10): ', grid(6, 5:10)
          !print*, 'expected_grid(6, 5:10): ', expected_grid(6, 5:10) 
  
      end subroutine
  
  end module
  