rayt2d.c

su 的源代码库

*****************************************************************************
Author: Zhenyue Liu, CSM 1995.
****************************************************************************/ 
{
	int na=raytr.na,nt=raytr.nt,nxo=geo2dt.nx,nzo=geo2dt.nz,
		nx=geo2dv.nx,nz=geo2dv.nz;
	float fxo=geo2dt.fx,fzo=geo2dt.fz,dxo=geo2dt.dx,dzo=geo2dt.dz,
		odxo=geo2dt.odx,odzo=geo2dt.odz;
	float fac=raytr.fac,dt=raytr.dt,fa=raytr.fa,da=raytr.da,xs=raytr.xs;
	float *vxx,*vxz,*vzz,*s,zs=raytr.zs;
	int nrs;
	Ray *ray;

 
 	int i,ia,ixo,izo;
	float a,xo,zo,exo,ezo;

/*	variables used for extrapolation	*/
	int it,jt,nxf,nxe,nzf,nze;
	float tc,xoc,zoc,xc,zc,r1,v0=0.0,norm2,terr,vd1,cuv,
      		sd,sind,cosd,vd,pxd,pzd,r2,t1,t2,r1x,r2x,t1x,t2x,t2xz,
      		p2d,q2d,gradv[2],d2t[3],dtvd,odt=1.0/dt,xcosd=0.0,*tvd;

	vxx = alloc1float(nx*nz);
	vxz = alloc1float(nx*nz);
	vzz = alloc1float(nx*nz);
	tvd = alloc1float(nt);
	s = alloc1float(nxo*nzo);
	ray = (Ray*)alloc1(nt,sizeof(Ray)); 

	/* compute second derivatives of velocity	*/
 	dv2(geo2dv,v,vxx,vxz,vzz); 


/*	maximum range of output points	*/
	exo = fxo+(nxo-1)*dxo;
	ezo = fzo+(nzo-1)*dzo;

/*	initialize traveltime and raypath length	*/
	for(ixo=0; ixo<nxo; ++ixo) 
		for(izo=0; izo<nzo; ++izo){
			i = izo+ixo*nzo;
			t[i] = INITIAL_T;
			s[i] = INITIAL_T;
		     if(npv) tv[i] = cs[i] = 0.0;
	}

/* 	loop over shooting angles at source point 	*/
	for(ia=0,a=fa; ia<na; ++ia,a+=da){  
	
/*		trace rays	*/
		makeRay(geo2dv,v,vxx,vxz,vzz,raytr,a,&nrs,ray,npv,pv);

/*		extropolate to grids near central rays	*/
		    v0 = ray[0].v;
		    r2 = 0.0;
		    xc = raytr.xs;
		    zc = raytr.zs;
		    sd = 0;
		    vd1 = v0;
		    if(npv) {
		      dtvd = 0.0;
		      for(it=1; it<nrs; ++it) {
			dtvd += 0.5*dt*(ray[it].dtv+ray[it-1].dtv);
			tvd[it] = ray[it].v*dtvd;
		      }
		    }
		    for (it=1; it<nrs; ++it) {
			xo = ray[it].x;
			zo = ray[it].z;
			vd = ray[it].v;
			sd = sd+0.5*dt*(vd+vd1);
			vd1 = vd;

/*		    	if the point is within the previous circle	*/
			if(r2 > (xc-xo)*(xc-xo)+(zc-zo)*(zc-zo)
      			  && it != nrs-1) continue;

/*		    	if the point is out of output range	*/
			if(xo<fxo || xo>exo || zo<fzo || zo>ezo) continue;
								
			xc = xo;
			zc = zo;
			q2d = ray[it].q2; 

/*			if caustics	*/
			if(q2d == 0.0) {
				r2 = 0.0;
				continue;
			}

			vd = ray[it].v;
			pxd = ray[it].px;
			pzd = ray[it].pz;

			p2d = ray[it].p2;
			sind = vd*pxd;
			cosd = vd*pzd;
			gradv[0] = ray[it].dvdx;
			gradv[1] = ray[it].dvdz;

/*			calculate second derivatives of traveltime*/
			ddt(p2d,q2d,cosd,sind,gradv,vd,d2t,&cuv);

/*			determine radius for extrapolation	*/
			tc = it*dt;
			terr = tc*fac;
			norm2 = sqrt(d2t[0]*d2t[0]+d2t[2]*d2t[2]+
      				2.0*d2t[1]*d2t[1]);

			r2 = terr/norm2;
			r1 = sqrt(r2);
/*			keep ray-centered coordinate system regular */
			if(r1*cuv > 0.1) r1 = 0.1/cuv;

/*			radius cannot be too large	*/
			if(r1 > 0.1*sd) r1 = 0.1*sd;

			r2 = r1*r1;

			nxf = (xc-r1-fxo)*odxo+0.9;
			if(nxf < 0) nxf = 0;
			nxe = (xc+r1-fxo)*odxo+0.1;
			if(nxe >= nxo) nxe = nxo-1;
/*       	fprintf(stderr,"x,z,t,r=%g %g %g %g\n",
			xc,zc,tc,r1);*/	 

			for(ixo=nxf; ixo<=nxe; ++ixo){
				xoc = fxo-xc+ixo*dxo;
				r2x = r2-xoc*xoc;
				if(r2x < 0) continue;

				r1x = sqrt(r2x);
				t1x = tc+xoc*pxd;
				t2x = tc*xoc*xoc*d2t[0];
				t2xz = 2.0*xoc*d2t[1];
				if(npv) xcosd = pzd+xoc*d2t[1];

				nzf = (zc-r1x-fzo)*odzo+0.9;
				if(nzf < 0) nzf = 0;
				nze = (zc+r1x-fzo)*odzo+0.1;
				if(nze>=nzo) nze = nzo-1;

				i = ixo*nzo;
				for(izo=nzf; izo<=nze; ++izo){ 
/*				    if ray path is shorter	*/
				    if(sd < s[i+izo] ) {
					zoc = fzo-zc+izo*dzo;
					t1 = t1x+zoc*pzd;
					t2 = t2x+zoc*(zoc*d2t[2]+t2xz);
					s[i+izo] = sd;
					t[i+izo] = t1*t1+tc*t2;
				      if(npv) {
					jt = (t1+0.5*t2)*odt+0.5;
					if(jt<0) jt = 0;
					if(jt>nrs-1) jt = nrs-1;
					tv[i+izo] = tvd[jt];
					cs[i+izo] = xcosd+d2t[2]*zoc;
				      }
				    }
 				}
 			}
 		}

	}

/*	square root of traveltimes */
	for (ixo=0; ixo<nxo; ++ixo){ 
		i = ixo*nzo;
		for (izo=0; izo<nzo; ++izo){  
			t[i+izo] = MAX(0.0,t[i+izo]);
			if(t[i+izo] < 999) t[i+izo] = sqrt(t[i+izo]);
  			if(npv) cs[i+izo] *= vo[i+izo];
		}
	}

/*	compute traveltime near source	*/
	nxf = (xs-fxo)/dxo-1.5;
	if(nxf<0) nxf = 0;
	nxe = (xs-fxo)/dxo+2.5;
	if(nxe>nxo-1) nxe = nxo-1;
	nzf = (zs-fzo)/dzo-0.5;
	if(nzf<0) nzf = 0;
	nze = (zs-fzo)/dzo+1.5;
	if(nze>=nzo) nze = nzo-1;
	for (ixo=nxf; ixo<=nxe; ++ixo){
		xo = fxo+ixo*dxo;
		i = ixo*nzo;
		for (izo=nzf; izo<=nze; ++izo){
			zo = fzo+izo*dzo;
			if(t[i+izo] > 999) 
			  t[i+izo] = sqrt((xo-xs)*(xo-xs)+(zo-zs)*(zo-zs))/v0;
 		}
	}

	free1float(vxx);
	free1float(vxz);
	free1float(vzz);
	free1float(s);
	free1((void*)ray);

}

void ddt(float p,float q,float c,float s, float *dv,float v,float *d2t,
	float *cuv)
/************************************************************************
ddt - calculate second derivatives of traveltime with respect to x,y
************************************************************************/
{
 	float ov2,m11,m12,m22;

	ov2 = 1.0/(v*v);
	m11 = p/q;

/*	calculate 2nd column of m	*/
	m12 = -(dv[0]*c-dv[1]*s)*ov2;
	m22 = -(dv[0]*s+dv[1]*c)*ov2;

/*	compute h*m*h^T	*/
	d2t[0] = m11*c*c+2.0*m12*c*s+m22*s*s;
	d2t[1] = (m22-m11)*c*s+m12*(c*c-s*s);
	d2t[2] = m11*s*s-2.0*m12*c*s+m22*c*c;

/*	compute the curvature of raypath	*/
	*cuv = fabs(m12)*v;

}


void trans(int nx,int nz,int nxt,int nx0,float *v,float *vt) 
{
	int ix,iz,i,i0;

	for(ix=0*nx; ix<nxt; ++ix){
		i = ix*nz;
		i0 = (ix+nx0)*nz;
		for(iz=0; iz<nz; ++iz)
			vt[i+iz] = v[i0+iz];
	}
}

void voint(Geo2d geo2dv,float *v,Geo2d geo2dt,float *ov2,int npv,float *vo)
{

	int nx=geo2dv.nx,nz=geo2dv.nz,nxo=geo2dt.nx,nzo=geo2dt.nz;
	int i,io,jx,jz,ixo,izo;
	float fx=geo2dv.fx,fz=geo2dv.fz,odx=geo2dv.odx,odz=geo2dv.odz;
	float fxo=geo2dt.fx,fzo=geo2dt.fz,dxo=geo2dt.dx,dzo=geo2dt.dz;
	float x,z,ax,az,sx,sz,temp;


	for(ixo=0,x=fxo; ixo<nxo; ++ixo,x+=dxo){
		ax = (x-fx)*odx;
		jx = ax;
		sx = ax-jx;
		if(jx < 0) jx = 0;
		if(jx > nx-2) jx = nx-2;
		for(izo=0,z=fzo; izo<nzo; ++izo,z+=dzo){
			az = (z-fz)*odz;
			jz = az;
			sz = az-jz;
			if(jz < 0) jz = 0;
			if(jz > nz-2) jz = nz-2;
			io = ixo*nzo+izo;
			i = jx*nz+jz;
		    	temp = (1.0-sx)*((1.0-sz)*v[i]+sz*v[i+1])
      			    	+sx*((1.0-sz)*v[i+nz]+sz*v[i+nz+1]);
			if(npv) vo[io] = temp;
			ov2[io] = 1.0/(temp*temp);
		}
	}
}

 
void dv2(Geo2d geo2dv,float *v,float *vxx,float *vxz,float *vzz) 
/*calculate second DeriVatives from a 2D VELocity grid by finite difference*/
/*  	input: 	velocity v 
  	output: vxx,vxz and vzz  */	 
{ 
  	int  nx=geo2dv.nx, nz=geo2dv.nz, ix, iz, i; 	 
 	float   odx=geo2dv.odx,odz=geo2dv.odz,odxx,odzz,odxz; 
   	
	odxx = odx*odx; 
	odzz = odz*odz; 
	odxz = 0.25*odx*odz; 

/*	initialize	*/
	for(ix=0; ix<nx; ++ix)    
 	    for(iz=0; iz<nz; ++iz){
 		i = ix*nz+iz;
		vxx[i] = 0.; 
		vxz[i] = 0.; 
		vzz[i] = 0.; 
	}
	
/*	finite difference 	*/
	for(ix=1; ix<nx-1; ++ix)    
 	    for(iz=1; iz<nz-1; ++iz){
 		i = ix*nz+iz;
 		vxx[i] = odxx*(v[i+nz]-2.0*v[i]+v[i-nz]); 
 		vxz[i] = odxz*(v[i+nz+1]-v[i+nz-1]-v[i-nz+1]+v[i-nz-1]); 
      		vzz[i] = odzz*(v[i+1]-2.0*v[i]+v[i-1]);
	}

}

void eiknl(Geo2d geo2dt,float *t,float *ov2,float tmax) 
/* 	compute traveltime in shadow zones by eikonal equation
  input
   t 		traveltimes form ray tracing
   ov2 		slowness squares
   tmax	 	upper limit of ordinary traveltime value
  output:
   t	traveltime (unchanged if t<=tmax)
*/
{
	int nx=geo2dt.nx,nz=geo2dt.nz,ix,iz;
	float dx=geo2dt.dx,dz=geo2dt.dz;
	float tx2,tz2,t0,tl,tr,temp,odx2,*tt1,*tt2;

	tt1 = alloc1float(nx);
	tt2 = alloc1float(nx);

	odx2 = 1.0/(dx*dx);

	for(ix=0; ix<nx; ++ix)
		tt1[ix] = t[ix*nz];

	for(iz=1; iz<nz; ++iz){ 
	    for(ix=0; ix<nx; ++ix) 
		tt2[ix] = t[ix*nz+iz];

	    for(ix=1; ix<nx; ++ix){
/*	if traveltime is abnormal and the upper is normal	*/
		t0 = tt1[ix];
		if(tt2[ix] > tmax && t0 <= tmax) {
			tl = tr = 0.0; 
			if(ix > 0) tl = tt1[ix-1];
			if(ix < nx-1) tr = tt1[ix+1];
                        tx2 = (tl-t0)*(tl-t0);
                        temp = (tr-t0)*(tr-t0);
                        if(tx2>temp) tx2 = temp;
                        temp = 0.25*(tl-tr)*(tl-tr);
                        if(tx2>temp) tx2 = temp;
                        tx2 *= odx2;

			tz2 = 0.5*(ov2[ix*nz+iz-1]+ov2[ix*nz+iz])-tx2;
			if(tz2 >= 0) tt2[ix] = t0+dz*sqrt(tz2);
		}
	    }
	    for(ix=0; ix<nx; ++ix){
		t[ix*nz+iz] = tt2[ix];
		tt1[ix] = tt2[ix];
	    }
	}

	free1float(tt1);
	free1float(tt2);

}

void vel2Interp(Geo2d geo2dv,float *v,float *vxx,float *vxz,float *vzz,
     float x,float z,
     float *u,float *ux,float *uz,float *uxx,float *uxz,float *uzz) 
/*************************************************************************
vel2Interp - Function to support interpolation of velocity and its
		derivatives	 
*************************************************************************
Input:
x,z		2D coordinate at which to interpolate 
v		velocity array(nz,nx) on unoform grids 
vxx,vxz,vzz	second derivaitve arrays(nz,nx) on uniform grids 

Output:  
u 	v(x,z)
ux  	dv/dx
uz  	dv/dz
uxx	ddv/dxdx
uxz 	ddv/dzdx
uzz	ddv/dzdz
*************************************************************************
Author: Zhenyue Liu, CSM 1995
*************************************************************************/
{
	int k0,k1,k2,k3,jx,jz,nx=geo2dv.nx,nz=geo2dv.nz;
	float fx=geo2dv.fx,dx=geo2dv.dx,odx=geo2dv.odx;
	float fz=geo2dv.fz,dz=geo2dv.dz,odz=geo2dv.odz;
	float ax,az,sx,sz,sxx,szz,a0,a1,a2,a3;
	float g0,g1,g2,g3,gx0,gx1,gx2,gx3,gz0,gz1,gz2,gz3;

/*	determine interpolate coefficients	*/
	    ax = (x-fx)*odx;
	    jx = ax;
	    if(jx<0) jx = 0;
	    if(jx>nx-2) jx = nx-2;
	    sx = ax-jx;
	    az = (z-fz)*odz;
	    jz = az;
	    if(jz<0) jz = 0;
	    if(jz>nz-2) jz = nz-2;
	    sz = az-jz;

	    sxx = 0.5*sx*(1.0-sx)*dx*dx;
	    szz = 0.5*sz*(1.0-sz)*dz*dz;

	    a0 = (1.0-sx)*(1.0-sz);
	    a1 = (1.0-sx)*sz;
	    a2 = sx*(1.0-sz);
	    a3 = sx*sz;

/*	    set the table of indices for interpolation	*/
	    k0 = nz*jx+jz;
	    k1 = k0+1;
 	    k2 = k0+nz;
	    k3 = k2+1;

	    g0 = v[k0];
	    g1 = v[k1];
	    g2 = v[k2];
	    g3 = v[k3];
	    gx0 = vxx[k0];
	    gx1 = vxx[k1];
	    gx2 = vxx[k2];
	    gx3 = vxx[k3];
	    gz0 = vzz[k0];
	    gz1 = vzz[k1];
	    gz2 = vzz[k2];
	    gz3 = vzz[k3];

/*	interpolation	*/
	    *uxx = a0*gx0+a1*gx1+a2*gx2+a3*gx3;
	    *uxz = a0*vxz[k0]+a1*vxz[k1]+a2*vxz[k2]+a3*vxz[k3];
	    *uzz = a0*gz0+a1*gz1+a2*gz2+a3*gz3;

	    *u = a0*g0+a1*g1+a2*g2+a3*g3-(sxx*(*uxx)+szz*(*uzz));
	    *ux = ((1.0-sz)*(g2-g0-sxx*(gx2-gx0)-szz*(gz2-gz0))
      		    +sz*(g3-g1-sxx*(gx3-gx1)-szz*(gz3-gz1)))*odx
      		    +(sx-0.5)*dx*(*uxx);
	    *uz = ((1.0-sx)*(g1-g0-sxx*(gx1-gx0)-szz*(gz1-gz0))
      		    +sx*(g3-g2-sxx*(gx3-gx2)-szz*(gz3-gz2)))*odz
      		    +(sz-0.5)*dz*(*uzz);
}

void vintp(Geo2d geo2dv,float *v,float x,float z,float *u) 
/*************************************************************************
 Function to support interpolation of a single fuction 
*************************************************************************
Input:
x,z		2D coordinate at which to interpolate
v		array(nz,nx) on unoform grids 
Output: 
u		v(x,z)
*************************************************************************
Author: Zhenyue Liu, CSM 1995
*************************************************************************/
{

	int k0,k1,k2,k3,jx,jz,nx=geo2dv.nx,nz=geo2dv.nz;
	float fx=geo2dv.fx,odx=geo2dv.odx;
	float fz=geo2dv.fz,odz=geo2dv.odz;
	float ax,az,sx,sz;

/*	determine interpolate coefficients	*/
	    ax = (x-fx)*odx;
	    jx = ax;
	    if(jx<0) jx = 0;
	    if(jx>nx-2) jx = nx-2;
	    sx = ax-jx;
	    az = (z-fz)*odz;
	    jz = az;
	    if(jz<0) jz = 0;
	    if(jz>nz-2) jz = nz-2;
	    sz = az-jz;

/*	    set the table of indices for interpolation	*/
	    k0 = nz*jx+jz;
	    k1 = k0+1;
 	    k2 = k0+nz;
	    k3 = k2+1;

/*	interpolation	*/
	    *u = (1.0-sx)*((1.0-sz)*v[k0]+sz*v[k1])
      			+sx*((1.0-sz)*v[k2]+sz*v[k3]);

}