rayt2d.c

su 的源代码库

/* Copyright (c) Colorado School of Mines, 2006.*/
/* All rights reserved.                       */

/* RAYT2D $Revision: 1.15 $; Date: 94/10/11 $	*/

#include "par.h"

/*********************** self documentation **********************/
char *sdoc[] = {
"									",
" RAYT2D - traveltime Tables calculated by 2D paraxial RAY tracing	",
"									",
"     rayt2d vfile= tfile= [optional parameters]			",
"									",
" Required parameters:							",
" vfile=stdin		file containning velocity v[nx][nz]		",
" tfile=stdout		file containning traveltime tables		",
"			t[nxs][nxo][nzo]				",
"									",
" Optional parameters							",
" dt=0.008		time sample interval in ray tracing		",
" nt=401		number of time samples in ray tracing		",
"									",
" fz=0			first depth sample in velocity			",
" nz=101		number of depth samples in velocity		",
" dz=100		depth interval in velocity			",
" fx=0			first lateral sample in velocity		",
" nx=101		number of lateral samples in velocity		",
" dx=100		lateral interval in velocity			",
"									",
" fzo=fz		first depth sample in traveltime table		",
" nzo=nz		number of depth samples in traveltime table	",
" dzo=dz		depth interval in traveltime table		",
" fxo=fx		first lateral sample in traveltime table	",
" nxo=nx		number of lateral samples in traveltime table	",
" dxo=dx		lateral interval in traveltime table		",
"									",
" surf=\"0,0;99999,0\"  Recording surface \"x1,z1;x2,z2;x3,z3;...\"	",
" fxs=fx		x-coordinate of first source			",
" nxs=1			number of sources				",
" dxs=2*dxo		x-coordinate increment of sources		",
" aperx=0.5*nx*dx  	ray tracing aperature in x-direction		",
"									",
" fa=-60		first take-off angle of rays (degrees)		",
" na=61			number of rays  				",
" da=2			increment of take-off angle  			",
" amin=0		minimum emergence angle 			",
" amax=90		maximum emergence angle 			",
"									",
" fac=0.01		factor to determine radius for extrapolation	",
" ek=1			flag of implementing eikonal in shadow zones 	",
" ms=10			print verbal information at every ms sources	",
" restart=n		job is restarted (=y yes; =n no)		",
" npv=0			flag of computing quantities for velocity analysis",
" if npv>0 specify the following three files				",
" pvfile=pvfile		input file of velocity variation pv[nxo][nzo]	",
" tvfile=tvfile		output file of traveltime variation tables  	",
"			tv[nxs][nxo][nzo]				",
" csfile=csfile		output file of cosine tables cs[nxs][nxo][nzo]	",
"									",
" Notes:								",
" 1. Each traveltime table is calculated by paraxial ray tracing; then 	",
"    traveltimes in shadow zones are compensated by solving eikonal	",
"    equation.								",
" 2. Input velocity is uniformly sampled and smooth one preferred.	",
" 3. Traveltime table and source ranges must be within velocity model.	",
" 4. Ray tracing aperature can be chosen as sum of migration aperature	",
"    plus half of maximum offset.					",
" 5. Memory requirement for this program is about			",
"      [nx*nz+4*mx*nz+3*nxo*nzo+na*(nx*nz+mx*nz+3*nxo*nzo)]*4 bytes	",
"    where mx = min(nx,2*(1+aperx/dx)).					",
"									",
NULL};
/*
 * Author:  Zhenyue Liu, 10/11/94,  Colorado School of Mines 
 *
 *          Trino Salinas, 01/01/96 included the option to handle nonflat
 *          reference surfaces.
 *          Subroutines from Dave Hale's modeling library were adapted in
 *          this code to define topography using cubic splines.
 *
 * References:
 *
 * Beydoun, W. B., and Keho, T. H., 1987, The paraxial ray method:
 *   Geophysics, vol. 52, 1639-1653.
 *
 * Cerveny, V., 1985, The application of ray tracing to the numerical
 *   modeling of seismic wavefields in complex structures, in Dohr, G.,
 *   ED., Seismic shear waves (part A: Theory): Geophysical Press,
 *   Number 15 in Handbook of Geophysical Exploration, 1-124.
 * 
 */
/**************** end self doc ********************************/

/* define initial value for traveltime as an arbitrary large number */
#define INITIAL_T 999999

/* Structures used internally */
/* one  ray */
typedef struct RayStruct {
	float t;		/* time */
	float px,pz;		/* slowness componets */
	float x,z;		/* x,z coordinates */
	float q2,p2;	/* Cerveny's dynamic ray tracing solution */
	float v,dvdx,dvdz;	/* velocity and its derivatives */
	float dtv;	/* traveltime variation */
} Ray;

/* ray tracing geometry parameters */
typedef struct RaytrStruct {
	float xs;		/* source lateral location */	 
	float zs;		/* source depth location */	 
	int nt;			/* number of time samples in ray tracing */
	float dt;		/* time sample interval in ray tracing  */	 
	int na;			/* number of rays */
	float fa;		/* first take-off angle of rays */
	float da;		/* increment of take-off angles */
	float amin;		/* minimum emergence angle  */	 
	float amax;		/* maximum emergence angle  */	 
	float fac;		/* factor to determine extrapolation radius*/ 
} Raytr;

/* 2D section geometry parameters */
typedef struct Geo2DStruct {
	int nx;			/* number of lateral samples  */
	float fx;		/* first lateral sample  */	 
	float dx;		/* lateral interval   */	 
	float odx;		/* 1 over dx   */	 
	int nz;			/* number of depth samples */
	float fz;		/* first depth sample  */	 
	float dz;		/* depth interval  */	
	float odz;		/* 1 over dz   */	 
}Geo2d;

/* Geometry of the recording surface */
typedef struct SurfaceSegmentStruct {
        float x;        /* x coordinate of segment midpoint */
        float z;        /* z coordinate of segment midpoint */
        float s;        /* x component of unit-normal-vector */
        float c;        /* z component of unit-normal-vector */
} SurfaceSegment;

typedef struct SurfaceStruct {
        int ns;               /* number of surface segments */
        float ds;             /* segment length */
        SurfaceSegment *ss;   /* array[ns] of surface segments */
} Surface;

/* Prototypes of functions used internally */
void makeRay (Geo2d geo2dv,float *v,float *vxx,float *vxz,float *vzz,
	Raytr raytr,float a0,int *nrs,Ray *ray,int npv,float *pv);

void raytime2d(Raytr raytr,Geo2d geo2dv,float *vt,Geo2d geo2dt,float *t,
	int npv,float *vo,float *pv,float *tv,float *cs);

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

void dv2(Geo2d geo2dv,float *v,float *vxx,float *vxz,float *vzz);

void eiknl(Geo2d geo2dt,float *t,float *ov2,float tmax);

void trans(int nx,int nz,int nxt,int nx0,float *v,float *vt);

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);

void vintp(Geo2d geo2dv,float *v,float x,float z,float *u);

void decodeSurfaces(int *nrPtr,int **nxzPtr,float ***xPtr,float ***zPtr);

int decodeSurface(char *string,int *nxzPtr,float **xPtr,float **zPtr);

void makesurf(float dsmax,int nr,int *nu,float **xu,float **zu,
             Surface **r);
void zcoorTopog(float fxs,float dxs,int nxs,Surface *srf,float *sz,
                float *nangl);


int
main(int argc, char **argv)
{
	int	na,nat,nt,nxs,nxo,nzo,nx,nz,nxt,nx0,mx,npv,nsrf,*nxzsrf;
	float	dt,xs,fxs,dxs,exs,fxo,fzo,dxo,dzo,exo,ezo,fa,ea,amin,eat,
		amax,da,fat,fac,tmax,aperx,temp,fx,fz,dx,dz,ex,ez,fxt,ext;
	float	*v,*vt,*t,*ov2,*pv=NULL,*pvt=NULL,*tv=NULL,*cs=NULL,*vo=NULL,
		**xsrf,**zsrf,*szi,*nangl;
	int	ixs,ixs0,nsize,isize,ek,ms;
	Geo2d geo2dv, geo2dt;
	Raytr raytr;
        Surface *srf;
	char *vfile="stdin",*tfile="stdout",*jpfile,*pvfile,*tvfile,*csfile;
	char *restart;
	FILE *vfp, *tfp, *jpfp, *pvfp=NULL, *tvfp=NULL, *csfp=NULL;

	/*hook up getpar to handle the parameters	*/
	initargs(argc,argv);
	requestdoc(1);

	/* get velocity      information */
	if( !getparstring("vfile",&vfile)) {
		vfp = stdin;
	} else  
		vfp = fopen(vfile,"r");
	if(!getparint("nx",&nx)) nx = 101;
	if(!getparint("nz",&nz)) nz = 101;
	if(nx<3 || nz<3) err("nx and nz must be not less than 3!\n");
	if(!getparfloat("fx",&fx)) fx = 0.;
	if(!getparfloat("fz",&fz)) fz = 0.;
	if(!getparfloat("dx",&dx)) dx = 100.;
	if(!getparfloat("dz",&dz)) dz = 100.;
	ex = fx+(nx-1)*dx;
	ez = fz+(nz-1)*dz;
	geo2dv.nx = nx;		geo2dv.fx = fx;	
	geo2dv.dx = dx;		geo2dv.odx = 1.0/dx;
	geo2dv.nz = nz;		geo2dv.fz = fz;	
	geo2dv.dz = dz;		geo2dv.odz = 1.0/dz;

	if(!getparint("nxo",&nxo)) nxo = nx;
	if(!getparint("nzo",&nzo)) nzo = nz;
	if(nxo<3 || nzo<3) err("nxo and nzo must be not less than 3!\n");
	if(!getparfloat("fxo",&fxo)) fxo = fx;
	if(!getparfloat("fzo",&fzo)) fzo = fz;
	if(!getparfloat("dxo",&dxo)) dxo = dx;
	if(!getparfloat("dzo",&dzo)) dzo = dz;
	exo = fxo+(nxo-1)*dxo;
	ezo = fzo+(nzo-1)*dzo;
	geo2dt.nx = nxo;	geo2dt.fx = fxo;	
	geo2dt.dx = dxo;	geo2dt.odx = 1.0/dxo;
	geo2dt.nz = nzo;	geo2dt.fz = fzo;	
	geo2dt.dz = dzo;	geo2dt.odz = 1.0/dzo;

	if(!getparint("nxs",&nxs)) nxs = 1;
	if(!getparfloat("fxs",&fxs)) fxs = fx;
	if(!getparfloat("dxs",&dxs)) dxs = 2*dxo;
	exs = fxs+(nxs-1)*dxs;
	if( !getparfloat("aperx",&aperx)) aperx = 0.5*(ex-fx);
	if(nxs>1) aperx += dxs;
	mx = 2.0*aperx/dx+2.99;
	if(mx>nx) mx = nx;

	if(!getparint("nt",&nt)) nt = 401;
	if(!getparfloat("dt",&dt)) dt = 0.008;
	tmax = (nt-1)*dt;
	if(!getparfloat("fa",&fa)) fa = -60.;
	if(!getparfloat("da",&da)) da = 2.;
	if(!getparint("na",&na)) na = 61;
	if(!getparfloat("amin",&amin)) amin = 0.;
	if(!getparfloat("amax",&amax)) amax = 90.;
	if(amax>180 || amin<0)  
	err("amin and amx must be within 0 to 180 degrees!\n");
	fa *= PI/180;
	da *= PI/180;
	ea = fa+(na-1)*da;
	amin *= PI/180;
	amax *= PI/180;
	if(!getparfloat("fac",&fac)) fac = 0.01;
	raytr.nt = nt;		raytr.dt = dt;
	raytr.da = da;		raytr.fac = fac;
	raytr.amin = amin;	raytr.amax = amax;

	if(!getparint("ek",&ek)) ek = 1;
	if( !getparstring("restart",&restart)) restart = "n";
	if(!getparint("ms",&ms)) ms = 10;
	if(!getparint("npv",&npv)) npv = 0;

	ixs0 = 0;
	isize = 0;
	nsize = nzo*nxo;
	if( !getparstring("tfile",&tfile)) {
		tfp = stdout;
	} else {
		if((tfp = fopen(tfile,"r"))!=NULL) {
			fclose(tfp);
			tfp = fopen(tfile,"r+");
		} else {
			tfp = fopen(tfile,"w");
		} 
		if(restart[0] =='y') {

			efseeko(tfp,(off_t) 0,SEEK_END);
			isize = (int) (eftello(tfp)/sizeof(float)/nsize);

			ixs0 = isize;
		}else {
			efclose(tfp);
			tfp = efopen(tfile,"w");
		} 
		efseeko(tfp,(off_t) (isize*nsize*sizeof(float)),SEEK_SET);
	}

	if(!getparstring("jpfile",&jpfile)) {
		jpfp = stderr;
	} else {
		jpfp = fopen(jpfile,"w");
	}
	
	if(npv){ 
	    if( !getparstring("pvfile",&pvfile)) 
		pvfile = "pvfile";
	    pvfp = fopen(pvfile,"r");
	    if( !getparstring("tvfile",&tvfile)) 
		tvfile = "tvfile";
	    tvfp = fopen(tvfile,"w");
	    if( !getparstring("csfile",&csfile)) 
		csfile = "csfile";
	    csfp = fopen(csfile,"w");
	}

	fprintf(jpfp,"\n");
	fprintf(jpfp," RAYT2D parameters\n");
	fprintf(jpfp," ================\n");
	fprintf(jpfp," vfile=%s \n",vfile);
	fprintf(jpfp," tfile=%s \n",tfile);
	fprintf(jpfp," one-way time: dt=%g nt=%d \n",dt,nt);
	fprintf(jpfp," nz=%d fz=%g dz=%g\n",nz,fz,dz);
	fprintf(jpfp," nx=%d fx=%g dx=%g\n",nx,fx,dx);
	fprintf(jpfp," nzo=%d fzo=%g dzo=%g\n",nzo,fzo,dzo);
	fprintf(jpfp," nxo=%d fxo=%g dxo=%g\n",nxo,fxo,dxo);
 	fprintf(jpfp," nxs=%d fxs=%g dxs=%g\n",nxs,fxs,dxs);
 	fprintf(jpfp," mx=%d aperx=%g \n",mx,aperx);
 	fprintf(jpfp," na=%d fa=%g da=%g\n",na,fa*180/PI,da*180/PI);
 	fprintf(jpfp," amin=%g amax=%g \n",amin*180/PI,amax*180/PI);
 	fprintf(jpfp," ek=%d ms=%d npv=%d restart=%s\n",ek,ms,npv,restart);
	if(npv)
 	  	fprintf(jpfp," pvfile=%s tvfile=%s csfile=%s\n",
			pvfile,tvfile,csfile);
	fflush(jpfp);

	/* ensure sources and output are in velocity grid	*/
	if(fx>fxs || ex<exs || fz>0) {
            warn("This condition must NOT be satisfied: fx>fxs || ex<exs || fz>0");
            warn("fx=%g fxs=%g ex=%g exs=%g fz=%g",
                    fx,fxo,ex,exo,fz);
		err("source lies outside of specified x grid  \n");
        }
	if(fx>fxo || ex<exo || fz>fzo || ez<ezo) {
            warn("This condition must NOT be satisfied: fx>fxo || ex<exo || fz>fzo || ez<ezo");
            warn("fx=%g fxo=%g ex=%g exo=%g fz=%g fzo=%g,ez=%g ezo=%g",
                    fx,fxo,ex,exo,fz,fzo,ez,ezo);
		err("output lies outside of specified x grid  \n");
        }

	/* allocate space */
	v = alloc1float(nx*nz);
	ov2 = alloc1float(nxo*nzo);
	t = alloc1float(nxo*nzo);
	if(npv) {
		pv = alloc1float(nx*nz);
		tv = alloc1float(nxo*nzo);
		vo = alloc1float(nxo*nzo);
		cs = alloc1float(nxo*nzo);
	}

	/* read velocity */
	if(fread(v,sizeof(float),nx*nz,vfp)!=nx*nz)
	    err("cannot read %d velocities from file %s",nx*nz,vfp);

	if(npv) 
	    if(fread(pv,sizeof(float),nx*nz,pvfp)!=nx*nz)
	    	err("cannot read %d values from file %s",nx*nz,pvfp);

	fprintf(jpfp," finish velocity input\n");


	/*interpolate velocity and compute slowness squares  */
	voint(geo2dv,v,geo2dt,ov2,npv,vo);

	fprintf(jpfp," begin traveltime calculation\n");
	fflush(jpfp);

	/* Estimation of the recording surface 	*/
	decodeSurfaces(&nsrf,&nxzsrf,&xsrf,&zsrf);
        makesurf(dxs*0.025,nsrf,nxzsrf,xsrf,zsrf,&srf);

        szi = alloc1float(nxs);
	nangl = ealloc1float(nxs);
        zcoorTopog(fxs,dxs,nxs,srf,szi,nangl);

	/* loop over sources */
	for(ixs=ixs0,xs=fxs+ixs0*dxs;ixs<nxs;ixs++,xs+=dxs){
		/* redefine the velocity model due to aperture	*/
		temp = fxo;
		if(xs<temp) temp = xs;
		if(xs-aperx>temp) temp = xs-aperx;
		nx0 = (temp-fx)/dx;
		fxt = fx+nx0*dx;
		temp = exo;
		if(xs>temp) temp = xs;
		if(xs+aperx<temp) temp = xs+aperx;
		nxt = 1+(int)((temp-fx)/dx+0.99)-nx0;
		ext = fxt+(nxt-1)*dx;

		vt = alloc1float(nxt*nz);
		if(npv) pvt = alloc1float(nxt*nz);

	        /* reducing velocity volume due to aperture	*/
		trans(nx,nz,nxt,nx0,v,vt);
		if(npv) trans(nx,nz,nxt,nx0,pv,pvt);

		/* determine range of take-off angles	*/
		fat = fa + nangl[ixs];
		eat = ea + nangl[ixs];
		nat = na;
		if (fat<-(PI/2)) fat = -PI/2;
		if (eat>(PI/2)) eat = PI/2;
		if(xs==fxt && fat<0) {
			fat = 0.;
			nat = eat/da+1.5;
		} else if(xs==ext && eat>0)
			nat = 1.5-fa/da;


		/* update geometry information	*/
		raytr.xs = xs;		raytr.zs = szi[ixs];
		raytr.na = nat;		raytr.fa = fat;		
		geo2dv.nx = nxt;	geo2dv.fx = fxt;
		
	        /* compute traveltime by paraxial ray tracing	*/
		raytime2d(raytr,geo2dv,vt,geo2dt,t,npv,vo,pvt,tv,cs);

		/*make up in shadow zones by eikonal equation	*/
		if(ek) eiknl(geo2dt,t,ov2,tmax);
	
		/*write traveltime	*/
		fwrite(t,sizeof(float),nxo*nzo,tfp);

		/*write quantities for velocity analysis	*/
		if(npv) {
			fwrite(tv,sizeof(float),nxo*nzo,tvfp);
			fwrite(cs,sizeof(float),nxo*nzo,csfp);
		}

		if(ixs%ms==0) {
		  fprintf(jpfp," traveltime computed at source xs=%g\n",xs);
		  fflush(jpfp);
		}
		free1float(vt);
		if(npv) free1float(pvt);
	}
 	fprintf(jpfp," finish program rayt2d\n\n");
	
	fclose(tfp);
	fclose(vfp);
	fclose(jpfp);
	
	free1float(v);
	free1float(t);
	if(npv){
		free1float(pv);
		free1float(tv);
		free1float(cs);
	}

	return(CWP_Exit());
}

/**********************************************************************
        Subroutines adapted from Dave Hale's modeling library

        decodeReflectors
        decodeReflector
        makeref                 Autor: Dave Hale, CSM, 09/17/91
**********************************************************************/

void decodeSurfaces(int *nrPtr,int **nxzPtr,float ***xPtr,float ***zPtr)
/*************************************************************************
decodeSurfaces - parse surfaces parameter string
**************************************************************************
Output:
nrPtr           pointer to nr an int specifying number of surfaces = 2
nxzPtr          pointer to nxz specifying number of (x,z) pairs defining the
                surfaces
xPtr            pointer to array[x][nr] of x values for each surface
zPtr            array[z][nr] of z values for each surface
**************************************************************************/
{
        int nr,*nxz,ir;
        float **x,**z;
	char t[6144],*s;

        /* count surfaces */
        nr = countparname("surf");
        if (nr==0) nr = 1;

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