pro piececub_interp, Xtab, Ytab, dYtab, Xint, Yint, dYint, err=err ;+ ; NAME: ; PIECECUB_INTERP ; PURPOSE: ; interpolate 1-d data and derivatives with piece-wise cubic functions ; EXPLANATION: ; Interpolation into a function with continuous values and ; derivatives, with tabulated values. Between grid points, the ; function is defined by a cubic polynomial. Unlike cubic ; splines, the second derivative is not necessarily continuous. ; ; Calculations are done in double precision. ; ; CALLING SEQUENCE: ; piececub_interp, Xtab, Ytab, dYtab, Xint, Yint, dYint ; ; INPUT PARAMETERS: ; Xtab - Vector containing the current independent variable grid. ; Must be monotonic increasing or decreasing ; Ytab - Vector containing the current dependent variable values at ; the XTAB grid points. ; dYtab - Vector containing the current derivative of the dependent ; variable values at the XTAB grid points. ; Xint - Scalar or vector containing the new independent variable grid ; points for which interpolated value(s) of the dependent ; variable are sought. Note that -- due to a limitation of the ; intrinsic INTERPOLATE() function -- Xint is always converted to ; floating point internally. ; ; OUTPUT PARAMETERS: ; Yint - Scalar or vector with the interpolated value(s) of the ; dependent variable at the XINT grid points. ; YINT is double precision if Ytab is, float otherwise. ; dYint - Scalar or vector with the interpolated value(s) of the ; derivative at the XINT grid points ; DYINT is double precision if DYtab is, float otherwise. ; ; ; OPTIONAL INPUT KEYWORD: ; ; EXAMPLE: ; ; PROCEDURES CALLED: ; TABINV, ZPARCHECK ; PROCEEDURE ; ; ________ DEFINITIONS ___________________ ; Start with some definitions: ; ; x the independent variable ; x0 x at node 0 (or node n) ; x1 x at node 1 (or node n+1) ; y the dependent variable ; y0 y at node 0 ; dy0 dy/dx at node 0 ; y1 y at node 0 ; dy1 dy/dx at node 0 ; ; Let's define some things to make our lives easier: ; d0 = x-x0 ; d1 = x-x1 ; del = x1-x0 ; so d1 and d0 are functions of x. ; ; d0(x0) = 0 d0(x1) = del ; d1(x0) = -del d1(x1) = 0 ; ; Using d0, d1, and del, define four functions: ; ; P0(x) = (3*d0-d1) * d1^2 / del^3 ; P1(x) = -(3*d1-d0) * d0^2 / del^3 ; Q0(x) = d0 * d1^2 / del^2 ; Q1(x) = d1 * d0^2 / del^2 ; ; with the following derivatives ; ; dP0(x)/dx = 6*d0*d1/del^3 ; dP1(x)/dx = -6*d0*d1/del^3 ; dQ0(x)/dx = d1*(d1+2*d0)/del^2 ; dQ1(x)/dx = d0*(d0+2*d1)/del^2 ; ; These functions have very tidy values and derivatives at the nodes: ; ; P0(x0) = 1 P0(x1) = 0 dP0(x0)/dx = 0 P0(x1)/dx = 0 ; P1(x0) = 0 P1(x1) = 1 dP1(x0)/dx = 0 P1(x1)/dx = 0 ; Q0(x0) = 0 Q0(x1) = 0 dQ0(x0)/dx = 1 Q0(x1)/dx = 0 ; Q1(x0) = 0 Q1(x1) = 0 dQ1(x0)/dx = 0 Q1(x1)/dx = 1 ; ; ________ CALCULATING THE CUBIC FUNCTION ___________________ ; This makes it trivial to define the cubic function that ; matches the nodes: ; ; y(x) = y0 * P0(x) + y1 * P1(x) + dy0 * Q0(x) + dy1 * Q1(x) [CUBIC] ; ; ; MODIFICATION HISTORY: ; ;- err = 1 ; set to true so we can just return after printing error message if N_params() LT 5 then begin print,'Syntax - PIECECUB_INTERP, Xtab, Ytab, dYtab, Xint, Yint, dYint' print,' Xtab, Ytab, dYtab - Input X, Y, dY/dX vectors' print,' Xint - Input X value (scalar or vector) at '$ +'which to interpolate' print,' Yint - Output interpolated Y value(s)' print,' dYint - Output interpolated dY/dX value(s)' return endif ; check parameters numeric = [indgen(5)+1,12,13,14,15] ;Numeric datatypes zparcheck, 'PIECECUB_INTERP', Xtab, 1, numeric, 1, 'Current X Vector' zparcheck, 'PIECECUB_INTERP', Ytab, 2, numeric, 1, 'Current Y Vector' zparcheck, 'PIECECUB_INTERP', dYtab, 3, numeric, 1, 'Current dY/dX Vector' zparcheck, 'PIECECUB_INTERP', Xint, 4, numeric, [0,1], $ 'New X Vector or Scalar' npts = n_elements(Xtab) if n_elements(Ytab) ne npts or n_elements(dYtab) ne npts then begin print,'PIECECUB_INTERP, Xtab, Ytab, dYtab, Xint, Yint, dYint' print,' Xtab, Ytab, and dYtab must be the same length' return end ; promote inputs (do all calculations in double precision) xt = double(xtab) yt = double(ytab) dyt = double(dytab) xi = double(xint) ; Determine index of data-points from which interpolation is made i0 = value_locate(xt, xi) i0 = 0 > i0 < (npts-2) i1 = i0 + 1 ; Perform cubic interpolation x0 = xt[i0] x1 = xt[i1] y0 = yt[i0] y1 = yt[i1] dy0 = dyt[i0] dy1 = dyt[i1] del = x1-x0 d0 = xi-x0 d1 = xi-x1 del = x1-x0 if min(abs(del) le 0.) then begin print,'PIECECUB_INTERP, Xtab, Ytab, dYtab, Xint, Yint, dYint' print,' Xtab cannot have duplicate entries' return end P0 = (3*d0-d1) * d1^2 / del^3 P1 = -(3*d1-d0) * d0^2 / del^3 Q0 = d0 * d1^2 / del^2 Q1 = d1 * d0^2 / del^2 dP0 = 6*d0*d1/del^3 dP1 = -6*d0*d1/del^3 dQ0 = d1*(d1+2*d0)/del^2 dQ1 = d0*(d0+2*d1)/del^2 yi = y0 * P0 + y1 * P1 + dy0 * Q0 + dy1 * Q1 if size(ytab,/type) eq 5 then yint = yi else yint = float(yi) if n_params() ge 5 then begin dyi = y0 * dP0 + y1 * dP1 + dy0 * dQ0 + dy1 * dQ1 if size(dytab,/type) eq 5 then dyint = dyi else dyint = float(dyi) end err = 0 ; all OK return end pro piececub_interp_test ; y = x^3 ; xtab = [0,1,2.d] ; ytab = xtab^3.d ; dytab = 3.*xtab^2.d ; xint = [-0.1,0.2,1.8,2.2] ; piececub_interp, Xtab, Ytab, dYtab, Xint, Yint, dYint ; print, xint, yint, xint^3.d, dyint, 3.*xint^2.d, for='(4F15.5)' ; y = x^2 - 1 (x < 1), y = 2x - 2 xtab = [-1.,1.,2.] ytab = [ 0.,0.,2.] dytab = [-2.,2.,2.] xint = [ 0.0, 0.5, 0.7, 1.2, 1.5] n = n_elements(xint) ycmp = fltarr(n) dycmp = fltarr(n) q = where(xint lt 1) ycmp[q] = xint[q]^2 - 1. dycmp[q] = 2*xint[q] l = where(xint ge 1) ycmp[l] = 2*xint[l] - 2. dycmp[l] = 2 piececub_interp, Xtab, Ytab, dYtab, Xint, Yint, dYint print, xint, yint, ycmp, dyint, dycmp, for='(5F15.5)' end