intt8c.c

seismic software,very useful

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

/*
FUNCTION:  interpolation of a uniformly-sampled complex function y(x)
via a table of 8-coefficient interpolators

PARAMETERS:
ntable		i number of tabulated interpolation operators; ntable>=2
table		i array of tabulated 8-point interpolation operators
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:
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.

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

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

#include "cwp.h"

void intt8c (int ntable, float table[][8],
	int nxin, float dxin, float fxin, complex yin[],
	complex yinl, complex yinr, int nxout, float xout[], complex yout[])
{
	int ioutb,ixout,ixoutn,kyin,ktable,itable,nxinm8;
	float xoutb,xoutf,xouts,xoutn,fntablem1,frac,
		yinlr,yinli,yinrr,yinri,yinir,yinii,sumr,sumi,
		*table00,*pyin,*ptable;
	complex *yin0;

	/* compute constants */
	ioutb = -3-8;
	xoutf = fxin;
	xouts = 1.0/dxin;
	xoutb = 8.0-xoutf*xouts;
	fntablem1 = (float)(ntable-1);
	nxinm8 = nxin-8;
	yin0 = &yin[0];
	table00 = &table[0][0];
	yinlr = yinl.r;  yinli = yinl.i;
	yinrr = yinr.r;  yinri = yinr.i;

	/* loop over output samples */
	for (ixout=0; ixout<nxout; ixout++) {

		/* determine pointers into table and yin */
		xoutn = xoutb+xout[ixout]*xouts;
		ixoutn = (int)xoutn;
		kyin = ioutb+ixoutn;
		pyin = (float*)(yin0+kyin);
		frac = xoutn-(float)ixoutn;
		ktable = frac>=0.0?frac*fntablem1+0.5:(frac+1.0)*fntablem1-0.5;
		ptable = table00+ktable*8;
		
		/* if totally within input array, use fast method */
		if (kyin>=0 && kyin<=nxinm8) {
			yout[ixout].r = 
				pyin[0]*ptable[0]+
				pyin[2]*ptable[1]+
				pyin[4]*ptable[2]+
				pyin[6]*ptable[3]+
				pyin[8]*ptable[4]+
				pyin[10]*ptable[5]+
				pyin[12]*ptable[6]+
				pyin[14]*ptable[7];
			yout[ixout].i = 
				pyin[1]*ptable[0]+
				pyin[3]*ptable[1]+
				pyin[5]*ptable[2]+
				pyin[7]*ptable[3]+
				pyin[9]*ptable[4]+
				pyin[11]*ptable[5]+
				pyin[13]*ptable[6]+
				pyin[15]*ptable[7];
		
		/* else handle end effects with care */
		} else {

			sumr = sumi = 0.0;
			for (itable=0; itable<8; itable++,kyin++) {
				if (kyin<0) {
					yinir = yinlr;
					yinii = yinli;
				} else if (kyin>=nxin) {
					yinir = yinrr;
					yinii = yinri;
				} else {
					yinir = yin[kyin].r;
					yinii = yin[kyin].i;
				}
				sumr += yinir*(*ptable);
				sumi += yinii*(*ptable++);
			}
			yout[ixout].r = sumr;
			yout[ixout].i = sumi;
		}
	}
}