sutifowler.c

su 的源代码库

	float *ttn;
	float *qtn;
	float *fold;
	float *big;
	float factor;
	float t;
	char *ccrt;

	nt=vnda->N[0];
	ntout=vnd->N[1]/nv;
	rt=(float *)VNDemalloc(ntout*sizeof(float),"sutifowler:cvstack: rt");
	buf=(float *)VNDemalloc(MAX(nt,ntout)*sizeof(float),"sutifowler:cvstack: buf");
	ttn=(float *)VNDemalloc(ntout*sizeof(float),"sutifowler:cvstack: ttn");
	qtn=(float *)VNDemalloc(ntout*sizeof(float),"sutifowler:cvstack: qtn");
	fold=(float *)VNDemalloc(ntout*sizeof(float),"sutifowler:cvstack: fold");
	big=(float *)VNDemalloc(noff*sizeof(float),"sutifowler:cvstack: big");
	ccrt=(char *)rt;
	
	/* check for completely dead traces so can skip them */
	for(ioff=0;ioff<noff;ioff++) {
		big[ioff]=0.;
		V2Dr0(vnda,ioff,(char *)buf,3);
		for(it=0;it<nt;it++)
			big[ioff]=MAX( big[ioff], fabs(buf[it]) );	
	}

	for(iv=0;iv<nv;iv++) {
		for(it=0;it<ntout;it++) rt[it]=0.;
		for(it=0;it<ntout;it++) fold[it]=0.;
		for(ioff=0;ioff<noff;ioff++) {
		    if(big[ioff]>0.) {
			V2Dr0(vnda,ioff,(char *)buf,4);

			factor=off[ioff]*off[ioff]*p2[iv];
			itmute=mute[ioff]/dtout;
			
			for(it=itmute;it<ntout;it++) {
				t=it*dtout;
				ttn[it]=sqrt(t*t+factor);
				}

			/* do nmo via 8-point sinc interpolation */
			ints8r(nt,dt,0.,buf,0.0,0.0,
				ntout-itmute,&ttn[itmute],&qtn[itmute]);
			
			/* apply linear ramp to taper mute */
			for (it=itmute; it<(itmute+lmute) && it<ntout; ++it)
				qtn[it] *= (float)(it-itmute+1)/(float)lmute;

			/* sum NMO corrected trace to stacked trace  */
			for(it=itmute;it<ntout;it++) rt[it]+=qtn[it];

			/* count fold information */
			for(it=itmute;it<ntout;it++) fold[it]+=1.;
		    }
		}
		for(it=0;it<ntout;it++) {
			if(fold[it]==0.) fold[it]=1.;
		}
		for(it=0;it<ntout;it++) rt[it]/=fold[it]; 
		key[0]=icmp;
		key[1]=0;
		VNDrw('w',0,vnd,1,key,0,ccrt,iv*ntout,1,ntout,1,"cv stacking");
	}

	VNDfree(rt,"cvstack: freeing rt");
	VNDfree(buf,"cvstack: freeing buf");
	VNDfree(ttn,"cvstack: freeing ttn");
	VNDfree(qtn,"cvstack: freeing qtn");
	VNDfree(fold,"cvstack: freeing fold");
	VNDfree(big,"cvstack: freeing big");
	if(icmp==20*(icmp/20)){
		fprintf(fpl,
		"sutifowler: completed stacking velocity analysis for cmp %d\n",icmp);
	}
	return;
}
VND *ptabledmo(int nv, float *v, float etamin, float deta, int neta,
	float d, float vsvp, int np, float dp, float dp2, char *file)
/*************************************************************************
	returns a VND table where dimension 0 is along v axis
	and dimension 1 is along p axis providing 

	1./(vnmo(p)*vnmo(p)*dp2)

	which will be the desired index in the stacking velocity
	analysis table for a given output velocity index and p
**************************************************************************
int nv		number of input velocities
float v[]	array of input velocities
float etamin	Alkhalifah's minimum eta value
float deta	Alkhalifah's eta increment
int neta	number of eta values
float d		Thomsen's delta
float vsvp	vs/vp ratio
int np		number of output slownesses
float dp	output slowness increment
float dp2	increment in slowness squared used for stacking
		velocity analysis
char *file	root name for temporary file to be created
***************************************************************************
Author: John Anderson (visitor to CSM from Mobil) Spring 1993
***************************************************************************/
{
	VND *vndvel;
	float a,b,dangle,angle,theta,dscale;
	float vel[5];
	int nangle,iangle,ip;
	long iv,newnv,newnp,ieta;
	float *pin;
	float *vin;
	float *pout;
	float *vout;
	float e;
	char *fname;

	newnv=nv*neta;
	newnp=np;
	fname=VNDtempname(file);
	dscale=1./sqrt(1+2*d);
	vndvel=V2Dop(2,125000,sizeof(float),
		fname,newnv,newnp);
	VNDfree(fname,"ptabledmo: freeing fname");

	nangle=10001;
	dangle = 89./(nangle-1);
	pin = (float *)VNDemalloc(nangle*sizeof(float),
			"sutifowler:ptable: allocating pin");
	vin = (float *)VNDemalloc(nangle*sizeof(float),
			"sutifowler:ptable: allocating vin");
	pout = (float *)VNDemalloc(np*sizeof(float),
			"sutifowler:ptable: allocating pout");
	vout = (float *)VNDemalloc(np*sizeof(float),
			"sutifowler:ptable: allocating vout");
	for(ip=0;ip<np;ip++) pout[ip]=ip*dp;
	for(ieta=0;ieta<neta;ieta++) {
	    e= (etamin + ieta*deta)*(1+2*d)+d;
	    for(iv=0;iv<nv;iv++) {
		a=v[iv]*dscale;	/* vertical p-wave phase velocity */
		b=vsvp*a;	/* vertical s-wave phase velocity */
		for(iangle=0;iangle<nangle;iangle++) {
			angle=iangle*dangle;
			theta=angle*PI/180.;
			vget(a,b,e,d,theta,vel);
			pin[iangle]=vel[4];
			vin[iangle]=vel[3];
		}
		intlin(nangle,pin,vin,vin[0],vin[nangle-1],np,pout,vout);
		for(ip=0;ip<np;ip++)
			vout[ip]=1./(vout[ip]*vout[ip]*dp2);
		V2Dw1(vndvel,iv+ieta*nv,(char *)vout,101);
	    }
	}
	VNDfree(pin,"ptabledmo: freeing pin");
	VNDfree(vin,"ptabledmo: freeing vin");
	VNDfree(pout,"ptabledmo: freeing pout");
	VNDfree(vout,"ptabledmo: freeing vout");
	return (vndvel);
}
VND *ptablemig(int nv, float *v, float etamin, float deta, int neta,
		float d, float vsvp, int np, char *file)
/***************************************************************************
	returns a VND table where dimension 0 is along k/(2*wmig) axis
	and dimension 1 is along v axis (note a=v/sqrt(1+2d)).

	The first element in dimension 0 will be dp, the mig slowness
	increment which corresponds to increments in tan(theta)/a. 
	Thereafter, the values will be

	factor = v/( a*cos(theta) )

	where p corresponds to (i-1)*dp and a is the vertical p-wave
	phase velocity.  A Stolt migration can use this factor to
	compute		

		wdmo = wmig*factor

	using k/(2*wmig) = tan(theta)/a to find the appropriate location
	in the table.
*****************************************************************************
int nv		number of input velocities
float v[]	array of input velocities
float etamin	Alkhalifah's minimum eta value
float deta	Alkhalifah's eta increment
int neta	number of eta values
float d	 Thomsen's delta
float vsvp      vs/vp ratio
int np		number of output slownesses
char *file	root name for temporary file to be created
***************************************************************************
Author: John Anderson (visitor to CSM from Mobil) Spring 1993
***************************************************************************/
{
	VND *vndvel;
	float a,b,dangle,angle,theta,pmax,dp,dscale;
	float vel[5];
	int nangle,iangle,ip;
	long iv,newnv,newnp,ieta;
	float *pin;
	float *vin;
	float *pout;
	float *vout;
	float e;
	char *fname;

	newnv=nv*neta;
	newnp=np;
	fname=VNDtempname(file);
	vndvel=V2Dop(2,125000,sizeof(float),
		fname,newnp,newnv);
	VNDfree(fname,"pmigtable: freeing fname");

	dscale=1./sqrt(1+2*d);
	nangle=10001;
	dangle = 89./(nangle-1);
	pin = (float *)VNDemalloc(nangle*sizeof(float),
			"sutifowler:ptable: allocating pin");
	vin = (float *)VNDemalloc(nangle*sizeof(float),
			"sutifowler:ptable: allocating vin");
	pout = (float *)VNDemalloc(np*sizeof(float),
			"sutifowler:ptable: allocating pout");
	vout = (float *)VNDemalloc(np*sizeof(float),
			"sutifowler:ptable: allocating vout");
	for(ieta=0;ieta<neta;ieta++) {
	    e= (etamin + ieta*deta)*(1+2*d)+d;
	    for(iv=0;iv<nv;iv++) {
		a=v[iv]*dscale;	/* vertical p-wave phase velocity */
		b=vsvp*a;	/* vertical s-wave phase velocity */
		for(iangle=0;iangle<nangle;iangle++) {
			angle=iangle*dangle;
			theta=angle*PI/180.;
			vget(a,b,e,d,theta,vel);
			pin[iangle]=tan(theta)/a;
			vin[iangle]=vel[0]/( a*cos(theta) );
		}
		pmax=pin[nangle-1];
		dp=pmax/(np-2);
		vout[0]=dp;
		for(ip=0;ip<(np-1);ip++) pout[ip]=ip*dp;
		intlin(nangle,pin,vin,vin[0],vin[iangle-1],np-1,pout,&vout[1]);
		V2Dw0(vndvel,iv+ieta*nv,(char *)vout,101);
	    }
	}
	VNDfree(pin,"ptablemig: freeing pin");
	VNDfree(vin,"ptablemig: freeing vin");
	VNDfree(pout,"ptablemig: freeing pout");
	VNDfree(vout,"ptablemig: freeing vout");
	return (vndvel);
}

static void vget( float a, float b, float e, float d, float theta, float *vel)
/***************************************************************************

This routine returns phase and NMO velocity information given
an input angle and Thomsen's parameters for a TI medium.
****************************************************************************
input parameters:
----------------
a      = alpha = vertical p-wave phase velocity = vrms/sqrt(1+2*d)   
b      = beta = vertical shear wave phase velocity
e      = epsilon = anisotropy factor for horizontally propagating p waves 
d      = delta = 0.5*(e + ds/(1-(b*b)/(a*a)))      
theta  = angle in radians

returned parameters:
-------------------
vel[0] = p-wave phase velocity
vel[1] = first derivative of p-wave phase velocity with respect to theta
vel[2] = second derivative of p-wave phase velocity with respect to theta
vel[3] = NMO velocity
vel[4] = ray parameter p = sin(theta)/vphase
***************************************************************************
Author: John Anderson (visitor to CSM from Mobil) Spring 1993
***************************************************************************/
{
	float p, psi, sint, cost, sin2t, cos2t, tant, eps;
	float vnmo,v,dv,d2v,gamma,dgamma,d2gamma,sqgamma;
	eps=1.0e-7;
	sint=sin(theta);
	cost=cos(theta);
	sin2t=sint*sint;
	cos2t=cost*cost;
	psi=1.-(b*b)/(a*a);
	gamma=0.25*psi*psi+(2*d-e)*psi*sin2t*cos2t+(psi+e)*e*sin2t*sin2t;
	dgamma=(2*d-e)*psi*2*sint*cost*cos2t +
		( 4*(psi+e)*e - 2*(2*d-e)*psi )*sint*sin2t*cost;
	d2gamma=(2*d-e)*psi*2*(cos2t*cos2t-3.*cos2t*sin2t) +
		( 4*(psi+e)*e - 2*(2*d-e)*psi )*(3.*cos2t*sin2t-sin2t*sin2t);	
	sqgamma=sqrt(gamma);
	v=a*sqrt(1.+e*sin2t-0.5*psi+sqgamma);
	dv=0.5*a*a*(2*e*sint*cost + 0.5*dgamma/sqgamma)/v;
	d2v=0.5*a*a*(2*e*(cos2t-sin2t)+0.5*d2gamma/sqgamma
		-0.25*dgamma/(sqgamma*gamma) )/v -
		0.5*a*a*dv*(2*e*sint*cost + 0.5*dgamma/sqgamma)/(v*v);
	if(cost<eps) cost=eps;
	tant=sint/cost;
	vnmo = ( v*sqrt( 1 + d2v/v) )/( cost * (1-tant*dv/v) );
	p=sint/v;
	vel[0]=v;
	vel[1]=dv;
	vel[2]=d2v;
	vel[3]=vnmo;
	vel[4]=p;
	return;
}

static void taper (int lxtaper, int lbtaper,
	int nx, int ix, int nt, float *trace)
/*****************************************************************************
Taper traces near left and right sides of trace window
******************************************************************************
Input:
lxtaper		length (in traces) of side taper
lbtaper		length (in samples) of bottom taper
nx		number of traces in window
ix		index of this trace (0 <= ix <= nx-1)
nt		number of time samples
trace		array[nt] containing trace

Output:
trace		array[nt] containing trace, tapered if within lxtaper of side
*****************************************************************************
Based on code by David Hale
*****************************************************************************/
{
	int it;
	float xtaper;

	/* if near left side */
	if (ix<lxtaper) {
		xtaper = 0.54+0.46*cos(PI*(lxtaper-ix)/lxtaper);
	
	/* else if near right side */
	} else if (ix>=nx-lxtaper) {
		xtaper = 0.54+0.46*cos(PI*(lxtaper+ix+1-nx)/lxtaper);
	
	/* else x tapering is unnecessary */
	} else {
		xtaper = 1.0;
	}

	/* if x tapering is necessary, apply it */
	if (xtaper!=1.0)
		for (it=0; it<nt; ++it)
			trace[it] *= xtaper;
	
	/* if requested, apply t tapering */
	for (it=MAX(0,nt-lbtaper); it<nt; ++it)
		trace[it] *= (0.54+0.46*cos(PI*(lbtaper+it+1-nt)/lbtaper));
}