sudmofk.c

来自「seismic software,very useful」· C语言 代码 · 共 529 行 · 第 1/2 页

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

	return EXIT_SUCCESS;
}

void dmooff (float offset, int nxfft, float dx, 
	int nt, float dt, float ft,
	float vmin, int ntvdmo, float *tdmo, float *vdmo,
	float **ptx)
/*****************************************************************************
perform dmo in w-k domain for one offset
******************************************************************************
Input:
offset		source receiver offset
nxfft		number of midpoints, zero-padded for prime factor fft
dx		midpoint sampling interval
nt		number of time samples
dt		time sampling interval
ft		first time
vmin		minimum velocity
ntvdmo		number of tdmo,vdmo pairs
tdmo		times at which rms velocities vdmo are specified
vdmo		rms velocities specified at times in tdmo
ptx		array[nxfft][nt] containing p(t,x) zero-padded for fft

Output:
ptx		array[nxfft][nt] containing p(t,x) after dmo
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 11/04/90
*****************************************************************************/
{
	int itmin,nu,nufft,nw,nk,ix,iu,iw,ik,it,iwn,iwnyq;
	float dw,dk,tcon,wwscl,scale,scales,amp,phase,fr,fi,pwr,pwi,fftscl,
		du,fu,w,k,
		*uoft,*tofu;
	complex czero=cmplx(0.0,0.0),**ptk,*pu,*pw;

	/* determine minimum time of first non-zero sample */
	for (ix=0,itmin=nt; ix<nxfft; ++ix) {
		for (it=0; it<itmin && ptx[ix][it]==0.0; ++it);
		itmin = it;
	}
	
	/* if all zeros, simply return */
	if (itmin>=nt) return;
	
	/* make stretch and compress functions t(u) and u(t) */
	maketu(offset,itmin,vmin,ntvdmo,tdmo,vdmo,nt,dt,ft,
		&uoft,&nu,&du,&fu,&tofu,&tcon);
	
	/* determine frequency sampling */
	nufft = npfao(2*nu,3*nu);
	nw = nufft;
	dw = 2.0*PI/(nufft*du);
	
	/* allocate workspace */
	pu = pw = ealloc1complex(nufft);
	
	/* determine wavenumber sampling and set pointers to complex p(t,k) */
	nk = nxfft/2+1;
	dk = 2.0*PI/(nxfft*dx);
	ptk = (complex**)ealloc1(nk,sizeof(complex*));
	for (ik=0; ik<nk; ++ik)
		ptk[ik] = (complex*)ptx[0]+ik*nt;
	
	/* compute fft scale factor */
	fftscl = 1.0/(nufft*nxfft);
	
	/* Fourier transform p(t,x) to p(t,k) */
	pfa2rc(-1,2,nt,nxfft,ptx[0],ptk[0]);

	/* loop over wavenumbers */
	for (ik=0,k=0.0; ik<nk; ++ik,k+=dk) {
		
		/* stretch p(t;k) to p(u) */
		ints8c(nt,dt,ft,ptk[ik],czero,czero,nu,tofu,pu);
		
		/* pad with zeros and Fourier transform p(u) to p(w) */
		for (iu=nu; iu<nufft; ++iu)
			pu[iu].r = pu[iu].i = 0.0;
		pfacc(1,nufft,pu);
		
		/* constants independent of w */
		iwnyq = nw/2;
		wwscl = pow(k*0.5*offset/tcon,2.0);
		
		/* zero p(w=0) (dc should be zero anyway) */
		pw[0].r = pw[0].i = 0.0;
		
		/* do dmo for frequencies between zero and Nyquist */
		for (iw=1,iwn=nw-1,w=dw; iw<iwnyq; ++iw,--iwn,w+=dw) {
			scales = 1.0+wwscl/(w*w);
			scale = sqrt(scales);
			amp = fftscl/scale;
			phase = w*tcon*(scale-1.0);
			fr = amp*cos(phase);
			fi = amp*sin(phase);
			pwr = pw[iw].r;
			pwi = pw[iw].i;
			pw[iw].r = pwr*fr-pwi*fi;
			pw[iw].i = pwr*fi+pwi*fr;
			pwr = pw[iwn].r;
			pwi = pw[iwn].i;
			pw[iwn].r = pwr*fr+pwi*fi;
			pw[iwn].i = pwi*fr-pwr*fi;
		}
	
		/* do dmo for the Nyquist frequency */
		w = iwnyq*dw;
		scales = 1.0+wwscl/(w*w);
		scale = sqrt(scales);
		amp = fftscl/scale;
		phase = w*tcon*(scale-1.0);
		fr = amp*cos(phase);
		fi = amp*sin(phase);
		pwr = pw[iwnyq].r;
		pwi = pw[iwnyq].i;
		pw[iwnyq].r = pwr*fr-pwi*fi;
		pw[iwnyq].i = pwr*fi+pwi*fr;
		
		/* Fourier transform p(w) to p(u) */
		pfacc(-1,nufft,pu);
		
		/* compress p(u) to p(t;k) */
		ints8c(nu,du,fu,pu,czero,czero,nt,uoft,ptk[ik]);
	}
	
	/* Fourier transform p(t,k) to p(t,x) */
	pfa2cr(1,2,nt,nxfft,ptk[0],ptx[0]);
	
	/* free workspace */
	free1float(tofu);
	free1float(uoft);
	free1complex(pu);
	free1(ptk);
}

void maketu (float offset, int itmin, float vmin, 
	int ntvdmo, float *tdmo, float *vdmo,
	int nt, float dt, float ft, float **uoftp,
	int *nup, float *dup, float *fup, float **tofup, float *tconp)
/*****************************************************************************
make stretch and compress functions t(u) and u(t)
******************************************************************************
Input:
offset		source receiver offset
itmin		index of minimum first non-zero sample for this offset
vmin		minimum velocity
ntvdmo		number of tdmo,vdmo pairs
tdmo		times at which rms velocities vdmo are specified
vdmo		rms velocities specified at times in tdmo
nt		number of time samples
dt		time sampling interval
ft		first time

Output:
uoftp		array[nt] of u(t)
nup		number of u (stretched t) samples
dup		u sampling interval
fup		first u
tofup		array[nu] of t(u)
tconp		time constant relating t(u) and u(t)
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 11/04/90
*****************************************************************************/
{
	int it,iu,numax,nu;
	float tmin,dumin,et,eu,t1,t2,vdmolo,vdmohi,
		v2,v22,v44,gamma,
		v2m,v22m,v44m,gammam,t,dv2,vi2,vi4,v24,
		*uoft,du,fu,*tofu,tcon;
	
	/* determine maximum number of u */
	numax = 500+log((float)nt)*(float)(nt-1);
		
	/* allocate space for u(t) */
	uoft = ealloc1float(nt);
	
	/* determine t1 and t2, rounded to nearest sampled times */
	tmin = ft+itmin*dt;
	et = ft+(nt-1)*dt;
	t1 = MIN(et,MAX(ft+dt,tmin));
	if (offset!=0.0)
		t2 = MAX(t1,1.0/(1.0/et+2.0*vmin*vmin*dt/(offset*offset)));
	else
		t2 = t1;
	t1 = ft+NINT(t1/dt)*dt;
	t2 = ft+NINT(t2/dt)*dt;
	
	/* compute u(t) */
	vdmolo = vdmo[0];
	vdmohi = vdmo[ntvdmo-1];
	intlin(ntvdmo,tdmo,vdmo,vdmolo,vdmohi,1,&ft,&v2);
	v22 = v2*v2;
	v44 = v22*v22;
	gamma = 1.0;
	for (it=0,t=ft; it<nt; ++it,t+=dt) {
		v2m = v2;
		v22m = v22;
		v44m = v44;
		gammam = gamma;
		if (t>0.0) {
			intlin(ntvdmo,tdmo,vdmo,vdmolo,vdmohi,1,&t,&v2);
			v22 = v2*v2;
			vi2 = (t*v22-(t-dt)*v22m)/dt;
			vi4 = vi2*vi2;
			v44 = (dt*vi4+(t-dt)*v44m)/t;
		} else {
			v2 = v2m;
			v22 = v22m;
			v44 = v44m;
		}
		dv2 = (v2-v2m)/dt;
		v24 = v22*v22;
		gamma = 1.5*v44/v24-t*dv2/v2-0.5;
		if (t<=t1) {
			uoft[it] = t-t1;
		} else if (t>t1 && t<=t2) {
			du = t1*(gamma*log(t/(t-0.5*dt)) -
				gammam*log((t-dt)/(t-0.5*dt)));
			dumin = 0.1*dt*t1/t;
			uoft[it] = uoft[it-1]+MAX(dumin,du);
		} else if (t>t2) {
			uoft[it] = 2.0*uoft[it-1]-uoft[it-2];
		}
	}
	
	/* determine minimum u(t)-u(t-dt) */
	dumin = uoft[1]-uoft[0];
	for (it=1; it<nt; ++it)
		dumin = MIN(dumin,uoft[it]-uoft[it-1]);
	
	/* determine u sampling for t(u) to avoid aliasing */
	fu = 0.0;
	eu = uoft[nt-1];
	du = dumin;
	nu = 1+NINT((eu-fu)/du);
	if (nu>numax) {
		nu = numax;
		du = (eu-fu)/(nu-1);
	}
	
	/* allocate space for t(u) */
	tofu = ealloc1float(nu);
	
	/* compute t(u) by inverse linear interpolation of u(t) */
	yxtoxy(nt,dt,ft,uoft,nu,du,fu,ft,et,tofu);
	
	/* set time constant */
	tcon = t1;
	
	/* set returned values */
	*uoftp = uoft;
	*nup = nu;
	*dup = du;
	*fup = fu;
	*tofup = tofu;
	*tconp = tcon;
}