mkdiff.c

seismic software,very useful

/* Copyright (c) Colorado School of Mines, 1990.
/* All rights reserved.                       */

/*
FUNCTION:  make discrete Taylor series approximation to n'th derivative

PARAMETERS:
n			i order of desired derivative (n>=0 && n<=m-l is required)
a			i fractional distance from integer sampling index (see notes)
h			i sampling interval
l			i sampling index of first coefficient (see notes below)
m			i sampling index of last coefficient (see notes below)
d			o array of coefficients for n'th order differentiator

NOTES:
The abscissae x of a smpled function f(x) can always be expressed as
x = (j+a)*h, where j is an integer, a is a fraction, and h is the
sampling interval.  To approximate the n'th order derivative fn(x)
of the sampled function f(x) at x = (j+a)*h, use the m-l+1 coefficients
in the output array d[] as follows:

	fn(x) = d[0]*f(j-l) + d[1]*f(j-l-1) +...+ d[m-l]*f(j-m)

i.e., convolve the coefficients in d with the samples in f.

m-l+1 (the number of coefficients) must not be greater than the
NCMAX parameter specified below.

For best approximations,
when n is even, use a = 0.0, l = -m
when n is odd, use a = 0.5, l = -m-1

AUTHOR:  Dave Hale, Colorado School of Mines, 06/02/89
*/

#include "cwp.h"

#define NCMAX 100

void mkdiff (int n, float a, float h, int l, int m, float d[])
{
	int nc,i;
	double fn,hn,v[NCMAX],b[NCMAX];

	/* determine number of coefficients */
	nc = m-l+1;

	/* build Vandermonde system of equations */
	for (i=0; i<nc; i++) {
		v[i] = -(l+i)-a;
		b[i] = 0.0;
	}
	for (i=1,fn=hn=1.0; i<=n; i++) {
		fn *= i;
		hn *= h;
	}
	b[n] = fn/hn;

	/* compute coefficients in double precision */
	vanded(nc,v,b,b);

	/* copy coefficients to output array */
	for (i=0; i<nc; i++)
		d[i] = b[i];
}