Fortran 95 source code to compute 3D Spline interpolation.

spline_interpolation.f95

!****************************************************
! 3次スプライン補間
!
! * 入力はテキストファイルをパイプ処理
!     1行目:     点の数
!     2行目以降: 1行に1点(x, y)
!
!   date          name            version
!   2018.12.20    mk-mode.com     1.00 新規作成
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
!****************************************************
!
module const
  ! SP: 単精度(4), DP: 倍精度(8)
  integer,     parameter :: SP = kind(1.0)
  integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
end module const

module spline
  use const
  implicit none
  private
  public :: interpolate

contains
  ! 3次スプライン補間
  !
  ! :param(in) integer(4)  :: n: 点の数
  ! :param(in) real(8)  p(:, :): 点の配列
  subroutine interpolate(n, p)
    use const
    implicit none
    integer(SP), intent(in) :: n
    real(DP),    intent(in) :: p(2, n)
    real(DP) :: h(n - 1), w(n - 2), mtx(n - 1, n - 2), v(n)
    real(DP) :: a(n - 1), b(n - 1), c(n - 1), d(n)
    real(DP)    :: xp, xp_2, xp_3, x, y
    integer(SP) :: idx, i

    ! 初期計算
    h = calc_h(n, p)
    w = calc_w(n, p, h)
    call gen_mtx(n, h, w, mtx)
    call gauss_jordan(n, mtx, v)
    v = cshift(v, -1)
    b = calc_b(n, v)
    a = calc_a(n, v, p)
    d = calc_d(n, p)
    c = calc_c(n, v, p)

    ! 補間
    do i = int(p(1, 1) * 10), int(p(1, n) * 10)
      x = i * 0.1_DP
      idx = find_idx(n, p, x)
      xp   = x - p(1, idx)
      xp_2 = xp * xp
      xp_3 = xp_2 * xp
      y = a(idx) * xp_3 + b(idx) * xp_2 + c(idx) * xp + d(idx)
      print '(F8.4, X, F8.4)', x, y
    end do
  end subroutine interpolate

  ! h 計算
  !
  ! :param(in) integer(4)     n: 点の数
  ! :param(in) real(8)  p(2, n): 点の配列
  ! :return    real(8) h(n - 1): h の配列
  function calc_h(n, p) result(h)
    implicit none
    integer(SP), intent(in) :: n
    real(DP),    intent(in) :: p(2, n)
    real(DP) :: h(n - 1)

    h = p(1, 2:n) - p(1, 1:n-1)
  end function calc_h

  ! w 計算
  !
  ! :param(in) integer(4)     n: 点の数
  ! :param(in) real(8)  p(2, n): 点の配列
  ! :param(in) real(8) h(n - 1): h の配列
  ! :return    real(8) w(n - 2): w の配列
  function calc_w(n, p, h) result(w)
    implicit none
    integer(SP), intent(in) :: n
    real(DP),    intent(in) :: p(2, n), h(n - 1)
    real(DP) :: w(n - 2)

    w = 6.0_DP * ((p(2, 3:n) - p(2, 2:n-1)) &
    & / h(2:n-1) - (p(2, 2:n-1) - p(2, 1:n-2)) / h(1:n-2))
  end function calc_w

  ! 行列生成
  !
  ! :param(in)  integer(4)             n: 点の数
  ! :param(in)  real(8)         h(n - 1): h の配列
  ! :param(in)  real(8)         w(n - 2): w の配列
  ! :param(out) real(8) mtx(n - 1, n - 2: 行列
  subroutine gen_mtx(n, h, w, mtx)
    implicit none
    integer(SP), intent(in)  :: n
    real(DP),    intent(in)  :: h(n - 1), w(n - 1)
    real(DP),    intent(out) :: mtx(n - 1, n - 2)
    integer(SP) :: i

    mtx(:, :) = 0.0_DP
    do i = 1, n - 2
      mtx(i, i)  = 2.0_DP * (h(i) + h(i + 1))
      mtx(n - 1, i) = w(i)
      if (i == 1) cycle
      mtx(i, i - 1) = h(i)
      mtx(i - 1, i) = h(i)
    end do
  end subroutine gen_mtx

  ! 連立一次方程式を解く（Gauss-Jordan 法）
  !
  ! :param(in)  integer(4)             n: 点の数
  ! :param(in)  real(8) mtx(n - 1, n - 2: 行列
  ! :param(out) real(8)             v(n): v の配列
  subroutine gauss_jordan(n, mtx, v)
    implicit none
    integer(SP), intent(in)  :: n
    real(DP),    intent(in)  :: mtx(n - 1, n - 2)
    real(DP),    intent(out) :: v(n)
    integer(SP) :: i, j
    real(DP)    :: mtx_w(n - 1, n - 2)  ! 作業用
    real(DP)    :: p, d

    mtx_w(:, :) = mtx(:, :)
    v(:) = 0.0_DP
    do j = 1, n - 2
      p = mtx_w(j, j)
      mtx_w(j:n-1, j) = mtx_w(j:n-1, j) / p
      do i = 1, n - 2
        if (i == j) cycle
        d = mtx_w(j, i)
        mtx_w(j:n-1, i) = mtx_w(j:n-1, i) - d * mtx_w(j:n-1, j)
      end do
    end do
    v(1:n-2) = mtx_w(n - 1, 1:n-2)
  end subroutine gauss_jordan

  ! a 計算
  !
  ! :param(in) integer(4)     n: 点の数
  ! :param(in) real(8)     v(n): v の配列
  ! :param(in) real(8)  p(2, n): 点の配列
  ! :return    real(8) a(n - 1): a の配列
  function calc_a(n, v, p) result(a)
    implicit none
    integer(SP), intent(in) :: n
    real(DP),    intent(in) :: v(n), p(2, n)
    real(DP) :: a(n - 1)

    a = (v(2:n) - v(1:n-1)) / (6.0_DP * (p(1, 2:n) - p(1, 1:n-1)))
  end function calc_a

  ! b 計算
  !
  ! :param(in) integer(4)     n: 点の数
  ! :param(in) real(8)     v(n): v の配列
  ! :return    real(8) b(n - 1): b の配列
  function calc_b(n, v) result(b)
    implicit none
    integer(SP), intent(in) :: n
    real(DP),    intent(in) :: v(n)
    real(DP) :: b(n - 1)

    b = v(1:n-1) / 2.0_DP
  end function calc_b

  ! c 計算
  !
  ! :param(in) integer(4)     n: 点の数
  ! :param(in) real(8)     v(n): v の配列
  ! :param(in) real(8)  p(2, n): 点の配列
  ! :return    real(8)     c(n): c の配列
  function calc_c(n, v, p) result(c)
    implicit none
    integer(SP), intent(in) :: n
    real(DP),    intent(in) :: v(n), p(2, n)
    real(DP) :: c(n - 1)

    c = (p(2, 2:n) - p(2, 1:n-1)) / (p(1, 2:n) - p(1, 1:n-1)) &
    & - (p(1, 2:n) - p(1, 1:n-1)) * (2.0_DP * v(1:n-1) + v(2:n)) &
    & / 6.0_DP
  end function calc_c

  ! d 計算
  !
  ! :param(in) integer(4)     n: 点の数
  ! :param(in) real(8)  p(2, n): 点の配列
  ! :return    real(8)     d(n): d の配列
  function calc_d(n, p) result(d)
    implicit none
    integer(SP), intent(in) :: n
    real(DP),    intent(in) :: p(2, n)
    real(DP) :: d(n)

    d = p(2, :)
  end function calc_d

  ! インデックス検索
  !
  ! :param(in) integer(4)     n: 点の数
  ! :param(in) real(8)  p(2, n): 点の配列
  ! :param(in) real(8)        x: x 値
  ! :return    integer(4)   idx: インデックス
  function find_idx(n, p, x) result(idx)
    implicit none
    integer(SP), intent(in) :: n
    real(DP),    intent(in) :: p(2, n), x
    integer(SP) :: idx
    integer(SP) :: i, j, k

    i = 1
    j = n
    do while (i < j)
      k = (i + j) / 2
      if (p(1, k) < x) then
        i = k + 1
      else
        j = k
      end if
    end do
    if (i > 1) i = i - 1
    idx = i
  end function find_idx
end module spline

program spline_interpolation
  use const
  use spline
  implicit none
  character(20), parameter :: F_SRC = "src.txt"
  integer(SP),   parameter :: UID_S = 10
  integer(SP) :: n                  ! 点の数
  real(DP), allocatable :: p(:, :)  ! 点の配列
  integer(SP) :: i, ios

  ! ファイル OPEN
  open(UID_S, file = F_SRC, status = 'old')

  ! 点の数の取得
  read(UID_S, *) n

  ! 点の配列メモリ確保
  allocate(p(2, n))

  ! 点の配列読み込み
  do i = 1, n
    read(UID_S, *) p(:, i)
  end do

  ! ファイル CLOSE
  close(UID_S)

  ! 補間
  call interpolate(n, p)

  ! 点の配列メモリ解放
  deallocate(p)
  
  stop
end program spline_interpolation