sudmotx.c

su 的源代码库

		if (cdp>=cdpmin && cdp<=cdpmax) {
		
			/* save trace header and update number of traces */
			efwrite(&tr,HDRBYTES,1,headerfp);
			ntrace++;
			
			/* determine output traces potentially modified
			   by input */
			icdp = cdp-cdpmin;
			jcdplo = MAX(0,icdp-0.5*ABS(offset/dxcdp));
			jcdphi = MIN(ncdp-1,icdp+0.5*ABS(offset/dxcdp));
			
			/* loop over potentially modified output traces */
			for (jcdp=jcdplo; jcdp<=jcdphi; ++jcdp) {
		
				/* do dmo for one output trace */
				dmotx(ns,ts,as,offset,(jcdp-icdp)*dxcdp,dxcdp,
				      0,nt,dt,ft,p,q[jcdp]);
			}

			/* remember offset */
			oldoffset = offset;
		}

		/* get next trace (if there is one) */
		if (!gettr(&tr)) gottrace = 0;
		
	} while (!done);

	/* clean up */
	efclose(headerfp);
	if (istmpdir) eremove(headerfile);
	return(CWP_Exit());
}

void maketa (float dx, float dt, float offmax, float tmute, float vrms,
	int nsmax, int *ns, float *ts, float *as)
/*****************************************************************************
make time shifts and amplitudes for non-aliased (t,x) domain DMO
******************************************************************************
Input:
dx		midpoint sampling interval (see notes)
dt		time sampling interval (see notes)
offmax		maximum offset (see notes)
tmute		mute time at maximum offset (see notes)
vrms		RMS velocity at mute time (see notes)
nsmax		maximum number of time shifts

Output:
ns		number of shifts
ts		array[ns] of time shifts (normalized by sampling interval dt)
as		array[ns] of amplitudes corresponding to time shifts
******************************************************************************
Notes:
dx, dt, offmax, tmute, and vrms must be greater than 0.0.

dx, dt, offmax, tmute, and vrms need not be specified precisely.
If these values are unknown, then one should overestimate dt and
offmax, and underestimate dx, tmute, and vrms.

The time shifts, ts, are computed to ensure that the DMO operator
is not aliased in time or space.

The number of time shifts, ns, is determined such that the steepest
reflection slope at the mute time (tmute) for the largest offset
is properly handled.  A typical value for ns is about ns=200, so
a reasonable value for nsmax is nsmax=400, just to be safe.  The
computation cost for (t,x) domain DMO is linearly proportional to ns.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 11/04/90
*****************************************************************************/
{
	int nsi,is,itaper,ltaper;
	float alpha,dalpha,dalphan,tsmax;
	
	/* ensure positive arguments */
	dx = ABS(dx);
	dt = ABS(dt);
	offmax = ABS(offmax);
	tmute = ABS(tmute);
	vrms = ABS(vrms);
	
	/* determine nominal increment in alpha = sqrt(ts) */
	dalphan = MIN(sqrt(2.0*(tmute/dt)*pow((dx/offmax),2.0)),1.0);
	
	/* determine maximum time shift for steepest reflection slope */
	tsmax = (tmute/dt)*(1.0-1.0/sqrt(1.0+pow(offmax/(vrms*tmute),2.0)));
	
	/* recursively compute time shifts and amplitudes */
	alpha = 0.0;
	dalpha = dalphan;
	for (nsi=0; alpha*alpha<=tsmax && nsi<nsmax; ++nsi) {
		ts[nsi] = alpha*alpha;
		as[nsi] = dalpha;
		alpha = alpha+dalpha;
		dalpha = sqrt(alpha*alpha+1.0)-alpha;
		if (dalpha>dalphan) dalpha = dalphan;
	}
	
	/* taper the amplitudes with a raised cosine */
	ltaper = nsi/3;
	for (is=nsi-ltaper,itaper=1; is<nsi; ++is,++itaper) 
		as[is] *= 0.54+0.46*cos(PI*itaper/ltaper);
	
	*ns = nsi;
}

void makeds (int ns, float *ts, float *as, int lds, int ifds, float *ds)
/*****************************************************************************
make shaping filter to correct DMO horizontal reflection response
******************************************************************************
Input:
ns		number of shifts
ts		array[ns] of time shifts (normalized by sampling interval)
as		array[ns] of amplitudes corresponding to time shifts
lds		length of shaping filter
ifds		index of first sample of shaping filter (see notes)

Output:
ds		array[lds] containing shaping filter
******************************************************************************
Notes:
Reasonable values for lds and ifds are lds=125 and ifds=-100.
The maximum permissible lds is the dimension of dd, di, and work below.

This function must be kept consistent with that used to perform DMO.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 11/04/90
*****************************************************************************/
#define LDMAX 400
#define LDSMAX 300
{
	int ld=LDMAX,ifd=1-ld,i,is;
	float one=1.0,d[LDMAX],dd[LDSMAX],di[LDSMAX],work[LDSMAX];
	
	/* compute dmo horizontal reflection response d(t) */
	for (i=0; i<ld; ++i)
		d[i] = 0.0;
	d[ld-1] = as[0];
	for (is=1; is<ns; ++is) {
		i = (int)ts[is];
		d[ld-i-1] += 2.0*as[is];
	}
	
	/* compute autocorrelation of d(t) */
	xcor(ld,ifd,d,ld,ifd,d,lds,0,dd);
	
	/* compute crosscorrelation of d(t) and desired impulse i(t) */
	xcor(ld,ifd,d,1,0,&one,lds,ifds,di);
	
	/* solve symmetric toeplitz system of equations for filter */
	stoepf(lds,dd,di,ds,work);
}

void dmotx (int ns, float *ts, float *as, float offset, float x, float dx,
	int itmute, int nt, float dt, float ft, float *p, float *q)
/*****************************************************************************
apply DMO in (t,x) domain for one input trace and one output trace
******************************************************************************
Input:
ns		number of shifts
ts		array[ns] of time shifts (normalized by sampling interval)
as		array[ns] of amplitudes corresponding to time shifts
offset		source-receiver offset of input trace
x		midpoint distance from input trace to output trace
dx		midpoint sampling interval
itmute		index of the first time sample that is not muted
nt		number of time samples
dt		time sampling interval
ft		first time sample
p		array[nt] containing input trace
q		array[nt] containing accumulated DMO output trace

Output:
q		array[nt] containing accumulated DMO output trace
******************************************************************************
Notes:
Time shifts and amplitudes (ts and as) must be computed by maketa().
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 11/04/90
*****************************************************************************/
{
	int it1,it2,its,it,is;
	float h,xb,ftodt,x1,x2,scale1,scale2,asi,tsi;
	
	/* compute half-offset */
	h = offset/2.0;
	
	/* round midpoint distance to nearest multiple of sampling interval */
	xb = ABS(NINT(x/dx)*dx);
	
	/* normalize the time of first sample */
	ftodt = ft/dt;
	
	/* if rounded midpoint distance is zero */
	if (xb==0.0) {
		
		/* compute midpoint distance at right edge of midpoint bin */
		x1 = xb+0.5*ABS(dx);
		
		/* compute scale factor used to determine index it1 below */
		scale1 = 1.0/(1.0-sqrt(MAX(0.0,1.0-x1*x1/MAX(h*h,x1*x1))));
		
		/* accumulate output for first (zero) time shift */
		asi = as[0];
		for (it=itmute; it<nt; ++it)
			q[it] += asi*p[it];
		
		/* loop over non-zero time shifts */
		for (is=1; is<ns; ++is) {
			tsi = ts[is];
			asi = as[is];
			
			/* round the normalized time shift to lower int */
			its = tsi;
			
			/* compute index of first time sample for this shift */
			it1 = scale1*tsi-ftodt;
			if (it1<itmute+its) it1 = itmute+its;
			
			/* if index is beyond last time sample, break */
			if (it1>=nt) break;
			
			/* accumulate shifted and weighted p(t) in q(t) */
			for (it=it1; it<nt; ++it)
				q[it-its] += 2.0*asi*p[it];
		}
	
	/* else if non-zero time shifts exist */
	} else if (ns>1) {
		
		/* compute midpoint distances at right and left edges of bin */
		x1 = xb+0.5*ABS(dx);
		x2 = xb-0.5*ABS(dx);
		
		/* if midpoint distance exceeds half-offset, return */
		if (x2*x2>h*h) return;
		
		/* compute scale factors for it1 and it2 below */
		scale1 = 1.0/(1.0-sqrt(MAX(0.0,1.0-x1*x1/MAX(h*h,x1*x1))));
		scale2 = 1.0/(1.0-sqrt(MAX(0.0,1.0-x2*x2/MAX(h*h,x2*x2))));
		
		/* loop over time shifts */
		for (is=1; is<ns; ++is) {
			tsi = ts[is];
			asi = as[is];
			
			/* round normalized time shift to lower int */
			its = tsi;
			
			/* compute indices of first and last samples */
			it1 = scale1*tsi-ftodt;
			if (it1<itmute+its) it1 = itmute+its;
			it2 = scale2*tsi-ftodt;
			if (it2>nt) it2 = nt;
			
			/* if first sample index is out of bounds, break */
			if (it1>=nt) break;
			
			/* compute contribution to q(t) for this shift */
			for (it=it1; it<it2; ++it)
				q[it-its] += asi*p[it];
		}
	}
}

/* for graceful interrupt termination */
static void closefiles(void)
{
	efclose(headerfp);
	eremove(headerfile);
	exit(EXIT_FAILURE);
}