sumiggbzo.c

su 的源代码库

				vu[mx+3][mz+1]*px[3]*pz[1]+
				vu[mx+3][mz+2]*px[3]*pz[2]+
				vu[mx+3][mz+3]*px[3]*pz[3]);
		}
		
	/* else handle end effects with constant extrapolation */
	} else {
		for (i=0; i<6; ++i) {
			px = &(tbl[kx[i]][0][0])+lx;
			pz = &(tbl[kz[i]][0][0])+lz;
			for (jx=0,jmx=mx,sum=0.0; jx<4; ++jx,++jmx) {
				mmx = jmx;
				if (mmx<0) mmx = 0;
				else if (mmx>=nx) mmx = nx-1;
				for (jz=0,jmz=mz; jz<4; ++jz,++jmz) {
					mmz = jmz;
					if (mmz<0) mmz = 0;
					else if (mmz>=nz) mmz = nz-1;
					sum += vu[mmx][mmz]*px[jx]*pz[jz];
				}
			}
			vd[i] = sx[kx[i]]*sz[kz[i]]*sum;
		}
	}

	/* set output variables */
	*v = vd[0];
	*vx = vd[1];
	*vz = vd[2];
	*vxx = vd[3];
	*vxz = vd[4];
	*vzz = vd[5];
}

/* hermite polynomials */
static float h00 (float x) {return 2.0*x*x*x-3.0*x*x+1.0;}
static float h01 (float x) {return 6.0*x*x-6.0*x;}
static float h02 (float x) {return 12.0*x-6.0;}
static float h10 (float x) {return -2.0*x*x*x+3.0*x*x;}
static float h11 (float x) {return -6.0*x*x+6.0*x;}
static float h12 (float x) {return -12.0*x+6.0;}
static float k00 (float x) {return x*x*x-2.0*x*x+x;}
static float k01 (float x) {return 3.0*x*x-4.0*x+1.0;}
static float k02 (float x) {return 6.0*x-4.0;}
static float k10 (float x) {return x*x*x-x*x;}
static float k11 (float x) {return 3.0*x*x-2.0*x;}
static float k12 (float x) {return 6.0*x-2.0;}

/* function to build interpolation tables */
static void buildTables(void)
{
	int i;
	float x;
	
	/* tabulate interpolator for 0th derivative */
	for (i=0; i<NTABLE; ++i) {
		x = (float)i/(NTABLE-1.0);
		tbl[0][i][0] = -0.5*k00(x);
		tbl[0][i][1] = h00(x)-0.5*k10(x);
		tbl[0][i][2] = h10(x)+0.5*k00(x);
		tbl[0][i][3] = 0.5*k10(x);
		tbl[1][i][0] = -0.5*k01(x);
		tbl[1][i][1] = h01(x)-0.5*k11(x);
		tbl[1][i][2] = h11(x)+0.5*k01(x);
		tbl[1][i][3] = 0.5*k11(x);
		tbl[2][i][0] = -0.5*k02(x);
		tbl[2][i][1] = h02(x)-0.5*k12(x);
		tbl[2][i][2] = h12(x)+0.5*k02(x);
		tbl[2][i][3] = 0.5*k12(x);
	}
	
	/* remember that tables have been built */
	tabled = 1;
}

#ifdef TEST2
#include "cwp.h"
#define NZ 2
#define NX 2
#define NXOUT 21
#define NZOUT 21
main()
{
	int nx=NX,nz=NZ,nxout=NXOUT,nzout=NZOUT,i,j;
	float dx=2.0,fx=-1.0,dxout=0.2,fxout=-2.0;
	float dz=4.0,fz=-2.0,dzout=0.4,fzout=-4.0;
	float x,z,v,vx,vz,vxx,vxz,vzz,**vu,**vo;
	void *vel2;

	vu = alloc2float(nz,nx);
	vo = alloc2float(nzout,nxout);

	vu[0][0] = 1.0;
	vu[1][0] = 2.0;
	vu[0][1] = 1.0;
	vu[1][1] = 2.0;

	vel2 = vel2Alloc(nx,dx,fx,nz,dz,fz,vu);

	for (i=0; i<nxout; ++i) {
		x = fxout+i*dxout;
		for (j=0; j<nzout; ++j) {
			z = fzout+j*dzout;
			vel2Interp(vel2,x,z,&v,&vx,&vz,&vxx,&vxz,&vzz);
			vo[i][j] = vz;
		}
	}
	vel2Free(vel2);
	fwrite(vo[0],sizeof(float),nxout*nzout,stdout);
}
#endif

/* Beam subroutines */
/* 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
*****************************************************************************/
{
	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=NULL,*gpad=NULL;
	complex **cf=NULL,**cg=NULL,*cfx=NULL,*cgpx=NULL;

	/* 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
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
*****************************************************************************/
{
	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)