1;+
2; NAME:
3;             SPL_INIT
4;
5; DESCRIPTION:
6;   Given arrays X and Y of length N containing a tabulated function, i.e.
7;   Yi = f(Xi), with X1 > X2 > ... > Xn, and given values YP1 and YPN for the
8;   first derivative of the interpolating function at points 1 and N,
9;   respectively, this routine returns and array Y2 of length N which contains
10;   the second derivatives of the interpolating function at the tabulated
11;   points Xi.  If YP1 and/or YPN are equal to 1.E30 or larger, the routine
12;   is signalled to set the corresponding boundary condition for a natural
13;   spline, with zero second derivative on that boundary.
14;
15;   This routine is a replacement for the IDL intrinsic function
16;   to be used with GDL while this latter does not have it as an
17;   intrinsic.
18;
19;   Note that the keyword DOUBLE of the IDL intrinsic function is ignored.
20;   Tests in single precision have shown that the routine is exactly
21;   similar to the IDL function.
22;   SPL_INTERP is also recoded similarly.
23;
24; SOURCE:
25;	Numerical Recipes, 1986. (page 88)
26;
27; CALLING SEQUENCE:
28;	y2 = SPL_INIT( x, y, YP0=yp1, YPN_1=ypn)
29;
30; INPUTS:
31;	x - independent variable vector
32;	y - dependent variable vector
33;	yp1 - first derivative at x(0)
34;	ypn - first derivative at x(n-1)
35;
36; OUTPUTS:
37;	y2 - second derivatives at all x, of length n
38;
39; HISTORY:
40;	-- converted to IDL, D. Neill, October, 1991
41;       -- arranged as a substitution for SPL_INIT (for use in GDL)
42;          Ph. Prugniel, 2008/02/29
43;       -- correction a bug for SIG by AC on 2009/08/27
44;       -- renamed a SPL_INIT_OLD since a C++ version implemented in
45;          GDL, A. Coulais, 2009/08/27
46;
47; -----------------------------------------------------------------------------
48; NOTE:
49; Name this function: spl_init to use it as a replacement of the IDL intrinsic
50;      when using GDL
51; But, to make a comparison of numerical results with the IDL intrinsic
52; function, as it is made in the attached program: test_splf, change its
53; name in "psplinf", so that IDL can execute either its intrinsic or
54; the substitute.
55; Anyway, this comparison has been made with success and if you just
56; want to use this function in GDL ... ignore this remark
57; -----------------------------------------------------------------------------
58
59;function PSPLINF, x, y, YP0=yp1, YPN_1=ypn, DOUBLE=double
60function SPL_INIT_OLD, x, y, YP0=yp1, YPN_1=ypn, DOUBLE=double, debug=debug
61
62print, 'GDL syntax, a C++ version is now available'
63
64n = N_ELEMENTS(x)
65
66; we should use the same type as for the input y !
67y2 = DBLARR(n)
68u = DBLARR(n)
69;
70; The lower boundary condition is set either to be "natural"
71;
72if (N_ELEMENTS(yp1) EQ 0) then begin
73    y2[0] = 0.
74    u[0] = 0.
75    ;;
76    ;; or else to have a specified first derivative
77    ;;
78endif else begin
79    y2[0] = -0.5
80    u[0] = ( 3. / ( x[1]-x[0] ) ) * ( ( y[1]-y[0] ) / ( x[1]-x[0] ) - yp1 )
81endelse
82
83; I suppose we can also take advantage here of the TRISOL function
84; from IDL... we can remove the for loops
85;
86; This is the decomposition loop of the tridiagonal algorithm.  Y2 and
87; U are used for temporary storage of the decomposed factors.
88;
89
90; AC 09/08/25 : bug which may appear with unsorted X data ...
91psig_old = DOUBLE((x - SHIFT(x, -1))) / (SHIFT(x, +1) - SHIFT(x, -1))
92psig=psig_old
93for ii=1, n-2 do psig[ii]=(x[ii]-x[ii-1])/((x[ii+1]-x[ii-1]))
94
95if KEYWORD_SET(debug) then begin
96   print, psig_old
97   print, psig
98endif
99
100pu = (DOUBLE(SHIFT(y,-1) - y) / (SHIFT(x,-1) - x) - $
101      (y - SHIFT(y,1)) / (x - SHIFT(x,1))) / (SHIFT(x,-1)- SHIFT(x,1))
102
103for i=1,n-2 do begin
104    p = psig[i] * y2[i-1] + 2.0D
105    y2[i] = ( psig[i]-1.0D ) / p
106    u[i]=( 6.0D * pu[i] - psig[i]*u[i-1] ) / p
107    ;;if KEYWORD_SET(debug) then print, i, psig[i], pu[i], p
108endfor
109
110;
111; The upper boundary condition is set either to be "natural"
112;
113if n_elements(ypn) eq 0 then begin
114    qn=0.
115    un=0.
116;
117; or else to have a specified first deriviative
118;
119endif else begin
120    qn=0.5
121    dx=x[n-1]-x[n-2]
122    un=(3./dx)*(ypn-(y[n-1]-y[n-2])/dx)
123endelse
124;
125y2[n-1] = ( un - qn * u[n-2] ) / ( qn * y2[n-2] + 1. )
126
127;
128; This is the backsubstitution loop of the tridiagonal algorithm
129;
130if KEYWORD_SET(debug) then begin
131   print, "psig:",psig
132   print, "pu:",pu
133   for ii=0, n-1 do print, ii, y2[ii] , u[ii]
134endif
135
136for k=n-2,0,-1 do begin
137    y2[k] = y2[k] * y2[k+1] + u[k]
138endfor
139
140;
141return, y2
142
143end	; spl_init.pro
144
145