taup.c

su 的源代码库

				/* maximum number that npfa can handle */
                                /* (see function npfa) */
        int it,iff; 		/* loop counters */
        int nfh,nfhp;  		/* half number of frequency coefficients */
	int ntfft;		/* time samples to pad */
	int nph=npoints/2;
        float f, df;            /* frequency and frequency sampling interval */
        complex *cx;            /* array of frequencies */

	/* oddness of npoints */
	if (2 * nph == npoints)   err ("npoints must be odd!\n");

	/* compare filter length with number of time samples */
	if (nph > nt)   
		err ("filter length larger than number of time samples!");

	/* compute padding factor */
	if (nt > maxnpfa)   
		err ("number of time samples too large for npfa!");
	ntfft = npfa(nt);

	/* allocate working space */
	cx = alloc1complex(ntfft);

	/* define constants */
	nfh = ntfft/2;
	nfhp = nfh+1;
	df = 1.0/(2*dt*nfh);

	/* compute filter coefficients */
	cx[0] = cmplx (0.0, 0.0);
	for (iff = 1, f = df; iff < nfhp; f = ++iff * df)
		cx[iff] = cx[ntfft-iff] = cmplx (f, 0.0);
  
	/* inverse Fourier transform from f to t */
	pfacc(-1,ntfft,cx);

	/* normalized output filter coefficients */
	rho[nph] = 1.0;
	for (it = 0; it < nph; it++)
		rho[nph-it-1] = rho[nph+it+1] = cx[it+1].r / cx[0].r;
  
	/* clean up */
	free1complex(cx);
}

/******************************************************************************
        Compute global forward slant stack in F-X domain
                        via radon transform

******************************************************************************/
void forward_p_transform(int nx, int nt, float dt, float pmax, float pmin,
        float dp, float depthref, float f1, float f2, float freq1, float freq2,
	int lagc, int lent, int lenx, int xopt, int ninterp, int nwin,
	float prewhite, float interoff, float offref, int igopt, float dx,
	float fx, float pmula, float pmulb, float **in_traces,
	float **out_traces)
/******************************************************************************
Input:
nx              number of input traces
nt              number of intput time samples
dt              sample rate in seconds
dx		offset sampling interval (distance between traces)
fx		first offset 
igopt           =1 parabolic trransform g(x)=offset**2
                =2 Foster/Mosher pseudo hyperbolic transform
                        g(x)=sqrt(depth**2+offset**2)
                =3 linear tau-p g(x)=offset
                =4 abs linear taup g(x)=abs(offset)
offref          reference maximum offset to which maximum and minimum moveout
                times are associated
interoff        intercept offset to which tau-p times are associated
                (usually zero)
pmax            maximum moveout in ms on reference offset
pmin            minimum moveout in ms on reference offset
dp              moveout sampling interval (ms/m)
depthref        reference depth for Foster/Mosher hyperbolic transform
f1              high-end  frequency before taper off
f2              high-end frequency
prewhite        prewhitening factor in percent (usually between 0.01 and 0.1)
nwin            number of windows to use through the mute zone (ignored)

Parameters with good suggested values:
freq1           low-end frequency for picking (usually 3 Hz)
freq1           high-end frequency for picking (usually 20 Hz)
lagc            length of AGC operator for picking (usually 400 ms)
lent            length of time smoother in samples for picking (usually 5)
lenx            length of space smoother in samples for picking (usually 1)
xopt            =1 use differences for spatial derivatives (works with
                irregular spacing)
                =0 use FFT derivative for spatial derivatives (more accurate
                but requires regular spacing and at least 16 input traces
                will switch to differences automatically is this is not met)
in_traces       2-D array[ntfft] of input t-x traces

Output:
out_traces      2-D array of output tau-p traces

offsets are not read from the headers but computed as offset[ix]=fx+ix*dx
*******************************************************************************
Credits:
	Adapted by Gabriel Alvarez (1995) from suradon.c written by
	John Anderson (1993)
******************************************************************************/
{
        int ix,ip,it,iw;                /* loop counters */
        int dummy; 	                /* dummy variable */
        int ntfft;                      /* length of time fft  */
        int nw;                         /* number of frequencies */
        int nxinterp;                   /* number of traces with interpol */
        long nmax;                       /* maximum number of traces */
        int np;                         /* number of sopes for taup transform*/
	int ipa,ipb;
	int ltaper;
        float w;                        /* frequency */
        float dw,df;                    /* frequency sampling intervals */
        float wa;                       /* some sort of frequencies ??? */
        float fac;                      /* scaling factor */
        float d,rsum;			/* auxiliary variables */
        float *rtr;			/* more auxiliary vectors */
        float *rt;			/* vector to hold trace */
        float *g,*offset; 		/* arrays of offsets */
        float *xin;  	                /* auxiliaty arrays */
        complex czero;                  /* complex zero number */
        complex *crt,*ccrt,*r,*rhs;     /* auxiliaty complex arrays */
        complex *wrk1,*wrk2,*wrk3,*wrk4;/* more complex arrays */
        float **interpolated_tr;        /* array of interpolated traces */
        char *fname;                    /* character pointer to file name */
        VND *vndresult;                 /* 2-D array of output (result) data */
        VND *vndinterp;                 /* 2-D array of input (interpol) data */

        /* define variables */
        ntfft=npfar(nt);
        fac = 1000.*gofx(igopt,offref,interoff,depthref);
        pmin /=fac;
        pmax /=fac;
        dp /=fac;
	pmula /=fac;
	pmulb /=fac;
	ipa = (pmula-pmin)/dp;
	ipb = (pmulb-pmin)/dp;
	if (ipa<0) ipa=0;
	ltaper=7;
	np =1+(pmax-pmin)/dp;
	nxinterp = (1+ninterp)*(nx-1)+1;
	nmax = MAX(nxinterp, np);
	nmax = MAX(nmax, ntfft+4);

	fprintf(stderr,"computing forward radon transform: pmin=%g pmax=%g"
		" dp=%g np=%d\n",pmin,pmax,dp,np);

	/* allocate 1-D VND arrays */
	offset =(float *)VNDemalloc(nx*sizeof(float),"fwdslant_taup:offset");
	rtr = (float *)VNDemalloc(nmax*sizeof(float),"fwdslant_taup:rtr");
	xin = (float *)VNDemalloc(nx*sizeof(float),"fwdslant_taup:trace");
	g = (float *)VNDemalloc(nxinterp*sizeof(float),"fwdslant_taup:gg");

	/* allocate 2-D VND arrays */
	fname = VNDtempname("radontem");
	vndinterp = V2Dop(2,2000000,sizeof(float),fname,nt,nxinterp);
	VNDfree(fname,"fwdslant_taup:fname 1");	
	fname = VNDtempname("radontmp");
	nmax = MAX(nxinterp,np);		/* max number of horiz samples*/
	vndresult = V2Dop(2,1000000,sizeof(float),fname,ntfft+2,nmax);
	VNDfree(fname,"fwdslant_taup:fname 2");

	/* compute offsets and offset function g(x) */
	for (ix=0; ix<nx; ix++) {

		/* compute offsets */
		offset[ix] = fx+ix*dx;	
	}

	/* interpolate traces */
	interpolated_tr=alloc2float(nt, nxinterp);
	ga_xinterpolate (in_traces, interpolated_tr, ninterp, nt, nx, freq1,
		freq2, lagc, lent, lenx, xopt, dt, 0);
	d = 1.0/(1+ninterp);

	/* copy interpolated traces to a VND array */
	for (ix=0; ix<nxinterp; ix++) {
		rtr[ix] = ix*d;	
		V2Dw0(vndinterp,ix,(char *)interpolated_tr[ix],1010);
	}
	free2float(interpolated_tr);
	
	/* interpolate offsets */
	for (ix=0; ix<nx; ix++)
		xin[ix] = ix;
	intlin (nx, xin, offset, offset[0], offset[nx-1], nxinterp, rtr, g);

	/* compute offset function depending on the type of transform */
	for (ix=0; ix<nxinterp; ix++) {
		g[ix] = gofx(igopt, g[ix], interoff, depthref);
	}
	
	/* define some variables */
	fac=1./ntfft;
	nw=1+ntfft/2;
	df=1./(ntfft*dt);
	dw=2.*PI*df;
	nmax=MAX(vndresult->N[0],vndresult->N[1]);
	czero.r=czero.i=0.;

	/* allocate working 1D complex space */
	crt=(complex *)VNDemalloc(nmax*sizeof(complex),
		"forward_transform:crt");
	ccrt=(complex *)VNDemalloc(MAX(2*np,vndresult->N[1])*
		sizeof(complex),"forward_transform:ccrt");
	r=(complex *)VNDemalloc(np*sizeof(complex),
		"forward_p_transform:r");
	rhs=(complex *)VNDemalloc(np*sizeof(complex),
		"forward_p_transform:rhs");
	wrk1=(complex *)VNDemalloc(np*sizeof(complex),
		"forward_p_transform:wrk1");
	wrk2=(complex *)VNDemalloc(np*sizeof(complex),
		"forward_p_transform:wrk2");
	wrk3=(complex *)VNDemalloc(np*sizeof(complex),
		"forward_p_transform:wrk3");
	wrk4=(complex *)VNDemalloc(np*sizeof(complex),
		"forward_p_transform:wrk4");

	/* pointers to complex arrays */
	rt=(float *)crt;

	/* do forward time to frequency fft */
	for(ix=0;ix<nxinterp;ix++) {
		V2Dr0(vndinterp,ix,(char *)rt,201);	/* read input data */
		for(it=0;it<nt;it++) rt[it]*=fac;	/* apply fft scale */
		for(it=nt;it<ntfft;it++) rt[it]=0.;	/* zero padding */
		pfarc(1,ntfft,rt,crt);
		V2Dw0(vndresult,ix,(char *)crt,202);
	}

	VNDr2c(vndresult); 		/* real to complex conversion */

	/* do radon transform, frequency by frequency, 
	for multiple spatial windows */

	for(iw=0;iw<nw;iw++) {
		wa=freqweight(iw,df,f1,f2);	/* compute frequency weights */
		if(wa>0.) {
		    	w=iw*dw;
		    	V2Dr1(vndresult,iw,(char *)crt,203);
		    	if(wa<1.)
         	       	   for(ix=0;ix<nxinterp;ix++) crt[ix]=crmul(crt[ix],wa);

			/* compute right hand side vector B+ and top row of
				Toeplitz matrix */
			compute_rhs(w,nxinterp,g,crt,np,pmin,dp,rhs);
			compute_r(w,nxinterp,g,np,dp,r);
			r[0].r *= (1.+prewhite);	/* apply prewhitening */
			for(rsum=0.,ip=1;ip<np;ip++) 
				rsum +=sqrt(r[ip].r*r[ip].r + r[ip].i*r[ip].i);
			rsum=rsum/r[0].r;  

			if (rsum>1.+np/5) { 

				/* conjugate gradient solution of Toeplitz eqn*/
				dummy=ctoephcg(np/7,np,r,&ccrt[0],rhs,
					wrk1,wrk2,wrk3,wrk4);
			} else {

				/* complex Hermitian solver of Toeplitz eqns */
				dummy=ctoep(np,r,&ccrt[0],rhs,wrk1,wrk2);
			}

			dummy=nwin|dummy;
		} else {
		    for(ip=0;ip<np;ip++) ccrt[ip]=czero;		
		}
		V2Dw1(vndresult,iw,(char *)ccrt,204);
	}

	/* do fourier transform from frequency to tau */
	for(ip=0; ip<np; ip++) {
		V2Dr0(vndresult,ip,(char *)crt,205);
		pfacr(-1,ntfft,crt,rt);

		/* do tau-p mute */	
		taupmute (ip,ipa,ipb,ntfft,ltaper,rt);	
	
		/* copy output data */
		for (it=0; it<nt; it++)
			out_traces[ip][it]=rt[it];
	}

	/* free allocated space */
	VNDcl(vndresult,1);
	VNDcl(vndinterp,1);
	VNDfree(crt,"forward_p_transform: crt");
	VNDfree(ccrt,"forward_p_transform: ccrt");
	VNDfree(r,"forward_p_transform: r");
	VNDfree(rhs,"forward_p_transform: rhs");
	VNDfree(wrk1,"forward_p_transform: wrk1");
	VNDfree(wrk2,"forward_p_transform: wrk2");
	VNDfree(wrk3,"forward_p_transform: wrk3");
	VNDfree(wrk4,"forward_p_transform: wrk4");	
	VNDfree(offset,"fwdslant_taup: offset");
	VNDfree(rtr,"fwdslant_taup: rtr");
	VNDfree(xin,"fwdslant_taup: xin");
	VNDfree(g,"fwdslant_taup: g");
	return;
}

/******************************************************************************

        Compute global inverse slant stack in frequency domain

******************************************************************************/
void inverse_p_transform(int nx, int nt, float dt, float pmax, float pmin,
	float dp, float depthref, float f1, float f2, float interoff,
	float offref, int igopt, float dx, float fx, float **in_traces,
	float **out_traces)
/******************************************************************************
Input:
nx              number of input traces
nt              number of intput time samples
dt              time sampling interval (seconds)
dx              spatial sampling interval (meters)
fx		first offset (meters)
igopt           =1 parabolic trransform g(x)=offset**2
                =2 Foster/Mosher pseudo hyperbolic transform
                        g(x)=sqrt(depth**2+offset**2)
                =3 linear tau-p g(x)=offset
                =4 abs linear taup g(x)=abs(offset)
offref          reference maximum offset to which maximum and minimum moveout
                times are associated
interoff        intercept offset to which tau-p times are associated
                (usually zero)
pmax            maximum moveout in ms on reference offset
pmin            minimum moveout in ms on reference offset
dp              moveout sampling interval (ms/m)
depthref        reference depth for Foster/Mosher hyperbolic transform
f1              high-end  frequency before taper off (hz)
f2              high-end frequency (hz)
in_traces       2-D array[np][ntfft]  of input taup traces

Output:
out_traces      2-D array[nx][nt] of output t-x traces

*******************************************************************************
Credits:
	Adapted by Gabriel Alvarez (1995) from suradon.c written by John
        Anderson (1993)
******************************************************************************/
{
	int ix,ip,it,iw;		/* loop counters */
	int ntfft;			/* length of time fft  */
	int nw;				/* number of frequencies */
	int np;				/* number of sopes for taup transform*/
	size_t nmax;			/* maximum space to allocate */
	float rsum,isum;		/* real and imaginary part of a sum */
	float w;			/* frequency */
	float p;			/* slope */
	float dw,df;			/* frequency sampling intervals */
	float wa;			/* frequency weight */
	float fac;			/* fft scaling factor */
	float dr,di,tr,ti;		/* r and im parts of aux complex num*/
	float *rt;			/* auxiliary vector to store a trace */
	float *g;			/* originally input arrays */
	float *offset;			/* array of offsets */
	complex czero;			/* complex zero number */
	complex *crt;			/*complex array for FFT's */
	complex *ctemp;			/* auxiliary complex array */
	char *fname;			/* character pointer to file name */
	VND *vnda;			/* 2-D VND scratch array */

	/* define variables */
	ntfft = npfar(nt);
	fac = 1000.*gofx(igopt,offref,interoff,depthref);
	pmin /= fac;
	pmax /= fac;
	dp /= fac;
	np = 1+(pmax-pmin)/dp;
	nmax = MAX(nx,np);
	nmax = MAX(nmax,ntfft+4);

	fprintf(stderr, "computing inverse radon transform: pmin=%g pmax=%g"
		" dp=%g np=%d\n",pmin,pmax,dp,np);

	/* allocate 1-D arrays */
	offset  = (float *)VNDemalloc(nx*sizeof(float),"invslant:offset");
	g 	= (float *)VNDemalloc(nx*sizeof(float),"invslant:g");