Path: utzoo!utgpu!news-server.csri.toronto.edu!helios.physics.utoronto.ca!ists!nereid!wagener
From: wagener@nereid.ists.ca (Richard Wagener)
Newsgroups: comp.lang.idl-pvwave
Subject: Hermite Cubic Spline function
Message-ID: <10443@ists.ists.ca>
Date: 15 Jul 90 00:11:13 GMT
Sender: news@ists.ists.ca
Lines: 171

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