-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathbicubic_module.f90
99 lines (79 loc) · 2.5 KB
/
bicubic_module.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
module bicubic_module
! bicubic interpolation
use kind_module, only : i4b, dp
implicit none
private
! call bcucof to calculate coefficents
! function bcuint returns the interpolated value
! Source: Numerical Recipes
! Author: T. Enomoto
! History: 3 March 2004
real(kind=dp), dimension(4,4), private :: c
public :: bcucof, bcuint, bcuintp
contains
subroutine bcucof(f, fx, fy, fxy, dx, dy)
implicit none
real(kind=dp), dimension(4), intent(in) :: f, fx, fy, fxy
real(kind=dp), intent(in) :: dx, dy
integer(kind=i4b) :: i, j, l
real(kind=dp), dimension(16) :: x, cl
real(kind=dp), dimension(16,16) :: wt
real(kind=dp) :: dxdy
data wt/&
1, 0,-3, 2, 0, 0, 0, 0,-3, 0, 9,-6, 2, 0,-6, 4, &
0, 0, 0, 0, 0, 0, 0, 0, 3, 0,-9, 6,-2, 0, 6,-4, &
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9,-6, 0, 0,-6, 4, &
0, 0, 3,-2, 0, 0, 0, 0, 0, 0,-9, 6, 0, 0, 6,-4, &
0, 0, 0, 0, 1, 0,-3, 2,-2, 0, 6,-4, 1, 0,-3, 2, &
0, 0, 0, 0, 0, 0, 0, 0,-1, 0, 3,-2, 1, 0,-3, 2, &
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 2, 0, 0, 3,-2, &
0, 0, 0, 0, 0, 0, 3,-2, 0, 0,-6, 4, 0, 0, 3,-2, &
0, 1,-2, 1, 0, 0, 0, 0, 0,-3, 6,-3, 0, 2,-4, 2, &
0, 0, 0, 0, 0, 0, 0, 0, 0, 3,-6, 3, 0,-2, 4,-2, &
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 2,-2, &
0, 0,-1, 1, 0, 0, 0, 0, 0, 0, 3,-3, 0, 0,-2, 2, &
0, 0, 0, 0, 0, 1,-2, 1, 0,-2, 4,-2, 0, 1,-2, 1, &
0, 0, 0, 0, 0, 0, 0, 0, 0,-1, 2,-1, 0, 1,-2, 1, &
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1, 0, 0,-1, 1, &
0, 0, 0, 0, 0, 0,-1, 1, 0, 0, 2,-2, 0, 0,-1, 1/
dxdy = dx*dy
do i=1, 4
x(i) = f(i)
x(i+4) = fx(i)*dx
x(i+8) = fy(i)*dy
x(i+12) = fxy(i)*dxdy
end do
do i=1, 16
cl(i) = sum(wt(i,:)*x)
end do
l = 0
do i=1, 4
do j=1, 4
l = l + 1
c(i,j) = cl(l)
end do
end do
end subroutine bcucof
function bcuint(t,u) result (fi)
implicit none
real(kind=dp), intent(in) :: t, u
real(kind=dp) :: fi
integer(kind=i4b) :: i
fi = 0.0_dp
do i=4,1,-1
fi = t*fi+((c(i,4)*u+c(i,3))*u+c(i,2))*u+c(i,1)
end do
end function bcuint
function bcuintp(t,u) result (fi)
implicit none
real(kind=dp), intent(in) :: t, u
real(kind=dp) :: fi
integer(kind=i4b) :: i
real(kind=dp), dimension(4) :: tt
tt = (/1.0_dp-t, 1.0_dp-t, t, t/)
fi = 0.0_dp
do i=4,1,-1
fi = tt(i)*fi+((c(i,4)*u+c(i,3))*u+c(i,2))*u+c(i,1)
end do
end function bcuintp
end module bicubic_module