1

I'm porting some MATLAB code to Fortran and need to replicate the functionality of scatteredInterpolant. scatteredInterpolant takes a set of sample points and returns what is essentially a function handle that can take a new point and return an interpolated value.

I could do this by returning a derived type with an "interpolate" function or something similar, but I'd like to keep the Fortran as close to the original MATLAB as possible. Is there a way to generate a function at runtime and return some way of accessing the function from another function?

mark2
  • 23
  • 5

2 Answers2

2

You could likely achieve the desired behavior the object-oriented way,

module Interpolation_mod

    use iso_fortran_env, only: RK => real64
    implicit none

    type :: Interpolation_type
        real(RK)              :: linearInterpSlope, linearInterpIntercept
        real(RK), allocatable :: X(:), Value(:)
    contains
        procedure, pass :: functionHandle
    end type Interpolation_type

    interface Interpolation_type
        module procedure :: scatteredInterpolant
    end interface Interpolation_type

contains

    ! creates an Interpolation object and sets up all the essentials of calling functionHandle
    function scatteredInterpolant(X,Value) result(InterpolationObject)
        implicit none
        real(RK), intent(in)        :: X(2), Value(2)
        type(Interpolation_type)    :: InterpolationObject
        InterpolationObject%X = X
        InterpolationObject%Value = Value
        InterpolationObject%linearInterpSlope = ( Value(2) - Value(1) ) / (X(2) - X(1))
        InterpolationObject%linearInterpIntercept = Value(1) - InterpolationObject%linearInterpSlope * X(1)
    end function scatteredInterpolant

    ! equivalent of MATLAB functionHandle
    function functionHandle(InterpolationObject,x) result(functionValue)
        implicit none
        class(Interpolation_type), intent(in)   :: InterpolationObject
        real(RK), intent(in)                    :: x
        real(RK)                                :: functionValue
        if ( x < InterpolationObject%X(1) .or. x > InterpolationObject%X(2) ) then
            write(*,"(A)") "This is interpolation, not extrapolation."
            error stop
        else
            functionValue = InterpolationObject%linearInterpSlope * x + InterpolationObject%linearInterpIntercept
        end if
    end function functionHandle

end module Interpolation_mod

program test_Interpolation
use Interpolation_mod, only: RK, Interpolation_type
implicit none
type(Interpolation_type) :: Interpolation
Interpolation = Interpolation_type( X = [0._RK,1._RK], Value = [0._RK,1._RK] )
write(*,"(*(g0,:,' '))") "Interpolation(", 0.5, ") =", Interpolation%functionHandle(0.5_RK)
end program test_Interpolation
C:\> ifort main.f90 -o main.exe

Intel(R) Visual Fortran Intel(R) 64 Compiler for applications running on Intel(R) 64, Version 19.0.4.245 Build 20190417
Copyright (C) 1985-2019 Intel Corporation.  All rights reserved.

Microsoft (R) Incremental Linker Version 14.22.27905.0
Copyright (C) Microsoft Corporation.  All rights reserved.

-out:main.exe
-subsystem:console
main.obj

C:\> main.exe
Interpolation( .5000000 ) = .5000000000000000
Scientist
  • 1,767
  • 2
  • 12
  • 20
0

What you would like is actually a closure of a certain sort. That is not possible in Fortran. Also, Fortran syntax does not allow calling a derived type directly. You will need one of the indirect approaches you propose. You can also use a module, but that only allows one instance of such closure/function.

Fortran does allow returning aa function pointer, but that is not enough for what you describe. See also related techniques in Fortran minimization of a function with additional arguments and in functions linked there under Related. It will provably not do exactly what you need, but it is related.