sutihaledmo.c

su 的源代码库

	float A;		/*	Hale's A			*/
	float B;		/*	Bleistein's scale factor	*/
	float phase;		/*	phase				*/
	float dw;		/*	angular freq inc		*/
	float tn;		/*	input event time after NMO	*/
	float eps=1.0e-20;	/*	zero offset*wavenumber thresh	*/
	float scale;		/*	forward/inverse fft scale factor*/
	int iw;			/*	freq index			*/
	int it;			/*	tn time index			*/
	int iwtest;		/*	test for pos or neg freq	*/
	

	hk = h*k;
	if(hk<eps) return;
	dw = 2*PI/(ntfft*dt);
	scale=1./ntfft;
	iwtest=ntfft/2+1;
	pwork[0]=cmplx(0.,0.);
	for(iw=1;iw<ntfft;iw++) {
		if(iw<iwtest)
			w=iw*dw;
		else
			w=(iw-ntfft)*dw;
		csum=cmplx(0.,0.);
		if( (k*v) < (2*fabs(w)) ) {
			hkpw=hk/(w);
			for(it=1;it<nt;it++) {
				tn=it*dt;
				hkpwt=hkpw/tn;
				A = sqrt(1. + hkpwt*hkpwt);
				B = sqrt(1. + 2*hkpwt*hkpwt);
				phase = w*tn*A;
				cfac = cmplx( B*cos(phase)/A, B*sin(phase)/A ); 		
				csum = cadd( csum, cmul(p[it],cfac) );
			}
		}
		pwork[iw]=csum;
	}
	pfacc(-1,ntfft,pwork);
	for(it=0;it<nt;it++) p[it]=crmul(pwork[it],scale);
	return;
}
void bleisteindmo2(float h,float k,complex *p,complex *pwork,
		int nt,int ntfft,float dt,float v)
/*  bleistein cos*cos weighted amplitude fk dmo

h	half offset
k	wavenumber
p	input NMO corrected data in wavenumber domain, p(tn,k,h)
pwork	complex array of lenght ntfft for work space
nt	number of time samples
ntfft	lenght of temporal fft
dt	temporal sample rate in seconds
v	velocity in m/sec (used to limit aperature of DMO response)

*/
{
	complex csum;		/* complex variable for sum 		*/
	complex cfac;	
	float w;		/* 	angular freq 			*/
	float hk;		/* 	h*k				*/
	float hkpw;		/*	h*k/w				*/
	float hkpwt;		/*	h*k/(w*tn)			*/
	float hkpwt2;		/*	square of above			*/
	float hpv;		/* 	2*h/v				*/
	float hpvt2;		/*	[(2*h)/(v*tn)]^2		*/
	float A;		/*	Hale's A			*/
	float B;		/*	Bleistein's scale factor	*/
	float phase;		/*	phase				*/
	float dw;		/*	angular freq inc		*/
	float tn;		/*	input event time after NMO	*/
	float scale;		/*	forward/inverse fft scale factor*/
	int iw;			/*	freq index			*/
	int it;			/*	tn time index			*/
	int iwtest;		/*	test for pos or neg freq	*/
	

	hk = h*k;
	dw = 2*PI/(ntfft*dt);
	scale=1./ntfft;
	iwtest=ntfft/2+1;
	pwork[0]=cmplx(0.,0.);
	hpv=2.*h/v;
	for(iw=1;iw<ntfft;iw++) {
		if(iw<iwtest)
			w=iw*dw;
		else
			w=(iw-ntfft)*dw;
		csum=cmplx(0.,0.);
		if( (k*v) < (2*fabs(w)) ) {
			hkpw=hk/(w);
			for(it=1;it<nt;it++) {
				tn=it*dt;
				hkpwt=hkpw/tn;
				hkpwt2=hkpwt*hkpwt;
				hpvt2=hpv/tn;
				hpvt2=hpvt2*hpvt2;
				A = sqrt(1. + hkpwt2);
				B = (1.+2*hkpwt2)*(1+hkpwt2)/(1+hpvt2);
				phase = w*tn*A;
				cfac = cmplx( B*cos(phase)/A, B*sin(phase)/A ); 		
				csum = cadd( csum, cmul(p[it],cfac) );
			}
		}
		pwork[iw]=csum;
	}
	pfacc(-1,ntfft,pwork);
	for(it=0;it<nt;it++) p[it]=crmul(pwork[it],scale);
	return;
}
void zhangdmo(float h,float k,complex *p,complex *pwork,
		int nt,int ntfft,float dt,float v)
/*  zhang fk dmo

h	half offset
k	wavenumber
p	input NMO corrected data in wavenumber domain, p(tn,k,h)
pwork	complex array of lenght ntfft for work space
nt	number of time samples
ntfft	lenght of temporal fft
dt	temporal sample rate in seconds
v	velocity in m/sec (used to limit aperature of DMO response)

*/
{
	complex csum;		/* complex variable for sum 		*/
	complex cfac;	
	float w;		/* 	angular freq 			*/
	float hk;		/* 	h*k				*/
	float hkpw;		/*	h*k/w				*/
	float hkpwt;		/*	h*k/(w*tn)			*/
	float hkpwt2;		/* 	square of above			*/
	float A;		/*	Hale's A			*/
	float B;		/*	Zhang's scale factor		*/
	float phase;		/*	phase				*/
	float dw;		/*	angular freq inc		*/
	float tn;		/*	input event time after NMO	*/
	float eps=1.0e-20;	/*	zero offset*wavenumber thresh	*/
	float scale;		/*	forward/inverse fft scale factor*/
	int iw;			/*	freq index			*/
	int it;			/*	tn time index			*/
	int iwtest;		/*	test for pos or neg freq	*/
	
	hk = h*k;
	if(hk<eps) return;
	dw = 2*PI/(ntfft*dt);
	scale=1./ntfft;
	iwtest=ntfft/2+1;
	pwork[0]=cmplx(0.,0.);
	for(iw=1;iw<ntfft;iw++) {
		if(iw<iwtest)
			w=iw*dw;
		else
			w=(iw-ntfft)*dw;
		csum=cmplx(0.,0.);
		if( (k*v) < (2*fabs(w)) ) {
			hkpw=hk/(w);
			for(it=1;it<nt;it++) {
				tn=it*dt;
				hkpwt=hkpw/tn;
				hkpwt2=hkpwt*hkpwt;
				A = sqrt(1. + hkpwt2);
				B = (1.+2*hkpwt2)/(1+hkpwt2);
				phase = w*tn*A;
				cfac = cmplx( B*cos(phase)/A, B*sin(phase)/A ); 		
				csum = cadd( csum, cmul(p[it],cfac) );
			}
		}
		pwork[iw]=csum;
	}
	pfacc(-1,ntfft,pwork);
	for(it=0;it<nt;it++) p[it]=crmul(pwork[it],scale);
	return;
}
void tsvankindmo(float h,float k,complex *p,complex *pwork,
		int nt,int ntfft,float dt,int np, float dp, float *vnmo)
/*  tsvankin fk dmo

h	half offset
k	wavenumber
p	input NMO corrected data in wavenumber domain, p(tn,k,h)
pwork	complex array of lenght ntfft for work space
nt	number of time samples
ntfft	lenght of temporal fft
dt	temporal sample rate in seconds
np	number of slowness values for which vnmo has been tabled
dp	slowness increment for vnmo table
vnmo	nmo velocity in m/sec as a function of p 

*/
{
	complex csum;		/* complex variable for sum 		*/
	complex cfac;	
	float w;		/* 	angular freq 			*/
	float hk;		/* 	h*k				*/
	double v;		/*	vnmo(p)				*/
	double v0;		/*	vnmo(0)				*/
	float A;		/*	Hale's A			*/
	float phase;		/*	phase				*/
	float dw;		/*	angular freq inc		*/
	float tn;		/*	input event time after NMO	*/
	float tn2;
	float eps=1.0e-20;	/*	zero offset*wavenumber thresh	*/
	float scale;		/*	forward/inverse fft scale factor*/
	float pp;		/*	ray parameter			*/
	double fac;
	float rp;
	float wgt;
	int iw;			/*	freq index			*/
	int it;			/*	tn time index			*/
	int it1;		/* 	first stable time index		*/
	int iwtest;		/*	test for pos or neg freq	*/
	int ip;
	int npm;

	hk = h*k;
	if(hk<eps) return;
	npm=np-2;
	dw = 2*PI/(ntfft*dt);
	scale=1./ntfft;
	iwtest=ntfft/2+1;
	pwork[0]=cmplx(0.,0.);
	for(iw=1;iw<ntfft;iw++) {

		if(iw<iwtest)
			w=iw*dw;
		else
			w=(iw-ntfft)*dw;

		csum=cmplx(0.,0.);
		pp=fabs(k/(2.*w));
		rp=pp/dp;
		ip=rp;
		if(ip<npm) {
			wgt=rp-ip;
			v=(1.-wgt)*vnmo[ip] + wgt*vnmo[ip+1];
			v0=vnmo[0];
			fac=4*h*h*(1./(v0*v0) - 1./(v*v) );
			if(fac<0.) {
				it1=2+sqrt(fabs(fac))/dt;
			}else{
				it1=1;
			}
			for(it=it1;it<nt;it++) {
				tn=it*dt;
				tn2=tn*tn;
				A=sqrt(1. + fac/tn2);
				phase = w*tn*A;
				cfac = cmplx( cos(phase)/A, sin(phase)/A ); 		
				csum = cadd( csum, cmul(p[it],cfac) );	
			}
		}
		pwork[iw]=csum;
	}
	pfacc(-1,ntfft,pwork);
	for(it=0;it<nt;it++) p[it]=crmul(pwork[it],scale);
	return;
}
void weakdmo(float h,float k,complex *p,complex *pwork,
		int nt,int ntfft,float dt,float e,float d, float v)
/*  tsvankin fk dmo weak approximation

h	half offset
k	wavenumber
p	input NMO corrected data in wavenumber domain, p(tn,k,h)
pwork	complex array of lenght ntfft for work space
nt	number of time samples
ntfft	lenght of temporal fft
dt	temporal sample rate in seconds
e	epsilon
d	delta
v	v0 velocity
*/
{
	complex csum;		/* complex variable for sum 		*/
	complex cfac;	
	float w;		/* 	angular freq 			*/
	float hk;		/* 	h*k				*/
	float A;		/*	Hale's A			*/
	float phase;		/*	phase				*/
	float dw;		/*	angular freq inc		*/
	float tn;		/*	input event time after NMO	*/
	float tn2;
	float eps=1.0e-20;	/*	zero offset*wavenumber thresh	*/
	float scale;		/*	forward/inverse fft scale factor*/
	float pp;		/*	ray parameter			*/
	float y;
	double fac;
	int iw;			/*	freq index			*/
	int it;			/*	tn time index			*/
	int it1;		/* 	first stable time index		*/
	int iwtest;		/*	test for pos or neg freq	*/

	hk = h*k;
	if(hk<eps) return;
	dw = 2*PI/(ntfft*dt);
	scale=1./ntfft;
	iwtest=ntfft/2+1;
	pwork[0]=cmplx(0.,0.);
	for(iw=1;iw<ntfft;iw++) {

		if(iw<iwtest)
			w=iw*dw;
		else
			w=(iw-ntfft)*dw;

		csum=cmplx(0.,0.);
		pp=fabs(k/(2.*w));
		if(pp<(1/v)) {
			y=pp*pp*v*v;
			fac=4*h*h*pp*pp*(1+2*(e-d)*(4*y*y-9*y+6));
			if(fac<0.) {
				it1=2+sqrt(fabs(fac))/dt;
			}else{
				it1=1;
			}
			for(it=it1;it<nt;it++) {
				tn=it*dt;
				tn2=tn*tn;
				A=sqrt(1. + fac/tn2);
				phase = w*tn*A;
				cfac = cmplx( cos(phase)/A, sin(phase)/A ); 		
				csum = cadd( csum, cmul(p[it],cfac) );	
			}
		}
		pwork[iw]=csum;
	}
	pfacc(-1,ntfft,pwork);
	for(it=0;it<nt;it++) p[it]=crmul(pwork[it],scale);
	return;
}