sumiggbzoan.c

su 的源代码库

	x0 = rs->x;
	z0 = rs->z;
	
	/* real and imaginary parts of v*v*m = v*v*p/q */
	v = rs->v;
	q1 = rs->q1;
	p1 = rs->p1;
	q2 = rs->q2;
	p2 = rs->p2;
	e = wmin*lmin*lmin;
	es = e*e;
	scale = v*v/(q2*q2+es*q1*q1);
	vvmr = (p2*q2+es*p1*q1)*scale;
	vvmi = -e*scale;
	
	/* components of slowness vector, px and pz */
	c = rs->c;
	s = rs->s;
	ov = 1/v;
	px = s*ov;
	pz = c*ov;
	pzpz = pz*pz;
	pxpx = px*px;
	pxpz = px*pz;
	
	/* velocity derivatives along tangent and normal */
	dvdx = rs->dvdx;
	dvdz = rs->dvdz;
	dvds = s*dvdx+c*dvdz;
	dvdn = c*dvdx-s*dvdz;
	
	/* real part of W matrix */
	wxxr = pzpz*vvmr-2.0*pxpz*dvdn-pxpx*dvds;
	wxzr = (pxpx-pzpz)*dvdn-pxpz*(vvmr+dvds);
	wzzr = pxpx*vvmr+2.0*pxpz*dvdn-pzpz*dvds;
	
	/* imaginary part of W matrix */
	wxxi = pzpz*vvmi;
	wxzi = -pxpz*vvmi;
	wzzi = pxpx*vvmi;
	
	/* vector from ray to cell */
	x = xc-x0;
	z = zc-z0;
	
	/* real and imaginary parts of complex time */
	tr = t0+(px+0.5*wxxr*x)*x+(pz+wxzr*x+0.5*wzzr*z)*z;
	ti = 0.5*wxxi*x*x+(wxzi*x+0.5*wzzi*z)*z;

	/* real and imaginary parts of complex amplitude */
	kmah = rs->kmah;
	scale = pow(es/(q2*q2+es*q1*q1),0.25);
	phase = 0.5*(atan2(q2,q1*e)+2.0*PI*((kmah+1)/2));
	ar = scale*cos(phase);
	ai = scale*sin(phase);
		/*fprintf(stderr,"tr=%f ti=%f ar=%f ai=%f\n",tr,ti,ar,ai);*/

	/* set cell parameters */
	cell[jx][jz].tr = tr;
	cell[jx][jz].ti = ti;
	cell[jx][jz].ar = ar;
	cell[jx][jz].ai = ai;

	/* return 1 if Gaussian amplitude is significant, 0 otherwise */
	return (wmin*ti>EXPMIN && tr<=tmax)?1:0;
}

static void accCell (Cells *cells, int jx, int jz)
/*****************************************************************************
Accumulate the contribution of a Gaussian beam in a cell and its
neighbors, recursively.
******************************************************************************
Input:
cells		pointer to cells
jx		x index of the cell to fill
jz		z index of the cell to fill
******************************************************************************
Notes:
To reduce the amount of memory required for recursion, the actual
accumulation is performed by function cellBeam(), 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-1 || jz<0 || jz>=cells->mz-1) return;

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

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

	/* if upper-left corner of cell is not live, return */
	if (cells->cell[jx][jz].live!=cells->live) return;

	/* if remaining three corners of cell are live */
	if (cells->cell[jx+1][jz].live==cells->live &&
		cells->cell[jx][jz+1].live==cells->live &&
		cells->cell[jx+1][jz+1].live==cells->live) {

		/* accumulate beam in cell */
		cellBeam(cells,jx,jz);
	}

	/* recursively accumulate in neighboring cells */
	accCell(cells,jx+1,jz);
	accCell(cells,jx-1,jz);
	accCell(cells,jx,jz+1);
	accCell(cells,jx,jz-1);
}

static void cellBeam (Cells *cells, int jx, int jz)
/*****************************************************************************
Accumulate Gaussian beam for one cell.
******************************************************************************
Input:
cells		pointer to cells
jx		x index of the cell in which to accumulate beam
jz		z index of the cell in which to accumulate beam
*****************************************************************************/
{
	int lx=cells->lx,lz=cells->lz,nx=cells->nx,nz=cells->nz;
	float **g=cells->g;
	Cell **cell=cells->cell;
	BeamData *bd=cells->bd;
	int ntr=bd->ntr,nti=bd->nti;
	float dtr=bd->dtr,ftr=bd->ftr,dti=bd->dti,fti=bd->fti;
	complex **cf=bd->cf;
	int kxlo,kxhi,kzlo,kzhi,kx,kz,itr,iti;
	float odtr,odti,t00r,t01r,t10r,t11r,t00i,t01i,t10i,t11i,
		a00r,a01r,a10r,a11r,a00i,a01i,a10i,a11i,
		tx0r,tx1r,tx0i,tx1i,ax0r,ax1r,ax0i,ax1i,
		txzr,txzi,axzr,axzi,
		dtx0r,dtx0i,dtx1r,dtx1i,dax0r,dax0i,dax1r,dax1i,
		dtxzr,dtxzi,daxzr,daxzi,xdelta,zdelta,
		trn,tin,trfrac,mtrfrac,tifrac,mtifrac,
		cf0r,cf0i,cf1r,cf1i,cfr,cfi;
	complex *cf0,*cf1;

	/* inverse of time sampling intervals */
	odtr = 1.0/dtr;
	odti = 1.0/dti;

	/* complex time and amplitude for each corner */
	t00r = cell[jx][jz].tr;
	t01r = cell[jx][jz+1].tr;
	t10r = cell[jx+1][jz].tr;
	t11r = cell[jx+1][jz+1].tr;
	t00i = cell[jx][jz].ti;
	t01i = cell[jx][jz+1].ti;
	t10i = cell[jx+1][jz].ti;
	t11i = cell[jx+1][jz+1].ti;
	a00r = cell[jx][jz].ar;
	a01r = cell[jx][jz+1].ar;
	a10r = cell[jx+1][jz].ar;
	a11r = cell[jx+1][jz+1].ar;
	a00i = cell[jx][jz].ai;
	a01i = cell[jx][jz+1].ai;
	a10i = cell[jx+1][jz].ai;
	a11i = cell[jx+1][jz+1].ai;

	/* x and z samples for cell */
	kxlo = jx*lx;
	kxhi = kxlo+lx;
	if (kxhi>nx) kxhi = nx;
	kzlo = jz*lz;
	kzhi = kzlo+lz;
	if (kzhi>nz) kzhi = nz;

	/* fractional increments for linear interpolation */
	xdelta = 1.0/lx;
	zdelta = 1.0/lz;

	/* increments for times and amplitudes at top and bottom of cell */
	dtx0r = (t10r-t00r)*xdelta;
	dtx1r = (t11r-t01r)*xdelta;
	dtx0i = (t10i-t00i)*xdelta;
	dtx1i = (t11i-t01i)*xdelta;
	dax0r = (a10r-a00r)*xdelta;
	dax1r = (a11r-a01r)*xdelta;
	dax0i = (a10i-a00i)*xdelta;
	dax1i = (a11i-a01i)*xdelta;

	/* times and amplitudes at top-left and bottom-left of cell */
	tx0r = t00r;
	tx1r = t01r;
	tx0i = t00i;
	tx1i = t01i;
	ax0r = a00r;
	ax1r = a01r;
	ax0i = a00i;
	ax1i = a01i;

	/* loop over x samples */
	for (kx=kxlo; kx<kxhi; ++kx) {

		/* increments for time and amplitude */
		dtxzr = (tx1r-tx0r)*zdelta;
		dtxzi = (tx1i-tx0i)*zdelta;
		daxzr = (ax1r-ax0r)*zdelta;
		daxzi = (ax1i-ax0i)*zdelta;

		/* time and amplitude at top of cell */
		txzr = tx0r;
		txzi = tx0i;
		axzr = ax0r;
		axzi = ax0i;

		/* loop over z samples */
		for (kz=kzlo; kz<kzhi; ++kz) {

			/* index of imaginary time */
			iti = tin = (txzi-fti)*odti;
			if (iti<0 || iti>=nti-1) continue;

			/* pointers to left and right imaginary time samples */
			cf0 = cf[iti];
			cf1 = cf[iti+1];

			/* imaginary time linear interpolation coefficients */
			tifrac = tin-iti;
			mtifrac = 1.0-tifrac;

			/* index of real time */
			itr = trn = (txzr-ftr)*odtr;
			if (itr<0 || itr>=ntr-1) continue;

			/* real time linear interpolation coefficients */
			trfrac = trn-itr;
			mtrfrac = 1.0-trfrac;

			/* real and imaginary parts of complex beam data */
			cf0r = mtrfrac*cf0[itr].r+trfrac*cf0[itr+1].r;
			cf1r = mtrfrac*cf1[itr].r+trfrac*cf1[itr+1].r;
			cfr = mtifrac*cf0r+tifrac*cf1r;
			cf0i = mtrfrac*cf0[itr].i+trfrac*cf0[itr+1].i;
			cf1i = mtrfrac*cf1[itr].i+trfrac*cf1[itr+1].i;
			cfi = mtifrac*cf0i+tifrac*cf1i;

			/* accumulate beam */
			g[kx][kz] += axzr*cfr-axzi*cfi;

			/* increment time and amplitude */
			txzr += dtxzr;
			txzi += dtxzi;
			axzr += daxzr;
			axzi += daxzi;
		}

		/* increment times and amplitudes at top and bottom of cell */
		tx0r += dtx0r;
		tx1r += dtx1r;
		tx0i += dtx0i;
		tx1i += dtx1i;
		ax0r += dax0r;
		ax1r += dax1r;
		ax0i += dax0i;
		ax1i += dax1i;
	}
}

/*****************************************************************
* Functions for Ray tracing					 *	
*****************************************************************/

/* circle for efficiently finding nearest ray step */
typedef struct CircleStruct {
	int irsf;               /* index of first raystep in circle */
	int irsl;               /* index of last raystep in circle */
	float x;                /* x coordinate of center of circle */
	float z;                /* z coordinate of center of circle */
	float r;                /* radius of circle */
} Circle;

/* functions defined and used internally */
static Circle *makeCircles (int nc, int nrs, RayStep *rs);

Ray *makeRay (float x0, float z0, float a0, int nt, float dt, float ft,
	float **a1111xz, float **a3333xz, float **a1133xz, float **a1313xz, 
	float **a1113xz, float **a3313xz, 
	int nx, float dx, float fx, int nz, float dz, float fz)
/*****************************************************************************
Trace a ray for uniformly sampled v(x,z).
******************************************************************************
Input:
x0		x coordinate of takeoff point
z0		z coordinate of takeoff point
a0		takeoff angle (radians)
nt		number of time samples
dt		time sampling interval
ft		first time sample
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
a1111		array[nx][nz] of uniformly sampled density normalized elastic coef.
a3333		array[nx][nz] of uniformly sampled density normalized elastic coef.
a1133           array[nx][nz] of uniformly sampled density normalized elastic coef.
a1313           array[nx][nz] of uniformly sampled density normalized elastic coef.
a1113           array[nx][nz] of uniformly sampled density normalized elastic coef.
a3313           array[nx][nz] of uniformly sampled density normalized elastic coef.
******************************************************************************
Returned:	pointer to ray parameters sampled at discrete ray steps
******************************************************************************
Notes:
The ray ends when it runs out of time (after nt steps) or with the first 
step that is out of the (x,z) bounds of the velocity function v(x,z).
*****************************************************************************
Technical Reference:

Alkhailfah, T., 1993, Gaussian beam migration for anisotropic media:
		submitted to Geophysics.

Hangya, A., 1986, Gaussian beams in anisotropic elastic media:
      Geophys. J. R. Astr. Soc., 85, 473--563.

Cerveny, V., 1972, Seismic rays and ray intensities 
	in inhomogeneous anisotropic media: 
	Geophys. J. R. Astr. Soc., 29, 1--13.
*****************************************************************************
 Credits: CWP: Tariq Alkhalifah
*****************************************************************************/
{
	int it,kmah;
	float t,x,z,c,s,p1,q1,p2,q2,px,pz,px2,pz2,pxz,
		lx,lz,cc,ss,
		a1111,da1111dx,da1111dz,dda1111dxdx,dda1111dxdz,dda1111dzdz,
		a3333,da3333dx,da3333dz,dda3333dxdx,dda3333dxdz,dda3333dzdz,
		a1133,da1133dx,da1133dz,dda1133dxdx,dda1133dxdz,dda1133dzdz,
		a1313,da1313dx,da1313dz,dda1313dxdx,dda1313dxdz,dda1313dzdz,
		a1113,da1113dx,da1113dz,dda1113dxdx,dda1113dxdz,dda1113dzdz,
		a3313,da3313dx,da3313dz,dda3313dxdx,dda3313dxdz,dda3313dzdz,
		da1111dn,dda1111dndn,da3333dn,dda3333dndn,da1133dn,dda1133dndn,
		da1313dn,dda1313dndn,da1113dn,dda1113dndn,da3313dn,dda3313dndn,
		gamm11,gamm13,gamm33,vp2,vp,ovp,sqr;
	Ray *ray;
	RayStep *rs;
	void *a11112;
	void *a33332;
	void *a11332;
	void *a13132;
	void *a11132;
	void *a33132;

	/* allocate and initialize velocities interpolator */
	a11112 = vel2Alloc(nx,dx,fx,nz,dz,fz,a1111xz);
	a33332 = vel2Alloc(nx,dx,fx,nz,dz,fz,a3333xz);
	a11332 = vel2Alloc(nx,dx,fx,nz,dz,fz,a1133xz);
	a13132 = vel2Alloc(nx,dx,fx,nz,dz,fz,a1313xz);
	a11132 = vel2Alloc(nx,dx,fx,nz,dz,fz,a1113xz);
	a33132 = vel2Alloc(nx,dx,fx,nz,dz,fz,a3313xz);

	/* last x and z in velocity model */
	lx = fx+(nx-1)*dx;
	lz = fz+(nz-1)*dz;

	/* ensure takeoff point is within model */
	if (x0<fx || x0>lx || z0<fz || z0>lz) return NULL;

	/* allocate space for ray and raysteps */
	ray = (Ray*)alloc1(1,sizeof(Ray));
	rs = (RayStep*)alloc1(nt,sizeof(RayStep));

	/* cosine and sine of takeoff angle */
	c = cos(a0);
	s = sin(a0);
	cc = c*c;
	ss = s*s;

	/* velocities and derivatives at takeoff point */
	vel2Interp(a11112,x0,z0,&a1111,&da1111dx,&da1111dz,&dda1111dxdx,
		&dda1111dxdz,&dda1111dzdz);
	da1111dn    = da1111dx*c-da1111dz*s;
	dda1111dndn = dda1111dxdx*cc-2.0*dda1111dxdz*s*c+dda1111dzdz*ss;

	vel2Interp(a33332,x0,z0,&a3333,&da3333dx,&da3333dz,&dda3333dxdx,
		&dda3333dxdz,&dda3333dzdz);
	da3333dn    = da3333dx*c-da3333dz*s;
	dda3333dndn = dda3333dxdx*cc-2.0*dda3333dxdz*s*c+dda3333dzdz*ss;
	
	vel2Interp(a11332,x0,z0,&a1133,&da1133dx,&da1133dz,&dda1133dxdx,
		&dda1133dxdz,&dda1133dzdz);
	da1133dn    = da1133dx*c-da1133dz*s;
	dda1133dndn = dda1133dxdx*cc-2.0*dda1133dxdz*s*c+dda1133dzdz*ss;

	vel2Interp(a13132,x0,z0,&a1313,&da1313dx,&da1313dz,&dda1313dxdx,
		&dda1313dxdz,&dda1313dzdz);
	da1313dn    = da1313dx*c-da1313dz*s;
	dda1313dndn = dda1313dxdx*cc-2.0*dda1313dxdz*s*c+dda1313dzdz*ss;

	vel2Interp(a11132,x0,z0,&a1113,&da1113dx,&da1113dz,&dda1113dxdx,
		&dda1113dxdz,&dda1113dzdz);
	da1113dn    = da1113dx*c-da1113dz*s;
	dda1113dndn = dda1113dxdx*cc-2.0*dda1113dxdz*s*c+dda1113dzdz*ss;

	vel2Interp(a33132,x0,z0,&a3313,&da3313dx,&da3313dz,&dda3313dxdx,
		&dda3313dxdz,&dda3313dzdz);
	da3313dn    = da3313dx*c-da3313dz*s;
	dda3313dndn = dda3313dxdx*cc-2.0*dda3313dxdz*s*c+dda3313dzdz*ss;

	/*computing the phase velocity for a0 angle */
	gamm11 = a1111*ss+ a1313*cc +2*a1113*s*c;
	gamm33 = a3333*cc + a1313*ss+2*a3313*s*c;
	gamm13 = (a1133+a1313)*s*c+ a1113*ss+ a3313*cc;