7

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

kyperros
  • 139
  • 1
  • 1
  • 5
  • There are many options. Search for interpolation at http://people.sc.fsu.edu/~%20jburkardt/f_src/f_src.html – stali Oct 08 '15 at 15:01

2 Answers2

10

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
Community
  • 1
Doug Lipinski
  • 4,591
  • 16
  • 25
  • 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. – Wolfgang Bangerth Oct 08 '15 at 12:07
  • 2
    @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. – Doug Lipinski Oct 08 '15 at 13:26
3

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.

cauchi
  • 897
  • 1
  • 9
  • 12
  • First link is broken. The official website is: http://www.numerical.recipes/ This includes some free ebooks - yay! Unfortunately, you need Adobe Flash or Adobe Acrobat to read them - boo! – Biggsy Nov 14 '19 at 10:47