wagener@nereid.ists.ca (Richard Wagener) (07/15/90)
Subject: Summary: Expires: Sender: Reply-To: wagener@nereid.ists.ca (Richard Wagener) Followup-To: Distribution: Distribution: Organization: Institute for Space and Terrestrial Science Keywords: Answering Tim Abbots' request for sources, although nothing as fancy as he had hoped for, here is an interpolation scheme using hermite cubic splines. It seems to be doing its job well, however, it is very slow for large arrays. I would appreciate any hints as to why this is so. Also included is function indx(x,ti) which returns the index of the point x(i) closest to ti, using the bisection rootfinding technique. ;********** HERMIT ****************************** function hermit,xi,y,t,ierr ;+ ; NAME: ; HERMIT ; PURPOSE: ; Cubic Spline Interpolation ; CATEGORY: ; Interpolation - E1 ; CALLING SEQUENCE: ; Result = HERMIT(X,Y,T) ... or: ; Result = HERMIT(X,Y,T,IERR) ; INPUTS: ; X = abcissa vector. ; Y = vector of ordinate values corresponding to X. ; T = vector of abcissae values for which ordinate is desired. ; ; OPTIONAL INPUT PARAMETERS: ; None. ; ; OUTPUTS: ; Result = vector of interpolated ordinates. Result(i) = value ; of function at T(i). ; ; OPTIONAL OUTPUT PARAMETERS: ; IERR = array of size(T): 1 where extrapolated, 0 otherwise ; ; COMMON BLOCKS: ; None. ; SIDE EFFECTS: ; None. ; RESTRICTIONS: ; Extrapolation is done by duplicating the endpoints, this is stable ; but doesn't give very satisfactory results in most cases. Future ; improvements to this might consider at least a linear extrapolation ; from the last two points at either end. ; Among the notable non-restrictions is the fact that the input ; x-array can be either monotonically increasing or decreasing. ; However, the program does not check the data for monotony. Also ; the vector t where the spline is supposed to be evaluated can be in ; any order, monotonic or not. This keeps the routine very general, ; however it is therefore not the most efficient for the case where t ; is in fact monotonic. ; PROCEDURE: ; Based on the FORTRAN routine published by Graham Hill "INTEP, an ; effective interpolation subroutine" in ; 'Publications of the Dominion Astrophysical Observatory', Vol. XVI, ; No.6, 67-70. The Hermite Spline is supposed to approximate best a ; curve that one would draw by hand, trying to connect the dots. ; MODIFICATION HISTORY: ; Author: Richard Wagener, SAL/ISTS, April, 1989. ; Example: ; 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 = HERMIT(X,Y,T) ;Do the interpolation. ; ;- ; doerr=n_params(0) ge 4 n = n_elements(xi) if n ne n_elements(y) then begin print,'ERROR in HERMIT: x-array size=',n,' ne y-size=',n_elements(y) return,0 endif nt = n_elements(t) p=t if doerr then ierr=t*0 ; if n le 1 then begin print,'Hermit, X and Y must be arrays' return,0 endif ; x = xi * 1. ;Make X values floating if not. idown=x(1) lt x(0) n1=n-1 if idown then begin il=where(t ge x(0)) ig=where(t le x(n1)) endif else begin il=where(t le x(0)) ig=where(t ge x(n1)) endelse s=size(il) ; Extrapolate by duplicating the endpoints if s(0) gt 0 then begin p(il)=y(0) ia=max(il)+1 if doerr then ierr(il)=1 endif else ia=0 s=size(ig) if s(0) gt 0 then begin p(ig)=y(n1) ie=min(ig)-1 if doerr then ierr(ig)=1 endif else ie=nt-1 for it=ia,ie do begin io=0 ti=t(it) i=indx(x,ti) if (x(i)-ti)*(x((i+1)<n1)-ti) gt 0. then i=i-1 if i ne io-1 then begin io=i+1 im=(i-1)>0 ip=(i+2)<n1 lp1=1./(x(i)-x(io)) lp2=-lp1 fp1=(y(io)-y(im))/(x(io)-x(im)) fp2=(y(ip)-y(i))/(x(ip)-x(i)) endif xpi1=ti-x(io) xpi=ti-x(i) l1=xpi1*lp1 l2=xpi*lp2 p(it)=y(i)*(1.-2.*lp1*xpi)*l1*l1 + y(io)*(1.-2.*lp2*xpi1)*l2*l2 $ + fp2*xpi1*l2*l2 + fp1*xpi*l1*l1 endfor return,p end ;****************** INDX ********************* function indx,x,t ; Using bisection, find the root of the function x-t thought to lie ; between x1 and x2. The root, returned as rtsec, is refined until its ; accuracy is +-1. (fashioned after FUNCTION RTSEC of the Numerical Recipes) func=x-t maxit=30 ; Maximum allowed number of iterations x1=0 x2=n_elements(x)-1 f=func(x1) fmid=func(x2) ; Return closes end point if t lies outside of x if f*fmid ge 0 then $ if abs(f) lt abs(fl) then return,x2 $ else return,x1 if f lt 0 then begin rtbis=x1 dx=x2-x1 endif else begin rtbis=x2 dx=x1-x2 endelse for j=1,maxit do begin dx=dx*0.5 xmid=long(rtbis+dx+0.5) fmid=func(xmid) if fmid le 0 then rtbis=xmid if abs(dx) lt 1 or fmid eq 0 then return,rtbis endfor print,'ERROR in INDX: Maximum number of iterations ',maxit,' reached.' end