added IDL lib files that may become obsolete or retired in future IDL releases

Annie Hughes
1 parent 86f1665a
Showing 3 changed files with 552 additions and 0 deletions Show diff stats
src/idl_special/int_tabulated.pro
src/idl_special/interpol.pro
src/idl_special/str_sep.pro
@@ -0,0 +1,171 @@
+; $Id: //depot/Release/ENVI52_IDL84/idl/idldir/lib/int_tabulated.pro#1 $
+;
+; Copyright (c) 1995-2014, Exelis Visual Information Solutions, Inc. All
+;       rights reserved. Unauthorized reproduction is prohibited.
+;+
+; NAME:
+;       INT_TABULATED
+;
+; PURPOSE:
+;       This function integrates a tabulated set of data { x(i) , f(i) },
+;       on the closed interval [min(X) , max(X)].
+;
+; CATEGORY:
+;       Numerical Analysis.
+;
+; CALLING SEQUENCE:
+;       Result = INT_TABULATED(X, F)
+;
+; INPUTS:
+;       X:  The tabulated X-value data. This data may be irregularly
+;           gridded and in random order. If the data is randomly ordered
+;	    you must set the SORT keyword to a nonzero value.
+;           Duplicate x values will result in a warning message.
+;       F:  The tabulated F-value data. Upon input to the function
+;           X(i) and F(i) must have corresponding indices for all
+;	    values of i. If X is reordered, F is also reordered.
+;
+;       X and F must be of floating point or double precision type.
+;
+; KEYWORD PARAMETERS:
+;       SORT:   A zero or non-zero scalar value.
+;               SORT = 0 (the default) The tabulated x-value data is
+;                        already in ascending order.
+;               SORT = 1 The tabulated x-value data is in random order
+;                        and requires sorting into ascending order. Both
+;			 input parameters X and F are returned sorted.
+;       DOUBLE: If set to a non-zero value, computations are done in
+;               double precision arithmetic.
+;
+; OUTPUTS:
+;       This fuction returns the integral of F computed from the tabulated
+;	data in the closed interval [min(X) , max(X)].
+;
+; RESTRICTIONS:
+;       Data that is highly oscillatory requires a sufficient number
+;       of samples for an accurate integral approximation.
+;
+; PROCEDURE:
+;       INT_TABULATED.PRO constructs a regularly gridded x-axis with a
+;	number of segments as an integer multiple of four. Segments
+;	are processed in groups of four using a 5-point Newton-Cotes
+;	integration formula.
+;       For 'sufficiently sampled' data, this algorithm is highly accurate.
+;
+; EXAMPLES:
+;       Example 1:
+;       Define 11 x-values on the closed interval [0.0 , 0.8].
+;         x = [0.0, .12, .22, .32, .36, .40, .44, .54, .64, .70, .80]
+;
+;       Define 11 f-values corresponding to x(i).
+;         f = [0.200000, 1.30973, 1.30524, 1.74339, 2.07490, 2.45600, $
+;              2.84299,  3.50730, 3.18194, 2.36302, 0.231964]
+;
+;       Compute the integral.
+;         result = INT_TABULATED(x, f)
+;
+;       In this example, the f-values are generated from a known function,
+;       (f = .2 + 25*x - 200*x^2 + 675*x^3 - 900*x^4 + 400*x^5)
+;
+;       The Multiple Application Trapazoid Method yields;  result = 1.5648
+;       The Multiple Application Simpson's Method yields;  result = 1.6036
+;				INT_TABULATED.PRO yields;  result = 1.6232
+;         The Exact Solution (4 decimal accuracy) yields;  result = 1.6405
+;
+;	Example 2: 
+;       Create 30 random points in the closed interval [-2 , 1].
+;         x = randomu(seed, 30) * 3.0 - 2.0
+;
+;       Explicitly define the interval's endpoints.
+;         x(0) = -2.0  &  x(29) = 1.0
+;
+;       Generate f(i) corresponding to x(i) from a given function.
+;         f = sin(2*x) * exp(cos(2*x))
+;
+;       Call INT_TABULATED with the SORT keyword.
+;         result = INT_TABULATED(x, f, /sort)
+;
+;       In this example, the f-values are generated from the function,
+;       f = sin(2*x) * exp(cos(2*x))
+;
+;       The result of this example will vary because the x(i) are random.
+;       Executing this example three times gave the following results:
+;		               INT_TABULATED.PRO yields;  result = -0.0702
+;		               INT_TABULATED.PRO yields;  result = -0.0731
+;		               INT_TABULATED.PRO yields;  result = -0.0698
+;        The Exact Solution (4 decimal accuracy) yields;  result = -0.0697
+;
+; MODIFICATION HISTORY:
+;           Written by:  GGS, RSI, September 1993
+;           Modified:    GGS, RSI, November  1993
+;                        Use Numerical Recipes cubic spline interpolation 
+;                        function NR_SPLINE/NR_SPLINT. Execution time is 
+;                        greatly reduced. Added DOUBLE keyword. The 'sigma' 
+;                        keyword is no longer supported.
+;           Modified:    GGS, RSI, April  1995
+;                        Changed cubic spline calls from NR_SPLINE/NR_SPLINT
+;                        to SPL_INIT/SPL_INTERP. Improved double-precision
+;                        accuracy.
+;           Modified:    GGS, RSI, April 1996
+;                        Replaced WHILE loop with vector operations. 
+;                        Check for duplicate points in x vector.  
+;                        Modified keyword checking and use of double precision.
+;-
+
+FUNCTION Int_Tabulated, X, F, Double = Double, Sort = Sort
+
+  ;Return to caller if an error occurs.
+  ON_ERROR, 2 
+
+  TypeX = SIZE(X)
+  TypeF = SIZE(F)
+
+  ;Check F data type.
+  if TypeF[TypeF[0]+1] ne 4 and TypeF[TypeF[0]+1] ne 5 then $
+    MESSAGE, "F values must be float or double."
+
+  ;Check length.
+  if TypeX[TypeX[0]+2] ne TypeF[TypeF[0]+2] then $
+    MESSAGE, "X and F arrays must have the same number of elements."
+
+  ;Check duplicate values.
+  if TypeX[TypeX[0]+2] ne N_ELEMENTS(UNIQ(X[SORT(X)])) then $
+    MESSAGE, "X array contains duplicate points."
+
+  ;If the DOUBLE keyword is not set then the internal precision and
+  ;result are identical to the type of input.
+  if N_ELEMENTS(Double) eq 0 then $
+    Double = (TypeX[TypeX[0]+1] eq 5 or TypeF[TypeF[0]+1] eq 5) 
+
+  Xsegments = TypeX[TypeX[0]+2] - 1L
+
+  ;Sort vectors into ascending order.
+  if KEYWORD_SET(Sort) ne 0 then begin
+    ii = SORT(x)
+    X = X[ii]
+    F = F[ii]
+  endif
+
+  while (Xsegments MOD 4L) ne 0L do $
+    Xsegments = Xsegments + 1L
+
+  Xmin = MIN(X)
+  Xmax = MAX(X)
+
+  ;Uniform step size.
+    h = (Xmax+0.0 - Xmin) / Xsegments
+  ;Compute the interpolates at Xgrid.
+    ;x values of interpolates >> Xgrid = h * FINDGEN(Xsegments + 1L) + Xmin
+    z = SPL_INTERP(X, F, SPL_INIT(X, F, Double = Double), $
+                   h * FINDGEN(Xsegments + 1L) + Xmin, Double = Double)
+  ;Compute the integral using the 5-point Newton-Cotes formula.
+    ii = (LINDGEN((N_ELEMENTS(z) - 1L)/4L)+1) * 4
+    if Double eq 0 then $
+      RETURN, FLOAT(TOTAL(2.0 * h * (7.0 * (z[ii-4] + z[ii]) + $
+                    32.0 * (z[ii-3] + z[ii-1]) + 12.0 * z[ii-2]) / 45.0)) $
+    else $
+      RETURN, TOTAL(2D * h * (7D * (z[ii-4] + z[ii]) + $
+                    32D * (z[ii-3] + z[ii-1]) + 12D * z[ii-2]) / 45D, /DOUBLE)
+
+END
+
@@ -0,0 +1,265 @@
+; $Id: //depot/Release/ENVI52_IDL84/idl/idldir/lib/interpol.pro#1 $
+;
+; Copyright (c) 1982-2014, Exelis Visual Information Solutions, Inc. All
+;       rights reserved. Unauthorized reproduction is prohibited.
+
+Function ls2fit, xx, y, xm
+
+COMPILE_OPT idl2, hidden
+
+x = xx - xx[0]                  ;Normalize to preserve significance.
+ndegree = 2L
+n = n_elements(xx)
+
+corrm = fltarr(ndegree+1, ndegree+1) ;Correlation matrix
+b = fltarr(ndegree+1)
+
+corrm[0,0] = n                  ;0 - Form the normal equations
+b[0] = total(y)
+z = x                           ;1
+b[1] = total(y*z)
+corrm[[0,1],[1,0]] = total(z)
+z = z * x                       ;2
+b[2] = total(y*z)
+corrm[[0,1,2], [2,1,0]] = total(z)
+z = z * x                       ;3
+corrm[[1,2],[2,1]] = total(z)
+corrm[2,2] = total(z*x)         ;4
+
+c = b # invert(corrm)
+xm0 = xm - xx[0]
+return, c[0] + c[1] * xm0 + c[2] * xm0^2
+end
+
+
+
+;+
+; NAME:
+;	INTERPOL
+;
+; PURPOSE:
+;	Linearly interpolate vectors with a regular or irregular grid.
+;	Quadratic or a 4 point least-square fit to a quadratic
+;	interpolation may be used as an option.
+;
+; CATEGORY:
+;	E1 - Interpolation
+;
+; CALLING SEQUENCE:
+;	Result = INTERPOL(V, N) 	;For regular grids.
+;
+;	Result = INTERPOL(V, X, XOUT)	;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.
+;
+;	Irregular grids:
+;	X:	The absicissae values for V.  This vector must have same # of
+;		elements as V.  The values MUST be monotonically ascending
+;		or descending.
+;
+;	XOUT:	The absicissae values for the result.  The result will have
+;		the same number of elements as XOUT.  XOUT does not need to be
+;		monotonic.  If XOUT is outside the range of X, then the
+;		closest two endpoints of (X,V) are linearly extrapolated.
+;
+; Keyword Input Parameters:
+; NAN = if set, then filter out NaN values before interpolating.
+;   The default behavior is to include the NaN values - by including NaN's
+;   the output will contain NaN's in locations where the interpolation
+;   result is undefined.
+;   
+;	LSQUADRATIC = if set, interpolate using a least squares
+;	  quadratic fit to the equation y = a + bx + cx^2, for each 4
+;	  point neighborhood (x[i-1], x[i], x[i+1], x[i+2]) surrounding
+;	  the interval, x[i] <= XOUT < x[i+1].
+;
+; 	QUADRATIC = if set, interpolate by fitting a quadratic
+; 	  y = a + bx + cx^2, to the three point neighborhood (x[i-1],
+; 	  x[i], x[i+1]) surrounding the interval x[i] <= XOUT < x[i+1].
+;
+;	SPLINE = if set, interpolate by fitting a cubic spline to the
+;	  4 point neighborhood (x[i-1], x[i], x[i+1], x[i+2]) surrounding
+;	  the interval, x[i] <= XOUT < x[i+1].
+;
+;	Note: if LSQUADRATIC or QUADRATIC or SPLINE is not set,  the
+;	default linear interpolation is used.
+;
+; OUTPUTS:
+;	INTERPOL returns a floating-point vector of N points determined
+;	by interpolating the input vector by the specified method.
+;
+;	If the input vector is double or complex, the result is double
+;	or complex.
+;
+; COMMON BLOCKS:
+;	None.
+;
+; SIDE EFFECTS:
+;	None.
+;
+; RESTRICTIONS:
+;	None.
+;
+; PROCEDURE:
+;	For linear interpolation,
+;	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 = XOUT(i).
+;		m = number of points of input vector.
+;
+;	  For QUADRATIC interpolation, the equation y=a+bx+cx^2 is
+;	solved explicitly for each three point interval, and is then
+;	evaluated at the interpolate.
+;	  For LSQUADRATIC interpolation, the coefficients a, b, and c,
+;	from the above equation are found, for the four point
+;	interval surrounding the interpolate using a least square
+;	fit.  Then the equation is evaluated at the interpolate.
+;	  For SPLINE interpolation, a cubic spline is fit over the 4
+;	point interval surrounding each interpolate, using the routine
+;	SPL_INTERP().
+;
+; 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.
+;	Modified, DMS, August, 1998.  Now use binary intervals which
+;		speed things up considerably when XOUT is random.
+;       DMS, May, 1999.  Use new VALUE_LOCATE function to find intervals,
+;       	which speeds things up by a factor of around 100, when
+;       	interpolating from large arrays.  Also added SPLINE,
+;       	QUADRATIC, and LSQUADRATIC keywords.
+;   CT, VIS, Feb 2009: CR54942: Automatically filter out NaN values.
+;     Clean up code.
+;   CT, VIS, Feb 2010: CR57403: Back out previous change.
+;     Add NAN keyword to control filtering out of NaN's
+;-
+;
+FUNCTION INTERPOL, VV, XX, XOUT, $
+  SPLINE=spline, LSQUADRATIC=ls2, QUADRATIC=quad, NAN=nan
+
+  COMPILE_OPT idl2
+  
+  on_error,2                      ;Return to caller if an error occurs
+  
+  regular = n_params(0) eq 2
+  
+  ; Make a copy so we don't overwrite the input arguments.
+  v = vv
+  x = xx
+  m = N_elements(v)               ;# of input pnts
+
+  if (regular) then nOut = LONG(x)
+  
+  ; Filter out NaN values in both the V and X arguments.
+  if (KEYWORD_SET(nan)) then begin
+    isNAN = FINITE(v, /NAN)
+    if (~regular) then isNAN or= FINITE(x, /NAN)
+    
+    if (~ARRAY_EQUAL(isNAN, 0)) then begin
+      good = WHERE(~isNAN, ngood)
+      if (ngood gt 0 && ngood lt m) then begin
+        v = v[good]
+        if (regular) then begin
+          ; We supposedly had a regular grid, but some of the values
+          ; were NaN (missing). So construct the irregular grid.
+          regular = 0b
+          x = LINDGEN(m)
+          xout = FINDGEN(nOut) * ((m-1.0) / ((nOut-1.0) > 1.0)) ;Grid points
+        endif
+        x = x[good]
+      endif
+    endif
+  endif
+
+  ; get the number of input points again, in case some NaN's got filtered
+  m = N_elements(v)
+  type = SIZE(v, /TYPE)
+
+  if regular && $                ;Simple regular case?
+    ((keyword_set(ls2) || keyword_set(quad) || keyword_set(spline)) eq 0) $
+    then begin
+    xout = findgen(nOut)*((m-1.0)/((nOut-1.0) > 1.0)) ;Grid points in V
+    xoutInt = long(xout)		;Cvt to integer
+    case (type) of
+    1: diff = v[1:*] - FIX(v)
+    12: diff = v[1:*] - LONG(v)
+    13: diff = v[1:*] - LONG64(v)
+    15: diff = LONG64(v[1:*]) - LONG64(v)
+    else: diff = v[1:*] - v
+    endcase
+    return, V[xoutInt] + (xout-xoutInt)*diff[xoutInt] ;interpolate
+  endif
+
+  if regular then begin           ;Regular intervals??
+    xout = findgen(nOut) * ((m-1.0) / ((nOut-1.0) > 1.0)) ;Grid points
+    s = long(xout)                 ;Subscripts
+  endif else begin                ;Irregular
+    if n_elements(x) ne m then $
+      message, 'V and X arrays must have same # of elements'
+    s = VALUE_LOCATE(x, xout) > 0L < (m-2) ;Subscript intervals.
+  endelse
+
+  ; Clip interval, which forces extrapolation.
+  ; XOUT[i] is between x[s[i]] and x[s[i]+1].
+  
+  CASE (1) OF
+  
+  KEYWORD_SET(ls2): BEGIN  ;Least square fit quadratic, 4 points
+    s = s > 1L < (m-3)          ;Make in range.
+    p = replicate(v[0]*1.0, n_elements(s)) ;Result
+    for i=0L, n_elements(s)-1 do begin
+        s0 = s[i]-1
+        p[i] = ls2fit(regular ? s0+findgen(4) : x[s0:s0+3], v[s0:s0+3], xout[i])
+    endfor
+    END
+    
+  KEYWORD_SET(quad): BEGIN ;Quadratic.
+    s = s > 1L < (m-2)          ;In range
+    x1 = regular ? float(s) : x[s]
+    x0 = regular ? x1-1.0 : x[s-1]
+    x2 = regular ? x1+1.0 : x[s+1]
+    p = v[s-1] * (xout-x1) * (xout-x2) / ((x0-x1) * (x0-x2)) + $
+        v[s] *   (xout-x0) * (xout-x2) / ((x1-x0) * (x1-x2)) + $
+        v[s+1] * (xout-x0) * (xout-x1) / ((x2-x0) * (x2-x1))
+    END
+    
+  KEYWORD_SET(spline): BEGIN
+    s = s > 1L < (m-3)          ;Make in range.
+    p = replicate(v[0], n_elements(s)) ;Result
+    sold = -1
+    for i=0L, n_elements(s)-1 do begin
+        s0 = s[i]-1
+        if sold ne s0 then begin
+            x0 = regular ? s0+findgen(4): x[s0: s0+3]
+            v0 = v[s0: s0+3]
+            q = spl_init(x0, v0)
+            sold = s0
+        endif
+        p[i] = spl_interp(x0, v0, q, xout[i])
+    endfor
+    END
+    
+  ELSE: begin              ;Linear, not regular
+    case (type) of
+    1: diff = v[s+1] - FIX(v[s])
+    12: diff = v[s+1] - LONG(v[s])
+    13: diff = v[s+1] - LONG64(v[s])
+    15: diff = LONG64(v[s+1]) - LONG64(v[s])
+    else: diff = v[s+1] - v[s]
+    endcase
+    p = (xout-x[s])*diff/(x[s+1] - x[s]) + v[s]
+    end
+    
+  ENDCASE
+
+  RETURN, p
+end
@@ -0,0 +1,116 @@
+; $Id: //depot/Release/ENVI52_IDL84/idl/idldir/lib/obsolete/str_sep.pro#1 $
+;
+; Copyright (c) 1992-2014, Exelis Visual Information Solutions, Inc and
+;       CreaSo Creative Software Systems GmbH.
+;       All rights reserved. Unauthorized reproduction prohibited.
+;+
+; NAME:
+;    STR_SEP
+;
+; PURPOSE:
+;    This routine cuts a string into pieces which are separated by the
+;    separator string.
+; CATEGORY:
+;    String processing.
+; CALLING SEQUENCE:
+;    arr = STR_SEP(str, separator)
+;
+; INPUTS:
+;    str - The string to be separated.
+;    separator - The separator.
+;
+; KEYWORDS:
+;    ESC = escape character.  Only valid if separator is a single character.
+;		Characters following the escape character are treated
+;		literally and not interpreted as separators.
+;		For example, if the separator is a comma,
+;		and the escape character is a backslash, the character
+;		sequence 'a\,b' is a single field containing the characters
+;		'a,b'.
+;    REMOVE_ALL = if set, remove all blanks from fields.
+;    TRIM = if set, remove only leading and trailing blanks from fields.
+;
+; OUTPUT:
+;    An array of strings as function value.
+;
+; COMMON BLOCKS:
+;    None
+;
+; SIDE EFFECTS:
+;    No known side effects.
+;
+; RESTRICTIONS:
+;    None.
+;
+; EXAMPLE:
+;    array = STR_SEP ("ulib.usca.test", ".")
+;
+; MODIFICATION HISTORY:
+;	July 1992, AH,	CreaSo		Created.
+;	December, 1994, DMS, RSI	Added TRIM and REMOVE_ALL.
+;-
+function STR_SEP, str, separator, REMOVE_ALL = remove_all, TRIM = trim, ESC=esc
+
+
+ON_ERROR, 2
+if n_params() ne 2 then message,'Wrong number of arguments.'
+
+spos = 0L
+if n_elements(esc) gt 0 then begin		;Check for escape character?
+  if strpos(str, esc) lt 0 then goto, no_esc	;None in string, use fast case
+  besc = (byte(esc))[0]
+  bsep = (byte(separator))[0]
+  new = bytarr(strlen(str)+1)
+  new[0] = byte(str)
+  j = 0L
+  for i=0L, n_elements(new)-2 do begin
+    if new[i] eq besc then begin
+	new[j] = new[i+1]
+	i = i + 1
+    endif else if new[i] eq bsep then new[j] = 1b $   ;Change seps to 1b char
+    else new[j] = new[i]
+    j = j + 1
+    endfor
+  new = string(new[0:j-1])
+  w = where(byte(new) eq 1b, count)  ;where seps are...
+  arr = strarr(count+1)
+  for i=0L, count-1 do begin
+	arr[i] = strmid(new, spos, w[i]-spos)
+	spos = w[i] + 1
+	endfor
+  arr[count] = strmid(new, spos, strlen(str))  ;Last element
+  goto, done
+  endif			;esc
+
+no_esc:
+if strlen(separator) le 1 then begin	;Single character separator?
+    w = where(byte(str) eq (byte(separator))[0], count)  ;where seps are...
+    arr = strarr(count+1)
+    for i=0, count-1 do begin
+	arr[i] = strmid(str, spos, w[i]-spos)
+	spos = w[i] + 1
+	endfor
+    arr[count] = strmid(str, spos, strlen(str))  ;Last element
+endif else begin		;Multi character separator....
+    n = 0L		   ; Determine number of seperators in string.
+    repeat begin
+	pos = strpos (str, separator, spos)
+	spos = pos + strlen(separator)
+	n = n+1
+    endrep until pos eq -1
+
+    arr = strarr(n)	   ; Create result array
+    spos = 0L
+    for i=0L, n-1 do begin   ; Separate substrings
+      pos = strpos (str, separator, spos)
+      if pos ge 0 then arr[i] = strmid (str, spos, pos-spos) $
+      else arr[i] = strmid(str, spos, strlen(str))
+      spos = pos+strlen(separator)
+   endfor
+endelse
+
+done:
+if keyword_set(trim) then arr = strtrim(arr,2) $
+else if keyword_set(remove_all) then arr = strcompress(arr, /REMOVE_ALL)
+return, arr
+end