; $Id: interpol.pro,v 1.4 1997/01/23 00:11:23 dave Exp $ ; ; Copyright (c) 1982-1997, Research Systems, Inc. All rights reserved. ; Unauthorized reproduction prohibited. FUNCTION INTERPOL, V, X, U ;+ ; NAME: ; INTERPOL ; ; PURPOSE: ; Linearly interpolate vectors with a regular or irregular grid. ; ; CATEGORY: ; E1 - Interpolation ; ; CALLING SEQUENCE: ; Result = INTERPOL(V, N) ;For regular grids. ; ; Result = INTERPOL(V, X, U) ;For irregular grids. ; ; INPUTS: ; V: The input vector can be any type except string. ; ; For regular grids: ; N: The number of points in the result when both input and ; output grids are regular. The output grid absicissa values ; equal FLOAT(i)/N_ELEMENTS(V), for i = 0, n-1. ; ; Irregular grids: ; X: The absicissae values for V. This vecotr must have same # of ; elements as V. The values MUST be monotonically ascending ; or descending. ; ; U: The absicissae values for the result. The result will have ; the same number of elements as U. U does not need to be ; monotonic. ; ; OPTIONAL INPUT PARAMETERS: ; None. ; ; OUTPUTS: ; INTERPOL returns a floating-point vector of N points determined ; by linearly interpolating the input vector. ; ; If the input vector is double or complex, the result is double ; or complex. ; ; COMMON BLOCKS: ; None. ; ; SIDE EFFECTS: ; None. ; ; RESTRICTIONS: ; None. ; ; PROCEDURE: ; Result(i) = V(x) + (x - FIX(x)) * (V(x+1) - V(x)) ; ; where x = i*(m-1)/(N-1) for regular grids. ; m = # of elements in V, i=0 to N-1. ; ; For irregular grids, x = U(i). ; m = number of points of input vector. ; ; MODIFICATION HISTORY: ; Written, DMS, October, 1982. ; Modified, Rob at NCAR, February, 1991. Made larger arrays possible ; and correct by using long indexes into the array instead of ; integers. ;- ; on_error,2 ;Return to caller if an error occurs m = N_elements(v) ;# of input pnts if N_params(0) eq 2 then begin ;Regular? r = findgen(x)*(m-1)/(x-1>1) ;Grid points in V rl = long(r) ;Cvt to integer s = size(v) if s(s(0)+1) eq 1 then dif = v[1:*]-fix(v) $ ;V[i+1]-v[i], signed for bytes else dif = v[1:*]-v ;Other types are already signed return, V[rl] + (r-rl)*dif[rl] ;interpolate endif ; if n_elements(x) ne m then $ stop,'INTERPOL - V and X must have same # of elements' n= n_elements(u) ;# of output points m2=m-2 ;last subs in v and x r= fltarr(n)+V[0] ;floating, dbl or cmplx result if x[1] - x[0] ge 0 then s1 = 1 else s1=-1 ;Incr or Decr X ; ix = 0L ;current point for i=0L,n-1 do begin ;point loop d = s1 * (u[i]-x[ix]) ;difference if d eq 0. then r[i]=v[ix] else begin ;at point if d gt 0 then while (s1*(u[i]-x[ix+1]) gt 0) and $ (ix lt m2) do ix=ix+1 else $ while (s1*(u[i]-x[ix]) lt 0) and (ix gt 0) do $ ix=ix-1 r[i] = v[ix] + (u[i]-x[ix])*(v[ix+1]-v[ix])/(x[ix+1]-x[ix]) endelse endfor return,r end