taup.c

su 的源代码库

	/*  Compute (er,ei) = 1./(vr,vi)   */
	er = vr/vsq;
	ei = -vi/vsq;
	for(k=0;k<j;k++) {
		f[k].r = er*f[k].r - ei*f[k].i;
		f[k].i = er*f[k].i + ei*f[k].r;	
	}
	return (j);
}

int ctoephcg( int niter, int n, complex *a, complex *x, complex *y, 
	complex *s, complex *ss, complex *g, complex *rr)

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

Hestenes and Stiefel conjugate gradient algorithm 
specialized for solving Hermitian Toeplitz
system.  a[] is input as a vector defining the only the
top row of A.  x[] is the solution vector returned.
y[] is input.  niter is the maximum number of conjugate 
gradient steps to compute.  The function returns as
the number of steps actually computed.  The other 
vectors provide workspace.

Complex Hermitian Toeplitz Solver for

N-1
Sum  A	     x  = y      for i=0,1,2,...,N-1
j=0   (i-j)   j    i

where A is Hermitian Toeplitz and x and y are complex.  For
an example 4 x 4 system,  x returns as the solution of


   A0  A1  A2  A3	x0	     y0

     *
   A1  A0  A1  A2	x1	     y1
				=    
     *   *
   A2  A1  A0  A1	x2	     y2

     *   *   *
   A3  A2  A1  A0	x3	     y3

*******************************************************************************
Author: John Anderson (visitor to CSM from Mobil) Spring 1993
******************************************************************************/
{
	int j, iter;
	complex czero;
	float alpha, beta, gamma, gammam, rsq, rp, test;
	float eps=1.0e-6;
	float rcdot(int n, complex *a, complex *b);
	rp   = rcdot(n,y,y);
	test = n*eps*eps*rp;
	czero.r=czero.i=0.;

	for(j=0;j<n;j++) {
		x[j]=czero;
		rr[j]=y[j];
	}
	htmul(n,a,rr,g);	   /*  adjoint matrix multiply */

	for(j=0;j<n;j++) s[j]=g[j];
	gammam=rcdot(n,g,g);

	for(iter=0;iter<niter;iter++) { /* forward matrix multiply  */
		htmul(n,a,s,ss);  
		alpha  = gammam/rcdot(n,ss,ss);
		for(j=0;j<n;j++) {
			x[j] =cadd(x[j],crmul(s[j],alpha));
			rr[j]=csub(rr[j],crmul(ss[j],alpha));
		}
		rsq = rcdot(n,rr,rr);
		if ( iter>0 && ( rsq==rp || rsq<test ) ) return(iter-1);
		rp = rsq;

		htmul(n,a,rr,g);   /*  adjoint matrix multiply  */
		gamma  = rcdot(n,g,g);
		if (gamma<eps) break;
		beta   = gamma/gammam;
		gammam = gamma;			

		for(j=0;j<n;j++) {
			s[j] =cadd(g[j],crmul(s[j],beta));
		}	
	}
return(iter);
}

void ga_xinterpolate(float **in_traces, float **out_traces, int ninterp, 
		int nt, int nx, float freq1, float freq2, int lagc, 
		int lent, int lenx, int xopt, float dt, int iopt)
/******************************************************************************
interpolate input data in space placing "ninterp" synthetic traces 
between each pair of original input traces
*******************************************************************************
Function parameters:
Input:
int ninterp		number of traces to interpolate between each input trace
int nt			number of time samples
int nx			number of input traces
float freq1		low-end frequency for picking (good default: 3 Hz)
float freq2		high-end frequency for picking (good default: 20 Hz)
int lagc		length of AGC operator for picking (good default:400 ms)
int lent		length of time smoother in samples for picker
                        (good default: 5 samples)
int lenx		length of space smoother in samples for picker
                        (good default: 1 sample)
int xopt		1 = use differences for spatial derivative
                            (works with irregular spacing)
                        0 = use FFT derivative for spatial derivatives
                            (more accurate but requires regular spacing and
                            at least 16 input tracs--will switch to differences
                            automatically if have less than 16 input traces)
float dt		sample rate in sec
int iopt		0 = interpolate: output 1+(nx-1)*(1+ninterp) traces
                           with ninterp traces between each pair of input traces
			1 = compute low-pass model: output nx traces
                            on original trace locations -- This is typically
                            used for Quality Control if the interpolator
                            is failing for any reason
			2 = compute dip picks in units of samples/trace: 
                            output nx traces on original trace locations
in_traces		2-D array of input data

Output:
out_traces		2-D array of interpolated traces


Notes:
This routine outputs 'ninterp' interpolated traces between each pair of 
input traces.  The values for lagc, freq1, and freq2 are only used for
event tracking. The output data will be full bandwidth with no agc.  The 
suggested default parameters typically will do a satisfactory job of 
interpolation for dips up to about 12 ms/trace.  Using a larger value for 
freq2 causes the algorithm to do a better job on the shallow dips, but to 
fail on the steep dips.  Only one dip is assumed at each time sample between 
each pair of input traces.  The original input traces are passed through
this routine without modification.

The key assumption used here is that the low frequency data are unaliased
and can be used for event tracking.  Those dip picks are used to interpolate
the original full-bandwidth data, giving some measure of interpolation
at higher frequencies which otherwise would be aliased.  Using iopt equal
to 1 allows you to visually check whether the low-pass picking model is aliased.
If you can't visually pick dips correctly on the low-pass picking 
model, this computer routine will fail.

The place this code is most likely to fail is on the first breaks.
*******************************************************************************
Credits:
Adapted by Gabriel Alvarez (1995) from suradon.c written by John 
	Anderson (1993)
******************************************************************************/
{
	int	ntfft,ntfftny,nxfft,nxfftny,j,k,ix,it,ixm;
	float	df,dff,wa,wb,dxx,eps=1.0e-30,f,fcl,fch;
	float 	*rt,*rrt,*a,*b,*p,*time,*aa,*bb,*save;
	complex	*crt,*ccrt;
	VND	*vnda,*vndb;
	char 	*fname;

	/* defensive programming */
	if(nx<2 || (iopt==0 && ninterp==0) ) {
		for (ix=0;ix<nx;ix++)
			for (it=0;it<nt;it++)
				out_traces[ix][it]=in_traces[ix][it];	
		return;
	}

	/* define useful variables */
	lent=1+2*(lent/2);
	lenx=1+2*(lenx/2);
	lagc=1 + lagc*0.001/dt;

	ntfft=npfar(nt);
	ntfftny=1+ntfft/2;
	nxfft=npfar(nx);
	nxfftny=1+nxfft/2;

	df=1./(ntfft*dt);

	/* allocate working VND space */
	crt = (complex *)VNDemalloc( MAX(ntfftny,nxfftny)*sizeof(complex),
		"ga_xinterpolate:allocating crt" );
	rt = (float *)crt;
	ccrt = (complex *)VNDemalloc( ntfftny*sizeof(complex),
		"ga_xinterpolate:allocating ccrt" );
	rrt = (float *)ccrt;
	a =  (float *)VNDemalloc( MAX(nx,nt)*sizeof(float),
		"ga_xinterpolate:allocating a" );
	b =  (float *)VNDemalloc( MAX(nx,nt)*sizeof(float),
		"ga_xinterpolate:allocating b" );
	p =  (float *)VNDemalloc( nt*sizeof(float),
		"ga_xinterpolate:allocating p" );
	time =  (float *)VNDemalloc( nt*sizeof(float),
		"ga_xinterpolate:allocating time" );
	aa =  (float *)VNDemalloc( MAX(nx,nt)*sizeof(float),
		"ga_xinterpolate:allocating aa" );
	bb =  (float *)VNDemalloc( MAX(nx,nt)*sizeof(float),
		"ga_xinterpolate:allocating bb" );

	fname = VNDtempname("ga_xinterpolate");
	vnda  = V2Dop(2,500000,sizeof(float),fname,nt,nx);
	VNDfree(fname,"ga_xinterpolate: fname");
	fname = VNDtempname("ga_xinterpolate");
	vndb  = V2Dop(2,500000,sizeof(float),fname,nt,nx);
	VNDfree(fname,"ga_xinterpolate: fname");

	/* loop computing filtered data for picking purposes in vnda */
	/* compute time derivative of filtered data in vndb */
	dff=2.*PI/ntfft;
	for(ix=0;ix<nx;ix++) {
		for(it=0;it<nt;it++) {
			rt[it]=in_traces[ix][it];
			a[it]=fabs(rt[it]);
		}
		runav(nt,lagc,a,b);
		runav(nt,lagc,b,a);
		for(it=0;it<nt;it++) rt[it]=rt[it]/(a[it]+eps);	
		for(it=nt;it<ntfft;it++) rt[it]=0.;
		pfarc(1,ntfft,rt,crt);
		for(it=0;it<ntfftny;it++){
			f=it*df;
			fcl=(f/freq1);
			fcl=fcl*fcl*fcl*fcl;
			fch=(f/freq2);
			fch=fch*fch*fch*fch;
			f=fcl/( (1.+fcl)*(1.+fch) );
			crt[it]=crmul(crt[it],f);
			ccrt[it]=cmul(crt[it],cmplx(0.,-it*dff));
		}
		pfacr(-1,ntfft,crt,rt); 
		V2Dw0(vnda,ix,(char *)rt,104);
		pfacr(-1,ntfft,ccrt,rrt); 
		V2Dw0(vndb,ix,(char *)rrt,105);
	} 

	if(iopt==1){
		for(ix=0;ix<nx;ix++){
			V2Dr0(vnda,ix,(char *)rt,104);
			for (it=0;it<nt;it++)
				out_traces[ix][it]=rt[it];
		}
		VNDcl(vnda,1);
		VNDcl(vndb,1);
		VNDfree(crt,"ga_xinterpolate: crt");
		VNDfree(ccrt,"ga_xinterpolate: ccrt");
		VNDfree(a,"ga_xinterpolate: a");
		VNDfree(b,"ga_xinterpolate: b");
		VNDfree(p,"ga_xinterpolate: p");
		VNDfree(time,"ga_xinterpolate: time");
		VNDfree(aa,"ga_xinterpolate: aa");
		VNDfree(bb,"ga_xinterpolate: bb");
		return;
	}

	/* loop computing spatial derivative of data for picking purposes*/
	nxfft=npfar(nx);
	nxfftny=1+nxfft/2;
	dxx=2.*PI/(nxfft*nxfft);
	if(nx<16) xopt=1;
	for(it=0;it<nt;it++) {
		V2Dr1(vnda,it,(char *)rt,106);
		if(xopt) {
			for(j=0;j<nx-1;j++) rt[j]=rt[j+1]-rt[j];
			rt[nx-1]=rt[nx-2];
		}else{
			for(j=nx;j<nxfft;j++) rt[j]=0.;
			pfarc(1,nxfft,rt,crt);
			for(j=0;j<nxfftny;j++){
				crt[j]=cmul(crt[j],cmplx(0.,-j*dxx));
			}
			pfacr(-1,nxfft,crt,rt); 
		}
		V2Dw1(vnda,it,(char *)rt,107);
	} 

	/* compute dot products and smooth over time */
	for(ix=0;ix<nx;ix++) {
		V2Dr0(vnda,ix,(char *)a,108);
		V2Dr0(vndb,ix,(char *)b,109);
		for(it=0;it<nt;it++) {
			aa[it]=a[it]*b[it];
			bb[it]=b[it]*b[it];
		}
		runav(nt,lent,aa,a);
		runav(nt,lent,a,aa);
		runav(nt,lent,bb,b);
		runav(nt,lent,b,bb);
		V2Dw0(vnda,ix,(char *)aa,110);
		V2Dw0(vndb,ix,(char *)bb,111);
	}

	/* smooth dot products in x */
	if(lenx>1){
	    for(it=0;it<nt;it++) {
		V2Dr1(vnda,it,(char *)a,112);
		V2Dr1(vndb,it,(char *)b,113);
		runav(nx,lenx,a,aa);
		runav(nx,lenx,aa,a);
		runav(nx,lenx,b,bb);
		runav(nx,lenx,bb,b);
		V2Dw1(vnda,it,(char *)a,114);
		V2Dw1(vndb,it,(char *)b,115);
	    }
	}

	/* loop computing p, interpolating, and outputting results */
	/* get first trace from input data */
	for (it=0; it<nt; it++)
		a[it]=in_traces[0][it];
	for(ix=1;ix<nx;ix++) {
		ixm=ix-1;
		V2Dr0(vnda,ixm,(char *)aa,117);
		V2Dr0(vndb,ixm,(char *)bb,118);
		for(it=0;it<nt;it++) {
			p[it] = - aa[it]/( bb[it] + eps );
		}
		
		/* get input traces one by one */
		for (it=0; it<nt; it++)
			b[it]=in_traces[ix][it];	
		if(iopt==2) {
	
			/* write to output array */
			for (it=0; it<nt; it++)
				out_traces[ixm][it]=p[it];
			/* don't output dip picks except on original traces */
		}else{
			/* write to output array */
			for (it=0; it<nt; it++)
				out_traces[ixm*(ninterp+1)][it]=a[it];
			for(k=0;k<ninterp;k++){
				wa=(1.+k)/(1+ninterp);
				wb=1.-wa;
				for(it=0;it<nt;it++) time[it] = it - p[it]*wa;		
				ints8r(nt,1.0,0.,a,0.0,0.0,nt,time,aa);
				for(it=0;it<nt;it++) time[it] = it + p[it]*wb;		
				ints8r(nt,1.0,0.,b,0.0,0.0,nt,time,bb);
				for(it=0;it<nt;it++) aa[it]=wb*aa[it]+wa*bb[it];

				/* write to output array */
				for (it=0; it<nt; it++)
					out_traces[k+1+ixm*(ninterp+1)][it]=
						aa[it];
			}
		}
		save=a;
		a=b;
		b=save;  
	} 
	if(iopt==2) {
		
		/* write to output array */
		for (it=0; it<nt; it++)
			out_traces[nx-1][it]=p[it];
	}else{
		/* write to output array */
		for (it=0; it<nt; it++)
			out_traces[(nx-1)*(ninterp+1)][it]=a[it];
	}


	/* close files, free temporary memory */
	VNDcl(vnda,1);
	VNDcl(vndb,1);
	VNDfree(crt,"ga_xinterpolate: crt");
	VNDfree(ccrt,"ga_xinterpolate: ccrt");
	VNDfree(a,"ga_xinterpolate: a");
	VNDfree(b,"ga_xinterpolate: b");
	VNDfree(p,"ga_xinterpolate: p");
	VNDfree(time,"ga_xinterpolate: time");
	VNDfree(aa,"ga_xinterpolate: aa");
	VNDfree(bb,"ga_xinterpolate: bb");

	return;
}

float rcdot(int n, complex *a, com