sudmovz.c

su 的源代码库

Output:
q	array[nx][nt] of output common offset section (DMO corrected)
****************************************************************************
Author: Craig Artley, Colorado School of Mines, 09/28/91
****************************************************************************/
{
	int it,ntn,itn,itnmin,itnmax,np,npdmo,ip,it0max,*it0min,*nt0;
	float h,lt,tn,ftn,dtn,dp2,tmin,**t0table,**tntable;

	/* compute minimum time based on stretch mute */
	tmin = tmute(vrms,nt,dt,ft,0.5*offset,smute);
	
	/* compute p sampling required to avoid aliasing */
	lt = (nt-1)*dt;
	if (tmin>=lt)
		err("tmin=%f >= lt=%f:  no data survive the mute!\n",tmin,lt);

	dtn = 10.0*dt;
	ftn = tmin;
	ntn = NINT((lt-ftn)/dtn)+1;
	dp2 = tmin/(fmax*0.25*offset*offset);
	dp2 *= speed;
	np = NINT(lp*lp/dp2)+1;

	/* allocate space for dmo mapping tables */
	it0min = ealloc1int(np);
	nt0 = ealloc1int(np);
	t0table = ealloc2float(ntn,np);
	tntable = ealloc2float(nt,np);

	/* zero the tables */
	for (ip=0; ip<np; ++ip)
		for (itn=0; itn<ntn; ++itn)
			t0table[ip][itn] = 0.0;
	for (ip=0; ip<np; ++ip)
		for (it=0; it<nt; ++it)
			tntable[ip][it] = 0.0;
	
	/* compute (unsigned) half-offset */
	h = ABS(offset/2.0);

	/* loop over nmo times tn */
	for (itn=0,tn=ftn; itn<ntn; ++itn,tn+=dtn) {

		/* compute dmo mapping t0(p) for this tn */
		dmomap(x,a,tau,nttab,dt,ft,nptab,dptab,np,dp2,fp,lp,
			itn,tn,h,vrms,verbose,t0table);
	}

	/* inverse interpolation to get tn(t0,p) */
	for (ip=0; ip<np; ++ip) {
		for (itn=0; itn<ntn; ++itn)
			if (t0table[ip][itn]>0.0) break;
		itnmin = itn;
		
		it0min[ip] = NINT(t0table[ip][itnmin]/dt);
		if (it0min[ip]>=nt) it0min[ip] = nt-1;

		for (itn=ntn-1; itn>=0; --itn)
			if (t0table[ip][itn]>0.0) break;
		itnmax = itn;
		if (itnmax <= itnmin) break;

		it0max = NINT(t0table[ip][itnmax]/dt);
		if (it0max>=nt) it0max = nt-1;

		nt0[ip] = it0max-it0min[ip]+1;

		/* check that table is monotonically increasing */
		for (itn=itnmin+1; itn<itnmax; ++itn)
			if (t0table[ip][itn]<=t0table[ip][itn-1] && verbose)
				fprintf(stderr,
					"t0table not monotonically "
					"increasing---extrapolating\n");

		yxtoxy(itnmax-itnmin+1,dtn,ftn+itnmin*dtn,t0table[ip]+itnmin,
			nt0[ip],dt,it0min[ip]*dt,
			ftn+itnmin*dtn,ftn+itnmax*dtn,tntable[ip]);

		/* extend tntable by linear extrapolation */
		for (it=nt0[ip]; it<nt; ++it) {
			tntable[ip][it] = 2.0*tntable[ip][it-1]
				-tntable[ip][it-2];
		}
		nt0[ip] = nt-it0min[ip];
	}

	/* reset the number of ray parameters for dmo */
	npdmo = MIN(ip,np);

	/* free space */
	free2float(t0table);

	/* dmo by dip decomposition */
	dmoapply(nx,dx,nt,dt,npdmo,dp2,h,tntable,it0min,nt0,tmin,fmax,q);

	/* free space */
	free1int(it0min);
	free1int(nt0);
	free2float(tntable);

}

void dmomap (float **x, float **a, float **tau, int nt, float dt, float ft,
	int nptab, float dptab, int np, float dp2, float fp, float lp,
	int itn, float tn, float h, float *vrms,
	int verbose, float **t0table)
/***************************************************************************
Compute DMO correction for one NMO time, tn, using Newton's method to solve
system of equations.
****************************************************************************
Input:
x	array[nptab][nt] of horizontal position of ray 
a	array[nptab][nt] of propagation angle along ray
tau	array[nptab][nt] of time depth of ray
nt	number of times in ray tables
dt	time sampling interval
ft	first time (must be zero for now)
nptab	number of ray parameters in ray tables
dptab	ray parameter sampling interval in ray tables
np	number of ray parameters for which to find DMO mapping
dp2	sampling interval in (ray parameter)^2 in DMO mapping table
fp	first ray parameter
lp	last ray parameter
itn	index of NMO time for which to compute DMO mapping
tn	NMO time for which to compute DMO mapping
h	source-receiver half offset
vrms	array[nt] of rms velocity as function of two-way vertical time
verbose	flag for diagnostic info (does nothing here)

Output:
t0table	array[np][ntn] containing DMO mapping t0(tn,p0)
	One call to this function fills in the itn'th row of this table.
****************************************************************************
Author: Craig Artley, Colorado School of Mines, 09/28/91
****************************************************************************/
{
	int ip,i,neqns=5,niter=NITER,iter,*ipvt;
	float p2,vtn,tsg,p0,tautp,dtaudt,dtaudp,error,rcond,
		lt,tol=TOL,*vars,*eqns,*work,t0,x0,**jacob;

	/* compute last time in ray tables */
	lt = ft+(nt-1)*dt + 0*verbose;

	/* allocate space for system of equations */
	vars = ealloc1float(neqns);		/* variables X */
	eqns = ealloc1float(neqns);		/* F(X) */
	jacob = ealloc2float(neqns,neqns);	/* Jacobian dF_i/dx_j */
	work = ealloc1float(neqns);		/* work array for sgeco() */
	ipvt = ealloc1int(neqns);	/* pivot index array for sgeco() */

	/* calculate tsg from tn via Dix approximation */
	vtn = vrms[(int)(tn/dt)];
	tsg = sqrt(tn*tn+4*h*h/(vtn*vtn));

	/* get initial values for parameters */
	getX(tsg,h,fp,vtn,vars);

	for (ip=0,p2=fp; ip<np; ++ip,p2+=dp2) {

		/* compute p0 */
		p0 = sqrt(p2);

		/* use Newton's method to improve the guess */
		for (iter=0; iter<niter; ++iter) {
	
			/* calculate F(X) and J(X) */
			getFJ(x,a,tau,nt,nptab,dt,dptab,tsg,h,p0,
				vars,eqns,jacob);

			/* compute total error */
			for (i=0,error=0.0; i<neqns; ++i)
				error += ABS(eqns[i]);

			/* if iteration has converged, break */
			if (error<=tol) break;

			/* solve system */
			sgeco(jacob,neqns,ipvt,&rcond,work);
			sgesl(jacob,neqns,ipvt,eqns,0);

			/* update X */
			if (rcond!=0.0) {
				for (i=0; i<neqns; ++i)
					vars[i] -= eqns[i];
			} else {
				/* system is ill-conditioned, break */
				break;
			}

			/* clip ps and pg values to pmax if necessary */
			if (ABS(vars[1])>lp) {
				vars[1] = SGN(vars[1])*lp;
			}
			if (ABS(vars[2])>lp) {
				vars[2] = SGN(vars[2])*lp;
			}

			/* clip tg to 0<tg<tsg */
			if (vars[3]<0.0) {
				vars[3] = 0.0;
			}
			if (vars[3]>tsg) {
				vars[3] = tsg;
			}

			/* clip t0 values to 0<t0<lt */
			if (vars[4]<0.0) {
				vars[4] = 0.0;
			}
			if (vars[4]>lt) {
				vars[4] = lt;
			}

		}

		/* vars[] now contains the solution to the system */

		/* compute depth of reflection point */
		intder(tau,nt,dt,ft,nptab,dptab,fp,vars[4],ABS(p0),
			&tautp,&dtaudt,&dtaudp);

		/* if reflector is at surface, break */
		if (tautp<dt) {
			break;
		}

		/* if system is ill-conditioned, break */
		if (rcond==0.0) {
			break;
		}

		/* now have t0(x0) for this tsg, p0, h -- store it */
		x0 = vars[0];
		t0 = vars[4];
		t0table[ip][itn] = t0+p0*x0;

	}

	/* free space */
	free1float(vars);
	free1float(eqns);
	free2float(jacob);
	free1float(work);
	free1int(ipvt);
}

void dmoapply (int nx, float dx, int nt, float dt, int np, float dp2,
	float h, float **tntable, int *it0min, int *nt0,
	float tmin, float fmax, float **qx)
/*****************************************************************************
Apply DMO correction to common offset section using previously computed
mapping tables.
******************************************************************************
Input:
nx		number of traces in common offset section
dx		trace sampling interval
nt		number of time samples per trace
dt		time sampling interval
np		number of ray parameters in DMO mapping tables
dp2		ray parameter (squared) sampling interval
h		source-receiver half-offset
tntable		array[np][nt0] of DMO mapping table tn(t0,p0)
it0min		array[np] of index of first entry in DMO mapping table 
nt0		array[np] of number of t0's in DMO mapping table
fmax		maximum frequency of data (frequencies above fmax are muted)
qx		array[nx][nt] of input NMO corrected data

Output:
qx		array[nx][nt] of output DMO corrected data
******************************************************************************
Author: Craig Artley, Colorado School of Mines, 08/29/91
******************************************************************************/
{
	int ix,nxfft,it,itmin,ik,nk;
	float k,dk,fk,ft,scale;
	complex **qk;

	/* determine wavenumber sampling (for real to complex FFT) */
	nxfft = npfar(nx);
	nk = nxfft/2+1;
	dk = 2.0*PI/(nxfft*dx);
	fk = 0.0;

	ft = 0.0;

	/* pointers to complex q(t,k) */
	qk = (complex**)ealloc1(nk,sizeof(complex*));
	for (ik=0; ik<nk; ++ik)
		qk[ik] = (complex*)qx[0]+ik*nt;
	
	/* Fourier transform x to k */
	pfa2rc(-1,2,nt,nxfft,qx[0],qk[0]);

	/* loop over wavenumber */
	for (ik=0,k=fk; ik<nk; ++ik,k+=dk) {

		/* apply DMO correction */
		jakub1k(k,h,np,dp2,nt,dt,ft,tntable,
			it0min,nt0,fmax,qk[ik]);
	}

	/* Fourier transform k to x (including FFT scaling and mute) */
	pfa2cr(1,2,nt,nxfft,qk[0],qx[0]);
	scale = 1.0/nxfft;
	itmin = MIN(nt,NINT(tmin/dt));
	for (ix=0; ix<nx; ++ix) {
		for (it=0; it<itmin; ++it)
			qx[ix][it] = 0.0;
		for (it=itmin; it<nt; ++it)
			qx[ix][it] *= scale;
	}

	/* free space */
	free1(qk);
}

void jakub1k (float k, float h, int np, float dp2, int nt, float dt, float ft,
	float **tntable, int *it0min, int *nt0, float fmax, complex *qq)
/*****************************************************************************
Jakubowicz's dip decomposition DMO for one wavenumber
******************************************************************************
Input:
k		wavenumber
h		half-offset (ignored)
np		number of ray paramaters
dp2		sampling interval in (ray parameter)^2
nt		number of time samples
dt		time sampling interval
ft		first time sample
qq		complex array[nt] of input NMO corrected data 

Output:
qq		complex array[nt] of output DMO corrected data
******************************************************************************
Author: Craig Artley, Colorado School of Mines, 08/29/91
******************************************************************************/
{
	int ntfft,nw,it,iw,iwlm,iwhm,iwlp,iwhp,ip;
	float p,ph,pl,ph2,dw,fw,wl,wh,scale;
	complex czero=cmplx(0.0,0.0),*qn,*q0;

	/* determine frequency sampling */
	ntfft = npfa(nt);
	nw = ntfft;
	dw = 2.0*PI/(ntfft*dt);
	fw = -PI/dt + 0.0*h;
	
	/* allocate workspace */
	qn = ealloc1complex(ntfft);
	q0 = ealloc1complex(ntfft);

	/* zero w0 domain accumulator */
	for (iw=0; iw<nw; ++iw) {
		qn[iw] = czero;
		q0[iw] = czero;
	}

	/* loop over zero-offset ray parameter */
	for (ip=0,p=pl=0.0,ph2=dp2; ip<np; ++ip,ph2+=dp2) {

		/* update ray parameter */
		ph = sqrt(ph2);

		/* determine limits of w0 for this wavenumber and slope */
		wl = k/ph;
		if (pl==0.0) {
			wh = 2.0*PI*fmax;
		} else {
			wh = MIN(k/pl,2.0*PI*fmax);
		}

		/* compute limits of sum for -p */
		iwlm = (int) ((-wl-fw)/dw);
		iwhm = (int) ((-wh-fw)/dw);
		if (iwhm<0) iwhm = 0;
		++iwhm;

		/* compute limits of sum for +p */
		iwlp = (int) ((wl-fw)/dw);
		iwhp = (int) ((wh-fw)/dw);
		if (iwhp>=nw) iwhp = nw-1;
		++iwlp;

		/* if this range of frequencies contributes to the sum */
		if ((iwhm<=iwlm) || (iwlp<=iwhp)) {
		
			/* interpolate, padding with zeros */
			for (it=0; it<it0min[ip]; ++it)
				qn[it] = czero;
			ints8c(nt,dt,ft,qq,czero,czero,nt0[ip],tntable[ip],
				qn+it0min[ip]);
			for (it=it0min[ip]+nt0[ip]; it<ntfft; ++it)
				qn[it] = czero;

			/* Fourier transform t to w, with w centered */
			for (it=1; it<ntfft; it+=2) {