taup.c

su 的源代码库

	/* allocate VND 2-D space */
	fname = VNDtempname("radontemp");
	vnda = V2Dop(2,1000000,sizeof(float),fname,ntfft+2,(long) nmax); 
	VNDfree(fname,"invslant:fname 1");

	/* compute offsets */
	for (ix=0; ix<nx; ix++) {

		/* compute offsets */
		offset[ix] = fx+ix*dx;	

		/* get g(x) (offset values) depending on type of transform */
		g[ix] = gofx (igopt, offset[ix], interoff, depthref);
	}
	/* copy input traces to ntfft VND array */
	for (ip=0; ip<np; ip++) {
		V2Dw0(vnda,ip,(char *)in_traces[ip],1002);
	} 
		
	nw=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)));

	/* allocate additional VND space */
	crt=(complex *)VNDemalloc(nmax,
		"inverse_p_transform:crt");
	rt=(float *)crt;
	ctemp=(complex *)VNDemalloc(np*sizeof(complex),
		"inverse_p_transform:ctemp");

	fac=1./ntfft;
	/* compute forward time (tau) to frequency Fourier transform */
	for(ip=0;ip<np;ip++) {

		V2Dr0(vnda,ip,(char *)rt,301);
		for (it=nt; it<ntfft; it++) rt[it]=0.0;
		for(it=0;it<ntfft;it++) rt[it]*=fac;
		pfarc(1,ntfft,rt,crt);
		V2Dw0(vnda,ip,(char *)crt,302);
	}
	VNDr2c(vnda);

	fac=1./np;
	for(iw=0;iw<nw;iw++) {
		wa=freqweight(iw,df,f1,f2);
		if(wa>0.) {
			w=iw*dw;
			V2Dr1(vnda,iw,(char *)crt,303);
			if(wa<1.) {
				for(ip=0;ip<np;ip++) crt[ip]=crmul(crt[ip],wa);
			}
			for(ip=0;ip<np;ip++) ctemp[ip]=crt[ip];
			for(ix=0;ix<nx;ix++) {
			    rsum = isum = 0.;
			    for(ip=0;ip<np;ip++) {
				p = pmin + ip*dp;
				tr = cos(w*p*g[ix]);
				ti = sin(w*p*g[ix]);
				dr = ctemp[ip].r;
				di = ctemp[ip].i;
				rsum += tr*dr - ti*di;
				isum += tr*di + ti*dr;
			    }
			    crt[ix].r   = fac*rsum;
			    crt[ix].i   = fac*isum;
			}
		}else{
			for(ix=0;ix<nx;ix++) crt[ix]=czero;
		}
		V2Dw1(vnda,iw,(char *)crt,304);
	}
	for(ix=0;ix<nx;ix++) {
		V2Dr0(vnda,ix,(char *)crt,305);
		pfacr(-1,ntfft,crt,rt);

		/* copy output data */
		for (it=0; it<nt; it++)
			out_traces[ix][it]=rt[it];
	}

	/* free allocated space */
	VNDcl(vnda,1);
	VNDfree(crt,"inverse_p_transform: crt");
	VNDfree(ctemp,"inverse_p_transform: ctemp");
	VNDfree(offset,"inverse_p_tarnsform: offset");
	VNDfree(g,"inverse_p_transform: g");
	return;
}

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 = modified Foster/Mosher pseudo hyperbolic option
			3 = linear tau-p
			4 = linear tau-p using absolute value of 
				offset
			5 = original Foster/Mosher pseudo hyperbolic option
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));
	}
	if(igopt==5) {
		return(sqrt(refdepth*refdepth+offset*offset)-refdepth);
	} else {
		return 0.0;
	}
}


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

void compute_r( float w, int nx, float *g, int np, float dp, complex *r)
/******************************************************************************
Compute the top row of the Hermitian Toeplitz Matrix
			+
		  R = B B

		  i w p g(x)
where B = (1/np) e	    for equal increments in p as

     +           -i w p g(x)
and B = (1/nx) e 

as used for the Discrete Radon Transform computation for
linear or parabolic tau-p.


		 nx-1	i w j dp g(x )
r[j] = 1/(nx*np) Sum	e	    k
		 k=0
						  2
g(x ) is initialized to  x  for linear tau-p or x   for the parabolic transform
   k		          k		         k
prior to calling this routine.  The use of g is intended to emphasize that the
spatial locations do not have to be equally spaced for either method.
In general, this routine can be called for g specified as any function
of spatial position only.  For a more general function of x, dp will
not correspond to an increment in slowness or slowness squared but
rather to a more general parameter.

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

float w	input as angular frequency component of interest
int   nx      number of spatial positions stored in g
float g[]     spatial function for this Radon Transform
int   np      number of slowness (or slowness squared) components
float dp      increment in slownes (or slowness squared)
float r[]     output vector of { real r0, imaginary r0, real r1, 
	      imaginary r1, ...}
*******************************************************************************
Author: John Anderson (visitor to CSM from Mobil) Spring 1993
******************************************************************************/
{
	int j,k;
	float rsum, isum, fac;
	fac= 1./(nx*np);
	for(j=0;j<np;j++) {
		rsum=0.;
		isum=0.;
		for(k=0;k<nx;k++) {
			rsum = rsum+cos( w*j*dp*g[k] );
			isum = isum+sin( w*j*dp*g[k] );	
		}
		r[j].r    = fac*rsum;
		r[j].i    = fac*isum;
	}
}

void compute_rhs( float w, int nx, float *g, complex *data, int np, 
		float pmin, float dp, complex *rhs)
/******************************************************************************
				     +
Compute the right-hand-side vector  B  data(x)

	+	    -i w p g(x)
where B   = (1/nx) e	        for equal increments in p as
used for the Discrete Radon Transform computation for
linear or parabolic tau-p.

Function parameters:

float w	input angular frequency of interest
int   nx	number of spatial positions ( defines length of g and data )
float g[]      spatial function corresponding to spatial locations of data
complex data[] data as a function of spatial position for a single
		angular frequency w as complex values 
int   np	number of output slownesses p (may be slowness squared
		or a more general function)
float pmin     starting value of output p
float dp	increment in output p
complex rhs[]  np complex values for the result 
*******************************************************************************
Author: John Anderson (visitor to CSM from Mobil) Spring 1993
******************************************************************************/
{
	int ip, ix;
	float p, rsum, isum, dr, di, tr, ti, fac;
	fac=1./nx;
	for(ip=0;ip<np;ip++) {
		p = pmin + ip*dp;
		rsum = isum = 0.;
		for(ix=0;ix<nx;ix++) {
			tr = cos(w*p*g[ix]);
			ti = -sin(w*p*g[ix]);
			dr = data[ix].r;
			di = data[ix].i;
			rsum += tr*dr - ti*di;
			isum += tr*di + ti*dr;
		}
		rhs[ip].r   = fac*rsum;
		rhs[ip].i   = fac*isum;
	}
}

int ctoep( int n, complex *r, complex *a, complex *b, complex *f, complex *g )
/*****************************************************************************	
Complex Hermitian Toeplitz Solver for

N-1
Sum  R	     A  = B      for i=0,1,2,...,N-1
j=0   (i-j)   j    i

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


   R0  R1  R2  R3	A0	     B0

     *
   R1  R0  R1  R2	A1	     B1
				=    
     *   *
   R2  R1  R0  R1	A2	     B2

     *   *   *
   R3  R2  R1  R0	A3	     B3

and


   R0  R1  R2  R3	F0	     1

     *
   R1  R0  R1  R2	F1	     0
				=    
     *   *
   R2  R1  R0  R1	F2	     0

     *   *   *
   R3  R2  R1  R0	F3	     0


*******************************************************************************	
where the function parameters are defined by

n     dimension of system
*r    provides the top row of the Hermitian Toeplitz matrix R 
*a    returns the complex solution vector A
*b    input as complex vector B (not changed during call)
*f    returns the complex spiking filter F
      (may be needed later for Simpson's sideways recursion
      if do search for optimum filter lag)
*g    work space of length n complex values

The function value returns as the number of successfully
computed complex filter coefficients (up to n) if successful or
0 if no coefficients could be computed.
*******************************************************************************	
Author: John Anderson (visitor to CSM from Mobil) Spring 1993
******************************************************************************/

{
	float er, ei, vr, vi, cr, ci, vsq;
	int j;	  	/*  for the jth iteration, j=0,n-1 	*/
	int k;			/*  for the kth component, k=0,j-1 	*/
	int jmk;		/*  j-k 				*/
	if (r[0].r==0.) return 0;

	f[0].r = 1.0/r[0].r;
	f[0].i = 0.;
	a[0].r = b[0].r/r[0].r;
	a[0].i = b[0].i/r[0].r;
	vr=1.;
	vi=0.;
	vsq=1.;

	for(j=1;j<n;j++) {     	/* iteration loop for iteration j	*/
	/*  	Compute spiking filter that outputs {v,0,0,...} 
		for this iteration step j			*/
		f[j].r=0.;
		f[j].i=0.;
		er=ei=0.;
		for(k=0;k<j;k++) {
			jmk=j-k;
			er+=r[jmk].r*f[k].r+r[jmk].i*f[k].i;
			ei+=r[jmk].r*f[k].i-r[jmk].i*f[k].r;
		}
		cr  = (er*vr - ei*vi)/vsq;
		ci  = (er*vi + ei*vr)/vsq;
		vr  = vr - (cr*er+ci*ei);
		vi  = vi + (cr*ei-ci*er);
		vsq =  vr*vr + vi*vi;
		if (vsq <= 0.) break;
		for(k=0;k<=j;k++) {
			jmk=j-k;
			g[k].r = f[k].r - cr*f[jmk].r - ci*f[jmk].i;
			g[k].i = f[k].i + cr*f[jmk].i - ci*f[jmk].r;		
		}
		for(k=0;k<=j;k++) {
			f[k]=g[k];
		}

		/*  Compute shaping filter for this iteration */
		a[j].r=0.;
		a[j].i=0.;
		er=ei=0.;
		for(k=0;k<j;k++) {
			jmk=j-k;
			er+=r[jmk].r*a[k].r+r[jmk].i*a[k].i;
			ei+=r[jmk].r*a[k].i-r[jmk].i*a[k].r;
		}
		er  = er-b[j].r;
		ei  = ei-b[j].i;
		cr  = (er*vr - ei*vi)/vsq;
		ci  = (er*vi + ei*vr)/vsq;
		for(k=0;k<=j;k++) {
			jmk=j-k;
			a[k].r += - cr*f[jmk].r - ci*f[jmk].i;
			a[k].i += + cr*f[jmk].i - ci*f[jmk].r;		
		}	
	}

	/* Properly normalize the spiking filter so that R F = {1,0,0,...} */
	/* instead of {v,0,0,...}.  To be accurate, recompute vr,vi,vsq */ 
	vr=vi=0.;
	for(k=0;k<j;k++) {
		vr+=r[k].r*f[k].r-r[k].i*f[k].i;
		vi+=r[k].r*f[k].i+r[k].i*f[k].r;
	}

	vsq = vr*vr + vi*vi;