diff --git a/src/idl_special/int_tabulated.pro b/src/idl_special/int_tabulated.pro new file mode 100644 index 0000000..7983bff --- /dev/null +++ b/src/idl_special/int_tabulated.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 + diff --git a/src/idl_special/interpol.pro b/src/idl_special/interpol.pro new file mode 100644 index 0000000..358313a --- /dev/null +++ b/src/idl_special/interpol.pro @@ -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 diff --git a/src/idl_special/str_sep.pro b/src/idl_special/str_sep.pro new file mode 100644 index 0000000..f59d2dc --- /dev/null +++ b/src/idl_special/str_sep.pro @@ -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 -- libgit2 0.21.2