V3FIT
cspline.f
1 c cspline -- dmc 15 Feb 1999
2 c
3 c a standard interface to the 1d spline setup routine
4 c modified dmc 3 Mar 2000 -- to use Wayne Houlberg's v_spline code.
5 c new BC options added.
6 c
7  subroutine cspline(x,nx,fspl,ibcxmin,bcxmin,ibcxmax,bcxmax,
8  > wk,iwk,ilinx,ier)
9 c
10  real x(nx) ! x axis (in)
11  real fspl(4,nx) ! spline data (in/out)
12  integer ibcxmin ! x(1) BC flag (in, see comments)
13  real bcxmin ! x(1) BC data (in, see comments)
14  integer ibcxmax ! x(nx) BC flag (in, see comments)
15  real bcxmax ! x(nx) BC data (in, see comments)
16  real wk(iwk) ! workspace of size at least nx
17  integer ilinx ! even spacing flag (out)
18  integer ier ! output, =0 means OK
19 c
20 c ** note wk(...) array is not used unless ibcxmin=-1 (periodic spline
21 c evaluation)
22 c
23 c this routine computes spline coefficients for a 1d spline --
24 c evaluation of the spline can be done by cspeval.for subroutines
25 c or directly by inline code.
26 c
27 c the input x axis x(1...nx) must be strictly ascending, i.e.
28 c x(i+1).gt.x(i) is required for i=1 to nx-1. This is checked and
29 c ier=1 is set and the routine exits if the test is not satisfied.
30 c
31 c on output, ilinx=1 is set if, to a reasonably close tolerance,
32 c all grid spacings x(i+1)-x(i) are equal. This allows a speedier
33 c grid lookup algorithm on evaluation of the spline. If on output
34 c ilinx=2, this means the spline x axis is not evenly spaced.
35 c
36 c the input data for the spline are given in f[j] = fspl(1,j). The
37 c output data are the spline coefficients fspl(2,j),fspl(3,j), and
38 c fspl(4,j), j=1 to nx. The result is a spline s(x) satisfying the
39 c boundary conditions and with the properties
40 c
41 c s(x(j)) = fspl(1,j)
42 c s'(x) is continuous even at the grid points x(j)
43 c s''(x) is continuous even at the grid points x(j)
44 c
45 c the formula for evaluation of s(x) is:
46 c
47 c let dx = x-x(i), where x(i).le.x.le.x(i+1). Then,
48 c s(x)=fspl(1,i) + dx*(fspl(2,i) +dx*(fspl(3,i) + dx*fspl(4,i)))
49 c
50 c ==>boundary conditions. Complete specification of a 1d spline
51 c requires specification of boundary conditions at x(1) and x(nx).
52 c
53 c this routine provides 4 options:
54 c
55 c -1 ***** PERIODIC BC
56 c ibcxmin=-1 -- periodic boundary condition. This means the
57 c boundary conditions s'(x(1))=s'(x(nx)) and s''(x(1))=s''(x(nx))
58 c are imposed. Note that s(x(1))=s(x(nx)) (i.e. fspl(1,1)=fspl(1,nx))
59 c is not required -- that is determined by the fspl array input data.
60 c The periodic boundary condition is to be preferred for periodic
61 c data. When splining periodic data f(x) with period P, the relation
62 c x(nx)=x(1)+n*P, n = the number of periods (usually 1), should hold.
63 c (ibcxmax, bcxmin, bcxmax are ignored).
64 c
65 c if a periodic boundary condition is set, this covers both boundaries.
66 c for the other types of boundary conditions, the type of condition
67 c chosen for the x(1) boundary need not be the same as the type chosen
68 c for the x(nx) boundary.
69 c
70 c 0 ***** NOT A KNOT BC
71 c ibcxmin=0 | ibcxmax=0 -- this specifies a "not a knot" boundary
72 c condition -- see cubsplb.for. This is a common way for inferring
73 c a "good" spline boundary condition automatically from data in the
74 c vicinity of the boundary. (bcxmin | bcxmax are ignored).
75 c
76 c 1 ***** BC: SPECIFIED SLOPE
77 c ibcxmin=1 | ibcxmax=1 -- boundary condition is to have s'(x(1)) |
78 c s'(x(nx)) match the passed value (bcxmin | bcxmax).
79 c
80 c 2 ***** BC: SPECIFIED 2nd DERIVATIVE
81 c ibcxmin=2 | ibcxmax=2 -- boundary condition is to have s''(x(1)) |
82 c s''(x(nx)) match the passed value (bcxmin | bcxmax).
83 c
84 c 3 ***** BC: SPECIFIED SLOPE = 0.0
85 c ibcxmin=3 | ibcxmax=3 -- boundary condition is to have s'(x(1)) |
86 c s'(x(nx)) equal to ZERO.
87 c
88 c 4 ***** BC: SPECIFIED 2nd DERIVATIVE = 0.0
89 c ibcxmin=4 | ibcxmax=4 -- boundary condition is to have s''(x(1)) |
90 c s''(x(nx)) equal to ZERO.
91 c
92 c 5 ***** BC: 1st DIVIDED DIFFERENCE
93 c ibcxmin=5 | ibcxmax=5 -- boundary condition is to have s'(x(1)) |
94 c s'(x(nx)) equal to the slope from the 1st|last 2 points
95 c
96 c 6 ***** BC: 2nd DIVIDED DIFFERENCE
97 c ibcxmin=6 | ibcxmax=6 -- boundary condition is to have s''(x(1)) |
98 c s''(x(nx)) equal to the 2nd derivative from the 1st|last 3 points
99 c
100 c 7 ***** BC: 3rd DIVIDED DIFFERENCE
101 c ibcxmin=7 | ibcxmax=7 -- boundary condition is to have s'''(x(1)) |
102 c s'''(x(nx)) equal to the 3rd derivative from the 1st|last 4 points
103 c
104 c---------------------------------------------------------------------
105  data half/0.5/
106  data sixth/0.166666666666666667/
107 c
108 c error checks
109 c
110  ier = 0
111  if(nx.lt.2) then
112  write(6,'('' ?cspline: at least 2 x points required.'')')
113  ier=1
114  endif
115  call ibc_ck(ibcxmin,'cspline','xmin',-1,7,ier)
116  if(ibcxmin.ge.0) call ibc_ck(ibcxmax,'cspline','xmax',0,7,ier)
117 c
118 c x axis check
119 c
120  call splinck(x,nx,ilinx,1.0e-3,ierx)
121  if(ierx.ne.0) ier=2
122 c
123  if(ier.eq.2) then
124  write(6,'('' ?cspline: x axis not strict ascending'')')
125  endif
126 c
127  if(ibcxmin.eq.-1) then
128  inum=nx
129  if(iwk.lt.inum) then
130  write(6,1009) inum,iwk,nx
131  1009 format(
132  > ' ?cspline: workspace too small. need: ',i6,' got: ',i6/
133  > ' (need = nx, nx=',i6)
134  ier=3
135  endif
136  endif
137 c
138  if(ier.ne.0) return
139 c
140 c OK -- evaluate spline
141 c
142  if(ibcxmin.eq.1) then
143  fspl(2,1)=bcxmin
144  else if(ibcxmin.eq.2) then
145  fspl(3,1)=bcxmin
146  endif
147 c
148  if(ibcxmax.eq.1) then
149  fspl(2,nx)=bcxmax
150  else if(ibcxmax.eq.2) then
151  fspl(3,nx)=bcxmax
152  endif
153 c
154  call v_spline(ibcxmin,ibcxmax,nx,x,fspl,wk)
155 c
156  do i=1,nx
157  fspl(3,i)=half*fspl(3,i)
158  fspl(4,i)=sixth*fspl(4,i)
159  enddo
160 c
161  return
162  end
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17