V3FIT
r8mkspline.f
1  subroutine r8mkspline(x,nx,fspl,ibcxmin,bcxmin,ibcxmax,bcxmax,
2  > ilinx,ier)
3 C
4 C make a 2-coefficient 1d spline
5 C
6 C only 2 coefficients, the data and its 2nd derivative, are needed to
7 C fully specify a spline. See e.g. Numerical Recipies in Fortran-77
8 C (2nd edition) chapter 3, section on cubic splines.
9 C
10 C input:
11 !============
12 ! idecl: explicitize implicit INTEGER declarations:
13  IMPLICIT NONE
14  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
15  INTEGER i,inwk
16 !============
17 ! idecl: explicitize implicit REAL declarations:
18  real*8 bcxmax
19 !============
20  integer nx ! no. of data points
21  real*8 x(nx) ! x axis data, strict ascending order
22 C
23 C input/output:
24  real*8 fspl(2,nx) ! f(1,*): data in; f(2,*): coeffs out
25 C
26 C f(1,j) = f(x(j)) on input (unchanged on output)
27 C f(2,j) = f''(x(j)) (of interpolating spline) (on output).
28 C
29 C ...boundary conditions...
30 C
31 C input:
32 C
33  integer ibcxmin ! b.c. flag @ x=xmin=x(1)
34  real*8 bcxmin ! b.c. data @xmin
35 C
36  integer ibcxmax ! b.c. flag @ x=xmax=x(nx)
37  real*8 bcxcmax ! b.c. data @xmax
38 C
39 C ibcxmin=-1 -- periodic boundary condition
40 C (bcxmin,ibcxmax,bcxmax are ignored)
41 C
42 C the output spline s satisfies
43 C s'(x(1))=s'(x(nx)) ..and.. s''(x(1))=s''(x(nx))
44 C
45 C if non-periodic boundary conditions are used, then the xmin and xmax
46 C boundary conditions can be specified independently:
47 C
48 C ibcxmin (ibcxmax) = 0 -- this specifies a "not a knot" boundary
49 C condition, see "cubsplb.for". This is a common way
50 C for inferring a "good" spline boundary condition
51 C automatically from data in the vicinity of the
52 C boundary. ... bcxmin (bcxmax) are ignored.
53 C
54 C ibcxmin (ibcxmax) = 1 -- boundary condition is to have s'(x(1))
55 C ( s'(x(nx)) ) match the passed value bcxmin (bcxmax).
56 C
57 C ibcxmin (ibcxmax) = 2 -- boundary condition is to have s''(x(1))
58 C ( s''(x(nx)) ) match the passed value bcxmin (bcxmax).
59 C
60 C ibcxmin (ibcxmax) = 3 -- boundary condition is to have s'(x(1))=0.0
61 C ( s'(x(nx))=0.0 )
62 C
63 C ibcxmin (ibcxmax) = 4 -- boundary condition is to have s''(x(1))=0.0
64 C ( s''(x(nx))=0.0 )
65 C
66 C ibcxmin (ibcxmax) = 5 -- boundary condition is to have s'(x(1))
67 C ( s'(x(nx)) ) match the 1st divided difference
68 C e.g. at x(1): d(1)/h(1), where
69 C d(j)=f(1,j+1)-f(1,j)
70 C h(j)=x(j+1)-x(j)
71 C
72 C ibcxmin (ibcxmax) = 6 -- BC is to have s''(x(1)) ( s''(x(nx)) )
73 C match the 2nd divided difference
74 C e.g. at x(1):
75 C e(1) = [d(2)/h(2) - d(1)/h(1)]/(0.5*(h(1)+h(2)))
76 C
77 C ibcxmin (ibcxmax) = 7 -- BC is to have s'''(x(1)) ( s'''(x(nx)) )
78 C match the 3rd divided difference
79 C e.g. at x(1): [e(2)-e(1)]/(0.33333*(h(1)+h(2)+h(3)))
80 C
81 C output:
82 C
83  integer ilinx ! =1: hint, x axis is ~evenly spaced
84 C
85 C let dx[avg] = (x(nx)-x(1))/(nx-1)
86 C let dx[j] = x(j+1)-x(j), for all j satisfying 1.le.j.lt.nx
87 C
88 C if for all such j, abs(dx[j]-dx[avg]).le.(1.0e-3*dx[avg]) then
89 C ilinx=1 is returned, indicating the data is (at least nearly)
90 C evenly spaced. Even spacing is useful, for speed of zone lookup,
91 C when evaluating a spline.
92 C
93 C if the even spacing condition is not satisfied, ilinx=2 is returned.
94 C
95  integer ier ! exit code, 0=OK
96 C
97 C an error code is returned if the x axis is not strict ascending,
98 C or if nx.lt.4, or if an invalid boundary condition specification was
99 C input.
100 C
101 C------------------------------------
102 C
103 C this routine calls traditional 4 coefficient spline software, and
104 C translates the result to 2 coefficient form.
105 C
106 C this could be done more efficiently but we decided out of conservatism
107 C to use the traditional software.
108 C
109 C------------------------------------
110 C workspaces -- f90 dynamically allocated
111 C
112  real*8, dimension(:,:), allocatable :: fspl4 ! traditional 4-spline
113  real*8, dimension(:), allocatable :: wk ! cspline workspace
114 C
115 C------------------------------------
116 C
117  allocate(fspl4(4,nx),wk(nx))
118 C
119 C make the traditional call
120 C
121  do i=1,nx
122  fspl4(1,i)=fspl(1,i)
123  fspl(2,i)=0.0_r8 ! for now
124  enddo
125 C
126  inwk=nx
127 C
128 C boundary conditions imposed by cspline...
129 C
130  call r8cspline(x,nx,fspl4,ibcxmin,bcxmin,ibcxmax,bcxmax,
131  > wk,inwk,ilinx,ier)
132 C
133  if(ier.eq.0) then
134 C
135 C copy the output -- careful of end point.
136 C
137  do i=1,nx-1
138  fspl(2,i)=2.0_r8*fspl4(3,i)
139  enddo
140  fspl(2,nx)=2.0_r8*fspl4(3,nx-1) +
141  > (x(nx)-x(nx-1))*6.0_r8*fspl4(4,nx-1)
142  endif
143 C
144  deallocate(fspl4,wk)
145 C
146  return
147  end
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17