suradon.c

su 的源代码库

				nxout= nx;
			}
						VNDmemchk(rt,"rt 01 e");
			erewind(headerfp);
			for(ix=0;ix<nxout;ix++) {
				if(ix<nx) efread(&tro,HDRBYTES,1,headerfp);
				V2Dr0(vndresult,ix,(char *)rt,1008);
				for(j=0;j< nt;j++) tro.data[j]=rt[j];
				icount++;
				tro.tracl = icount;
				tro.tracr = ix+1;
				if( choose==0) {
					tro.f2=1000.*( pmin+ ix* dp)*
						gofx( igopt, offref,
						     intercept_off, depthref);
					tro.d2=1000.*dp*gofx( igopt, offref,
							     intercept_off,
							     depthref);
				}
				puttr(&tro);
			}
		   }else{  /* pass input data unmodified */
			erewind(headerfp);
			efread(&tro,HDRBYTES,1,headerfp);
			V2Dr0(vndorig,0,(char *)rt,1009);
			for(j=0;j< nt;j++) tro.data[j]=rt[j];
			icount++;
			tro.tracl = icount;
			tro.tracr = 1;
			puttr(&tro);
		   }
		   oldcmp=hdrwd.i;
		   erewind(headerfp);
		   nx=0;
		}
		/* save trace and header */
		efwrite(&tr,HDRBYTES,1,headerfp);
		gethval(&tr,offindex,&hdrwd);
		offset[nx]=hdrwd.i;
		g[nx]=gofx(igopt,offset[nx],intercept_off,depthref);
		mutetime[nx]=tr.muts;
		V2Dw0(vndorig,nx,(char *)tr.data,1010);
		nx++;
		if(ieod) {
			iend++;
		}else{
			if(gettr(&tr)) {
				iend=0;
			}else{
				iend=1;
				ieod=1;
			}
		}
	}while(iend<2);

/* free memory and temporary file space */
    	VNDfree(offset,"suradon_main: offset"); 
	VNDfree(rt,"suradon_main: rt");
	VNDfree(trace,"suradon_main: trace");  
	VNDfree(xin,"suradon_main: xin");
	VNDfree(g,"suradon_main: g");
	VNDfree(gg,"suradon_main: gg");
	VNDfree(mutetime,"suradon_main: mutetime");
	fclose(headerfp);
	remove(headerfile);
	VNDfree(headerfile,"suradon_main: headerfile");
	VNDcl(vndorig,1);
	VNDcl(vndinterp,1);
	VNDcl(vndresult,1);
	if(VNDtotalmem()!=0) {
		fprintf(stderr,"Warning, not all of the VND memory \n");
		fprintf(stderr,"has been freed and checked for overruns\n");
		fprintf(stderr,"Total VND memory at end of job = %ld\n",
		VNDtotalmem());
	}
	return(CWP_Exit());
}

static void forward_p_transform(VND *vnda,VND *vndb,int nx, int nt, float *g,
	float dt, int ntfft, int np, float pmin, float dp,
	float *mutetime, float *offset, int nk,float f1, float f2,
	float prewhite)
/*******************************************************************
do forward generalized radon transform

******************************************************************
Function parameters:

VND *vnda	input data in time-space domain
VND *vndb	output data in tau-p domain
int nx		number of input traces
int nt		number of intput time samples
float *g	g[ix]=offset**2 for parabolic transform
float dt	sample rate in seconds
int ntfft	length of time fft and output tau-p data in samples
int np		number of generalized ray parameters
float pmin	minimum generalized ray parameter
float dp	generalized ray parameter increment
float *mutetime array of mute times in ms
float *offset   array of offsets (ignored)
int nk          number of offset ranges to transform separately
		through the mute zone
float f1        max freq without taper
float f2        max non-zero freq component
prewhite	0.01 means prewhiten 1 percent

key assumption: offsets are sorted to increase with index
*******************************************************************
Author: John Anderson (visitor to CSM from Mobil) Spring 1993
*******************************************************************/
{
	int ix,ip,it,iw,j,ntfftny,k,nxx,nxxinc,ik,ik2;
	size_t nmax;
	float *rt,*rrt,*kindex=NULL,*tindex=NULL,w,dw,rsum,fac,wa,wb,rk[2],rit[2],df;
	complex czero,*crt,*ccrt,*r,*rhs,*wrk1,*wrk2,*wrk3,*wrk4;
	VND *vndc;
	char *fname;
	
	fac=1./ntfft;
	ntfftny=1+ntfft/2;
	df=1./(ntfft*dt);
	dw=2.*PI*df;
	nmax=MAX(vndb->N[0],vndb->N[1]);
	nxxinc=1+(nx-1)/nk;
	czero.r=czero.i=0.;

	if(nk>1) {
/* allocate file space and build a set of (mute time, group index) pairs */
		fname=VNDtempname("radontmp");
		vndc = V2Dop(2,1000000,sizeof(complex),fname,
				ntfftny,np*nk);
		VNDfree(fname,"forward_p_transform: fname");
		kindex = offset; /* dummy assignment */
		kindex=(float *)VNDemalloc(nk*sizeof(float),
			"forward_p_transform:kindex");
		tindex=(float *)VNDemalloc(nk*sizeof(float),
			"forward_p_transform:tindex");
		for(k=0;k<nk;k++) {
			kindex[k]=k;
			nxx=MIN(nx,(k+0.75)*nxxinc);
			tindex[k]=0.001*mutetime[nxx]/dt;
		}	
	}else{
		vndc = vndb;
	}

	crt=(complex *)VNDemalloc(nmax*sizeof(complex),
		"forward_transform:crt");
	rt=(float *)crt;
	ccrt=(complex *)VNDemalloc(MAX((nk+1)*np,vndb->N[1])*sizeof(complex),
		"forward_transform:ccrt");
	rrt=(float *)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");


/* do forward time to frequency fft */
	for(ix=0;ix<nx;ix++) {
		V2Dr0(vnda,ix,(char *)rt,201);
		for(it=0;it<nt;it++) rt[it]*=fac;
		for(j=nt;j<ntfft;j++) rt[j]=0.;
		pfarc(1,ntfft,rt,crt);
		V2Dw0(vndb,ix,(char *)crt,202);
	}
	VNDr2c(vndb);

/* do radon transform, frequency by frequency, for multiple spatial windows */
	for(iw=0;iw<ntfftny;iw++) {
		wa=freqweight(iw,df,f1,f2);
		if(wa>0.) {
		    w=iw*dw;
		    V2Dr1(vndb,iw,(char *)crt,203);
		    if(wa<1.) {
		    	for(ix=0;ix<nx;ix++) crt[ix]=crmul(crt[ix],wa);
		    }
		    for(k=0;k<nk;k++) {
			nxx=MIN(nx,(k+1)*nxxinc);
			compute_rhs(w,nxx,g,crt,np,pmin,dp,rhs);
			compute_r(w,nxx,g,np,dp,r);
			r[0].r *= (1.+prewhite);
			for(rsum=0.,j=1;j<np;j++) 
				rsum +=sqrt(r[j].r*r[j].r + r[j].i*r[j].i);
			rsum=rsum/r[0].r;  
			if (rsum>1.+np/5) { 
				j=ctoephcg(np/7,np,r,&ccrt[k*np],rhs,
					wrk1,wrk2,wrk3,wrk4);
			}else{
				j=ctoep(np,r,&ccrt[k*np],rhs,wrk1,wrk2);
			}
		    }
		}else{
		    for(ip=0;ip<np*nk;ip++) ccrt[ip]=czero;		
		}
		V2Dw1(vndc,iw,(char *)ccrt,204);
	}

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

/* merge appropriate tau-p transform for each window using mute zone
   information */
	if(nk>1) {
		VNDc2r(vndb);
		for(it=0;it<ntfft;it++) {
			rit[0]=it;
			intlin(nk,tindex,kindex,kindex[0],kindex[nk-1],1,rit,rk);
			ik=rk[0];
			ik2=MIN(ik+1,nk-1);
			wb=rk[0]-ik;
			wa=1-wb;
			V2Dr1(vndc,it,(char *)rrt,207);
			for(ip=0;ip<np;ip++) {
				rt[ip]=wa*rrt[ik*np+ip]+ wb*rrt[ik2*np+ip];
			}
			V2Dw1(vndb,it,(char *)rt,208);		
		}
		VNDcl(vndc,1);
		VNDfree(tindex,"forward_p_transform: tindex");
		VNDfree(kindex,"forward_p_transform: kindex");
	}
	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");	
	return;
}
static float gofx(int igopt, float offset, float intercept_off,float refdepth)
/*******************************************************************
return g(x) for various options
******************************************************************
Function parameters:

int igopt		1 = parabolic transform
			2 = Foster/Mosher pseudo hyperbolic option
			3 = linear tau-p
			4 = linear tau-p using absolute value of 
				offset
float offset		offset in m
float intercept_off	offset corresponding to intercept time
float refdepth		reference depth in m for igopt=2
*******************************************************************
Author: John Anderson (visitor to CSM from Mobil) Spring 1993
*******************************************************************/
{
	offset=offset-intercept_off;
	if(igopt==1) {
		return(offset*offset);
	}
	if(igopt==2) {
		return( sqrt(refdepth*refdepth + offset*offset) );
	}
	if(igopt==3) {
		return(offset);
	}
	if(igopt==4) {
		return(fabs(offset));
	}
	return(offset);
}
static float freqweight(int j, float df, float f1, float f2)
/*******************************************************************
return weight for each frequency
******************************************************************
Function parameters:

int j		freq index
float df	freq increment
float f1	taper off freq
float f2	freq beyond which all components are zero
*******************************************************************
Author: John Anderson (visitor to CSM from Mobil) Spring 1993
*******************************************************************/
{
	float w;
	float f=j*df;
	if(f<=f1) return (1.);
	if(f>=f2) return (0.);
	w = (f2-f)/(f2-f1);
	return (w);
}

static void taupmute(int ip,int ipa,int ipb,int nt, int ltap, float *rt)
/*******************************************************************
do simple tau-p mute to elliminate multiples
******************************************************************
Function parameters:

int ip		current ray parameter index
int ipa		max ray parameter primary  pick at maximum time
int ipb		max ray parmater primary pick at minimum time
int nt		number of time samples
int ltap	length of mute taper in samples
float rt[nt]	tau-p data for all tau values

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

{
	int j,k;
	float w;
	if(ip>=ipb) return;
	if(ip<=ipa) {
		for(k=0;k<nt;k++) rt[k]=0;
		return;
	}
	w=MAX(ipb-ipa,1);
	w=(ipb-ip)/w;
	j=w*nt;
	for(k=0;k<j;k++) rt[k]=0.;
	for(k=j;k<MIN(nt,j+ltap);k++) rt[k]*=(k-j)/ltap;
}

static void inverse_p_transform(VND *vnda,int nx, float *g, float dt, 
	int ntfft, int np, float pmin, float dp, int ip1, 
	float f1, float f2)
/*******************************************************************
do inverse generalized radon transform

******************************************************************
Function parameters:

VND *vnda	output data in tau-p domain and data in time-space domain
int nx		number of output traces
int nt		number of output time samples
float *g	g[ix]=offset**2 for parabolic transform
float dt	sample rate in seconds
int ntfft	length of time fft and input tau-p data in samples
int np		number of generalized ray parameters
float pmin	minimum generalized ray parameter
float dp	generalized ray parameter increment
int ip1		starting generalized ray parameter index to compute
		(use as 0 for full inversion, positive integer to
		just invert multiples)
float f1        max freq without taper
float f2        max non-zero freq component
*******************************************************************
Author: John Anderson (visitor to CSM from Mobil) Spring 1993
*******************************************************************/
{
	int ip,iw,ntfftny,ix,it;
	size_t nmax;
	float w,dw,p,rsum,isum,dr,di,tr,ti,fac,wa,df;
	float *rt;
	complex *crt, *ctemp,czero;
	
	ntfftny=1+ntfft/2;
	df=1./(ntfft*dt);
	dw=2.*PI*df;
	czero.r=czero.i=0.;

	nmax=MAX(vnda->N[0],2*vnda->N[1])*vnda->NumBytesPerNode;
	nmax=MAX(nmax,nx*sizeof(complex));