6
$\begingroup$

Is there a Fortran subroutine which performs linear interpolation in one-dimenional data? I need something similar to MATLAB function interp1.

$\endgroup$
8
$\begingroup$

There is no built-in Fortran functionality to do linear interpolation. You could either use a library or write your own routine. I haven't tried compiling or testing and my fortran may be a bit rusty, but something like the following should work.

subroutine interp1( xData, yData, xVal, yVal )
! Inputs: xData = a vector of the x-values of the data to be interpolated
!         yData = a vector of the y-values of the data to be interpolated
!         xVal  = a vector of the x-values where interpolation should be performed
! Output: yVal  = a vector of the resulting interpolated values

  implicit none

  real, intent(in) :: xData(:), yData(:), xVal(:)
  real, intent(out) :: yVal(:)
  integer :: inputIndex, dataIndex
  real :: minXdata, minYdata, xRange, weight

  ! Possible checks on inputs could go here
  ! Things you may want to check:
  !   monotonically increasing xData
  !   size(xData) == size(yData)
  !   size(xVal) == size(yVal)

  minXData = xData(1)
  maxXData = xData(size(xData))
  xRange = maxXData - minXData

  do inputIndex = 1, size(xVal)
      ! possible checks for out of range xVal could go here

      ! this will work if x is uniformly spaced, otherwise increment
      ! dataIndex until xData(dataIndex+1)>xVal(inputIndex)
      dataIndex = floor((xVal(inputIndex)-minXData)/xRange);

      weight = (xVal - xData(dataIndex))/(xData(dataIndex+1)-xData(dataIndex));
      yVal(inputIndex) = (1.0-weight)*yData(dataIndex) + ...
                         weight*yData(dataIndex+1);
  end do
end subroutine
$\endgroup$
  • $\begingroup$ In your loop, I think you want to have an if statement whether the point xVal is within the current interval. Your loop should also only go to size(xVal)-1 to avoid the out-of-range access. $\endgroup$ – Wolfgang Bangerth Oct 8 '15 at 12:07
  • 2
    $\begingroup$ @WolfgangBangerth 1) That's why I put the comment in the loop about adding out of range checks, the OP can add those for him/herself 2) No, you want to loop over all x values to be interpolated and fortran indices start from one so the loop is correct. $\endgroup$ – Doug Lipinski Oct 8 '15 at 13:26
2
$\begingroup$

Take a look at the Numerical Recipes book, a classic, found here. Look in the interpolation chapter. You can find an explanation, plus the code. It's in pdf. A simple google search also produced a github library by Thomas Robitaille, where a routine called interp1d exists.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.