suinterp.c

su 的源代码库

Author: John Anderson (visitor to CSM from Mobil) Spring 1993
*******************************************************************/

{
	int	ntfft,ntfftny,n2fft,n2fftny,j,k,ixm;
	long	ix,it,ntlong,n2long;
	float	df,dff,wa,wb,dxx,eps=1.0e-30,f,fcl,fch;
	float 	*rt,*rrt,*a,*b,*p,*time,*aa,*bb,*save,*tin=NULL;
	complex	*crt,*ccrt;
	VND	*vnda,*vndb;
	char 	*file;

	ntlong=nt;
	n2long=n2;
	lent=1+2*(lent/2);
	len2=1+2*(len2/2);
	lagc=1 + lagc*0.001/dt;

	ntfft=npfar(nt);
	ntfftny=1+ntfft/2;
	n2fft=npfar(n2);
	n2fftny=1+n2fft/2;

	df=1./(ntfft*dt);

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

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

	ccrt = (complex *)VNDemalloc( ntfftny*sizeof(complex),
		"jea_xinterpolate:allocating ccrt" );
	rrt = (float *)ccrt;
	a =  (float *)VNDemalloc( MAX(n2,nt)*sizeof(float),
		"jea_xinterpolate:allocating a" );
	b =  (float *)VNDemalloc( MAX(n2,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(n2,nt)*sizeof(float),
		"jea_xinterpolate:allocating aa" );
	bb =  (float *)VNDemalloc( MAX(n2,nt)*sizeof(float),
		"jea_xinterpolate:allocating bb" );
	if(linear) {
		tin=(float *)VNDemalloc(nt*sizeof(float),
			"jea_xinterpolate: allocating tin");
		for(j=0;j<nt;j++) tin[j]=j;
	}

	file=VNDtempname("suinterp");
	vnda  = V2Dop(2,500000,sizeof(float),file,ntlong,n2long);
	VNDfree(file,"file");
	file=VNDtempname("suinterp");
	vndb  = V2Dop(2,500000,sizeof(float),file,ntlong,n2long);
	VNDfree(file,"file");

	/* 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<n2;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.;
		if(deriv) {
			for(j=nt-1;j>0;j--){
				rt[j]=rt[j]-rt[j-1];
			}
			rt[0]=rt[1];
		}
		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<n2;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");
		if(linear) VNDfree(tin,"jea_xinterpolate: tin");
		return;
	}

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

	/* compute dot products and smooth over time */
	for(ix=0;ix<n2;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(len2>1){
	    for(it=0;it<nt;it++) {
		V2Dr1(vnda,it,(char *)a,112);
		V2Dr1(vndb,it,(char *)b,113);
		runav(n2,len2,a,aa);
		runav(n2,len2,aa,a);
		runav(n2,len2,b,bb);
		runav(n2,len2,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<n2;ix++) {
		ixm=((int) (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] = ((float ) it) - p[it]*wa;
				if(linear){
					intlin(nt,tin,a,a[0],a[nt-1],nt,
					       time,aa);
				}else{		
					ints8r(nt,1.0,0.,a,0.0,0.0,nt,time,aa);
				}
				for(it=0;it<nt;it++) time[it] = ((float ) it) + p[it]*wb;
				if(linear){
					intlin(nt,tin,a,a[0],a[nt-1],nt,
					       time,bb);
				}else{		
					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,n2-1,(char *)p,122);
	}else{
		V2Dw0(vndinterp,(n2-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");
	if(linear) VNDfree(tin,"jea_xinterpolate: tin");
	return;
}
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;
}

/* for graceful interrupt termination */
static void closefiles(void)
{
	efclose(headerfp);
	eremove(headerfile);
	exit(EXIT_FAILURE);
}