anchor.diff:  
       end subroutine
   
  +    !
  +    ! Suppose input_x, input_y are some data with input_x being sorted and decreasing
  +    ! We are given output_x, and wish to get output_y by linearly interpolating the
  +    ! original series.
  +    pure subroutine linear_interpolation_input_x_decreasing(n_input, input_x, input_y, n_output, output_x, output_y)
  +        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(inout):: output_y(n_output)
  +
  +        integer(ip):: i, j
  +        real(dp):: gradient
  +
  +        do i = 1, n_output
  +            ! set j = nearest index to output_x(i) in input_x
  +            call nearest_index_sorted_x_decreasing(n_input, input_x, output_x(i), j)
  +            ! interpolate
  +            if (input_x(j) < output_x(i)) then
  +                if(j > 1) then
  +                    gradient = (input_y(j) - input_y(j-1))/(input_x(j) - input_x(j-1))
  +                    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
  +                    gradient = (input_y(j+1) - input_y(j))/(input_x(j+1) - input_x(j))
  +                    output_y(i) = input_y(j) +  gradient * (output_x(i) - input_x(j))
  +                else
  +                    output_y(i) = input_y(n_input)
  +                end if
  +            end if
  +        end do
  +
  +    end subroutine
  +
   end module
  
positive.code: 
      end subroutine
  
      !
      ! 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.
      pure subroutine linear_interpolation(n_input, input_x, input_y, n_output, output_x, output_y)
          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(inout):: output_y(n_output)
  
negative.code: 
          do i = 1, n_output
              ! set j = nearest index to output_x(i) in input_x
              call nearest_index_sorted_x_increasing(n_input, input_x, output_x(i), j)
              ! interpolate
              if (input_x(j) > output_x(i)) then
                  if(j > 1) then
  