7
$\begingroup$

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

Improve this question
$\endgroup$
1

2 Answers 2

Reset to default
9
$\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
Improve this answer
$\endgroup$
2
  • $\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, 2015 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, 2015 at 13:26
3
$\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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ First link is broken. The official website is: numerical.recipes This includes some free ebooks - yay! Unfortunately, you need Adobe Flash or Adobe Acrobat to read them - boo! $\endgroup$
    – Biggsy
    Nov 14, 2019 at 10:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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