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