; $Id: spline.pro,v 1.15 2004/01/21 15:55:03 scottm Exp $ ; ; Copyright (c) 1983-2004, Research Systems, Inc. All rights reserved. ; Unauthorized reproduction prohibited. ;+ ; NAME: ; SPLINE ; ; PURPOSE: ; This function performs cubic spline interpolation. ; ; CATEGORY: ; Interpolation - E1. ; ; CALLING SEQUENCE: ; Result = SPLINE(X, Y, T [, Sigma]) ; ; INPUTS: ; X: The abcissa vector. Values MUST be monotonically increasing. ; ; Y: The vector of ordinate values corresponding to X. ; ; T: The vector of abcissae values for which the ordinate is ; desired. The values of T MUST be monotonically increasing. ; ; OPTIONAL INPUT PARAMETERS: ; Sigma: The amount of "tension" that is applied to the curve. The ; default value is 1.0. If sigma is close to 0, (e.g., .01), ; then effectively there is a cubic spline fit. If sigma ; is large, (e.g., greater than 10), then the fit will be like ; a polynomial interpolation. ; ; KEYWORDS: ; DOUBLE: Set this keyword to force computations to be done ; using double-precision arithmetic. ; ; OUTPUTS: ; SPLINE returns a vector of interpolated ordinates. ; Result(i) = value of the function at T(i). ; ; RESTRICTIONS: ; Abcissa values must be monotonically increasing. ; ; EXAMPLE: ; The commands below show a typical use of SPLINE: ; ; X = [2.,3.,4.] ;X values of original function ; Y = (X-3)^2 ;Make a quadratic ; T = FINDGEN(20)/10.+2 ;Values for interpolated points. ; ;twenty values from 2 to 3.9. ; Z = SPLINE(X,Y,T) ;Do the interpolation. ; ; MODIFICATION HISTORY: ; Author: Walter W. Jones, Naval Research Laboratory, Sept 26, 1976. ; Reviewer: Sidney Prahl, Texas Instruments. ; Adapted for IDL: DMS, Research Systems, March, 1983. ; CT, RSI, July 2003: Added double precision support and DOUBLE keyword, ; use vector math to speed up the loops. ; CT, RSI, August 2003: Must have at least 3 points. ; ;- function spline,x,y,t,sigmaIn, DOUBLE=double compile_opt idl2 on_error,2 ;Return to caller if an error occurs n = N_ELEMENTS(x) < N_ELEMENTS(y) if (n le 2) then $ MESSAGE, 'X and Y must be arrays of 3 or more elements.' dbl = (N_ELEMENTS(double) gt 0) ? KEYWORD_SET(double) : $ (SIZE(x,/TYPE) eq 5) || (SIZE(y,/TYPE) eq 5) || (SIZE(t,/TYPE) eq 5) sigma = (n_params(0) lt 4) ? (dbl ? 1d : 1.0) : $ (sigmaIn > (dbl ? 0.001d : 0.001)) ;in range? xx = dbl ? DOUBLE(x) : FLOAT(x) yy = dbl ? DOUBLE(y) : FLOAT(y) tt = dbl ? DOUBLE(t) : FLOAT(t) yp = dbl ? DBLARR(2*n) : FLTARR(2*n) ;temp storage delx1 = xx[1] - xx[0] ;1st incr dx1 = (yy[1] - yy[0])/delx1 nm1 = n - 1L np1 = n + 1L delx2 = xx[2]-xx[1] delx12 = xx[2]-xx[0] c1 = -(delx12+delx1)/delx12/delx1 c2 = delx12/delx1/delx2 c3 = -delx1/delx12/delx2 slpp1 = c1*yy[0]+c2*yy[1]+c3*yy[2] deln = xx[nm1]-xx[nm1-1] delnm1 = xx[nm1-1]-xx[nm1-2] delnn = xx[nm1]-xx[nm1-2] c1 = (delnn+deln)/delnn/deln c2 = -delnn/deln/delnm1 c3 = deln/delnn/delnm1 slppn = c3*yy[nm1-2]+c2*yy[nm1-1]+c1*yy[nm1] sigmap = sigma*nm1/(xx[nm1]-xx[0]) dels = sigmap*delx1 exps = exp(dels) sinhs = 0.5d*(exps-1./exps) sinhin = 1./(delx1*sinhs) diag1 = sinhin*(dels*0.5d*(exps+1./exps)-sinhs) diagin = 1./diag1 yp[0] = diagin*(dx1-slpp1) spdiag = sinhin*(sinhs-dels) yp[n] = diagin*spdiag ; Do as much work using vectors as possible. delx2 = xx[1:*] - xx dx2 = (yy[1:*] - yy)/delx2 dels = sigmap*delx2 exps = exp(dels) sinhs = 0.5d *(exps-1./exps) sinhin = 1./(delx2*sinhs) diag2 = sinhin*(dels*(0.5d*(exps+1./exps))-sinhs) diag2 = [0, diag2 + diag2[1:*]] dx2nm1 = dx2[nm1-1] ; need to save this to calc yp[nm1] dx2 = [0, dx2[1:*] - dx2] spdiag = sinhin*(sinhs-dels) ; Need to do an iterative loop for this part. for i=1L,nm1-1 do begin diagin = 1./(diag2[i] - spdiag[i-1]*yp[i+n-1]) yp[i] = diagin*(dx2[i] - spdiag[i-1]*yp[i-1]) yp[i+n] = diagin*spdiag[i] endfor diagin = 1./(diag1-spdiag[nm1-1]*yp[n+nm1-1]) yp[nm1] = diagin*(slppn-dx2nm1-spdiag[nm1-1]*yp[nm1-1]) for i=n-2,0,-1 do $ yp[i] = yp[i] - yp[i+n]*yp[i+1] m = n_elements(t) subs = replicate(long(nm1),m) ;subscripts s = xx[nm1]-xx[0] sigmap = sigma*nm1/s j = 0L for i=1L,nm1 do begin ;find subscript where xx[subs] > t(j) > xx[subs-1] while tt[j] lt xx[i] do begin subs[j]=i j++ if j eq m then goto,done endwhile endfor done: subs1 = subs - 1 del1 = tt-xx[subs1] del2 = xx[subs] - tt dels = xx[subs]-xx[subs1] exps1 = exp(sigmap*del1) sinhd1 = 0.5*(exps1-1./exps1) exps = exp(sigmap*del2) sinhd2 = 0.5*(exps-1./exps) exps = exps1*exps sinhs = 0.5*(exps-1./exps) spl = (yp[subs]*sinhd1+yp[subs1]*sinhd2)/sinhs + $ ((yy[subs]-yp[subs])*del1+(yy[subs1]-yp[subs1])*del2)/dels return, (m eq 1) ? spl[0] : spl end