taup.c

su 的源代码库

	/* determine frequency and wavenumber sampling */
	nw = ntfft/2+1;
	dw = 2.0*PI/(ntfft*dt);
	fw = 0.0;
	nk = nxfft;
	dk = 2.0*PI/(nxfft*dx);
	fk = -PI/dx;
	/* lk = PI/dx; */
	fka = fk-lwrap*dk;
	nka = lwrap+nk+lwrap;

	/* allocate working space for FFT's */
	tr_fft = alloc1float(ntfft);
	ctr = alloc2complex(nw,nk);
	ctr_p = alloc2complex(nw,np);
	tr_ka = alloc1complex(nka);
	tr_x = tr_k=tr_ka+lwrap;        /* pointers to tr_ka[lwrap] */
	hp = alloc1complex(np);
	kp = alloc1float(np);

	/* loop over traces */
	for (ix=0; ix<nx; ix++) {

		/* pad time axis with zeros */
		for (it=0; it<nt; it++)
		tr_fft[it]=traces[ix][it];
		for (it=nt; it<ntfft; it++)
		tr_fft[it]=0.0;

		/* Fourier transform time to frequency */
		pfarc(1,ntfft,tr_fft,ctr[ix]);
	}

	/* loop over w */
	for (iw=0, w=fw; iw<nw; ++iw, w+=dw) {

	/* scale tr_x by x sampling interval */
	for (ix=0; ix<nx; ix++) {
		tr_x[ix].r = ctr[ix][iw].r*dx*fftscl;
		tr_x[ix].i = ctr[ix][iw].i*dx*fftscl;
	}

	/* pad tr_x with zeros */
	for (ix=nx; ix<nxfft; ix++)
		tr_x[ix].r = tr_x[ix].i = 0.0;

		/* negate every other sample for k-axis centered */
		for (ix=1; ix<nx; ix+=2) {
			tr_x[ix].r = -tr_x[ix].r;
			tr_x[ix].i = -tr_x[ix].i;
		}

		/* Fourier transform tr_x to tr_k */
		pfacc (-1, nxfft, tr_x);
				
		/* wrap-around tr_k to avoid interpol end effects */
		for (ik=0; ik<lwrap; ik++)
			tr_ka[ik] = tr_k[ik+nk-lwrap];
		for (ik=lwrap+nk; ik<lwrap+nk+lwrap; ik++)
			tr_ka[ik] = tr_k[ik-lwrap-nk];

		/* phase shift to account for non-centered x-axis */
		xshift = 0.5*(nx-1)*dx;
		for (ik=0, k=fka; ik<nka; ik++, k+=dk) {
			phase = k*xshift;
			c = cos(phase);
			s = sin(phase);
			temp = tr_ka[ik].r*c-tr_ka[ik].i*s;
			tr_ka[ik].i = tr_ka[ik].r*s+tr_ka[ik].i*c;
			tr_ka[ik].r=temp;
		}

		/* compute k values at which to interpolate tr_k */
		for (ip=0, p=fp; ip<np; ip++, p+=dp) {
			kp[ip] = w*p;

			/* if outside Nyq bounds do not interpolate 
			if (kp[ip]<fk && kp[ip]<lk)
				kp[ip] = fk-1000.0*dk;
			else if (kp[ip]>fk && kp[ip]>lk)
				kp[ip] = lk+1000.0*dk; */
		}

		/* 8-point sinc interpolation of tr_k to obtain h(p) */
		ints8c (nka, dk, fka, tr_ka, czero, czero, np, kp, hp);

		/* phase shift to account for non-centered x-axis */
		xshift = -fx - 0.5*(nx-1)*dx;
		for (ip=0; ip<np; ip++) {
			phase = kp[ip]*xshift;
			c = cos(phase);
			s = sin(phase);
			ctr_p[ip][iw].r = hp[ip].r*c-hp[ip].i*s;
			ctr_p[ip][iw].i = hp[ip].r*s+hp[ip].i*c;
		}
	}

	/* loop over p */
	for (ip=0; ip<np; ip++) {

		/* Fourier transform frequency to time */
		pfacr(-1, ntfft, ctr_p[ip], tr_fft);

		/* copy to output array */
		for (itau=0; itau<ntau; itau++)
		out_traces[ip][itau] = tr_fft[itau];
	}

	/* clean up */
	free1float(tr_fft);
	free2complex(ctr);
	free2complex(ctr_p);
	free1complex(tr_ka);
	free1complex(hp);
	free1float(kp);
}

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

	Subroutine to compute a forward slant stack (taup transform)
				in t-x domain

******************************************************************************/
void fwd_tx_sstack (float dt, int nt, int nx, float xmin, float dx, int np,
	float pmin, float dp, float **traces, float **out_traces)
/******************************************************************************
Input:
dt              time sampling interval
nt              number of time samples
nx              number of horizontal samples (traces)
np              number of slopes
pmin            minimum slope for tau-p transform
dp		slope sampling interval
traces          2-D array of input traces in t-x domain

Output:
traces          2-D array of output traces in tau-p domain

Credits:
written by Gabriel Alvarez, CWP
modified by Bjoern E. Rommel, IKU, Petroleumsforskning
******************************************************************************/
{
	int ip,ix,it;		/* loop counters */
	float x,p;		/* auxilairy variables */
	int fit,lit;		/* first and last time samples */
	int id;                 /* auxiliary variables */
	float dfrac,delay;	/* more auxiliary variables */

	/* loop over slopes */
	for (ip=0, p=pmin; ip<np; p=pmin+(++ip)*dp) {

		/* initialize output array */
		for (it=0; it<nt; it++)
			out_traces[ip][it]=0.0;

		/* loop over traces */
		for (ix=0, x=xmin; ix<nx; x=xmin+(++ix)*dx) {

			/* compute two point interpolator */
			delay=p*x/dt;
			if (delay>=0) {
				id = (int)delay;
				fit = id+1;
				lit = nt-1;
			} else {
				id = (int)delay-1;
				fit = 1;
				lit = nt+id;
			}	
			dfrac = delay-id;

			/* compute the actual slant stack */
			for (it=fit; it<lit; it++) {
				out_traces[ip][it-id]+=fabs(dx)*
				    (traces[ix][it]+
				     dfrac*(traces[ix][it+1]-traces[ix][it]));
			}
		}
        }
}

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

        Compute inverse global slant stack in FK domain

******************************************************************************/
void inv_FK_sstack (float dt, int nt, int nx, float xmin, float dx, int np,
	float pmin, float dp, float fmin, float **traces, float **out_traces) 
/******************************************************************************
Input Parameters:
dt              time sampling interval
nt              number of time samples
nx              number of horizontal samples
xmin		minimum horizontal distance (ignored)
np              number of slopes
pmin            minimum slope for inverse tau-p transform
dp		slope sampling interval
fmin            minimum frequency of interest (ignored)
traces          2-D array of input traces in tau-p domain

Output Parameters:
out_traces      2-D array of output traces in t-x domain
******************************************************************************/

{
        int ix,ip,itau,iw,ik,it;/* loop counters */
        float fw;               /* first frequency */
        float dw,dk;            /* frequency and wavenumber sampling interval */
        int nw,nk;              /* number of frequencies and wavenumbers */
        int ntaufft;            /* number of padded time samples */
        int nkfft;              /* number of padded wavenumber samples */
	int nkmax;
	int ntau=nt;
	float pmax;
	float dtau=dt;		/* time intercept sampling ingerval */
        float fp=0.,fk;            /* first slope and wavenumber */
	float w,k;		/* angular frequency, slope and wavenumber */
	float fftscl;  		/* scale factor for FFT */
        float *pk;		/* k-interpolated w trace */
        float *tr_fft;          /* padded trace for FFT */
	complex *tr_p;
        complex **ctr;          /* */
        complex **ctr_x;        /* */
        complex *hk;            /* slant stacked single trace */
	complex czero;		/* complex zero number */

        /* compute slope sampling interval */
	czero = cmplx(0.0*fmin,0.0*xmin);
	pmax=pmin+(np-1)*dp;

	/* determine lengths and scale factors for FFT's */
	fp=pmin+0.5*(pmax-pmin-(np-1)*dp);
	pmax=fp+(np-1)*dp;	
	ntaufft = npfar(2*ntau); 
	fftscl = 1.0/ntaufft;
	nkfft = npfa(4*np);

	/* determine frequency and wavenumber sampling */
	nw = ntaufft/2+1;
	dw = 2.0*PI/(ntaufft*dtau);
	fw = fk = 0.0;
	dk = 2.0*PI/(nkfft*dx);
	nkmax = PI*pmax/(dt*dk) + 1;

	/* allocate working space for FFT's */
	tr_fft = alloc1float(ntaufft);
	ctr = alloc2complex(nw,np);
	ctr_x = alloc2complex(nw,nx);
	tr_p = alloc1complex(np); 
	hk = alloc1complex(nkmax);

	/* loop over traces in tau-p domain */
	for (ip=0; ip<np; ip++) {

		/* pad tau axis with zeros */
		for (itau=0; itau<ntau; itau++)
		tr_fft[itau]=traces[ip][itau];
		for (itau=ntau; itau<ntaufft; itau++)
			tr_fft[itau]=0.0;

		/* Fourier transform tau to frequency */
		pfarc(1,ntaufft,tr_fft,ctr[ip]);
	}

	/* loop over w */
	for (iw=0, w=fw; iw<nw; ++iw, w+=dw) {

	/* scale ctr by p sampling interval */
	for (ip=0; ip<np; ip++)
		tr_p[ip] = crmul(ctr[ip][iw],dp*fftscl);

		/* compute number of k's for interpolation */
		nk = pmax*w/dk + 1;
               	pk = alloc1float(nk);

		/* compute p values at which to interpolate tr_p */
		for (ik=0, k=fk; ik<nk; ik++, k+=dk) {
			if (w==0.0) pk[ik]=0.0;  
			else pk[ik] = k/w;
		}
	
		/* interpolate tr_pa to obtain tr_k */ 
		ints8c (np, dp, fp, tr_p, czero, czero, nk, pk, hk); 
	
		/* free temporary space */
		free1float(pk);	

		/*pad hk with zeros  */
		for (ik=nk; ik<nkfft; ik++) 
			hk[ik]=czero;

		/* inverse Fourier transform from k to x */
		pfacc (1, nkfft, hk);

		/* copy 1-D interpolated array to 2-D array */
		for (ix=0; ix<nx; ix++)
			ctr_x[ix][iw]=hk[ix];
	}

	/* inverse Fourier transform from frequency to time */
	for (ix=0; ix<nx; ix++) {
		pfacr(-1,ntaufft,ctr_x[ix],tr_fft);
		for (it=0; it<nt; it++)
			out_traces[ix][it]=tr_fft[it];	
	} 	

	/* clean up */
	free1float(tr_fft);
	/* free1complex(hk);  */
	/* free1complex(tr_p); */
	free2complex(ctr);
	free2complex(ctr_x);
}

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

	Subroutine to compute an inverse slant stack (taup transform)
				in t-x domain

******************************************************************************/
void inv_tx_sstack (float dt, int nt, int nx, int npoints, float xmin, 
	float dx, int np, float pmin, float dp, float **traces,
	float **out_traces) 
/******************************************************************************
Input Parameters:
dt              time sampling interval
nt              number of time samples
nx              number of horizontal samples
np              number of slopes
pmin            minimum slope for inverse tau-p transform
dp		slope sampling interval
traces          2-D array of input traces in tau-p domain

Output Parameters:
out_traces      2-D array of output traces in t-x domain
******************************************************************************/
{
	int ip;			/* loop counter */
	int np2=npoints/2;	/* half number of points in rho filter */
	float *rho;		/* auxiliary array for rho filter */
	float **rho_traces;	/* aux array for traces convolved with rho */

	/* allocate space */
	rho_traces = alloc2float(nt,np);
	rho = alloc1float(npoints);

	/* compute rho filter */
	rho_filter (npoints, nt, dt, rho);

       	/* convolve input traces with rho filter */
        for (ip=0; ip<np; ip++)
       	        conv(nt,-np2,traces[ip],npoints,0,rho,nt,0,rho_traces[ip]);

	/* inverse transform == forward transform with negative slopes */
	fwd_tx_sstack (dt, nt, np, -pmin, -dp, nx, xmin, dx, rho_traces,
		out_traces);

	/* clean up */
	free2float(rho_traces);
	free1float(rho);
}

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

        Subroutine to compute the rho filter in frequenccy
         domain for the time domain inverse slant stack

******************************************************************************/
void rho_filter (int npoints, int nt, float dt, float *rho)
/******************************************************************************
Input Parameters:
npoints         number of point for the rho filter
nt              number of time samples
dt              time sampling interval

Output Parameters:
rho             1-D array of filter points

Credits:
written by Gabriel Alvarez, CWP
corrected by Bjoern E. Rommel, IKU, Petroleumsforskning
******************************************************************************/
{
	const int maxnpfa = 720720;