cwplib.doc

seismic software,very useful

cwp/lib:  	TABLE OF CONTENTS	(07/17/90)


 1. alloc   - package of routines that allocate and free arrays

 2. complex - package for complex number manipulation

 3. bm      - package to handle big matrices in virtual memory

 4. Prime Factor FFT routines
	a. npfa   - find next allowed fft size for pfacc
	b. npfao  - find optimal fft size for pfacc
	c. npfar  - find next allowed fft size for pfarc
	d. npfaro - find optimal fft size for pfarc
	e. pfacc  - prime factor complex-to-complex fft
	f. pfarc  - prime factor real-to-complex fft
	g. pfacr  - prime factor complex-to-real fft
	h. pfa2cc - 2d prime factor complex-to-complex fft
	i. pfa2rc - 2d prime factor real-to-complex fft
	j. pfa2cr - 2d prime factor complex-to-real fft
	k. pfamcc - multiple dimension prime factor complex-to-complex fft

 5. Interpolation routines
	 a. fsinc  - single precision sinc function
	 b. dsinc  - double precision sinc function
	 c. mksinc - compute least-squares optimal sinc coefficients
	 d. ints8r - uniformly-sampled real function by sinc
	 e. ints8c - uniformly-sampled complex function by sinc
	 f. intt8r - uniformly-sampled real function by table
	 g. intt8c - uniformly-sampled complex function by table
         h. ress8r - resample uniformly-sampled real function by sinc
         i. ress8c - resample uniformly-sampled complex function by sinc
         j. shfs8r - shift a uniformly-sampled real-valued function by sinc
	 k. xindex - determine index of x with respect to array of x values
	 l. intl2b - bilinear interpolation of a 2-D array of bytes
	 m. intlin - eval. y(x) by linear interp of y(x[0]), y(x[1]), ...
	 n. intcub - eval. y(x), y'(x), y''(x), ... by piecewise cubic interp
	 o. cakima - compute cubic interpolation coefficients by Akima's method
	 p. cmonot - compute cubic interpolation coefficients by Fritsch-Carlson
	 q. csplin - cubic spline interp coeffs with continuous 2nd derivatives
	 r. yxtoxy - compute x(y) from monotone increasing y(x)

 6. bf      - package of Butterworth filter routines

 7. Differentiator approximations
	a. mkdiff - discrete Taylor series approximation to n'th derivative
	b. mkhdiff - filter approximating the bandlimited half-differentiator

 8. General Signal Processing
	a. conv - compute z = x convolved with y
	b. xcor - compute z = x cross-correlated with y
	c. hilbert - compute Hilbert transform y of x
	d. antialias - anti-alias filter for use before sub-sampling

 9. Discrete Abel transform package

10. Discrete Hankel transform package

11. Sorting and Searching routines
	a. hpsort - sort an array so that a[0] <= a[1] <= ... <= a[n-1]
	b. qksort - sort functions based on Hoare's quicksort algorithm
	c. qkisort - index sort functions based on Hoare's quicksort algorithm

12. Statistical routines
	a. quest - return an estimate of a specified quantile 
	b. questinit - alloc,init,return pointer to state of quantile estimator 
	c. questupdate - update and return a quantile estimate 

13. Single precision BLAS (basic linear algebra subroutines)
	a. isamax - return index of element with maximum absolute value
	b. sasum - return sum of absolute values
	c. saxpy - compute y[i] = a*x[i]+y[i]
	d. scopy - copy x[i] to y[i] (i.e., set y[i] = x[i])
	e. sdot - return sum of x[i]*y[i] (i.e., dot product of x and y)
	f. snrm2 - return square root of sum of squares of x[i]
	g. sscal - compute x[i] = a*x[i]
	h. sswap - swap x[i] and y[i]

14. Double precision BLAS (basic linear algebra subroutines)
	a. idamax - return index of element with maximum absolute value
	b. dasum - return sum of absolute values
	c. daxpy - compute y[i] = a*x[i]+y[i]
	d. dcopy - copy x[i] to y[i] (i.e., set y[i] = x[i])
	e. ddot - return sum of x[i]*y[i] (i.e., dot product of x and y)
	f. dnrm2 - return square root of sum of squares of x[i]
	g. dscal - compute x[i] = a*x[i]
	h. dswap - swap x[i] and y[i]

15. Single precision LINPACK general matrix functions
	a. sgefa - LU factorization
	b. sgeco - LU factorization and condition number
	c. sgesl - Solve linear system Ax = b or A'x = b after LU factorization.

16. Double precision LINPACK general matrix functions
	a. dgefa - LU factorization
	b. dgeco - LU factorization and condition number
	c. dgesl - Solve linear system Ax = b or A'x = b after LU factorization.

17. Single precision LINPACK QR decomposition functions
	a. sqrdc - QR matrix decomposition
	b. sqrsl - Least squares solutions, etc. by QR
	c. sqrst - Under-determined or over-determined least squares

18. Other Linear System Solving routines
	a. stoepd - double precision symmetric Toeplitz linear systems
	b. stoepf - single precision symmetric Toeplitz linear systems
	c. vanded - double precision Vandermonde systems
	d. vandef - single precision Vandermonde systems
	e. tridif - single precision tridiagonal linear systems

19. Graphics Utility routines
	a. rfwtva - rasterize a float array as wiggle-trace-variable-area
	b. scaxis - compute a readable scale for use in plotting axes
	c. yclip  - clip a uniformly sampled function y(x)

20. Special Function routines
	a. airya - return approximation to the Airy function Ai(x)
	b. airyb - return approximation to the Airy function Bi(x)

21. Timing routines
	a. cpusec  - return cpu time (UNIX user time) in seconds
	b. wallsec - return elapsed time (wall clock time) in seconds
	c. cputime  - same as a using ANSI C built-ins
	d. walltime - same as b using ANSI C built-ins

22. Random Number Generating routines
	a. franuni - return a uniformly distributed random float
	b. sranuni - seed random number generator for uniform distribution
	c. frannor - return a normally distributed random float
	d. srannor - seed random number generator for normal distribution

23. Miscellaneous routines
	a. pp1d   - printer plot of a 1-dimensional array
	b. pplot1 - printer plot of a 1-dimensional array
	c. filestat  - determine type of file from file descriptor
 	d. printstat - print the filetype as a string

24. Dip filtering routines
	a. dpf.f  - 2d dip filtering 
------------------------------------------------------------------------

1. alloc.c - package of routines that allocate and free arrays

Contents:

	alloc1, alloc2, alloc3, alloc4 - generic 1D,2D,3D,4D allocation
	free1,  free2,  free3,  free4  - ... their free routines

	alloc1int,     alloc2int,     alloc3int     - 1D,2D,3D int arrays
	alloc1float,   alloc2float,   alloc3float   - ... float arrays
	alloc1double,  alloc2double,  alloc3double  - ... double arrays
	alloc1complex, alloc2complex, alloc3complex - ... complex arrays

	free1int,     free2int,     free3int     - int alloc free routines
	free1float,   free2float,   free3float   - ... float
	free1double,  free2double,  free3double  - ... double
	free1complex, free2complex, free3complex - ... complex

	realloc1        - reallocation routine for alloc1
	realloc1int     - ... for alloc1int
	realloc1float   - ... for alloc1float
	realloc1double  - ... for alloc1double
	realloc1complex - ... for alloc1complex

Notes:
	These functions simplify manipulation of multi-dimensional arrays
	in C.  These functions are useful because true multi-dimensional
	arrays in C cannot have variable dimensions (as in FORTRAN).
	For example, the following function IS NOT valid in C:
		void badFunc(a,n1,n2)
		float a[n2][n1];
		{
			a[n2-1][n1-1] = 1.0;
		}
	However, the following function IS valid in C:
		void goodFunc(a,n1,n2)
		float **a;
		{
			a[n2-1][n1-1] = 1.0;
		}
	Therefore, the functions defined below do not allocate true
	multi-dimensional arrays, as described in the C specification.
	Instead, they allocate and initialize pointers (and pointers to 
	pointers) so that, for example, a[i2][i1] behaves like a 2-D array.

	The array dimensions are numbered, which makes it easy to add 
	functions for arrays of higher dimensions.  In particular,
	the 1st dimension of length n1 is always the fastest dimension,
	the 2nd dimension of length n2 is the next fastest dimension,
	and so on.  Note that the 1st (fastest) dimension n1 is the 
	first argument to the allocation functions defined below, but 
	that the 1st dimension is the last subscript in a[i2][i1].
	(This is another important difference between C and Fortran.)

	The allocation of pointers to pointers implies that more storage
	is required than is necessary to hold a true multi-dimensional array.
	The fraction of the total storage allocated that is used to hold 
	pointers is approximately 1/(n1+1).  This extra storage is unlikely
	to represent a significant waste for large n1.

	The functions defined below are significantly different from similar 
	functions described by Press et al, 1988, Numerical Recipes in C.
	In particular, the functions defined below:
		(1) Allocate arrays of arbitrary size elements.
		(2) Allocate contiguous storage for arrays.
		(3) Return NULL if allocation fails (just like malloc).
		(4) Do not provide arbitrary lower and upper bounds for arrays.

	Contiguous storage enables an allocated multi-dimensional array to
	be passed to a C function that expects a one-dimensional array.
	For example, to allocate and zero an n1 by n2 two-dimensional array
	of floats, one could use
		a = alloc2(n1,n2,sizeof(float));
		zeroFloatArray(n1*n2,a[0]);
	where zeroFloatArray is a function defined as
		void zeroFloatArray(int n, float *a)
		{
			int i;
			for (i=0; i<n; i++)
				a[i] = 0.0;
		}

	Internal error handling and arbitrary array bounds, if desired,
	should be implemented in functions that call the functions defined 
	below, with the understanding that these enhancements may limit 
	portability.  (See .../par/lib/ealloc.c for routines that gracefully
	abort if the allocation fails.)


Prototypes:

void *alloc1 (size_t n1, size_t size);
void *realloc1 (void *v, size_t n1, size_t size);
void **alloc2 (size_t n1, size_t n2, size_t size);
void ***alloc3 (size_t n1, size_t n2, size_t n3, size_t size);
void ****alloc4 (size_t n1, size_t n2, size_t n3, size_t n4, size_t size);
void free1 (void *p);
void free2 (void **p);
void free3 (void ***p);
void free4 (void ****p);
int *alloc1int (size_t n1);
int *realloc1int (int *v, size_t n1);
int **alloc2int (size_t n1, size_t n2);
int ***alloc3int (size_t n1, size_t n2, size_t n3);
float *alloc1float (size_t n1);
float *realloc1float (float *v, size_t n1);
float **alloc2float (size_t n1, size_t n2);
float ***alloc3float (size_t n1, size_t n2, size_t n3);
double *alloc1double (size_t n1);
double *realloc1double (double *v, size_t n1);
double **alloc2double (size_t n1, size_t n2);
double ***alloc3double (size_t n1, size_t n2, size_t n3);
complex *alloc1complex (size_t n1);
complex *realloc1complex (complex *v, size_t n1);
complex **alloc2complex (size_t n1, size_t n2);
complex ***alloc3complex (size_t n1, size_t n2, size_t n3);
void free1int (int *p);
void free2int (int **p);
void free3int (int ***p);
void free1float (float *p);
void free2float (float **p);
void free3float (float ***p);
void free1double (double *p);
void free2double (double **p);
void free3double (double ***p);
void free1complex (complex *p);
void free2complex (complex **p);
void free3complex (complex ***p);

------------------------------------------------------------------------
 2. complex.c - complex number manipulation

Contents:

	cadd  - add two complex numbers
	csub  - subtract two complex numbers
	cmul  - multiply two complex numbers
	cdiv  - divide two complex numbers
	cmplx - make a complex number from two real numbers
	conjg - complex conjugate of a complex number 
	cneg  - negate a complex number
	cinv  - invert a complex number
	csqrt - complex square root of a complex number
	cexp  - complex exponential of a complex number
	crmul - multiply a complex number by a real number 
	fcabs	real magnitude of a complex number

Notes:
	From Numerical Recipes in C, Press, et al. (cinv was added)

Prototypes:

complex cadd (complex a, complex b);
complex csub (complex a, complex b);
complex cmul (complex a, complex b);
complex cdiv (complex a, complex b);
float fcabs (complex z);
complex cmplx (float re, float im);
complex cneg (complex z);
complex cinv (complex z);
complex csqrt (complex z);
complex cexp (complex z);
complex crmul (complex a, float x);

------------------------------------------------------------------------
 3. bm.c - handle big matrices in virtual memory

Contents:

	bmalloc  allocate a big matrix
	bmfree   free a big matrix 
	bmread   read a vector from a big matrix
	bmwrite  write a vector to a big matrix

Notes:

	The bm functions provide access to a big 2-dimensional matrix
	along either the 1st or 2nd dimensions.  Although the matrix
	must be small enough to fit in virtual memory, it can be larger
	than available real memory.  These functions provide equally
	efficient (or equally inefficient) access to vectors in a big
	matrix along either the 1st or 2nd dimensions.

	For example, the following algorithm will efficiently transpose an
	n1 by n2 array of (n1*n2) floats stored in a file:

	bmstate *bm;
	float *v;
	bm = (bmstate*)bmalloc(sizeof(float),n1,n2);
	for (i2=0; i2<n2; i2++) {
		(read n1 floats from input file into array v);
		bmwrite(bm,1,0,i2,n1,(char*)v);
	}
	for (i1=0; i1<n1; i1++) {
		bmread(bm,2,i1,0,n1,(char*)v);
		(write n2 floats in array v to output file);
	}
	bmfree(bm);

Prototypes:

bmstate *bmalloc (int nbpe, int n1, int n2);
void bmfree (bmstate *state);
void bmread (bmstate *state, int dir, int k1, int k2, int n, void *v);
void bmwrite (bmstate *state, int dir, int k1, int k2, int n, void *v);

------------------------------------------------------------------------
4. Prime Factor FFT routines:

	a. npfa - find next allowed fft size for pfacc

Notes:
	npfa returns the smallest valid n not less than nmin for prime
	factor complex-to-complex fft.

	npfa will be composed of mutually prime factors from
	the set {2,3,4,5,7,8,9,11,13,16}.  Because npfa cannot exceed
	720720 = 5*7*9*11*13*16, npfa = 720720 is returned if nmin
	exceeds 720720.


Prototype and parameters:

int npfa (int nmin);
nmin	i lower bound on returned value

-----------------------------------------
	b. npfao  - find optimal fft size for pfacc

Notes:
	npfao returns the optimal n between nmin and nmax
	for prime factor complex-to-complex fft.

	npfao will be composed of mutually prime factors from the set
	{2,3,4,5,7,8,9,11,13,16}.  Because npfao cannot exceed
	720720 = 5*7*9*11*13*16, npfao = 720720 is returned if nmin
	exceeds 720720.  If nmin does not exceed 720720, then the
	returned npfao will not be less than nmin.  The optimal npfao