sumiggbzoan.c

su 的源代码库

dt		time sampling interval
ft		first time sample
f		array[nt] containing data for one ray f(t)		
nx		number of x samples
dx		x sampling interval
fx		first x sample
nz		number of z samples
dz		z sampling interval
fz		first z sample
g		array[nx][nz] in which to accumulate beam

*****************************************************************************
Output:
g		array[nx][nz] after accumulating beam
*****************************************************************************
Technical Reference:
 Hale, D., 1992, Migration by the Kirchhoff,slant stack,
 and Gaussian beam methods:
 CWP,1992 Report 121, Colorado School of Mines.
 
 Hale, D., 1992, Computational Aspects of Gaussian Beam migration:
      	CWP,1992 Report 121, Colorado School of Mines.
*****************************************************************************
Credits: CWP: Dave Hale
*****************************************************************************/
{
	int lx,lz,mx,mz,jx,jz,live,dead;
	float wmin;
	RayStep *rs=ray->rs;
	Cells *cells;
	Cell **cell;
	BeamData *bd;

	/* frequency in radians per unit time */
	wmin = 2.0*PI*fmin;

	/* random numbers used to denote live and dead cells */
	live = 1+(int)(1.0e7*franuni());
	dead = 1+(int)(1.0e7*franuni());
	
	/* number of samples per cell */
	lx = CELLSIZE;
	lz = CELLSIZE;

	/* number of cells */
	mx = 2+(nx-1)/lx;
	mz = 2+(nz-1)/lz;

	/* compute complex beam data */
	bd = beamData(wmin,nt,dt,ft,f);

	/* allocate cells */
	cells = (Cells*)alloc1(1,sizeof(Cells));
	cell = (Cell**)alloc2(mz,mx,sizeof(Cell));

	/* set information needed to set and fill cells */
	cells->nt = nt;
	cells->dt = dt;
	cells->ft = ft;
	cells->lx = lx;
	cells->mx = mx;
	cells->nx = nx;
	cells->dx = dx;
	cells->fx = fx;
	cells->lz = lz;
	cells->mz = mz;
	cells->nz = nz;
	cells->dz = dz;
	cells->fz = fz;
	cells->live = live;
	cells->dead = dead;
	cells->wmin = wmin;
	cells->lmin = lmin;
	cells->cell = cell;
	cells->ray = ray;
	cells->bd = bd;
	cells->g = g;

	/* cell closest to initial point on ray will be first live cell */
	jx = NINT((rs[0].x-fx)/dx/lx);
	jz = NINT((rs[0].z-fz)/dz/lz);
	
	/* set first live cell and its neighbors recursively */
	setCell(cells,jx,jz);

	/* accumulate beam in first live cell and its neighbors recursively */
	accCell(cells,jx,jz);

	/* free complex beam data */
	free2complex(bd->cf);
	free1((void*)bd);

	/* free cells */
	free2((void**)cells->cell);
	free1((void*)cells);
}

/* functions for internal use only */

static void xtop (float w,
	int nx, float dx, float fx, complex *g,
	int np, float dp, float fp, complex *h)
/*****************************************************************************
Slant stack for one frequency w, where slant stack is defined by

           fx+(nx-1)*dx
    h(p) =   integral   exp(-sqrt(-1)*w*p*x) * g(x) * dx
                fx
******************************************************************************
Input:
w		frequency (radians per unit time)
nx		number of x samples
dx		x sampling interval
fx		first x sample
g		array[nx] containing g(x)
np		number of p samples
dp		p sampling interval
fp		first p sample
******************************************************************************
Output:
h		array[np] containing h(p)
******************************************************************************
Notes:
The units of w, x, and p must be consistent.

Slant stack is performed via FFT and 8-point (tapered-sinc) interpolation.

The Fourier transform over time (t) is assumed to have been performed with
positive sqrt(-1)*w*t in the exponent;  if negative sqrt(-1)*w*t was used
instead, call this function with negative w.
*****************************************************************************/
{
	int nxfft,nk,nka,ix,ik,ip,lwrap;
	float dk,fk,ek,fka,k,p,phase,c,s,x,xshift,temp,*kp;
	complex czero=cmplx(0.0,0.0),*gx,*gk,*gka,*hp;

	/* number of samples required to make wavenumber k periodic */
	lwrap = 8;
	
	/* wavenumber k sampling */
	nxfft = npfa((nx+lwrap)*2);
	nk = nxfft;
	dk = 2.0*PI/(nxfft*dx);
	fk = -PI/dx;
	ek = PI/dx;
	fka = fk-lwrap*dk;
	nka = lwrap+nk+lwrap;
	
	/* allocate workspace */
	gka = alloc1complex(nka);
	gx = gk = gka+lwrap;
	hp = alloc1complex(np);
	kp = alloc1float(np);
	
	/* scale g(x) by x sampling interval dx */
	for (ix=0; ix<nx; ++ix,x+=dx) {
		gx[ix].r = dx*g[ix].r;
		gx[ix].i = dx*g[ix].i;
	}
	
	/* pad g(x) with zeros */
	for (ix=nx; ix<nxfft; ++ix)
		gx[ix].r = gx[ix].i = 0.0;
	
	/* negate every other sample so k-axis will be centered */
	for (ix=1; ix<nx; ix+=2) {
		gx[ix].r = -gx[ix].r;
		gx[ix].i = -gx[ix].i;
	}
	
	/* Fourier transform g(x) to g(k) */
	pfacc(-1,nxfft,gx);
	
	/* wrap-around g(k) to avoid interpolation end effects */
	for (ik=0; ik<lwrap; ++ik)
		gka[ik] = gk[ik+nk-lwrap];
	for (ik=lwrap+nk; ik<lwrap+nk+lwrap; ++ik)
		gka[ik] = gk[ik-lwrap-nk];
	
	/* phase shift to account for non-centered x-axis */
	xshift = 0.5*(nx-1)*dx;
	for (ik=0,k=fka; ik<nka; ++ik,k+=dk) {
		phase = k*xshift;
		c = cos(phase);
		s = sin(phase);
		temp = gka[ik].r*c-gka[ik].i*s;
		gka[ik].i = gka[ik].r*s+gka[ik].i*c;
		gka[ik].r = temp;
	}
	
	/* compute k values at which to interpolate g(k) */
	for (ip=0,p=fp; ip<np; ++ip,p+=dp) {
		kp[ip] = w*p;
		
		/* if outside Nyquist bounds, do not interpolate */
		if (kp[ip]<fk && kp[ip]<ek)
			kp[ip] = fk-1000.0*dk;
		else if (kp[ip]>fk && kp[ip]>ek)
			kp[ip] = ek+1000.0*dk;
	}
		
	/* interpolate g(k) to obtain h(p) */
	ints8c(nka,dk,fka,gka,czero,czero,np,kp,hp);
	
	/* phase shift to account for non-centered x-axis and non-zero fx */
	xshift = -fx-0.5*(nx-1)*dx;
	for (ip=0; ip<np; ++ip) {
		phase = kp[ip]*xshift;
		c = cos(phase);
		s = sin(phase);
		h[ip].r = hp[ip].r*c-hp[ip].i*s;
		h[ip].i = hp[ip].r*s+hp[ip].i*c;
	}

	/* free workspace */
	free1complex(gka);
	free1complex(hp);
	free1float(kp);
}

static BeamData* beamData (float wmin, int nt, float dt, float ft, float *f)
/*****************************************************************************
Compute filtered complex beam data as a function of real and imaginary time.
******************************************************************************
Input:
wmin		minimum frequency (in radians per unit time)
nt		number of time samples
dt		time sampling interval
ft		first time sample
f		array[nt] containing data to be filtered
*****************************************************************************
Returned:	pointer to beam data
*****************************************************************************/
{
	int ntpad,ntfft,nw,iwnyq,ntrfft,ntr,nti,nwr,it,itr,iti,iw;
	float dw,fw,dtr,ftr,dti,fti,w,ti,scale,*fa;
	complex *ca,*cb,*cfi,**cf;
	BeamData *bd;

	/* pad to avoid wraparound in Hilbert transform */
	ntpad = 25;

	/* fft sampling */
	ntfft = npfaro(nt+ntpad,2*(nt+ntpad));
	nw = ntfft/2+1;
	dw = 2.0*PI/(ntfft*dt);
	fw = 0.0;
	iwnyq = nw-1;

	/* real time sampling (oversample for future linear interpolation) */
	ntrfft = nwr = npfao(NOVERSAMPLE*ntfft,NOVERSAMPLE*ntfft+ntfft);
	dtr = dt*ntfft/ntrfft;
	ftr = ft;
	ntr = 1+(nt+ntpad-1)*dt/dtr;

	/* imaginary time sampling (exponential decay filters) */
	nti = NFILTER;
	dti = EXPMIN/(wmin*(nti-1));
	fti = 0.0;

	/* allocate space for filtered data */
	cf = alloc2complex(ntr,nti);

	/* allocate workspace */
	fa = alloc1float(ntfft);
	ca = alloc1complex(nw);
	cb = alloc1complex(ntrfft);

	/* pad data with zeros */
	for (it=0; it<nt; ++it)
		fa[it] = f[it];
	for (it=nt; it<ntfft; ++it)
		fa[it] = 0.0;

	/* Fourier transform and scale to make complex analytic signal */
	pfarc(1,ntfft,fa,ca);
	for (iw=1; iw<iwnyq; ++iw) {
		ca[iw].r *= 2.0;
		ca[iw].i *= 2.0;
	}

	/* loop over imaginary time */
	for (iti=0,ti=fti; iti<nti; ++iti,ti+=dti) {

		/* apply exponential decay filter */
		for (iw=0,w=fw; iw<nw; ++iw,w+=dw) {
			scale = exp(w*ti);
			cb[iw].r = ca[iw].r*scale;
			cb[iw].i = ca[iw].i*scale;
		}

		/* pad with zeros */
		for (iw=nw; iw<nwr; ++iw)
			cb[iw].r = cb[iw].i = 0.0;

		/* inverse Fourier transform and scale */
		pfacc(-1,ntrfft,cb);
		cfi = cf[iti];
		scale = 1.0/ntfft;
		for (itr=0; itr<ntr; ++itr) {
			cfi[itr].r = scale*cb[itr].r;
			cfi[itr].i = scale*cb[itr].i;
		}
	}

	/* free workspace */
	free1float(fa);
	free1complex(ca);
	free1complex(cb);

	/* return beam data */
	bd = (BeamData*)alloc1(1,sizeof(BeamData));
	bd->ntr = ntr;
	bd->dtr = dtr;
	bd->ftr = ftr;
	bd->nti = nti;
	bd->dti = dti;
	bd->fti = fti;
	bd->cf = cf;
	return bd;
}

static void setCell (Cells *cells, int jx, int jz)
/*****************************************************************************
Set a cell by computing its Gaussian beam complex time and amplitude.  
If the amplitude is non-zero, set neighboring cells recursively.
******************************************************************************
Input:
cells		pointer to cells
jx		x index of the cell to set
jz		z index of the cell to set
******************************************************************************
Notes:
To reduce the amount of memory required for recursion, the actual 
computation of complex time and amplitude is performed by the cellTimeAmp() 
function, so that no local variables are required in this function, except
for the input arguments themselves.
*****************************************************************************/
{
	/* if cell is out of bounds, return */
	if (jx<0 || jx>=cells->mx || jz<0 || jz>=cells->mz) return;

	/* if cell is live, return */
	if (cells->cell[jx][jz].live==cells->live) return;

	/* make cell live */
	cells->cell[jx][jz].live = cells->live;

	/* compute complex time and amplitude.  If amplitude is
	 * big enough, recursively set neighboring cells. */
	if (cellTimeAmp(cells,jx,jz)) {
		setCell(cells,jx+1,jz);
		setCell(cells,jx-1,jz);
		setCell(cells,jx,jz+1);
		setCell(cells,jx,jz-1);
	}
}

static int cellTimeAmp (Cells *cells, int jx, int jz)
/*****************************************************************************
Compute complex and time and amplitude for a cell.
******************************************************************************
Input:
cells		pointer to cells
jx		x index of the cell to set
jz		z index of the cell to set

*****************************************************************************
Returned:	1 if Gaussian amplitude is significant, 0 otherwise
*****************************************************************************
Technical Reference:
	Hale, D., 1992, Migration by the Kirchhoff, 
 	slant stack, and Gaussian beam methods:
       CWP,1992 Report 121, Colorado School of Mines.
 
 	Cerveny, V., 1972, Seismic rays and ray intensities 
 	in inhomogeneous anisotropic media: 
 	Geophys. J. R. Astr. Soc., 29, 1--13.
*****************************************************************************
Credits: CWP: Tariq Alkhalifah,  based on MIGGBZO by Dave Hale
*****************************************************************************/
{
	int lx=cells->lx,lz=cells->lz;
	float dx=cells->dx,fx=cells->fx,dz=cells->dz,fz=cells->fz,
		wmin=cells->wmin,lmin=cells->lmin;
	Ray *ray=cells->ray;
	Cell **cell=cells->cell;
	int irs,kmah;
	float tmax,xc,zc,t0,x0,z0,x,z,tr,ti,ar,ai,
		v,q1,p1,q2,p2,e,es,scale,phase,vvmr,vvmi,c,s,ov,
		px,pz,pzpz,pxpx,pxpz,dvdx,dvdz,dvds,dvdn,
		wxxr,wxzr,wzzr,wxxi,wxzi,wzzi;
	RayStep *rs;

	/* maximum time */
	tmax = ray->rs[ray->nrs-1].t;

	/* cell coordinates */
	xc = fx+jx*lx*dx;
	zc = fz+jz*lz*dz;

	/* ray step nearest to cell */
	irs = nearestRayStep(ray,xc,zc);
	rs = &(ray->rs[irs]);
		
	/* ray time and coordinates */
	t0 = rs->t;