cwplib.doc

seismic software,very useful

nxin		i number of x values at which y(x) is input
dxin		i x sampling interval for input y(x)
fxin		i x value of first sample input
yin		i array of input y(x) values:  yin[0] = y(fxin), etc.
yinl		i value used to extrapolate yin values to left of yin[0]
yinr		i value used to extrapolate yin values to right of yin[nxin-1]
nxout		i number of x values a which y(x) is output
xout		i array of x values at which y(x) is output
yout		o array of output y(x) values:  yout[0] = y(xout[0]), etc.


Notes on parameters:

	ntable must not be less than 2.

	The table of interpolation operators must be as follows:

	Let d be the distance, expressed as a fraction of dxin, from a
	particular xout value to the sampled location xin just to the
	left of xout.  Then, for d = 0.0,

	table[0][0:7] = 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0

	are the weights applied to the 8 input samples nearest xout.
	Likewise, for d = 1.0,

	table[ntable-1][0:7] = 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0

	are the weights applied to the 8 input samples nearest xout.
	In general, for d = (float)itable/(float)(ntable-1),
	table[itable][0:7] are the weights applied to the 8 input
	samples nearest xout.  If the actual sample distance d does not
	exactly equal one of the values for which interpolators are
	tabulated, then the interpolator corresponding to the nearest
	value of d is used.

-----------------------------------------
    h. ress8r.c  - resample a uniformly-sampled real function y(x) by sinc
    i. ress8c.c  - resample a uniformly-sampled complex function y(x) by sinc

Notes:

	Resample a uniformly-sampled real function y(x) via 8-coefficient sinc
	approximations; maximum error for frequiencies less than 0.6 nyquist is 
	less than one percent.

	Because extrapolation of the input function y(x) is defined by the
	left and right values yinl and yinr, the output x values defined
	by nxout, dxout, and fxout are not restricted to lie within the range 
	of input x values defined by nxin, dxin, and fxin.


Prototypes and parameters:

void ress8r (int nxin, float dxin, float fxin, float yin[], 
	float yinl, float yinr, 
	int nxout, float dxout, float fxout, float yout[]);
void ress8c (int nxin, float dxin, float fxin, complex yin[], 
	complex yinl, complex yinr, 
	int nxout, float dxout, float fxout, complex yout[]);

nxin		i number of x values at which y(x) is input
dxin		i x sampling interval for input y(x)
fxin		i x value of first sample input
yin		i array of input y(x) values:  yin[0] = y(fxin), etc.
yinl		i value used to extrapolate yin values to left of yin[0]
yinr		i value used to extrapolate yin values to right of yin[nxin-1]
nxout		i number of x values at which y(x) is output
dxout		i x sampling interval for output y(x)
fxout		i x value of first sample output
yout		o array of output y(x) values:  yout[0] = y(xout[0]), etc.

-----------------------------------------
  j. shfs8r - shift a uniformly-sampled real-valued function by sinc

Notes:

	This routine shifts a uniformly-sampled real-valued function
	y(x) via a table of 8-coefficient sinc approximations; maximum
	error for frequencies less than 0.6*nyquist is less than one
	percent.

	Because extrapolation of the input function y(x) is defined by
	the left and right values yinl and yinr, the output samples
	defined by dx, nxout, and fxout are not restricted to lie
	within the range of input sample locations defined by dx, nxin,
	and fxin.


Prototype and parameters:

void shfs8r (float dx, int nxin, float fxin, float yin[], 
	float yinl, float yinr, int nxout, float fxout, float yout[]);

dx			i x sampling interval for both input and output y(x)
nxin		i number of x values at which y(x) is input
fxin		i x value of first sample input
yin			i array of input y(x) values:  yin[0] = y(fxin), etc.
yinl		i value used to extrapolate yin values to left of yin[0]
yinr		i value used to extrapolate yin values to right of yin[nxin-1]
nxout		i number of x values a which y(x) is output
fxout		i x value of first sample output
yout		o array of output y(x) values:  yout[0] = y(fxout), etc.

-----------------------------------------
   k. xindex - determine index of x with respect to an array of x values

Notes:

	This function is designed to be particularly efficient when
	called repeatedly for slightly changing x values; in such
	cases, the index returned from one call should be used in the
	next.

Prototype and parameters:
void xindex (int nx, float ax[], float x, int *index);

nx	i number of x values in array ax
ax	i array[nx] of monotonically increasing or decreasing x values
x	i the value for which index is to be determined
index	i index determined previously (used to begin search)

index	o for monotonically increasing ax values, the largest index
	  for which ax[index]<=x, except index=0 if ax[0]>x;
	  for monotonically decreasing ax values, the largest index
	  for which ax[index]>=x, except index=0 if ax[0]<x

-----------------------------------------
  l. intl2b - bilinear interpolation of a 2-D array of bytes

Notes:
	The arrays zin and zout must passed as pointers to the first
	element of a two-dimensional contiguous array of unsigned char
	values.

	Constant extrapolation of zin is used to compute zout for
	output x and y outside the range of input x and y.
 
	For efficiency, this function builds a table of interpolation
	coefficents pre-multiplied by byte values.  To keep the table
	reasonably small, the interpolation does not distinguish
	between x and y values that differ by less than dxin/ICMAX and
	dyin/ICMAX, respectively, where ICMAX is a parameter #defined
	in the code (now 100).


Prototype and parameters:

void intl2b(int nxin, float dxin, float fxin,
	int nyin, float dyin, float fyin, unsigned char *zin,
	int nxout, float dxout, float fxout,
	int nyout, float dyout, float fyout, unsigned char *zout);

nxin	i number of x samples input (fast dimension of zin)
dxin	i x sampling interval input
fxin	i first x sample input
nyin	i number of y samples input (slow dimension of zin)
dyin	i y sampling interval input
fyin	i first y sample input
zin	i array[nyin][nxin] of input samples (see notes)
nxout	i number of x samples output (fast dimension of zout)
dxout	i x sampling interval output
fxout	i first x sample output
nyout	i number of y samples output (slow dimension of zout)
dyout	i y sampling interval output
fyout	i first y sample output

zout	o array[nyout][nxout] of output samples (see notes)

-----------------------------------------
   m. intlin - evaluate y(x) via linear interpolation of y(x[0]), y(x[1]), ...

Notes:
	xin values must be monotonically increasing or decreasing.

	Extrapolation of the function y(x) for xout values outside the
	range spanned by the xin values is performed as follows:

	For monotonically increasing xin values,
		yout=yinl if xout<xin[0], and yout=yinr if xout>xin[nin-1].

	For monotonically decreasing xin values, 
		yout=yinl if xout>xin[0], and yout=yinr if xout<xin[nin-1].

	If nin==1, then the monotonically increasing case is used.


Prototype and parameters:

void intlin (int nin, float xin[], float yin[], float yinl, float yinr, 
	int nout, float xout[], float yout[]);

nin	i length of xin and yin arrays
xin	i array[nin] of monotonically increasing or decreasing x values
yin	i array[nin] of input y(x) values
yinl	i value used to extrapolate y(x) to left of input yin values
yinr	i value used to extrapolate y(x) to right of input yin values
nout	i length of xout and yout arrays
xout	i array[nout] of x values at which to evaluate y(x)

yout	o array[nout] of linearly interpolated y(x) values

-----------------------------------------
   n. intcub - eval. y(x), y'(x), y''(x), ... via piecewise cubic interpolation

Notes:

	xin values must be monotonically increasing or decreasing.

	Extrapolation of the function y(x) for xout values outside the
	range spanned by the xin values is performed using the
	derivatives in ydin[0][0:3] or ydin[nin-1][0:3], depending on
	whether xout is closest to xin[0] or xin[nin-1], respectively.


Prototype and parameters:

void intcub (int ideriv, int nin, float xin[], float ydin[][4], 
	int nout, float xout[], float yout[]);

ideriv	i =0 if y(x) desired; =1 if y'(x) desired, ...
nin	i length of xin and ydin arrays; nin>=2 is required
xin	i array[nin] of monotonically increasing or decreasing x values
ydin	i array[nin][4] of y(x), y'(x), y''(x), and y'''(x)
nout	i length of xout and yout arrays
xout	i array[nout] of x values at which to evaluate y(x), y'(x), ...

yout	o array[nout] of y(x), y'(x), ... values

-----------------------------------------
   o. cakima - compute cubic interpolation coefficients via Akima's method

Notes:

	The computed cubic spline coefficients are as follows:
	yd[i][0] = y(x[i])    (the value of y at x = x[i])
	yd[i][1] = y'(x[i])   (the 1st derivative of y at x = x[i])
	yd[i][2] = y''(x[i])  (the 2nd derivative of y at x = x[i])
	yd[i][3] = y'''(x[i]) (the 3rd derivative of y at x = x[i])

	To evaluate y(x) for x between x[i] and x[i+1] and h = x-x[i],
	use the computed coefficients as follows:
	y(x) = yd[i][0]+h*(yd[i][1]+h*(yd[i][2]/2.0+h*yd[i][3]/6.0))

	Akima's method provides continuous 1st derivatives, but 2nd and
	3rd derivatives are discontinuous.  Akima's method is not linear, 
	in that the interpolation of the sum of two functions is not the 
	same as the sum of the interpolations.

	For more information, see Akima, H., 1970, A new method for 
	interpolation and smooth curve fitting based on local procedures,
	Journal of the ACM, v. 17, n. 4, p. 589-602.


Prototype and parameters:

void cakima (int n, float x[], float y[], float yd[][4]);

n	i number of samples (n>=2 is required and assumed)
x  	i array of monotonically increasing or decreasing abscissae
y	i array of ordinates
yd	o array of cubic interpolation coefficients

-----------------------------------------
   p. cmonot - compute cubic interpolation coefficients by Fritsch-Carlson

Notes:
	The Fritsch-Carlson method preserves monotonicity.

	The computed cubic spline coefficients are as follows:
	yd[i][0] = y(x[i])    (the value of y at x = x[i])
	yd[i][1] = y'(x[i])   (the 1st derivative of y at x = x[i])
	yd[i][2] = y''(x[i])  (the 2nd derivative of y at x = x[i])
	yd[i][3] = y'''(x[i]) (the 3rd derivative of y at x = x[i])

	To evaluate y(x) for x between x[i] and x[i+1] and h = x-x[i],
	use the computed coefficients as follows:
	y(x) = yd[i][0]+h*(yd[i][1]+h*(yd[i][2]/2.0+h*yd[i][3]/6.0))

	The Fritsch-Carlson method yields continuous 1st derivatives,
	but 2nd and 3rd derivatives are discontinuous.  The method will
	yield a monotonic interpolant for monotonic data.  1st
	derivatives are set to zero wherever first divided differences
	change sign.

	For more information, see Fritsch, F. N., and Carlson, R. E., 1980, 
	Monotone piecewise cubic interpolation:  SIAM J. Numer. Anal., v. 17,
	n. 2, p. 238-246.

	Also, see the book by Kahaner, D., Moler, C., and Nash, S., 1989, 
	Numerical Methods and Software, Prentice Hall.  This function was 
	derived from SUBROUTINE PCHEZ contained on the diskette that comes 
	with the book.


Prototype and parameters:

void cmonot (int n, float x[], float y[], float yd[][4]);

n	i number of samples (n>=2 is required and assumed)
x  	i array[n] of monotonically increasing or decreasing abscissae
y	i array[n] of ordinates

yd	o array[n][4] of cubic interpolation coefficients (see notes)

-----------------------------------------
   q. csplin - cubic spline interpolation coeffs with continuous 2nd derivatives

Notes:
	The computed cubic spline coefficients are as follows:
	yd[i][0] = y(x[i])    (the value of y at x = x[i])
	yd[i][1] = y'(x[i])   (the 1st derivative of y at x = x[i])
	yd[i][2] = y''(x[i])  (the 2nd derivative of y at x = x[i])
	yd[i][3] = y'''(x[i]) (the 3rd derivative of y at x = x[i])

	To evaluate y(x) for x between x[i] and x[i+1] and h = x-x[i],
	use the computed coefficients as follows:
	y(x) = yd[i][0]+h*(yd[i][1]+h*(yd[i][2]/2.0+h*yd[i][3]/6.0))

Prototype and parameters:

void csplin (int n, float x[], float y[], float yd[][4]);

n	i number of samples (n>=2 is required and assumed)
x  	i array[n] of monotonically increasing or decreasing abscissae
y	i array[n] of ordinates

yd	o array[n][4] of cubic interpolation coefficients (see notes)

-----------------------------------------
   r. yxtoxy - compute x(y) from monotone increasing y(x)

Notes:

	This routine computes a regularly-sampled, monotonically
	increasing function x(y) from a regularly-sampled,
	monotonically increasing function y(x) by inverse linear
	interpolation.

	User must ensure that:
	(1) dx>0.0 && dy>0.0
	(2) y[0] < y[1] < ... < y[nx-1]


Prototype and parameters:

void yxtoxy (int nx, float dx, float fx, float y[], 
	int ny, float dy, float fy, float xylo, float xyhi, float x[]);

nx	i number of samples of y(x)
dx	i x sampling interval; dx>0.0 is required
fx	i first x
y	i array of y(x) values; y[0] < y[1] < ... < y[nx-1] required
ny	i number of samples of x(y)
dy	i y sampling interval; dy>0.0 is required
fy	i first y
xylo	i x value assigned to x(y) when y is less than smallest y(x)
xyhi	i x value assigned to x(y) when y is greater than largest y(x)
x	o array of x(y) values


------------------------------------------------------------------------
6. bf      - package of Butterworth filter routines

Contents:

	bfdesign	design a Butterworth filter
	bfhighpass	apply a high-pass Butterworth filter 
	bflowpass	apply a low-pass Butterworth filter 

Prototypes:

void bfhighpass (int npoles, float f3db, int n, float p[], float q[]);
void bflowpass (int npoles, float f3db, int n, float p[], float q[]);
void bfdesign (float fpass, float apass, float fstop, float astop,
	int *npoles, float *f3db);

------------------------------------------------------------------------
 7. Differentiator approximations