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
|
! Copyright (C) 2018 J. K. Dewhurst, S. Sharma and E. K. U. Gross.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.
real(8) function splint(n,x,f)
implicit none
! arguments
integer, intent(in) :: n
real(8), intent(in) :: x(n),f(n)
! local variables
integer i
real(8) x0,x1,x2,x3,y0,y1,y2,y3
real(8) t0,t1,t2,t3,t4,t5,t6,t7
if (n.lt.5) then
write(*,*)
write(*,'("Error(splint): n < 5 : ",I8)') n
write(*,*)
stop
end if
! fit piecewise cubic spline to data and integrate
y0=f(1)
y1=f(2)-y0; y2=f(3)-y0; y3=f(4)-y0
x0=x(1)
x1=x(2)-x0; x2=x(3)-x0; x3=x(4)-x0
t1=x1*x2*y3; t2=x2*x3*y1; t3=x3*x1*y2
t4=x1-x2; t5=x1-x3; t6=x2-x3
t0=0.5d0/(x1*x3*t4*t5*t6)
t7=t1*t4+t2*t6-t3*t5
t4=x1**2; t5=x2**2; t6=x3**2
y1=t3*t6-t1*t5; y3=t2*t5-t3*t4; y2=t1*t4-t2*t6
t1=x1*y1+x2*y2+x3*y3
t2=y1+y2+y3
splint=x2*(y0+t0*(t1+x2*(0.5d0*t7*x2-0.6666666666666666667d0*t2)))
do i=3,n-3
y0=f(i)
y1=f(i-1)-y0; y2=f(i+1)-y0; y3=f(i+2)-y0
x0=x(i)
x1=x(i-1)-x0; x2=x(i+1)-x0; x3=x(i+2)-x0
t1=x1*x2*y3; t2=x2*x3*y1; t3=x3*x1
t4=x1-x2; t5=x1-x3; t6=x2-x3
t0=0.5d0/(t3*t4*t5*t6)
t3=t3*y2
t7=t1*t4+t2*t6-t3*t5
t4=x1**2; t5=x2**2; t6=x3**2
y1=t3*t6-t1*t5; y2=t1*t4-t2*t6; y3=t2*t5-t3*t4
t1=x1*y1+x2*y2+x3*y3
t2=y1+y2+y3
splint=splint+x2*(y0+t0*(t1+x2*(0.5d0*t7*x2-0.6666666666666666667d0*t2)))
end do
y0=f(n-2)
y1=f(n-3)-y0; y2=f(n-1)-y0; y3=f(n)-y0
x0=x(n-2)
x1=x(n-3)-x0; x2=x(n-1)-x0; x3=x(n)-x0
t1=x1*x2*y3; t2=x2*x3*y1; t3=x3*x1*y2
t4=x1-x2; t5=x1-x3; t6=x2-x3
t0=0.5d0/(x1*x2*t4*t5*t6)
t7=t1*t4+t2*t6-t3*t5
t4=x1**2; t5=x2**2; t6=x3**2
y1=t3*t6-t1*t5; y2=t1*t4-t2*t6; y3=t2*t5-t3*t4
t1=x1*y1+x2*y2+x3*y3
t2=y1+y2+y3
splint=splint+x3*(y0+t0*(t1+x3*(0.5d0*t7*x3-0.6666666666666666667d0*t2)))
return
end function
|
