Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Fortran 95 source code to compute 3D Spline interpolation.
!****************************************************
! 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.