anchor.diff:      ! Suppose input_x, input_y are some data with input_x being sorted and increasing 
       ! We are given output_x, and wish to get output_y by linearly interpolating the
       ! original series.
  -    SUBROUTINE linear_interpolation(n_input, input_x, input_y, n_output, output_x, output_y)
  +    PURE SUBROUTINE linear_interpolation(n_input, input_x, input_y, n_output, output_x, output_y, precomputed_gradient)
           INTEGER(ip), INTENT(IN):: n_input, n_output
           REAL(dp), INTENT(IN):: input_x(n_input), input_y(n_input), output_x(n_output)
           REAL(dp), INTENT(OUT):: output_y(n_output)
  +        !REAL(dp), OPTIONAL, INTENT(IN):: precomputed_gradient(n_input-1)
  +        REAL(dp), INTENT(IN):: precomputed_gradient(n_input-1)
   
           INTEGER(ip):: i, j
           REAL(dp):: gradient
  
positive.code:         
          INTEGER(ip), PARAMETER:: ninterp_points = 601
          LOGICAL, PARAMETER:: log_transform_interpolation = .FALSE.
          REAL(dp):: kj_xs(ninterp_points), kj_ys(ninterp_points)
  
          ! Make the kajiura_G_interpolator
          DO i = 1, ninterp_points
  
negative.code:             ! interpolate
              if (input_x(j) > output_x(i)) then
                  if(j > 1) then
                      !if(present(precomputed_gradient)) then
                          gradient = precomputed_gradient(j-1)
                      !else
                      !    gradient = (input_y(j) - input_y(j-1))/(input_x(j) - input_x(j-1))
                      !end if
                          
                      output_y(i) = input_y(j-1) +  gradient * (output_x(i) - input_x(j-1))
                  else
                      output_y(i) = input_y(1)
                  end if
              else
                  if(j < n_input) then
                      !if(present(precomputed_gradient)) then
                          gradient = precomputed_gradient(j)
                      !else
                      !    gradient = (input_y(j+1) - input_y(j))/(input_x(j+1) - input_x(j))
                      !end if
                      output_y(i) = input_y(j) +  gradient * (output_x(i) - input_x(j))
                  else
                      output_y(i) = input_y(n_input)
  