fortran实现一维线性插值：与matlab的interp1的线性插值功能一致

program main
        implicit none
        integer :: i
        real*8  :: x(9), y(9), t, z
        
        !// x与y的数据来源于matlab的interp1帮助文档, t为被插值节点，见帮助文档的变量xq
        x = [0.d0,0.785398163397448d0,1.57079632679490d0,2.35619449019235d0,3.14159265358979d0,3.92699081698724d0,4.71238898038469d0,5.49778714378214d0,6.28318530717959d0]
        y = [0.d0,0.707106781186548d0,1.d0,0.707106781186548d0,0.d0,-0.707106781186548d0,-1.d0,-0.707106781186548d0,0.d0]
        
        open( 101, file = 'interp_m.dat' )
        open( 102, file = 'interp_f.dat' )
        do i = 1, 33
                read(101,*) t
                z = 0.d0
                call interp1( x, y, size(x), t, z )
                write(102,*) t, z
        end do
        close( 101 )
        close( 102 )
        
end program main
        
        
subroutine interp1( x, y, n, t, tmp )
        implicit none
        integer :: i, j, k1, k2
        integer :: n
        real*8  :: x(n), y(n), t, tmp, k

        !// 如果插值节点与被插值节点相同，则被插值节点处的值等于原节点处的值
        do i = 1, n
                if ( abs( t - x(i) ) < 1d-15 ) then
                        tmp = y(i)
                        return
                end if
        end do
        
        !// 寻找被插值节点属于的区间[X_i, X_i+1]
        i = 1
        do
                if ( t < x(i) ) then
                        exit
                else
                        i = i + 1
                end if
        end do
        
        k2 = i; k1 = k2 - 1

        k = ( y(k2) - y(k1) ) / ( x(k2) - x(k1) )
        tmp = y(k1) + k*(t-x(k1))

End subroutine interp1