rayt2d.c

su 的源代码库

        z = ealloc1(nr,sizeof(float*));

        /* get surfaces */
        for (ir=0; ir<nr; ++ir) {
                if (!getnparstring(ir+1,"surf",&s)) s = "0,0;99999,0";
                strcpy(t,s);
                if (!decodeSurface(t,&nxz[ir],&x[ir],&z[ir]))
                        err("Surface number %d specified "
                                "incorrectly!\n",ir+1);
        }

        /* set output parameters before returning */
        *nrPtr = nr;
        *nxzPtr = nxz;
        *xPtr = x;
        *zPtr = z;
}

/* parse one surface specification; return 1 if valid, 0 otherwise */

int decodeSurface (char *string,int *nxzPtr, float **xPtr, float **zPtr)
/**************************************************************************
decodeSurface - parse one surface specification
***************************************************************************
Input:
string          string representing surface

Output:
nxzPtr          pointer to number of x,z pairs
xPtr            array of x values for one surface
zPtr            array of z values for one surface
**************************************************************************/
{
        int nxz,ixz;
        float *x,*z;
        char *s,*t;

        s = string;
        s = strtok(s,",;\0");
        /* count x and z values, while splitting string into tokens */
        for (t=s,nxz=0; t!=NULL; ++nxz)
                t = strtok(NULL,",;\0");

        /* total number of values must be even */
        if (nxz%2) return 0;

        /* number of (x,z) pairs */
        nxz /= 2;

        /* 2 or more (x,z) pairs are required */
        if (nxz<2) return 0;

        /* allocate space */
        x = ealloc1(nxz,sizeof(float));
        z = ealloc1(nxz,sizeof(float));

        /* convert (x,z) values */
        for (ixz=0; ixz<nxz; ++ixz) {
                x[ixz] = atof(s);
                s += strlen(s)+1;
                z[ixz] = atof(s);
                s += strlen(s)+1;
        }

        /* set output parameters before returning */
        *nxzPtr = nxz;
        *xPtr = x;
        *zPtr = z;
        return 1;
}

void makesurf(float dsmax,int nr,int *nu,float **xu,float **zu,
        Surface **r)
/*****************************************************************************
Make piecewise cubic surfaces
******************************************************************************
Input:
dsmax           maximum length of surface segment
nr              number of surfaces = 2
nu              array[nr] of numbers of (x,z) pairs; u = 0, 1, ..., nu[ir]
xu              array[nr][nu[ir]] of surface x coordinates x(u)
zu              array[nr][nu[ir]] of surface z coordinates z(u)

Output:
r               array[nr] of surfaces
******************************************************************************
Notes:
Space for the nu, xu, and zu arrays is freed by this function, since
they are only used to construct the surfaces.
*****************************************************************************/
{
        int ir,iu,nuu,iuu,ns,is;
        float x,z,xlast,zlast,dx,dz,duu,uu,ds,fs,rsx,rsz,rsxd,rszd,
                *u,*s,(*xud)[4],(*zud)[4],*us;
        SurfaceSegment *ss;
        Surface *rr;

        /* allocate space for surfaces */
        *r = rr = ealloc1(nr,sizeof(Surface));

        /* loop over surfaces */
        for (ir=0; ir<nr; ++ir) {

                /* compute cubic spline coefficients for uniformly sampled u */
                u = ealloc1float(nu[ir]);
                for (iu=0; iu<nu[ir]; ++iu)
                        u[iu] = iu;
                xud = (float(*)[4])ealloc1float(4*nu[ir]);
                csplin(nu[ir],u,xu[ir],xud);
                zud = (float(*)[4])ealloc1float(4*nu[ir]);
                csplin(nu[ir],u,zu[ir],zud);

                /* finely sample x(u) and z(u) and compute length s(u) */
                nuu = 20*nu[ir];
                duu = (u[nu[ir]-1]-u[0])/(nuu-1);
                s = ealloc1float(nuu);
                s[0] = 0.0;
                xlast = xu[ir][0];
                zlast = zu[ir][0];
                for (iuu=1,uu=duu; iuu<nuu; ++iuu,uu+=duu) {
                        intcub(0,nu[ir],u,xud,1,&uu,&x);
                        intcub(0,nu[ir],u,zud,1,&uu,&z);
                        dx = x-xlast;
                        dz = z-zlast;
                        s[iuu] = s[iuu-1]+sqrt(dx*dx+dz*dz);
                        xlast = x;
                        zlast = z;
                }

                /* compute u(s) from s(u) */
                ns = 1+s[nuu-1]/dsmax;
                ds = s[nuu-1]/ns;
                fs = 0.5*ds;
                us = ealloc1float(ns);
                yxtoxy(nuu,duu,0.0,s,ns,ds,fs,0.0,(float)(nu[ir]-1),us);

                /* compute surface segments uniformly sampled in s */
                ss = ealloc1(ns,sizeof(SurfaceSegment));
                for (is=0; is<ns; ++is) {
                        intcub(0,nu[ir],u,xud,1,&us[is],&rsx);
                        intcub(0,nu[ir],u,zud,1,&us[is],&rsz);
                        intcub(1,nu[ir],u,xud,1,&us[is],&rsxd);
                        intcub(1,nu[ir],u,zud,1,&us[is],&rszd);
                        ss[is].x = rsx;
                        ss[is].z = rsz;
                        ss[is].c = rsxd/sqrt(rsxd*rsxd+rszd*rszd);
                        ss[is].s = -rszd/sqrt(rsxd*rsxd+rszd*rszd);
                }

                /* fill in surface structure */
                rr[ir].ns = ns;
                rr[ir].ds = ds;
                rr[ir].ss = ss;

                /* free workspace */
                free1float(us);
                free1float(s);
                free1float(u);
                free1float((float*)xud);
                free1float((float*)zud);

                /* free space replaced by surface segments */
                free1(xu[ir]);
                free1(zu[ir]);
        }

        /* free space replaced by surface segments */
        free1(nu);
        free1(xu);
        free1(zu);
}

/* Estimation of the Z coordinate of the topography and its normal vector */
  void zcoorTopog(float fxs,float dxs,int nxs,Surface *srf,float *sz,
		  float *nangl)
/*****************************************************************************
From the cubic spline calculation, the Z coor. of the surface are estimated
******************************************************************************
Input:
fxs             x-coordinate of the first shot (travel-time tables)
dxs             x-coordinate increment of shots (travel-time tables)
nxs             number of shots (travel-time tables)
srf             surface structure

Output:
sz              z[] coordinates of the current surface (shot positions)
nangl           array of angles that the normal form with the vertical
*****************************************************************************
Author: Trino Salinas, CSM, 1996.
*****************************************************************************/
{
        int   ik,jx,ns,ixi,is,ix;
        float x,ax,ax0,az,temp;
        float *ssx,*ssz,*sss;
        SurfaceSegment *ss;

        ns = srf[0].ns;
        ss = srf[0].ss;
        ssx = alloc1float(ns);
        ssz = alloc1float(ns);
        sss = alloc1float(ns);

        /* number of segments, segment length */
        for (jx=0; jx<ns ; jx++) {
            ssx[jx] = ss[jx].x;
            ssz[jx] = ss[jx].z;
            sss[jx] = ss[jx].s;
        }
        ixi = 0;
        for (ix=0; ix<nxs; ++ix) {
            x = fxs + ix*dxs;
            for (ik=ixi; ik<ns-1; ++ik) {
                if (ssx[ik]<=x && ssx[ik+1]>=x) {
                   az = ssz[ik] - ssz[ik+1];
                   ax0 = ssx[ik+1] - x;
                   ax = ssx[ik+1] - ssx[ik];
                   sz[ix] = ax0*az/ax+ssz[ik+1];
                   az = sss[ik] - sss[ik+1];
                   temp = ax0*az/ax + sss[ik+1];
                   nangl[ix] = asin(temp);
                   ixi = ik;
                }
            }
        }
        sz[0] = 2*sz[1] - sz[2];
        nangl[0] = 2*nangl[1] - nangl[2];
        sz[nxs-1] = 2*sz[nxs-2] - sz[nxs-3];
        nangl[nxs-1] = 2*nangl[nxs-2] - nangl[nxs-3];

         for(is=0;is<nxs;is++)
           warn("sz[%d]=%g,nangl[%d]=%g",is,sz[is],is,nangl[is]);

        free1float(ssx);
        free1float(ssz);
        free1float(sss);
}

/* Prototypes of functions used interally */
 void dfrungk(float v,float vx,float vz,float vxx,float vxz,float vzz,
     	float x,float z,float px,float pz,float p2,float q2,
      	float *dx,float *dz,float *dpx,float *dpz,float *dp2,float *dq2);
 void sum2(float h,float a1,float a2,float a3,float a4,float a5,float a6,
	float b1,float b2,float b3,float b4,float b5,float b6,
	float *c1,float *c2,float *c3,float *c4,float *c5,float *c6);

void makeRay(Geo2d geo2dv,float *vel,float *vxx,float *vxz,float *vzz,
    Raytr raytr,float a0,int *nrs,Ray *ray,int npv,float *pv) 
/*****************************************************************************
 Trace rays in uniformly sampled v(nz,nx)
*****************************************************************************
Input:
geo2dv		grid parameters of velocity 
v2		sampled velocity array
raytr		geometry parameters of ray tracing 
vxx,vxz,vzz	sampled second derivatives of velocity  
a0  		shooting angles at source point
npv		flag of velocity analysis
pv		velocity variation if npv>0

Output:
ray		pointer to ray parameters sampled at discrete ray steps
nrs		number of points on ray 
*****************************************************************************
Note:
The ray ends when it runs out of time or with the first step that is out
of the velocity boudary or out of emergence angle boundary
*****************************************************************************
Author: Zhenyue Liu, CSM, 1995.
****************************************************************************/ 
{
	int nx=geo2dv.nx,nz=geo2dv.nz,nt=raytr.nt;
	float fx=geo2dv.fx,fz=geo2dv.fz;
	float dt=raytr.dt,amin=raytr.amin,amax=raytr.amax;

	int it;
	float tzmin,tzmax,ov,lx,lz,h,h2,h6;

	float v,dvdx,dvdz,uxx,uxz,uzz,tzt,dtv;
	float x,z,px,pz,p2,q2,
	      dx,dz,dpx,dpz,dp2,dq2;
	float xt,zt,pxt,pzt,p2t,q2t,
	    dxt,dzt,dpxt,dpzt,dp2t,dq2t;


	/* velocity and derivatives at source	*/
	x = raytr.xs;
	z = raytr.zs;
	vel2Interp(geo2dv,vel,vxx,vxz,vzz,x,z,&v,&dvdx,&dvdz,&uxx,&uxz,&uzz);
	ov = 1.0/v;

	p2 = 1.0;
	q2 = 0.0;
	*nrs = nt;

/*	compute slowness at source	*/
	px = ov*sin(a0);
	pz = ov*cos(a0);

/*	first ray step output	*/
	tzt = v*pz;
	ray[0].x = x;
	ray[0].z = z;
	ray[0].px = px;
	ray[0].pz = pz;
	ray[0].q2 = q2;
	ray[0].p2 = p2;
	ray[0].v = v;
	ray[0].dvdx = dvdx;
	ray[0].dvdz = dvdz;
	if(npv) ray[0].dtv = 0;

/*	compute minimum and maximum z-component of unit ray vector	*/
	tzmin = cos(amax)-0.01;
	tzmax = cos(amin)+0.01;

/*	determine fraction steps for Rung-Kuta	*/
	h = dt;		h2 = h*0.5;	h6 = h/6.0;

	lx = fx+(nx-1)*geo2dv.dx;
	lz = fz+(nz-1)*geo2dv.dz;

/*	loop over time steps	*/
	for(it=1; it<nt; ++it){
		if(x<fx || x>lx || z<fz || z>lz ||
      		   tzt<tzmin || tzt>tzmax) {
			it -= 1; 
			break;
		}

/*	step 1 of 4th-order Rung-Kuta	*/

/*	  compute k1 = f(y0)	*/
	    dfrungk(v,dvdx,dvdz,uxx,uxz,uzz,
      		x,z,px,pz,p2,q2,
      		&dx,&dz,&dpx,&dpz,&dp2,&dq2);
          
/*	  compute y1 = y0+0.5*h*k1	*/
	    sum2(h2,x,z,px,pz,p2,q2,
      		dx,dz,dpx,dpz,dp2,dq2,
      		&xt,&zt,&pxt,&pzt,&p2t,&q2t);

/*	velocity interpolation	*/
	vel2Interp(geo2dv,vel,vxx,vxz,vzz,
      	  xt,zt,&v,&dvdx,&dvdz,&uxx,&uxz,&uzz);

/*	step 2 of 4th-order Rung-Kuta	*/

/*	  compute k2 = f(y1)	*/
	    dfrungk(v,dvdx,dvdz,uxx,uxz,uzz,
      		xt,zt,pxt,pzt,p2t,q2t,
      		&dxt,&dzt,&dpxt,&dpzt,&dp2t,&dq2t);

/*	  compute y2 = y0+0.5*h*k2	*/
	    sum2(h2,x,z,px,pz,p2,q2,
      		dxt,dzt,dpxt,dpzt,dp2t,dq2t,
      		&xt,&zt,&pxt,&pzt,&p2t,&q2t);

/*	  compute k = k1+2.0*k2	*/
	    sum2(2.0,dx,dz,dpx,dpz,dp2,dq2,
      		dxt,dzt,dpxt,dpzt,dp2t,dq2t,
      		&dx,&dz,&dpx,&dpz,&dp2,&dq2);
 
/*	velocity interpolation	*/
	vel2Interp(geo2dv,vel,vxx,vxz,vzz,
      	  xt,zt,&v,&dvdx,&dvdz,&uxx,&uxz,&uzz);

/*	step 3 of 4th-order Rung-Kuta	*/

/*	  compute k3 = f(y2)	*/
	    dfrungk(v,dvdx,dvdz,uxx,uxz,uzz,
      		xt,zt,pxt,pzt,p2t,q2t,
      		&dxt,&dzt,&dpxt,&dpzt,&dp2t,&dq2t);

/*	  compute y3 = y0+h*k2	*/
 	    sum2(h,x,z,px,pz,p2,q2,
      		dxt,dzt,dpxt,dpzt,dp2t,dq2t,
      		&xt,&zt,&pxt,&pzt,&p2t,&q2t);

/*	  compute k = k1+2.0*k2+2.0*k3	*/
	    sum2(2.0,dx,dz,dpx,dpz,dp2,dq2,
      		dxt,dzt,dpxt,dpzt,dp2t,dq2t,
      		&dx,&dz,&dpx,&dpz,&dp2,&dq2);
 
/*	velocity interpolation	*/
	vel2Interp(geo2dv,vel,vxx,vxz,vzz,
      	  xt,zt,&v,&dvdx,&dvdz,&uxx,&uxz,&uzz);

/*	step 4 of 4th-order Rung-Kuta	*/

/*	  compute k4 = f(y3)	*/
	    dfrungk(v,dvdx,dvdz,uxx,uxz,uzz,
      		xt,zt,pxt,pzt,p2t,q2t,
      		&dxt,&dzt,&dpxt,&dpzt,&dp2t,&dq2t);

/*	  compute k = k1+2.0*k2+2.0*k3+k4	*/
	    sum2(1.0,dx,dz,dpx,dpz,dp2,dq2,
      		dxt,dzt,dpxt,dpzt,dp2t,dq2t,
      		&dx,&dz,&dpx,&dpz,&dp2,&dq2);
 
/*	  compute y4 = y0+h*k/6	*/
 	    sum2(h6,x,z,px,pz,p2,q2,
      		dx,dz,dpx,dpz,dp2,dq2,
      		&x,&z,&px,&pz,&p2,&q2);
         
/*	velocity interpolation	*/
	    vel2Interp(geo2dv,vel,vxx,vxz,vzz,
      	  	x,z,&v,&dvdx,&dvdz,&uxx,&uxz,&uzz);

/*	velocity itself interpolation	*/
	if(npv)	vintp(geo2dv,pv,x,z,&dtv);

		tzt = v*pz;

/*	save ray parameters */
		ray[it].x = x;
		ray[it].z = z;
		ray[it].px = px;
		ray[it].pz = pz;
		ray[it].p2 = p2;
		ray[it].q2 = q2;
		ray[it].v = v;
		ray[it].dvdx = dvdx;
		ray[it].dvdz = dvdz;
		if(npv)  ray[it].dtv = dtv/v;
	}
	*nrs = it; 

}

/***********************************************************************/
  void dfrungk(float v,float vx,float vz,float vxx,float vxz,float vzz,
     	float x,float z,float px,float pz,float p2,float q2,
      	float *dx,float *dz,float *dpx,float *dpz,float *dp2,float *dq2)
{

	float ov,vv,c,s,v11;

		ov = 0.0*z+ 0.0*x + 1.0/v;
		vv = v*v;
		*dx = vv*px;
		*dz = vv*pz;
		*dpx = -ov*vx;
		*dpz = -ov*vz;

		c = v*pz;
		s = v*px;
		v11 = vxx*c*c-2.0*vxz*c*s+vzz*s*s;
		*dp2 = -ov*v11*q2;
		*dq2 = vv*p2;

}


/************************ c = a+h*b ***********************************/
 void sum2(float h,float a1,float a2,float a3,float a4,float a5,float a6,
	float b1,float b2,float b3,float b4,float b5,float b6,
	float *c1,float *c2,float *c3,float *c4,float *c5,float *c6)
{

	*c1 = a1+h*b1;
	*c2 = a2+h*b2;
	*c3 = a3+h*b3;
	*c4 = a4+h*b4;
	*c5 = a5+h*b5;
	*c6 = a6+h*b6;

}

/* Prototype of function used internally */
void ddt(float p,float q,float c,float s, float *dv,float v,float *d2t,
	float *cuv);

void raytime2d(Raytr raytr,Geo2d geo2dv,float *v,Geo2d geo2dt,float *t,
	int npv,float *vo,float *pv,float *tv,float *cs) 
/****************************************************************************
 raytime2d - calculate traveltimes by paraxial ray tracing
*****************************************************************************
Input:
geo2dt		grid parameters of traveltime 
geo2dv		grid parameters of velocity 
v2		sampled velocity array
raytr		geometry parameters of ray tarcing 

Output:
t		traveltime 
*****************************************************************************
Note: when traveltime has multiple values, the one has the shortest 
 ray path is chosen.