sumiggbzoan.c

su 的源代码库

		a0,x0,z0,bwhc,**b;
	float f1,f2,f3,f4,f5,f6,alpha,beta,gamma,det,
		rad,signbeta,eps,pz2,q,vp,px2,pz;
	float aa,bb,ccc,dd,ee,fxx,dfx,f7,f8,f9,f10,f11,dpz;
	float gamm11,gamm33,gamm13,vpmin,vpmax,vp2,s,c,ss,cc;
	Ray *ray;

	/* first t, x, and z assumed to be zero */
	ft = fx = fz = 0.0;

	/* convert minimum and maximum angles to radians */
	amin *= PI/180.0;
	amax *= PI/180.0;
	if (amin>amax) {
		float atemp=amin;
		amin = amax;
		amax = atemp;
	}
	
	/* window half-width */
	xwh = 3.0*bwh;

	/* beam center sampling */
	dxb = NINT((bwh*sqrt(2.0*fmin/fmax))/dx)*dx;
	nxb = 1+(nx-1)*dx/dxb;
	fxb = fx+0.5*((nx-1)*dx-(nxb-1)*dxb);

	fprintf(stderr,"nxb=%d dxb=%g fxb=%g bwh=%g\n",nxb,dxb,fxb,bwh);

	/* minimum velocity at surface z=0 */
	ixm = 0;
	for (ix=1,vmin=a3333[0][0]; ix<nx; ++ix)
		if (a3333[ix][0]<vmin) { 
			ixm=ix;
			vmin = a3333[ix][0];
		}

	/* phase velocity at min. and max. angles */
	c = cos(amin);
	s = sin(amin);
	cc = c*c;
	ss = s*s;
	gamm11 = a1111[ixm][0]*ss+ a1313[ixm][0]*cc +2*a1113[ixm][0]*s*c;
	gamm33 = a3333[ixm][0]*cc + a1313[ixm][0]*ss+2*a3313[ixm][0]*s*c;
	gamm13 = (a1133[ixm][0]+a1313[ixm][0])*s*c+ a1113[ixm][0]*ss+ 
		a3313[ixm][0]*cc;
	vp2    = gamm11+gamm33+sqrt((gamm11+gamm33)*(gamm11+gamm33)-
			4*(gamm11*gamm33-gamm13*gamm13));
	vpmin     = sqrt(vp2*.5);
	c = cos(amax);
	s = sin(amax);
	cc = c*c;
	ss = s*s;
	gamm11 = a1111[ixm][0]*ss+ a1313[ixm][0]*cc +2*a1113[ixm][0]*s*c;
	gamm33 = a3333[ixm][0]*cc + a1313[ixm][0]*ss+2*a3313[ixm][0]*s*c;
	gamm13 = (a1133[ixm][0]+a1313[ixm][0])*s*c+ a1113[ixm][0]*ss+ 
		a3313[ixm][0]*cc;
	vp2    = gamm11+gamm33+sqrt((gamm11+gamm33)*(gamm11+gamm33)-
			4*(gamm11*gamm33-gamm13*gamm13));
	vpmax     = sqrt(vp2*.5);

	/* beam sampling */
	pxmin = sin(amin)/vpmin;
	pxmax = sin(amax)/vpmax;
	dpx = sqrt(vmin)*1.0/(2.0*xwh*sqrt(fmin*fmax)*.5*(vpmin+vpmax));
	npx = 1+(pxmax-pxmin)/dpx;
	fpx = pxmin+0.5*(pxmax-pxmin-(npx-1)*dpx);
	taupad = MAX(ABS(pxmin),ABS(pxmax))*xwh;
	taupad = NINT(taupad/dt)*dt;
	ftau = ft-taupad;
	dtau = dt;
	ntau = nt+2.0*taupad/dtau;

	fprintf(stderr,"npx=%d dpx=%g fpx=%g\n",npx,dpx,fpx);

	/* zero migrated image */
	for (ix=0; ix<nx; ++ix)
		for (iz=0; iz<nz; ++iz)
			g[ix][iz] = 0.0;

	/* loop over beam centers */
	for (ixb=0,xb=fxb; ixb<nxb; ++ixb,xb+=dxb) {

		/* horizontal window */
		ix = NINT((xb-fx)/dx);
		ixlo = MAX(ix-NINT(xwh/dx),0);
		ixhi = MIN(ix+NINT(xwh/dx),nx-1);
		nxw = 1+ixhi-ixlo;
		fxw = fx+(ixlo-ix)*dx;

		fprintf(stderr,"ixb/nxb = %d/%d  ix = %d\n",ixb,nxb,ix);

		/* allocate space for beams */
		b = alloc2float(ntau,npx);

		/* form beams at surface */
		formBeams(bwh,dxb,fmin,
			nt,dt,ft,nxw,dx,fxw,&f[ixlo],
			ntau,dtau,ftau,npx,dpx,fpx,b);
		
		/* loop over beams */
		for (ipx=0,px=fpx; ipx<npx; ++ipx,px+=dpx) {

			/*calculating emergence angle from ray parameter*/
			f1 = a1111[ix][0]+a1313[ix][0];
			f2 = a3333[ix][0]+a1313[ix][0];
			f3 = a1111[ix][0]-a1313[ix][0];
			f4 = a3333[ix][0]-a1313[ix][0];
			f5 = 2.*(a1133[ix][0]+a1313[ix][0])*
				(a1133[ix][0]+a1313[ix][0]);
			f6 = 2;

			eps   = .0001;
			px2   = px*px;
	    		alpha = f2*f2-f4*f4;
	    		beta  = 2*((f1*f2+f3*f4-f5)*px2-f2*f6);
	    		gamma = f6*f6-(2.*f1*f6-(f1*f1-f3*f3)*px2)*px2;
	    		det   = beta*beta-4.*alpha*gamma;

			if (det<0) continue;
			rad = sqrt(det);
			if(ABS(beta)>eps)   signbeta = ABS(beta)/beta;
	    		else                signbeta = 1.;
			q    = -.5*(beta+signbeta*rad);
			pz2  = gamma/q;
			if(pz2<0) continue;
			pz   = sqrt(pz2);

			f7   = 2*(a1113[ix][0]+a3313[ix][0]);
			f8   = 2*(a1113[ix][0]+a3313[ix][0]);		
			f9   = a3313[ix][0];
			f10  = a1113[ix][0];
			f11  = 2*(a1133[ix][0]+a1313[ix][0]);
			aa   = f2*f2-f4*f4-f9*f9;
			bb   = (2*f2*f7+2*f4*f8-2*f11*f9)*px;
			ccc  = f7*f7*px2-2*(f6-f1*px2)*f2-f8*f8*px2+
				2*f3*f4*px2-f11*f11*px2-2*f10*f9*px2;
			dd   = -2*(f6-f1*px2)*f7*px-2*f3*f8*px2*px-
				2*f11*f10*px*px2;
			ee   = (f6-f1*px2)*(f6-f1*px2)-f3*f3*px2*px2-
				f10*f10*px2*px2;
			for (i=0; i<20; i++) {
				fxx = aa*pz2*pz2+bb*pz2*pz+ccc*pz2+
					dd*pz+ee;
				dfx = 4*aa*pz*pz2+3*bb*pz2+2*ccc*pz+dd;
				if(dfx==0) break;
				dpz = -fxx/dfx;
				if(ABS(dpz)<.000001) break;
				pz += dpz;
				pz2 = pz*pz;
			}
			if (pz<0) continue;
			vp   = 1/sqrt(px2+pz2);

			/* sine of emergence angle; skip if out of bounds */
			if (px*vp>sin(amax)+0.01) continue;
			if (px*vp<sin(amin)-0.01) continue;

			/* emergence angle and location */
			a0 = -asin(px*vp);
			x0 = fx+ix*dx;
			z0 = fz;

			/* beam half-width adjusted for cosine of angle */
			bwhc = bwh*cos(a0);

			/* trace ray */
			ray = makeRay(x0,z0,a0,nt,dt,ft,
			a1111,a3333,a1133,a1313,a1113,a3313,
			nx,dx,fx,nz,dz,fz);

			/* accumulate contribution of beam in migrated image */
			accBeam(ray,fmin,bwhc,
				ntau,dtau,ftau,b[ipx],
				nx,dx,fx,nz,dz,fz,g);
			
			/* free ray */
			freeRay(ray);
		}

		/* free space for beams */
		free2float(b);
	}
}

/************************************************************************
* Functions for Gaussian Beam computation		                *
*************************************************************************
*  Functions:							        *
*	     formBeams: Form beams (filtered slant stacks) for later    *
*		        superposition of Gaussian beams.	        *
*	     accBeam:   Accumulate contribution of one Gaussian beam.   *
*************************************************************************/


/* size of cells in which to linearly interpolate complex time and amplitude */
#define CELLSIZE 8

/* factor by which to oversample time for linear interpolation of traces */
#define NOVERSAMPLE 4

/* number of exponential decay filters */
#define NFILTER 6

/* exp(EXPMIN) is assumed to be negligible */
#define EXPMIN (-5.0)

/* filtered complex beam data as a function of real and imaginary time */
typedef struct BeamDataStruct {
	int ntr;		/* number of real time samples */
	float dtr;		/* real time sampling interval */
	float ftr;		/* first real time sample */
	int nti;		/* number of imaginary time samples */
	float dti;		/* imaginary time sampling interval */
	float fti;		/* first imaginary time sample */
	complex **cf;		/* array[nti][ntr] of complex data */
} BeamData;

/* one cell in which to linearly interpolate complex time and amplitude */
typedef struct CellStruct {
	int live;	/* random number used to denote a live cell */
	int dead;	/* random number used to denote a dead cell */
	float tr;	/* real part of traveltime */
	float ti;	/* imaginary part of traveltime */
	float ar;	/* real part of amplitude */
	float ai;	/* imaginary part of amplitude */
} Cell;

/* structure containing information used to set and fill cells */
typedef struct CellsStruct {
	int nt;		/* number of time samples */
	float dt;	/* time sampling interval */
	float ft;	/* first time sample */
	int lx;		/* number of x samples per cell */
	int mx;		/* number of x cells */
	int nx;		/* number of x samples */
	float dx;	/* x sampling interval */
	float fx;	/* first x sample */
	int lz;		/* number of z samples per cell */
	int mz;		/* number of z cells */
	int nz;		/* number of z samples */
	float dz;	/* z sampling interval */
	float fz;	/* first z sample */
	int live;	/* random number used to denote a live cell */
	int dead;	/* random number used to denote a dead cell */
	float wmin;	/* minimum (reference) frequency */
	float lmin;	/* minimum beamwidth for frequency wmin */
	Cell **cell;	/* cell array[mx][mz] */
	Ray *ray;	/* ray */
	BeamData *bd;	/* complex beam data as a function of complex time */
	float **g;	/* array[nx][nz] containing g(x,z) */
} Cells;

/* functions defined and used internally */
static void xtop (float w,
	int nx, float dx, float fx, complex *g,
	int np, float dp, float fp, complex *h);
static BeamData* beamData (float wmin, int nt, float dt, float ft, float *f);
static void setCell (Cells *cells, int jx, int jz);
static void accCell (Cells *cells, int jx, int jz);
static int cellTimeAmp (Cells *cells, int jx, int jz);
static void cellBeam (Cells *cells, int jx, int jz);


/* functions for external use */

void formBeams (float bwh, float dxb, float fmin,
	int nt, float dt, float ft, 
	int nx, float dx, float fx, float **f,
	int ntau, float dtau, float ftau, 
	int npx, float dpx, float fpx, float **g)
/*****************************************************************************
Form beams (filtered slant stacks) for later superposition of Gaussian beams.
******************************************************************************
Input:
bwh		horizontal beam half-width
dxb		horizontal distance between beam centers
fmin		minimum frequency (cycles per unit time)
nt		number of input time samples
dt		input time sampling interval
ft		first input time sample
nx		number of horizontal samples
dx		horizontal sampling interval
fx		first horizontal sample
f		array[nx][nt] of data to be slant stacked into beams
ntau		number of output time samples
dtau		output time sampling interval (currently must equal dt)
ftau		first output time sample
npx		number of horizontal slownesses
dpx		horizontal slowness sampling interval
fpx		first horizontal slowness

*****************************************************************************
Output:
g		array[npx][ntau] containing beams
*****************************************************************************
Technical Reference:
  Hale, D., 1992, Migration by the Kirchhoff, 
 	slant stack, and Gaussian beam methods:
       CWP,1992 Report 121, Colorado School of Mines.
*****************************************************************************
Credits: CWP: Dave Hale
*****************************************************************************/
{
	int ntpad,ntfft,nw,ix,iw,ipx,it,itau;
	float wmin,pxmax,xmax,x,dw,fw,w,fftscl,
		amp,phase,scale,a,b,as,bs,es,cfr,cfi,
		*fpad,*gpad;
	complex **cf,**cg,*cfx,*cgpx;

	/* minimum frequency in radians */
	wmin = 2.0*PI*fmin;

	/* pad time axis to avoid wraparound */
	pxmax = (dpx<0.0)?fpx:fpx+(npx-1)*dpx;
	xmax = (dx<0.0)?fx:fx+(nx-1)*dx;
	ntpad = ABS(pxmax*xmax)/dt;

	/* fft sampling */
	ntfft = npfar(MAX(nt+ntpad,ntau));
	nw = ntfft/2+1;
	dw = 2.0*PI/(ntfft*dt);
	fw = 0.0;
	fftscl = 1.0/ntfft;

	/* allocate workspace */
	fpad = alloc1float(ntfft);
	gpad = alloc1float(ntfft);
	cf = alloc2complex(nw,nx);
	cg = alloc2complex(nw,npx);
	cfx = alloc1complex(nx);
	cgpx = alloc1complex(npx);

	/* loop over x */
	for (ix=0; ix<nx; ++ix) {

		/* pad time with zeros */
		for (it=0; it<nt; ++it)
			fpad[it] = f[ix][it];
		for (it=nt; it<ntfft; ++it)
			fpad[it] = 0.0;
		
		/* Fourier transform time to frequency */
		pfarc(1,ntfft,fpad,cf[ix]);
	}

	/* loop over w */
	for (iw=0,w=fw; iw<nw; ++iw,w+=dw) {
		
		/* frequency-dependent amplitude scale factors */
		scale = -0.5*w/(wmin*bwh*bwh);
		amp = fftscl*dpx*w/(2.0*PI)*sqrt(w/(PI*wmin))*dxb/bwh;

		/* phase shift to account for ft not equal to ftau */
		phase = w*(ft-ftau);

		/* apply complex filter */
		a = amp*cos(phase);
		b = amp*sin(phase);
		for (ix=0,x=fx; ix<nx; ++ix,x+=dx) {
			es = exp(scale*x*x);
			as = a*es;
			bs = b*es;
			cfr = cf[ix][iw].r;
			cfi = cf[ix][iw].i;
			cfx[ix].r = as*cfr-bs*cfi;
			cfx[ix].i = bs*cfr+as*cfi;
		}

		/* transform x to p */
		xtop(w,nx,dx,fx,cfx,npx,dpx,fpx,cgpx);
		for (ipx=0; ipx<npx; ++ipx) {
			cg[ipx][iw].r = cgpx[ipx].r;
			cg[ipx][iw].i = cgpx[ipx].i;
		}
	}

	/* loop over px */
	for (ipx=0; ipx<npx; ++ipx) {
		
		/* Fourier transform frequency to time */
		pfacr(-1,ntfft,cg[ipx],gpad);

		/* copy to output array */
		for (itau=0; itau<ntau; ++itau)
			g[ipx][itau] = gpad[itau];
	}

	free1float(fpad);
	free1float(gpad);
	free2complex(cf);
	free2complex(cg);
	free1complex(cfx);
	free1complex(cgpx);
}

void accBeam (Ray *ray, float fmin, float lmin,
	int nt, float dt, float ft, float *f,
	int nx, float dx, float fx, int nz, float dz, float fz, float **g)
/*****************************************************************************
Accumulate contribution of one Gaussian beam.
******************************************************************************
Input:
ray		ray parameters sampled at discrete ray steps
fmin		minimum frequency (cycles per unit time)
lmin		initial beam width for frequency wmin
nt		number of time samples