velpertan.c

su 的源代码库

}

/**************************************************************************/
/****** the numerical algorithm to calculate the right takeoff angle ******/
/**************************************************************************/
float zbrentou(float da,float da_2,float tol,float off,float x, float z,
		float a_normal,int nt,float dt,float **a1111,
		float **a3333,float **a1133, float **a1313,float **a1113,
		float **a3313,int nx,float dx,float fx,int nz,float dz,
		float fz,float amax,float amin)
{
 int NRS,iter,ITMAX;
 float a=da,b=da_2,c=da_2,d=0,e=0,min1,min2,X1,X2;
 float fa,fb,fc,p,q,r,s,tol1,xm,EPS;
 Ray *ray;

ITMAX=1000; EPS=3e-8;

ray = makeRay(x,z,a_normal+a,nt,dt,a1111,a3333,a1133,a1313,
                        a1113,a3313,nx,dx,fx,nz,dz,fz,amax,amin);
NRS=ray->nrs;
	X1=-((ray->rs[NRS-1].x-ray->rs[NRS-2].x)/(ray->rs[NRS-1].z-
		ray->rs[NRS-2].z))*ray->rs[NRS-2].z+ray->rs[NRS-2].x;
freeRay(ray);
ray = makeRay(x,z,a_normal-a,nt,dt,a1111,a3333,a1133,a1313,
                        a1113,a3313,nx,dx,fx,nz,dz,fz,amax,amin);
NRS=ray->nrs;
	X2=-((ray->rs[NRS-1].x-ray->rs[NRS-2].x)/(ray->rs[NRS-1].z-
		ray->rs[NRS-2].z))*ray->rs[NRS-2].z+ray->rs[NRS-2].x;
freeRay(ray);
fa=X2-X1-off;
/*printf("%f \t",X2-X1);*/
ray = makeRay(x,z,a_normal+b,nt,dt,a1111,a3333,a1133,a1313,
                        a1113,a3313,nx,dx,fx,nz,dz,fz,amax,amin);
NRS=ray->nrs;
	X1=-((ray->rs[NRS-1].x-ray->rs[NRS-2].x)/(ray->rs[NRS-1].z-
		ray->rs[NRS-2].z))*ray->rs[NRS-2].z+ray->rs[NRS-2].x;
freeRay(ray);
ray = makeRay(x,z,a_normal-b,nt,dt,a1111,a3333,a1133,a1313,
                        a1113,a3313,nx,dx,fx,nz,dz,fz,amax,amin);
NRS=ray->nrs;
	X2=-((ray->rs[NRS-1].x-ray->rs[NRS-2].x)/(ray->rs[NRS-1].z-
		ray->rs[NRS-2].z))*ray->rs[NRS-2].z+ray->rs[NRS-2].x;
freeRay(ray);
fb=X2-X1-off;
/*printf("%f \n",X2-X1);*/
if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) 
    printf("Root must be bracketed in zbrent");

fc=fb;
for (iter = 1;iter<=ITMAX;iter++){
   if ((fb > 0.0 && fc > 0.0)||(fb < 0.0 && fc < 0.0)) {
       c=a;
       fc=fa;
       e=d=b-a;
}
if (fabs(fc) < fabs(fb)) {
       a=b;
       b=c;
       c=a;
       fa=fb;
       fb=fc;
       fc=fa;
}

tol1=2.0*EPS*fabs(b)+0.5*tol;
xm=0.5*(c-b);
if ((fabs(xm)<=tol1)||(fb==0) ){
return b;
}
if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
       s=fb/fa;
       if (a==c) {
          p=2.0*xm*s;
          q=1.0-s;
       } else {
          q=fa/fc;
          r=fb/fc;
          p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1));
          q=(q-1.0)*(r-1.0)*(s-1.0);
}
if (p>0.0) q = -q;
p=fabs(p);
min1=3.0*xm*q-fabs(tol1*q);
min2=fabs(e*q);
if(2.0*p < (min1-min2 ? min1 : min2)){
     e=d;
     d=p/q;
} else {
     d=xm;
     e=d;
}
}else{
   d=xm;
   e=d;
}
a=b;
fa=fb;
if (fabs(d) > tol1)
    b+=d;
else
    b += fabs(tol1)*xm/fabs(xm);
ray = makeRay(x,z,a_normal+b,nt,dt,a1111,a3333,a1133,a1313,
                        a1113,a3313,nx,dx,fx,nz,dz,fz,amax,amin);
NRS=ray->nrs;
	X1=-((ray->rs[NRS-1].x-ray->rs[NRS-2].x)/(ray->rs[NRS-1].z-
		ray->rs[NRS-2].z))*ray->rs[NRS-2].z+ray->rs[NRS-2].x;
freeRay(ray);
ray = makeRay(x,z,a_normal-b,nt,dt,a1111,a3333,a1133,a1313,
                        a1113,a3313,nx,dx,fx,nz,dz,fz,amax,amin);
NRS=ray->nrs;
	X2=-((ray->rs[NRS-1].x-ray->rs[NRS-2].x)/(ray->rs[NRS-1].z-
		ray->rs[NRS-2].z))*ray->rs[NRS-2].z+ray->rs[NRS-2].x;
freeRay(ray);
fb=X2-X1-off;
}
printf("Maximum number of iterations exceeded in zbrent\n");
return 0.0;
}

/****************************************************************************/
/***************************  Tariq's ray tracer ****************************/
/****************************************************************************/
Ray *makeRay (float x0, float z0, float a0, int nt, float dt,
	float **a1111xz, float **a3333xz, float **a1133xz, float **a1313xz, 
	float **a1113xz, float **a3313xz, 
	int nx, float dx, float fx, int nz, float dz, float fz, float amax, 
	float amin)
/*****************************************************************************
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
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
amax            maximum emergence angle
amin            minimum emergence angle
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:

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;

	/*Convert degrees to radians*/
	a0=a0*3.1415/180; amax=amax*3.1415/180; amin=amin*3.1415/180;
	/* 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 ){
	warn("The CRP lies outside the defined grid");
	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;
	sqr    = sqrt((gamm11+gamm33)*(gamm11+gamm33)-
			4*(gamm11*gamm33-gamm13*gamm13));
	vp2    = gamm11+gamm33+sqr;
	vp     = sqrt(vp2*.5);
	ovp    = 1/vp;
	px     = s*ovp;
	pz     = c*ovp;

	/* first ray step */
	rs[0].t = t = 0;
	rs[0].x = x = x0;
	rs[0].z = z = z0;
	rs[0].q1 = q1 = 1.0;
	rs[0].p1 = p1 = 0.0;
	rs[0].q2 = q2 = 0.0;
	rs[0].p2 = p2 = 1.0;
	rs[0].kmah = kmah = 0;
	rs[0].c = c;
	rs[0].s = s;
	rs[0].v = vp;
	rs[0].dvdx = .5*da3333dx*vp/a3333;
	rs[0].dvdz = .5*da3333dz*vp/a3333;

	/* loop over time steps */
	for (it=1; it<nt; ++it) {

		/* variables used for Runge-Kutta integration */
		float h=dt,hhalf=dt/2.0,hsixth=dt/6.0,
			q2old,xt,zt,p1t,q1t,p2t,q2t,
			dx,dz,dp1,dq1,dp2,dq2,
			dxt,dzt,dp1t,dq1t,dp2t,dq2t,
			dxm,dzm,dp1m,dq1m,dp2m,dq2m,
			gamma11,gamma33,gamma13,g11,g13,g33,den,
			sxfactor,szfactor,snfact,dpx,dpz,pxt,pzt,dpxt,
			dpzt,dpxm,dpzm,dxdn,dzdn,snfactor,
			dxx,dzz,dcdp1,dcdp3,dcdp13,ddcdnn,ddcdqq,
			ddcdpn,dgamma11dpx,dgamma11dpz,dgamma33dpx,
			dgamma33dpz,dgamma13dpx,dgamma13dpz,dg11dpx,
			dg11dpz,dg33dpx,dg33dpz,dg13dpx,dg13dpz,ddxdpx,
			ddzdpz,ddxdpz,dgamma11dn,dgamma33dn,dgamma13dn,
			dg11dn,dg33dn,dg13dn,dsnfactordn,ddxdn,ddzdn;

		/* if ray is out of bounds, break */
		if (x<fx || x>lx || z<fz || z>lz || c>(cos(amin)+0.01) || c<(cos(amax))-0.01) break;

		/* remember old q2 */
		q2old = q2;
		
	        /* step 1 of 4th-order Runge-Kutta */
		px2   = px*px;
		pz2   = pz*pz;
		pxz   = px*pz;

		/*anisotropy parameters*/
		gamma11 = a1111*px2+ a1313*pz2 +2*a1113*pxz;
		gamma33 = a3333*pz2 + a1313*px2+2*a3313*pxz;
		gamma13 = (a1133+a1313)*pxz+ a1113*px2+ a3313*pz2;
		den     = 1/(gamma11+gamma33-2);
		g11     = (gamma33-1)*den;
		g33     = (gamma11-1)*den;
		g13     = -gamma13*den;
		sxfactor = da1111dx*px2*g11+da3333dx*pz2*g33+
			2*(da1133dx+da1313dx)*pxz*g13+da1313dx*(px2*g33+pz2*g11)+
			2*da3313dx*(pz2*g13+pxz*g33)+2*da1113dx*(pxz*g11+px2*g13);
		szfactor = da1111dz*px2*g11+da3333dz*pz2*g33+
			2*(da1133dz+da1313dz)*pxz*g13+da1313dz*(px2*g33+pz2*g11)+
			2*da3313dz*(pz2*g13+pxz*g33)+2*da1113dz*(pxz*g11+px2*g13);
		snfact = sxfactor*c-szfactor*s;
		
		/*computing ray velocities and derivatives*/
		dx =  (a1111*px*g11+(a1133+a1313)*pz*g13+a3313*pz*g33
			+a1113*(pz*g11+2*px*g13)+a1313*g33*px);
		dz =  (a3333*pz*g33+(a1133+a1313)*px*g13+a1113*px*g11
			+a3313*(px*g33+2*pz*g13)+a1313*g11*pz);

		dgamma11dpx = 2*a1111*px+2*a1113*pz;
		dgamma11dpz = 2*a1313*pz+2*a1113*px;
		dgamma33dpx = 2*a1313*px+2*a3313*pz;
		dgamma33dpz = 2*a3333*pz+2*a3313*px;
		dgamma13dpx= (a1133+a1313)*pz+2*a1113*px;
		dgamma13dpz= (a1133+a1313)*px+2*a3313*pz;
		dgamma11dn = da1111dn*px2+ da1313dn*pz2 +2*da1113dn*pxz;
		dgamma33dn = da3333dn*pz2 + da1313dn*px2+2*da3313dn*pxz;
		dgamma13dn = (da1133dn+da1313dn)*pxz+ da1113dn*px2+ da3313dn*pz2;
		dg11dpx = -(gamma33-1)*(dgamma11dpx+dgamma33dpx-4*dx)*den*den+
			(dgamma33dpx-2*dx)*den;
		dg11dpz = -(gamma33-1)*(dgamma11dpz+dgamma33dpz-4*dz)*den*den+
			(dgamma33dpz-2*dz)*den;
		dg33dpx = -(gamma11-1)*(dgamma11dpx+dgamma33dpx-4*dx)*den*den+
			(dgamma11dpx-2*dx)*den;
		dg33dpz = -(gamma11-1)*(dgamma11dpz+dgamma33dpz-4*dz)*den*den+
			(dgamma11dpz-2*dz)*den;
		dg13dpx = gamma13*(dgamma11dpx+dgamma33dpx-4*dx)*den*den-
			dgamma13dpx*den;
		dg13dpz = gamma13*(dgamma11dpz+dgamma33dpz-4*dz)*den*den-
			dgamma13dpz*den;
		dg11dn = -(gamma33-1)*(dgamma11dn+dgamma33dn-2*snfact)*den*den+
			(dgamma33dn-snfact)*den;
		dg33dn = -(gamma11-1)*(dgamma11dn+dgamma33dn-2*snfact)*den*den+
			(dgamma11dn-snfact)*den;