sudmotivz.c

su 的源代码库

		qq[it].i = q0[it].i*scale;
	}
	for (it=1; it<nt; it+=2) {
		qq[it].r = -q0[it].r*scale;
		qq[it].i = -q0[it].i*scale;
	}

	/* free workspace */
	free1complex(qn);	
	free1complex(q0);	
}
	
float tmute (float *vrms, int nt, float dt, float ft, float h, float smute)
/***************************************************************************
Returns MUTE Time for given rms velocity function and half-offset.
****************************************************************************
Input:
vrms	array[nt] of rms velocity as a function of NMO time
nt	number of times
dt	time sampling interval
ft	first time
h	source-receiver half-offset
smute	maximum allowable NMO stretch t/tn [(100+stretch_percentage)/100]

Return:
	gretest time that is subject to given NMO stretch mute
****************************************************************************
Author: Craig Artley, Colorado School of Mines, 09/28/91
****************************************************************************/
{
	int it;
	float tn,*t;

	/* allocate space */
	t = ealloc1float(nt);

	for (it=0,tn=ft; it<nt; ++it,tn+=dt)
		t[it] = sqrt(tn*tn+4.0*h*h/(vrms[it]*vrms[it]));

	/* find greatest nmo time that is subject to stretch mute */
	for (it=nt-1; it>0; --it)
		if (dt/(t[it]-t[it-1])>smute) break;

	/* free space */
	free1float(t);

	return ft+it*dt;
}

void getX (float tsg, float h, float p0, float v, float *vars)
/***************************************************************************
GET trial solution X to system of equations assuming constant velocity.
****************************************************************************
Input:
tsg	recording time (from source to receiver)
h	source-receiver half-offset
p0	ray parameter of zero-offset ray
v	constant velocity

Output:
vars	array[5] of variables that satisfies the system for const. v
****************************************************************************
Author: Craig Artley, Colorado School of Mines, 09/21/91
****************************************************************************/
{
	float x,z,a2,b2,h2,x0,ps,pg,tg,t0,sine,tangent;

	/* compute constants that determine prestack migration ellipse */
	h2 = h*h;
	a2 = v*v*tsg*tsg/4.0;
	b2 = MAX(a2-h2,0.0);

	/* compute sine and tangent of propagation angle */
	sine = 0.5*p0*v;
	tangent = sine/sqrt(MAX(1.0-sine*sine,0.0));
	
	/* find reflection point (x,z) on prestack migration ellipse */
	x = tangent*sqrt(a2/(b2+tangent*tangent));
	z = sqrt(b2*MAX(1.0-x*x/a2,0.0));
	
	/* find zero-offset location */
	x0 = h2/a2*x;

	/* find times along receiver and zero-offset rays */
	tg = sqrt((h-x)*(h-x)+z*z)/v;
	t0 = sqrt((x0-x)*(x0-x)+z*z)*2.0/v;

	/* find ray parameters for source and receiver rays */
	ps = (x+h)/(v*v*(tsg-tg));
	pg = (x-h)/(v*v*tg);

	/* load trial solution into the variables array */
	vars[0] = x0;
	vars[1] = ps;
	vars[2] = pg;	
	vars[3] = tg;
	vars[4] = t0;
}

void getFJ (float **x, float **a, float **tau, float **s, int nt, int np,
	float dt, float dp, float tsg, float h, float p0,
	float *vars, float *eqns, float **jacob)
/***************************************************************************
Compute equation array and Jacobian matrix for system of equations for
a given set of parameters.
****************************************************************************
Input:
x	array[np][nt] of horizontal location of tabled rays
a	array[np][nt] of propagation angle of tabled rays
tau	array[np][nt] of travel-time depth of tabled rays
s	array[np][nt] of slownesses of tabled rays
nt	number of times in ray tables
np	number of ray parameters in ray tables
dt	time sampling interval in tables
dp	ray parameter sampling interval in tables
tsg	recording time (from source to receiver)
h	source-receiver half offset
p0	zero-offset ray parameter
vars	array[5] of variables of system of equations

Output:
eqns	array[5] of right-hand-sides of equations
jacob	pointer to array[5][5] of Jacobian partial derivatives
****************************************************************************
Author: Craig Artley and Tariq Alkhalifah, Colorado School of Mines, 94
****************************************************************************/
{
	int i,j,neqns=5;
	float x0,ps,pg,ts,tg,t0,ft=0.0,fp=0.0,ftp,dfdt,dfdp,sign,f0tp,
		df0dt,df0dp,stp,dsdt,dsdp,signp0;

	/* extract parameters from variables array */
	x0 = vars[0];
	ps = vars[1];
	pg = vars[2];
	tg = vars[3];
	t0 = vars[4];
	ts = MAX(tsg-tg,0.0);

	/* zero out F and J */
	for (i=0; i<neqns; ++i)
		eqns[i] = 0.0;
	for (j=0; j<neqns; ++j)
		for (i=0; i<neqns; ++i)
			jacob[j][i] = 0.0;

	/* insert constant terms into equations */
	eqns[0] += 2.0*h;
	eqns[1] += h-x0;
	jacob[0][1] -= 1.0;

	/* get x(ps,2*ts) and partial derivatives, load into equations */
	intder(x,nt,dt,ft,np,dp,fp,2.0*ts,ABS(ps),&ftp,&dfdt,&dfdp);
	sign = SGN(ps);
	eqns[0] -= sign*ftp;
	jacob[1][0] -= dfdp;
	jacob[3][0] += 2.0*sign*dfdt;

	/* get x(pg,2*tg) and partial derivatives, load into equations  */
	intder(x,nt,dt,ft,np,dp,fp,2.0*tg,ABS(pg),&ftp,&dfdt,&dfdp);
	sign = SGN(pg);
	eqns[0] += sign*ftp;
	eqns[1] += sign*ftp;
	jacob[2][0] += dfdp;
	jacob[2][1] += dfdp;
	jacob[3][0] += 2.0*sign*dfdt;
	jacob[3][1] += 2.0*sign*dfdt;

	/* get x(p0,t0) and partial derivatives, load into equations */
	intder(x,nt,dt,ft,np,dp,fp,t0,ABS(p0),&ftp,&dfdt,&dfdp);
	sign = SGN(p0);
	eqns[1] -= sign*ftp;
	jacob[4][1] -= sign*dfdt;

	/* get tau(ps,2*ts) and partial derivatives, load into equations */
	intder(tau,nt,dt,ft,np,dp,fp,2.0*ts,ABS(ps),&ftp,&dfdt,&dfdp);
	sign = SGN(ps);
	eqns[2] -= ftp;
	jacob[1][2] -= sign*dfdp;
	jacob[3][2] += 2.0*dfdt;

	/* get tau(pg,2*tg) and partial derivatives, load into equations */
	intder(tau,nt,dt,ft,np,dp,fp,2.0*tg,ABS(pg),&ftp,&dfdt,&dfdp);
	sign = SGN(pg);
	eqns[2] += ftp;
	eqns[3] += ftp;
	jacob[2][2] += sign*dfdp;
	jacob[2][3] += sign*dfdp;
	jacob[3][2] += 2.0*dfdt;
	jacob[3][3] += 2.0*dfdt;

	/* get tau(p0,t0) and partial derivatives, load into equations */
	intder(tau,nt,dt,ft,np,dp,fp,t0,ABS(p0),&ftp,&dfdt,&dfdp);
	sign = SGN(p0);
	eqns[3] -= ftp;
	jacob[4][3] -= dfdt;

	/* get a(ps,2*ts) and partial derivatives, load into equations */
	intder(a,nt,dt,ft,np,dp,fp,2.0*ts,ABS(ps),&ftp,&dfdt,&dfdp);
	intder(a,nt,dt,ft,np,dp,fp,t0,ABS(p0),&f0tp,&df0dt,&df0dp);
	intder(s,nt,dt,ft,np,dp,fp,2.0*ts,ABS(ps),&stp,&dsdt,&dsdp);
	sign = SGN(ps);
	signp0 = SGN(p0);
	eqns[4] -= stp*sin(sign*ftp-signp0*f0tp);
	jacob[1][4] += -sign*dsdp*sin(sign*ftp-signp0*f0tp)-
			stp*cos(sign*ftp-signp0*f0tp)*dfdp;
	jacob[3][4] += 2.0*(dsdt*sin(sign*ftp-signp0*f0tp)+
			sign*stp*cos(sign*ftp-signp0*f0tp)*dfdt);
	jacob[4][4] +=  2.*signp0*stp*cos(sign*ftp-signp0*f0tp)*df0dt;

	/* get a(pg,2*tg) and partial derivatives, load into equations */
	intder(a,nt,dt,ft,np,dp,fp,2.0*tg,ABS(pg),&ftp,&dfdt,&dfdp);
	intder(s,nt,dt,ft,np,dp,fp,2.0*tg,ABS(pg),&stp,&dsdt,&dsdp);
	sign = SGN(pg);
	eqns[4] -= stp*sin(sign*ftp-signp0*f0tp);
	jacob[2][4] += -sign*dsdp*sin(sign*ftp-signp0*f0tp)-
			stp*cos(sign*ftp-signp0*f0tp)*dfdp;
	jacob[3][4] -= 2.0*(dsdt*sin(sign*ftp-signp0*f0tp)+
			sign*stp*cos(sign*ftp-signp0*f0tp)*dfdt);
	jacob[4][4] +=  2.*signp0*stp*cos(sign*ftp-signp0*f0tp)*df0dt;
}

void intder (float **f, int nx, float dx, float fx,
	int ny, float dy, float fy, float x, float y,
	float *fxy, float *dfdx, float *dfdy)
/***************************************************************************
Finds f(x,y), df/dx, df/dy given function values on a mesh f[iy][ix],
via linear interpolation and finite differencing.
****************************************************************************
Input:
f	array[ny][nx] of sampled function values
nx	number of samples in x direction (fast dimension)
dx	x sampling interval
fx	first sample in x direction
ny	number of samples in y direction (slow dimension)
dy	y sampling interval
fy	first sample in y direction
x	x value at which to evaluate functions
y	y value at which to evaluate functions

Output:
fxy	interpolated function value f(x,y)
dfdx	partial derivative df/dx evaluated at (x,y)
dfdy	partial derivative df/dy evaluated at (x,y)
****************************************************************************
Note: Constant extrapolation of function and zero derivatives for points
	outside the table.  This is desirable, but somewhat awkward.
****************************************************************************
Author: Craig Artley, Colorado School of Mines, 08/31/91
****************************************************************************/
{
	int	ix0,ix1,iy0,iy1;
	float	x0i,y0i,fracx,fracy,lx,ly,f00,f10,f01,f11;

	/* compute x and y extent of table */
	lx = fx+(nx-1)*dx;
	ly = fy+(ny-1)*dy;

	if (x<=fx) {
		ix0 = ix1 = 0;
		fracx = 0.0;
	} else if (x>=lx) {
		ix0 = ix1 = nx-1;
		fracx = 1.0;
	} else {
		x0i = (x-fx)/dx;
		ix0 = x0i;
		ix1 = ix0+1;
		fracx = x0i-ix0;
	}

	if (y<=fy) {
		iy0 = iy1 = 0;
		fracy = 0;
	} else if (y>=ly) {
		iy0 = iy1 = ny-1;
		fracy = 1.0;
	} else {
		y0i = (y-fy)/dy;
		iy0 = y0i;
		iy1 = iy0+1;
		fracy = y0i-iy0;
	}

	/* get tabled function values nearest (x,y) */
	f00 = f[iy0][ix0];	f01 = f[iy1][ix0];
	f10 = f[iy0][ix1];	f11 = f[iy1][ix1];

	/* linear interpolation to find f(x,y) */
	*fxy = (1-fracx)*(1-fracy)*f00 + (1-fracx)*fracy*f01
		+ fracx*(1-fracy)*f10 + fracx*fracy*f11;

	/* interpolate and finite difference to find partial derivatives */
	*dfdx = ((1-fracy)*(f10-f00) + fracy*(f11-f01))/dx;
	*dfdy = ((1-fracx)*(f01-f00) + fracx*(f11-f10))/dy;
}

void rayvt (float p, int nt, float dt,
	float a1111[], float a3333[], float a1313[], float a1133[],
	float tau[], float x[], float a[], float s[]) 
/*****************************************************************************
trace ray for the given elastic coeffecient as a function of verical time, 
tau=ONE-way vertical time
******************************************************************************
Input:
p		ray parameter
nt		number of times (including t=0.0)
dt		time sampling interval
a1111		array[nt] of elastic coeffecient a1111 (a1111 = a1111(tau))
a3333		array[nt] of elastic coeffecient a3333 (a3333 = a3333(tau))
a1313		array[nt] of elastic coeffecient a1313 (a1313 = a1313(tau))
a1133		array[nt] of elastic coeffecient a1133 (a1133 = a1133(tau))

Output:
tau		array[nt] of vertical times tau(t)
x		array[nt] of horizontal distances x(t)
a		array[nt] of propagation angles a(t)
s  		array[nt] of slownesses s(t)
******************************************************************************
Notes:
Ray tracing begins at tau=x=t=0.0.  If the ray turns and reaches tau=0.0
at some time index it, then tau[it+1:nt-1] = tau[it], x[it+1:nt-1] = x[it],
and a[it+1:nt-1] = a[it].  In other words, constant extrapolation is used
to fill the remaining values in the output arrays.
******************************************************************************
Author:  Tariq Alkhalifah, Colorado School of Mines, 94
******************************************************************************/
{
	int it,jt,itau,turn;
	float at,taut,xt,taui,frac;
	float a1111t,a3333t,a1313t,a1133t;
	float f1,f2,f3,f4,f5,f6,eps,pz,pxz,px2,pz2,alpha,beta,gamma,det;
	float rad,signbeta,q,vp,gamma11,gamma33,gamma13,den,g11,g33,g13;
	float fac,dx,dz,dtau,sp;

	/* initialize ray tracing parameters */
	a1111t = a1111[0];
	a3333t = a3333[0];
	a1313t = a1313[0];
	a1133t = a1133[0];

	f1 = a1111t+a1313t;
	f2 = a3333t+a1313t;
	f3 = a1111t-a1313t;
	f4 = a3333t-a1313t;
	f5 = 2.*(a1133t+a1313t)*(a1133t+a1313t);
	f6 = 2;


	eps   = .00001;
	px2   = p*p;
	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) return;
	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) return;
	sp   = sqrt(px2+pz2);
	vp   = 1/sp;
	pz   = sqrt(pz2);
	pxz  = p*pz;

	gamma11 = a1111t*px2+ a1313t*pz2;
	gamma33 = a3333t*pz2 + a1313t*px2;
	gamma13 = (a1133t+a1313t)*pxz;
	den     = 1/(gamma11+gamma33-2);
	g11     = (gamma33-1)*den;
	g33     = (gamma11-1)*den;
	g13     = -gamma13*den;
	fac=1/sqrt(a3333t);

	dx =  a1111t*p*g11+(a1133t+a1313t)*pz*g13+a1313t*g33*p;
	dz = a3333t*pz*g33+(a1133t+a1313t)*p*g13+a1313t*g11*pz;
	dtau =  MIN(1.0,dz*fac);

	at = acos(MIN(1.0,pz*vp));

	taut = xt = 0.0;
	tau[0] = taut;
	x[0] = xt;
	a[0] = at;
	s[0] = sp;
	turn=1;

	/* loop over times t */
	for (it=1; it<nt; ++it) {
		
		/* update vertical time */
		tau[it] = tau[it-1]+dt*dtau;

		/* if ray has returned to surface, break */
		if (tau[it]<=0.0) break;

		/* update horizontal distance */
		x[it] = x[it-1]+dx*dt;

		/* update velocity and derivative */
		if (it<nt-1) {
			taui = tau[it]/dt;
			itau = taui;
			frac = taui-itau;
			a1111t = (1.0-frac)*a1111[itau] + frac*a1111[itau+1];
			a3333t = (1.0-frac)*a3333[itau] + frac*a3333[itau+1];
			a1313t = (1.0-frac)*a1313[itau] + frac*a1313[itau+1];
			a1133t = (1.0-frac)*a1133[itau] + frac*a1133[itau+1];
		}

		f1 = a1111t+a1313t;
		f2 = a3333t+a1313t;
		f3 = a1111t-a1313t;
		f4 = a3333t-a1313t;
		f5 = 2.*(a1133t+a1313t)*(a1133t+a1313t);

		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){
		fprintf(stderr,"det=%f\n",det);
		return;
		}
		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)
			turn *=-1;
		pz2 = ABS(pz2);
		pz = turn*sqrt(pz2);

		sp   = sqrt(px2+pz2);
		vp   = 1/sp;
		
		pxz  = p*pz;
		gamma11 = a1111t*px2+ a1313t*pz2;
		gamma33 = a3333t*pz2 + a1313t*px2;
		gamma13 = (a1133t+a1313t)*pxz;
		den     = 1/(gamma11+gamma33-2);
		g11     = (gamma33-1)*den;
		g33     = (gamma11-1)*den;
		g13     = -gamma13*den;
		fac=1/sqrt(a3333t);

		dx =  a1111t*p*g11+(a1133t+a1313t)*pz*g13+a1313t*g33*p;
		dz = a3333t*pz*g33+(a1133t+a1313t)*p*g13+a1313t*g11*pz;
		dtau =  dz*fac;

		at = acos(MIN(1.0,pz*vp));
		a[it] = at;
		s[it] = sp;
	}

	/* if ray returned to surface, extrapolate with surface values */
	for (jt=it; jt<nt; ++jt) {
		tau[jt] = tau[jt-1];
		x[jt] = x[jt-1];
		a[jt] = a[jt-1];
		s[jt] = s[jt-1];
	}
}