suradon.c

su 的源代码库

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

static float rcdot(int n, complex *a, complex *b)
/********************************************************************  
return the real part of a complex dot product where
    the first vector is the one complex conjugated
*********************************************************************
Author: John Anderson (visitor to CSM from Mobil) Spring 1993
*********************************************************************/
{
	int j;
	float sum=0.;
	for(j=0;j<n;j++) sum += a[j].r * b[j].r + a[j].i * b[j].i;
	return(sum);
}

static void htmul(int n, complex *a, complex *x, complex *y)

/*******************************************************************
   Hermitian Toeplitz matrix multiply

     solve for y = A x   where A is Hermitian Toeplitz

     and defined by the vector a giving the top row of A.
     x is input.  y is output. 
*******************************************************************
Author: John Anderson (visitor to CSM from Mobil) Spring 1993
*******************************************************************/
{
	int j,irow;
	complex czero;
	czero.r=czero.i=0.;

	for(irow=0;irow<n;irow++) {
		y[irow]=czero;
		for(j=0;j<irow;j++)
			y[irow] = cadd(cmul(conjg(a[irow-j]),x[j]),y[irow]);
		for(j=irow;j<n;j++)
			y[irow] = cadd(cmul(a[j-irow],x[j]),y[irow]);
	}
}

static void jea_xinterpolate(VND *vndorig, VND *vndinterp, 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:

VND *vndorig		VND file with input data
VND *vndinterp		VND file with output original plus interpolated data
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 in Hz for picking
						(good default: 3 Hz)
float freq2		high-end frequency in Hz 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

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.

This routine assumes that the input and output files hav been allocated in
the calling routine as

char *fname;
fname=VNDtempname("main_prog");
vndorig = V2Dop(2,1000000,sizeof(float),fname,nt,nxmax);
VNDfree(fname,"jea_xinterpolate: fname");
fname=VNDtempname("main_prog");
vndinterp = V2Dop(2,1000000,sizeof(float),fname,
			nt,1+(nxmax-1)*(ninterp+1));
VNDfree(fname,"main: fname");

where nxmax is the maximum number of input traces and nt is the number 
of time samples.
*******************************************************************
Author: John Anderson (visitor to CSM from Mobil) Spring 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;

	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);

	crt = (complex *)VNDemalloc( MAX(ntfftny,nxfftny)*sizeof(complex),
		"jea_xinterpolate:allocating crt" );
	rt = (float *)crt;

	if(nx<2 || (iopt==0 && ninterp==0) ) {
	    	for(ix=0;ix<nx;ix++) {
			V2Dr0(vndorig,ix,(char *)rt,101);	
			V2Dw0(vndinterp,ix,(char *)rt,102);
		}
		VNDfree(crt,"crt");
		return;
	}

	ccrt = (complex *)VNDemalloc( ntfftny*sizeof(complex),
		"jea_xinterpolate:allocating ccrt" );
	rrt = (float *)ccrt;
	a =  (float *)VNDemalloc( MAX(nx,nt)*sizeof(float),
		"jea_xinterpolate:allocating a" );
	b =  (float *)VNDemalloc( MAX(nx,nt)*sizeof(float),
		"jea_xinterpolate:allocating b" );
	p =  (float *)VNDemalloc( nt*sizeof(float),
		"jea_xinterpolate:allocating p" );
	time =  (float *)VNDemalloc( nt*sizeof(float),
		"jea_xinterpolate:allocating time" );
	aa =  (float *)VNDemalloc( MAX(nx,nt)*sizeof(float),
		"jea_xinterpolate:allocating aa" );
	bb =  (float *)VNDemalloc( MAX(nx,nt)*sizeof(float),
		"jea_xinterpolate:allocating bb" );

	fname = VNDtempname("jea_xinterpolate");
	vnda  = V2Dop(2,500000,sizeof(float),fname,nt,nx);
	VNDfree(fname,"jea_xinterpolate: fname");
	fname = VNDtempname("jea_xinterpolate");
	vndb  = V2Dop(2,500000,sizeof(float),fname,nt,nx);
	VNDfree(fname,"jea_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++) {
		V2Dr0(vndorig,ix,(char *)rt,103);
		for(j=0;j<nt;j++) a[j]=fabs(rt[j]);
		runav(nt,lagc,a,b);
		runav(nt,lagc,b,a);
		for(j=0;j<nt;j++) rt[j]=rt[j]/(a[j]+eps);	
		for(j=nt;j<ntfft;j++) rt[j]=0.;
		pfarc(1,ntfft,rt,crt);
		if(freq1>0.){
			for(j=0;j<ntfftny;j++){
				f=j*df;
				fcl=(f/freq1);
				fcl=fcl*fcl*fcl*fcl;
				fch=(f/freq2);
				fch=fch*fch*fch*fch;
				f=fcl/( (1.+fcl)*(1.+fch) );
				crt[j]=crmul(crt[j],f);
				ccrt[j]=cmul(crt[j],cmplx(0.,-j*dff));
			}
		}else{
			for(j=0;j<ntfftny;j++){
				f=j*df;
				fch=(f/freq2);
				f=1./(1.+fch*fch*fch*fch);
				crt[j]=crmul(crt[j],f);
				ccrt[j]=cmul(crt[j],cmplx(0.,-j*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);
			V2Dw0(vndinterp,ix,(char *)rt,104);
		}
		VNDcl(vnda,1);
		VNDcl(vndb,1);
		VNDfree(crt,"jea_xinterpolate: crt");
		VNDfree(ccrt,"jea_xinterpolate: ccrt");
		VNDfree(a,"jea_xinterpolate: a");
		VNDfree(b,"jea_xinterpolate: b");
		VNDfree(p,"jea_xinterpolate: p");
		VNDfree(time,"jea_xinterpolate: time");
		VNDfree(aa,"jea_xinterpolate: aa");
		VNDfree(bb,"jea_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 */
	V2Dr0(vndorig,0,(char *)a,116);
	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 );
		}
		V2Dr0(vndorig,ix,(char *)b,119);
		if(iopt==2) {
			V2Dw0(vndinterp,ixm,(char *)p,120);
			/* don't output dip picks except on original traces */
		}else{
			V2Dw0(vndinterp,ixm*(ninterp+1),(char *)a,120);
			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];
				V2Dw0(vndinterp,k+1+ixm*(ninterp+1),
				      (char *)aa,121);
			}
		}
		save=a;
		a=b;
		b=save;  
	} 
	if(iopt==2) {
		V2Dw0(vndinterp,nx-1,(char *)p,122);
	}else{
		V2Dw0(vndinterp,(nx-1)*(ninterp+1),(char *)a,122);
	}


/* close files, free temporary memory, and return results in file vndinterp */
	VNDcl(vnda,1);
	VNDcl(vndb,1);
	VNDfree(crt,"jea_xinterpolate: crt");
	VNDfree(ccrt,"jea_xinterpolate: ccrt");
	VNDfree(a,"jea_xinterpolate: a");
	VNDfree(b,"jea_xinterpolate: b");
	VNDfree(p,"jea_xinterpolate: p");
	VNDfree(time,"jea_xinterpolate: time");
	VNDfree(aa,"jea_xinterpolate: aa");
	VNDfree(bb,"jea_xinterpolate: bb");

	return;
}
static void runav(int n,int len,float *a,float *b)
/******************************************************************
compute a boxcar running average filter
*******************************************************************
int n   	number of samples in a[] and b[]
int len 	length of running average in samples
float a[n]	input array
float b[n]	output array
*******************************************************************
Author: John Anderson (visitor to CSM from Mobil) Spring 1993
*******************************************************************/
{
	float sum=0.;
	int j,lenh=len/2;
	if(len<=1) {
		for(j=0;j<n;j++) b[j]=a[j];
		return;
	}
	for(j=0;j<MIN(len,n);j++) sum+=a[j];
	for(j=0;j<MIN(lenh+1,n);j++) b[j]=sum;
	for(j=lenh+1;j<n-lenh;j++) {
		sum=sum+a[j+lenh]-a[j-lenh-1];
		b[j]=sum;
	}
	for(j=MAX(0,n-lenh);j<n;j++) b[j]=sum;
	sum=1./len;
	for(j=0;j<n;j++) b[j]*=sum;
	return;
}