V3FIT
v_spline.f
1  SUBROUTINE v_spline(k_bc1,k_bcn,n,x,f,wk)
2 !***********************************************************************
3 !V_SPLINE evaluates the coefficients for a 1d cubic interpolating spline
4 !References:
5 ! Forsythe, Malcolm, Moler, Computer Methods for Mathematical
6 ! Computations, Prentice-Hall, 1977, p.76
7 ! Engeln-Muellges, Uhlig, Numerical Algorithms with Fortran, Springer,
8 ! 1996, p.251
9 ! W.A.Houlberg, D.McCune 3/2000
10 !Input:
11 ! k_bc1-option for BC at x1 = x(1)
12 ! =-1 periodic, ignore k_bcn
13 ! =0 not-a-knot
14 ! =1 s'(x1) = input value of f(2,1)
15 ! =2 s''(x1) = input value of f(3,1)
16 ! =3 s'(x1) = 0.0
17 ! =4 s''(x1) = 0.0
18 ! =5 match first derivative to first 2 points
19 ! =6 match second derivative to first 3 points
20 ! =7 match third derivative to first 4 points
21 ! =else use not-a-knot
22 ! k_bcn-option for boundary condition at xn = x(n)
23 ! =0 not-a-knot
24 ! =1 s'(xn) = input value of f(2,n)
25 ! =2 s''(xn) = input value of f(3,n)
26 ! =3 s'(xn) = 0.0
27 ! =4 s''(xn) = 0.0
28 ! =5 match first derivative to last 2 points
29 ! =6 match second derivative to lasst 3 points
30 ! =7 match third derivative to last 4 points
31 ! =else use knot-a-knot
32 ! n-number of data points or knots-(n.ge.2)
33 ! x(n)-abscissas of the knots in strictly increasing order
34 ! f(1,i)-ordinates of the knots
35 ! f(2,1)-input value of s'(x1) for k_bc1=1
36 ! f(2,n)-input value of s'(xn) for k_bcn=1
37 ! f(3,1)-input value of s''(x1) for k_bc1=2
38 ! f(3,n)-input value of s''(xn) for k_bcn=2
39 ! wk(n)-scratch work area for periodic BC
40 !Output:
41 ! f(2,i)=s'(x(i))
42 ! f(3,i)=s''(x(i))
43 ! f(4,i)=s'''(x(i))
44 !Comments:
45 ! s(x)=f(1,i)+f(2,i)*(x-x(i))+f(3,i)*(x-x(i))**2/2!
46 ! +f(4,i)*(x-x(i))**3/3! for x(i).le.x.le.x(i+1)
47 ! W_SPLINE can be used to evaluate the spline and its derivatives
48 ! The cubic spline is twice differentiable (C2)
49 !
50 ! modifications -- dmc 18 Feb 2010:
51 ! Deal with n.lt.4 -- the general tridiagonal spline method
52 ! does not have the right formulation for n.eq.3 "not a knot" or periodic
53 ! boundary conditions, nor for n.eq.2 with any boundary conditions.
54 !
55 ! Apply boundary conditions even for n=2, when the "spline" is really
56 ! just a single cubic polynomial.
57 ! In this case, several boundary condition (BC) options are mapped to
58 ! BC option 5, 1st divided difference. If 5 is used on both sides
59 ! of an n=2 "spline" you get a linear piece which is what the old
60 ! code always gave, regardless of BC option settings. The following
61 ! BC controls are mapped to 5:
62 ! periodic (-1)
63 ! not a knot (0) (for n=2 no grid point exists for knot location).
64 ! option (5) is preserved
65 ! options 6 and 7 -- mapped to (5); higher divided differences
66 ! need n>2; in fact 7 needs n>3; for n=3 option 6 is substituted.
67 !
68 ! The boundary condition screening is done at the start of the code;
69 ! passed controls k_bc1 and k_bcn are mapped to i_bc1 and i_bcn.
70 !
71 ! ALSO: for n=3, "not a knot" from both left and right needs special
72 ! interpretation, since the 2 boundary conditions overlap. The
73 ! chosen interpretation is a parabolic fit to the 3 given data points.
74 ! and so f''' = 0 and f'' = constant. If "not a knot" is used on
75 ! one side only, the solution is a single cubic piece and special
76 ! code is also needed.
77 ! ALSO: for n=3, "periodic" boundary condition needs special code; this
78 ! is added.
79 !
80 ! bugfixes -- dmc 24 Feb 2004:
81 ! (a) fixed logic for not-a-knot:
82 ! ! Set f(3,1) for not-a-knot
83 ! IF(k_bc1.le.0.or.k_bc1.gt.7) THEN ...
84 ! instead of
85 ! ! Set f(3,1) for not-a-knot
86 ! IF(k_bc1.le.0.or.k_bc1.gt.5) THEN ...
87 ! and similarly for logic after cmt
88 ! ! Set f(3,n) for not-a-knot
89 ! as required since k_bc*=6 and k_bc*=7 are NOT not-a-knot BCs.
90 !
91 ! (b) the BCs to fix 2nd derivative at end points did not work if that
92 ! 2nd derivative were non-zero. The reason is that in those cases
93 ! the off-diagonal matrix elements nearest the corners are not
94 ! symmetric; i.e. elem(1,2).ne.elem(2,1) and
95 ! elem(n-1,n).ne.elem(n,n-1) where I use "elem" to refer to
96 ! the tridiagonal matrix elements. The correct values for the
97 ! elements is: elem(1,2)=0, elem(2,1)=x(2)-x(1)
98 ! elem(n,n-1)=0, elem(n-1,n)=x(n)-x(n-1)
99 ! the old code in effect had these as all zeroes. Since this
100 ! meant the wrong set of linear equations was solved, the
101 ! resulting spline had a discontinuity in its 1st derivative
102 ! at x(2) and x(n-1). Fixed by introducing elem21 and elemnn1
103 ! to represent the non-symmetric lower-diagonal values. Since
104 ! elem21 & elemnn1 are both on the lower diagonals, logic to
105 ! use them occurs in the non-periodic forward elimination loop
106 ! only. DMC 24 Feb 2004.
107 !***********************************************************************
108  IMPLICIT NONE
109 !Declaration of input variables
110  INTEGER k_bc1, k_bcn,
111  & n
112  REAL x(*), wk(*),
113  & f(4,*)
114 !Declaration in local variables
115  INTEGER i, ib,
116  & imax, imin
117  REAL a1, an,
118  & b1, bn,
119  & q, t,
120  & hn
121  REAL elem21, elemnn1 ! (dmc)
122 
123  integer :: i_bc1,i_bcn ! screened BC controls
124 
125  integer :: iord1,iord2 ! used for n=2 only
126  real :: h,f0,fh ! used for n=2,3 only
127  real :: h1,h2,h3,dels ! used for n=3 special cases
128  real :: f1,f2,f3,aa,bb ! used for n=3
129 
130  integer :: i3knots ! for n=3, number of not-a-knot BCs
131  integer :: i3perio ! for n=3, periodic BC
132 
133 !------------------------------------------------------------
134 ! screen the BC options (DMC Feb. 2010...)
135 
136  i_bc1=k_bc1
137  i_bcn=k_bcn
138 
139  if((i_bc1.lt.-1).or.(i_bc1.gt.7)) i_bc1=0 ! outside [-1:7] -> not-a-knot
140  if((i_bcn.lt.0).or.(i_bcn.gt.7)) i_bcn=0 ! outside [0:7] -> not-a-knot
141 
142  if(i_bc1.eq.-1) i_bcn=-1 ! periodic BC
143 
144  i3knots=0
145  i3perio=0
146  if(n.eq.3) then
147  i_bc1=min(6,i_bc1)
148  i_bcn=min(6,i_bcn)
149  if(i_bc1.eq.0) i3knots = i3knots + 1
150  if(i_bcn.eq.0) i3knots = i3knots + 1
151  if(i_bc1.eq.-1) i3perio = 1
152  endif
153 
154  if(n.eq.2) then
155  if(i_bc1.eq.-1) then
156  i_bc1=5
157  i_bcn=5
158  endif
159  if((i_bc1.eq.0).or.(i_bc1.gt.5)) i_bc1=5
160  if((i_bcn.eq.0).or.(i_bcn.gt.5)) i_bcn=5
161 
162  if((i_bc1.eq.1).or.(i_bc1.eq.3).or.(i_bc1.eq.5)) then
163  iord1=1 ! 1st derivative match on LHS
164  else
165  iord1=2 ! 2nd derivative match on LHS
166  endif
167 
168  if((i_bcn.eq.1).or.(i_bcn.eq.3).or.(i_bcn.eq.5)) then
169  iord2=1 ! 1st derivative match on RHS
170  else
171  iord2=2 ! 2nd derivative match on RHS
172  endif
173  endif
174 
175 !Set default range
176  imin=1
177  imax=n
178 !Set first and second BC values
179  a1=0.0
180  b1=0.0
181  an=0.0
182  bn=0.0
183  IF(i_bc1.eq.1) THEN
184  a1=f(2,1)
185  ELSEIF(i_bc1.eq.2) THEN
186  b1=f(3,1)
187  ELSEIF(i_bc1.eq.5) THEN
188  a1=(f(1,2)-f(1,1))/(x(2)-x(1))
189  ELSEIF(i_bc1.eq.6) THEN
190  b1=2.0*((f(1,3)-f(1,2))/(x(3)-x(2))
191  & -(f(1,2)-f(1,1))/(x(2)-x(1)))/(x(3)-x(1))
192  ENDIF
193  IF(i_bcn.eq.1) THEN
194  an=f(2,n)
195  ELSEIF(i_bcn.eq.2) THEN
196  bn=f(3,n)
197  ELSEIF(i_bcn.eq.5) THEN
198  an=(f(1,n)-f(1,n-1))/(x(n)-x(n-1))
199  ELSEIF(i_bcn.eq.6) THEN
200  bn=2.0*((f(1,n)-f(1,n-1))/(x(n)-x(n-1))
201  & -(f(1,n-1)-f(1,n-2))/(x(n-1)-x(n-2)))/(x(n)-x(n-2))
202  ENDIF
203 !Clear f(2:4,n)
204  f(2,n)=0.0
205  f(3,n)=0.0
206  f(4,n)=0.0
207  IF(n.eq.2) THEN
208  if((i_bc1.eq.5).and.(i_bcn.eq.5)) then
209 !Coefficients for n=2 (this was the original code)
210  f(2,1)=(f(1,2)-f(1,1))/(x(2)-x(1))
211  f(3,1)=0.0
212  f(4,1)=0.0
213  f(2,2)=f(2,1)
214  f(3,2)=0.0
215  f(4,2)=0.0
216  else if((iord1.eq.1).and.(iord2.eq.1)) then
217  ! LHS: match a1 for 1st deriv; RHS: match an for 1st deriv.
218  f(2,1)=a1
219  f(2,2)=an
220  h = (x(2)-x(1))
221  f0 = f(1,1)
222  fh = f(1,2)
223 
224  ! setting xx = x-x(1),
225  ! f = c1*xx**3 + c2*xx**2 + a1*xx + f0
226  ! --> c1*h**3 + c2*h**2 = fh - f0 - a1*h
227  ! and 3*c1*h**2 + 2*c2*h = an - a1
228  ! f' = 3*c1*xx*2 + 2*c2*xx + a1
229  ! f'' = 6*c1*xx + 2*c2
230  ! f''' = 6*c1
231 
232  ! solve 2x2 system for c1 -> f(4,1)/6 and c2 -> f(3,1)/2
233 
234  f(3,1) = (3*(fh-f0)/(h*h) - (2*a1 + an)/h)*2 ! 2*c2
235  f(4,1) = (-2*(fh-f0)/(h*h*h) + (a1 + an)/(h*h))*6 ! 6*c1
236 
237  f(4,2) = f(4,1)
238  f(3,2) = f(4,1)*h + f(3,1)
239 
240  else if((iord1.eq.1).and.(iord2.eq.2)) then
241  ! LHS: match a1 for 1st deriv; RHS: match bn for 2nd deriv.
242  f(2,1)=a1
243  f(3,2)=bn
244  h = (x(2)-x(1))
245  f0 = f(1,1)
246  fh = f(1,2)
247 
248  ! setting xx = x-x(1),
249  ! f = c1*xx**3 + c2*xx**2 + a1*xx + f0
250  ! --> c1*h**3 + c2*h**2 = fh - f0 - a1*h
251  ! and 6*c1*h + 2*c2 = bn
252  ! f' = 3*c1*xx*2 + 2*c2*xx + a1
253  ! f'' = 6*c1*xx + 2*c2
254  ! f''' = 6*c1
255 
256  ! solve 2x2 system for c1 -> f(4,1)/6 and c2 -> f(3,1)/2
257 
258  f(3,1) = (-bn/4 + 3*(fh-f0)/(2*h*h) - 3*a1/(2*h))*2 ! 2*c2
259  f(4,1) = (bn/(4*h) - (fh-f0)/(2*h*h*h) + a1/(2*h*h))*6 ! 6*c1
260 
261  f(4,2) = f(4,1)
262  f(2,2) = f(4,1)*h*h/2 + f(3,1)*h + a1
263  else if((iord1.eq.2).and.(iord2.eq.1)) then
264  ! LHS: match b1 for 2nd deriv; RHS: match an for 1st deriv.
265  f(3,1)=b1
266  f(2,2)=an
267  h = (x(2)-x(1))
268  f0 = f(1,1)
269  fh = f(1,2)
270 
271  ! setting xx = x-x(1),
272  ! f = c1*xx**3 + (b1/2)*xx**2 + c3*xx + f0
273  ! --> c1*h**3 + c3*h = fh - f0 - b1*h**2/2
274  ! and 3*c1*h**2 + c3 = an - b1*h
275  ! f' = 3*c1*xx*2 + b1*xx + c3
276  ! f'' = 6*c1*xx + b1
277  ! f''' = 6*c1
278 
279  ! solve 2x2 system for c1 -> f(4,1)/6 and c2 -> f(3,1)/2
280 
281  f(2,1) = 3*(fh-f0)/(2*h) - b1*h/4 - an/2 ! c3
282  f(4,1) = (an/(2*h*h) - (fh-f0)/(2*h*h*h) - b1/(4*h))*6 ! 6*c1
283 
284  f(4,2) = f(4,1)
285  f(3,2) = f(4,1)*h + f(3,1)
286  else if((iord1.eq.2).and.(iord2.eq.2)) then
287  ! LHS: match b1 for 2nd deriv; RHS: match bn for 2nd deriv.
288  f(3,1)=b1
289  f(3,2)=bn
290  h = (x(2)-x(1))
291  f0 = f(1,1)
292  fh = f(1,2)
293 
294  ! setting xx = x-x(1),
295  ! f = c1*xx**3 + (b1/2)*xx**2 + c3*xx + f0
296  ! --> c1*h**3 + c3*h = fh - f0 - b1*h**2/2
297  ! and 6*c1*h = bn - b1
298  ! f' = 3*c1*xx*2 + b1*xx + c3
299  ! f'' = 6*c1*xx + b1
300  ! f''' = 6*c1
301 
302  ! solve 2x2 system for c1 -> f(4,1)/6 and c2 -> f(3,1)/2
303 
304  f(2,1) = (fh-f0)/h -b1*h/3 -bn*h/6 ! c3
305  f(4,1) = (bn-b1)/h ! 6*c1
306 
307  f(4,2) = f(4,1)
308  f(2,2) = f(4,1)*h*h/2 + b1*h + f(2,1)
309  endif
310 
311  ELSE IF(i3perio.eq.1) then
312 !Special case: nx=3 periodic spline
313  h1=x(2)-x(1)
314  h2=x(3)-x(2)
315  h=h1+h2
316 
317  dels=(f(1,3)-f(1,2))/h2 - (f(1,2)-f(1,1))/h1
318 
319  f(2,1)= (f(1,2)-f(1,1))/h1 + (h1*dels)/h
320  f(3,1)= -6*dels/h
321  f(4,1)= 12*dels/(h1*h)
322 
323  f(2,2)= (f(1,3)-f(1,2))/h2 - (h2*dels)/h
324  f(3,2)= 6*dels/h
325  f(4,2)= -12*dels/(h2*h)
326 
327  f(2,3)=f(2,1)
328  f(3,3)=f(3,1)
329  f(4,3)=f(4,1)
330 
331 
332  ELSE IF(i3knots.eq.2) then
333 !Special case: nx=3, not-a-knot on both sides
334  h1=x(2)-x(1)
335  h2=x(3)-x(2)
336  h=h1+h2
337  ! this is just a quadratic fit through 3 pts
338  f1=f(1,1)-f(1,2)
339  f2=f(1,3)-f(1,2)
340 
341 ! quadratic around origin at (x(2),f(1,2))
342 ! aa*h1**2 - bb*h1 = f1
343 ! aa*h2**2 + bb*h2 = f2
344 
345  aa = (f2*h1 + f1*h2)/(h1*h2*h)
346  bb = (f2*h1*h1 - f1*h2*h2)/(h1*h2*h)
347 
348  f(4,1:3)=0.0 ! f''' = 0 (quadratic polynomial)
349  f(3,1:3)=2*aa ! f'' = const
350 
351  f(2,1)=bb-2*aa*h1
352  f(2,2)=bb
353  f(2,3)=bb+2*aa*h2
354 
355  ELSE IF(i3knots.eq.1) then
356 !Special cases: nx=3, not-a-knot on single side
357  if((i_bc1.eq.1).or.(i_bc1.eq.3).or.(i_bc1.eq.5)) then
358  ! f' LHS condition; not-a-knot RHS
359 ! a1 = f' LHS BC
360  h2=x(2)-x(1)
361  h3=x(3)-x(1)
362 
363  f2=f(1,2)-f(1,1)
364  f3=f(1,3)-f(1,1)
365 
366 ! find cubic aa*xx**3 + bb*xx**2 + a1*xx
367 ! satisfying aa*h2**3 + bb*h2**2 + a1*h2 = f2
368 ! and aa*h3**3 + bb*h3**2 + a1*h3 = f3
369 
370  aa=a1/(h2*h3) + f3/(h3*h3*(h3-h2)) - f2/(h2*h2*(h3-h2))
371  bb=-a1*(h3*h3-h2*h2)/(h2*h3*(h3-h2))
372  > + f2*h3/(h2*h2*(h3-h2)) - f3*h2/(h3*h3*(h3-h2))
373 
374  f(2,1)=a1
375  f(3,1)=2*bb
376  f(4,1)=6*aa
377 
378  f(2,2)=3*aa*h2*h2 + 2*bb*h2 + a1
379  f(3,2)=6*aa*h2 + 2*bb
380  f(4,2)=6*aa
381 
382  f(2,3)=3*aa*h3*h3 + 2*bb*h3 + a1
383  f(3,3)=6*aa*h3 + 2*bb
384  f(4,3)=6*aa
385 
386  else if((i_bc1.eq.2).or.(i_bc1.eq.4).or.(i_bc1.eq.6)) then
387  ! f'' LHS condition; not-a-knot RHS
388 ! b1 = f'' LHS BC
389  h2=x(2)-x(1)
390  h3=x(3)-x(1)
391 
392  f2=f(1,2)-f(1,1)
393  f3=f(1,3)-f(1,1)
394 
395 ! find cubic aa*xx**3 + (b1/2)*xx**2 + bb*xx
396 ! satisfying aa*h2**3 + bb*h2 = f2 -(b1/2)*h2**2
397 ! and aa*h3**3 + bb*h3 = f3 -(b1/2)*h3**2
398 
399  aa= -(b1/2)*(h3-h2)/(h3*h3-h2*h2)
400  > -f2/(h2*(h3*h3-h2*h2)) + f3/(h3*(h3*h3-h2*h2))
401  bb= -(b1/2)*h2*h3*(h3-h2)/(h3*h3-h2*h2)
402  > +f2*h3*h3/(h2*(h3*h3-h2*h2))
403  > -f3*h2*h2/(h3*(h3*h3-h2*h2))
404 
405  f(2,1)=bb
406  f(3,1)=b1
407  f(4,1)=6*aa
408 
409  f(2,2)=3*aa*h2*h2 + b1*h2 + bb
410  f(3,2)=6*aa*h2 + b1
411  f(4,2)=6*aa
412 
413  f(2,3)=3*aa*h3*h3 + b1*h3 + bb
414  f(3,3)=6*aa*h3 + b1
415  f(4,3)=6*aa
416 
417  else if((i_bcn.eq.1).or.(i_bcn.eq.3).or.(i_bcn.eq.5)) then
418  ! f' RHS condition; not-a-knot LHS
419 ! an = f' RHS BC
420  h2=x(2)-x(3)
421  h3=x(1)-x(3)
422 
423  f2=f(1,2)-f(1,3)
424  f3=f(1,1)-f(1,3)
425 
426 ! find cubic aa*xx**3 + bb*xx**2 + an*xx
427 ! satisfying aa*h2**3 + bb*h2**2 + an*h2 = f2
428 ! and aa*h3**3 + bb*h3**2 + an*h3 = f3
429 
430  aa=an/(h2*h3) + f3/(h3*h3*(h3-h2)) - f2/(h2*h2*(h3-h2))
431  bb=-an*(h3*h3-h2*h2)/(h2*h3*(h3-h2))
432  > + f2*h3/(h2*h2*(h3-h2)) - f3*h2/(h3*h3*(h3-h2))
433 
434  f(2,3)=an
435  f(3,3)=2*bb
436  f(4,3)=6*aa
437 
438  f(2,2)=3*aa*h2*h2 + 2*bb*h2 + an
439  f(3,2)=6*aa*h2 + 2*bb
440  f(4,2)=6*aa
441 
442  f(2,1)=3*aa*h3*h3 + 2*bb*h3 + an
443  f(3,1)=6*aa*h3 + 2*bb
444  f(4,1)=6*aa
445 
446  else if((i_bcn.eq.2).or.(i_bcn.eq.4).or.(i_bcn.eq.6)) then
447  ! f'' RHS condition; not-a-knot LHS
448 ! bn = f'' RHS BC
449  h2=x(2)-x(3)
450  h3=x(1)-x(3)
451 
452  f2=f(1,2)-f(1,3)
453  f3=f(1,1)-f(1,3)
454 
455 ! find cubic aa*xx**3 + (bn/2)*xx**2 + bb*xx
456 ! satisfying aa*h2**3 + bb*h2 = f2 -(bn/2)*h2**2
457 ! and aa*h3**3 + bb*h3 = f3 -(bn/2)*h3**2
458 
459  aa= -(bn/2)*(h3-h2)/(h3*h3-h2*h2)
460  > -f2/(h2*(h3*h3-h2*h2)) + f3/(h3*(h3*h3-h2*h2))
461  bb= -(bn/2)*h2*h3*(h3-h2)/(h3*h3-h2*h2)
462  > +f2*h3*h3/(h2*(h3*h3-h2*h2))
463  > -f3*h2*h2/(h3*(h3*h3-h2*h2))
464 
465  f(2,3)=bb
466  f(3,3)=bn
467  f(4,3)=6*aa
468 
469  f(2,2)=3*aa*h2*h2 + bn*h2 + bb
470  f(3,2)=6*aa*h2 + bn
471  f(4,2)=6*aa
472 
473  f(2,1)=3*aa*h3*h3 + bn*h3 + bb
474  f(3,1)=6*aa*h3 + bn
475  f(4,1)=6*aa
476 
477  endif
478  ELSE IF(n.gt.2) THEN
479 !Set up tridiagonal system for A*y=B where y(i) are the second
480 ! derivatives at the knots
481 ! f(2,i) are the diagonal elements of A
482 ! f(4,i) are the off-diagonal elements of A
483 ! f(3,i) are the B elements/3, and will become c/3 upon solution
484  f(4,1)=x(2)-x(1)
485  f(3,2)=(f(1,2)-f(1,1))/f(4,1)
486  DO i=2,n-1
487  f(4,i)=x(i+1)-x(i)
488  f(2,i)=2.0*(f(4,i-1)+f(4,i))
489  f(3,i+1)=(f(1,i+1)-f(1,i))/f(4,i)
490  f(3,i)=f(3,i+1)-f(3,i)
491  ENDDO
492 !
493 ! (dmc): save now:
494 !
495  elem21=f(4,1)
496  elemnn1=f(4,n-1)
497 !
498 ! BC's
499 ! Left
500  IF(i_bc1.eq.-1) THEN
501  f(2,1)=2.0*(f(4,1)+f(4,n-1))
502  f(3,1)=(f(1,2)-f(1,1))/f(4,1)-(f(1,n)-f(1,n-1))/f(4,n-1)
503  wk(1)=f(4,n-1)
504  DO i=2,n-3
505  wk(i)=0.0
506  ENDDO
507  wk(n-2)=f(4,n-2)
508  wk(n-1)=f(4,n-1)
509  ELSEIF(i_bc1.eq.1.or.i_bc1.eq.3.or.i_bc1.eq.5) THEN
510  f(2,1)=2.0*f(4,1)
511  f(3,1)=(f(1,2)-f(1,1))/f(4,1)-a1
512  ELSEIF(i_bc1.eq.2.or.i_bc1.eq.4.or.i_bc1.eq.6) THEN
513  f(2,1)=2.0*f(4,1)
514  f(3,1)=f(4,1)*b1/3.0
515  f(4,1)=0.0 ! upper diagonal only (dmc: cf elem21)
516  ELSEIF(i_bc1.eq.7) THEN
517  f(2,1)=-f(4,1)
518  f(3,1)=f(3,3)/(x(4)-x(2))-f(3,2)/(x(3)-x(1))
519  f(3,1)=f(3,1)*f(4,1)**2/(x(4)-x(1))
520  ELSE ! not a knot:
521  imin=2
522  f(2,2)=f(4,1)+2.0*f(4,2)
523  f(3,2)=f(3,2)*f(4,2)/(f(4,1)+f(4,2))
524  ENDIF
525 ! Right
526  IF(i_bcn.eq.1.or.i_bcn.eq.3.or.i_bcn.eq.5) THEN
527  f(2,n)=2.0*f(4,n-1)
528  f(3,n)=-(f(1,n)-f(1,n-1))/f(4,n-1)+an
529  ELSEIF(i_bcn.eq.2.or.i_bcn.eq.4.or.i_bcn.eq.6) THEN
530  f(2,n)=2.0*f(4,n-1)
531  f(3,n)=f(4,n-1)*bn/3.0
532 !xxx f(4,n-1)=0.0 ! dmc: preserve f(4,n-1) for back subst.
533  elemnn1=0.0 ! lower diaganol only (dmc)
534  ELSEIF(i_bcn.eq.7) THEN
535  f(2,n)=-f(4,n-1)
536  f(3,n)=f(3,n-1)/(x(n)-x(n-2))-f(3,n-2)/(x(n-1)-x(n-3))
537  f(3,n)=-f(3,n)*f(4,n-1)**2/(x(n)-x(n-3))
538  ELSEIF(i_bc1.ne.-1) THEN ! not a knot:
539  imax=n-1
540  f(2,n-1)=2.0*f(4,n-2)+f(4,n-1)
541  f(3,n-1)=f(3,n-1)*f(4,n-2)/(f(4,n-1)+f(4,n-2))
542  ENDIF
543  IF(i_bc1.eq.-1) THEN
544 !Solve system of equations for second derivatives at the knots
545 ! Periodic BC
546 ! Forward elimination
547  DO i=2,n-2
548  t=f(4,i-1)/f(2,i-1)
549  f(2,i)=f(2,i)-t*f(4,i-1)
550  f(3,i)=f(3,i)-t*f(3,i-1)
551  wk(i)=wk(i)-t*wk(i-1)
552  q=wk(n-1)/f(2,i-1)
553  wk(n-1)=-q*f(4,i-1)
554  f(2,n-1)=f(2,n-1)-q*wk(i-1)
555  f(3,n-1)=f(3,n-1)-q*f(3,i-1)
556  ENDDO
557 ! Correct the n-1 element
558  wk(n-1)=wk(n-1)+f(4,n-2)
559 ! Complete the forward elimination
560 ! wk(n-1) and wk(n-2) are the off-diag elements of the lower corner
561  t=wk(n-1)/f(2,n-2)
562  f(2,n-1)=f(2,n-1)-t*wk(n-2)
563  f(3,n-1)=f(3,n-1)-t*f(3,n-2)
564 ! Back substitution
565  f(3,n-1)=f(3,n-1)/f(2,n-1)
566  f(3,n-2)=(f(3,n-2)-wk(n-2)*f(3,n-1))/f(2,n-2)
567  DO ib=3,n-1
568  i=n-ib
569  f(3,i)=(f(3,i)-f(4,i)*f(3,i+1)-wk(i)*f(3,n-1))/f(2,i)
570  ENDDO
571  f(3,n)=f(3,1)
572  ELSE
573 ! Non-periodic BC
574 ! Forward elimination
575 ! For Not-A-Knot BC the off-diagonal end elements are not equal
576  DO i=imin+1,imax
577  IF((i.eq.n-1).and.(imax.eq.n-1)) THEN
578  t=(f(4,i-1)-f(4,i))/f(2,i-1)
579  ELSE
580  if(i.eq.2) then
581  t=elem21/f(2,i-1)
582  else if(i.eq.n) then
583  t=elemnn1/f(2,i-1)
584  else
585  t=f(4,i-1)/f(2,i-1)
586  endif
587  ENDIF
588  IF((i.eq.imin+1).and.(imin.eq.2)) THEN
589  f(2,i)=f(2,i)-t*(f(4,i-1)-f(4,i-2))
590  ELSE
591  f(2,i)=f(2,i)-t*f(4,i-1)
592  ENDIF
593  f(3,i)=f(3,i)-t*f(3,i-1)
594  ENDDO
595 ! Back substitution
596  f(3,imax)=f(3,imax)/f(2,imax)
597  DO ib=1,imax-imin
598  i=imax-ib
599  IF((i.eq.2).and.(imin.eq.2)) THEN
600  f(3,i)=(f(3,i)-(f(4,i)-f(4,i-1))*f(3,i+1))/f(2,i)
601  ELSE
602  f(3,i)=(f(3,i)-f(4,i)*f(3,i+1))/f(2,i)
603  ENDIF
604  ENDDO
605 ! Reset d array to step size
606  f(4,1)=x(2)-x(1)
607  f(4,n-1)=x(n)-x(n-1)
608 ! Set f(3,1) for not-a-knot
609  IF(i_bc1.le.0.or.i_bc1.gt.7) THEN
610  f(3,1)=(f(3,2)*(f(4,1)+f(4,2))-f(3,3)*f(4,1))/f(4,2)
611  ENDIF
612 ! Set f(3,n) for not-a-knot
613  IF(i_bcn.le.0.or.i_bcn.gt.7) THEN
614  f(3,n)=f(3,n-1)+(f(3,n-1)-f(3,n-2))*f(4,n-1)/f(4,n-2)
615  ENDIF
616  ENDIF
617 !f(3,i) is now the sigma(i) of the text and f(4,i) is the step size
618 !Compute polynomial coefficients
619  DO i=1,n-1
620  f(2,i)=
621  > (f(1,i+1)-f(1,i))/f(4,i)-f(4,i)*(f(3,i+1)+2.0*f(3,i))
622  f(4,i)=(f(3,i+1)-f(3,i))/f(4,i)
623  f(3,i)=6.0*f(3,i)
624  f(4,i)=6.0*f(4,i)
625  ENDDO
626  IF(i_bc1.eq.-1) THEN
627  f(2,n)=f(2,1)
628  f(3,n)=f(3,1)
629  f(4,n)=f(4,1)
630  ELSE
631  hn=x(n)-x(n-1)
632  f(2,n)=f(2,n-1)+hn*(f(3,n-1)+0.5*hn*f(4,n-1))
633  f(3,n)=f(3,n-1)+hn*f(4,n-1)
634  f(4,n)=f(4,n-1)
635  IF(i_bcn.eq.1.or.i_bcn.eq.3.or.i_bcn.eq.5) THEN
636  f(2,n)=an
637  ELSE IF(i_bcn.eq.2.or.i_bcn.eq.4.or.i_bcn.eq.6) THEN
638  f(3,n)=bn
639  ENDIF
640  ENDIF
641  ENDIF
642  RETURN
643  END
Generated on Thu Mar 5 2020 15:49:25 for V3FIT by   doxygen 1.8.17