rayt2dan.c

su 的源代码库

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

/* RAYT2DAN -- RAYTracing for P- and SV-waves in a 2D ANisotropic media*/

#include "par.h"

/*********************** self documentation **********************/
char *sdoc[] = {
"									",
" RAYT2DAN -- P- and SV-wave raytracing in 2D anisotropic media		",
"									",
" rayt2dan > ttime parameterfiles= nt= nx= nz= [optional parameters]	",
"									",
" Required Parameters:							",
" VP0file=		 name of file containing VP0(x,z)		",
" nt=                    number of time samples	for each ray		",
" nx=                    number of samples (x) for the parameter fields	",
" nz=                    number of samples (z) for the parameter fields	",
" 									",
" Optional Parameters:							",
" SV=0			 for P-waves, SV=1 for Shear waves		",
"									",
" Parameters defining the velocity field				",
" dt=0.008                 time sampling interval			",	
" fx=0                     first lateral sample (x) in parameter field  ",
" fz=0                     first lateral sample (z) in parameter field  ",
" dx=100.0                 sample spacing (x) for the parameter fields	",
" dz=100.0                 sample spacing (z) for the parameter fields	",
" 									",	
" Parameters defining the takeoff angle of a ray at a source position	",	
" 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; must be > -90 degrees	",
" amax=90                maximum emergence angle; must be < 90 degrees	",
" 									",
" Parameters defining the output traveltime table			",
" fxo=fx                 first lateral sample in traveltime table	",
" nxo=nx                 number of later samples in traveltime table	",
" dxo=dx                 lateral interval in traveltime table		",
" fzo=fz                 first depth sample in traveltime table		",
" nzo=nz                 number of depth samples in traveltime table	",
" dzo=dz                 depth interval in traveltime table		",
" fac=0.01               factor to determine the radius of extrap. 	",
" 									",
" Parameters defining the source positions			        ",
" fsx=fx                 x-coordinate of first source                   ",
" nsx=1                  number of sources                              ",
" dsx=2*dxo              x-coordinate increment of sources              ",
" aperx=0.5*nx*dx        ray tracing aperature in x-direction           ",
"									",
" Files for general anisotropic parameters confined to a vertical plane:",
" a1111file		name of file containing a1111(x,z)		",
" a1133file          	name of file containing a1133(x,z)		",
" VS0file	       	name of file containing VS0(x,z)		",
" a1113file          	name of file containing a1113(x,z)		",
" a3313file          	name of file containing a3313(x,z)		",
"									",
" For transversely isotropic media Thomsen's parameters could be used:	",
" deltafile		name of file containing delta(x,z)		",
" epsilonfile		name of file containing epsilon(x,z)		",
"									",
" if anisotropy parameters are not given the program considers		", 
" isotropic media.							", 
"									",
NULL};

/* Credits:
 *      Debashish Sarkar                  
 *		 
 *		
 *   Technical Reference:
 *
 *	Cerveny, V., 1972, Seismic rays and ray intensities 
 *	in inhomogeneous anisotropic media: 
 *	Geophys. J. R. Astr. Soc., 29, 1--13.
 *
 * 	Hangya, A., 1986, Gaussian beams in anisotropic elastic media:
 *      Geophys. J. R. Astr. Soc., 85, 473--563.
 *	
 *	Gajewski, D. and Psencik, I., 1987, Computation of high frequency 
 *	seismic wavefields in 3-D  laterally inhomogeneous anisotropic 
 * 	media: Geophys. J. R. Astr. Soc., 91, 383-411.
 *
 */

/**************** end self doc ********************************/

/* initial value for traveltime as an arbitary large number */
# define INITIAL 99999999

/* functions defined and used internally */
/* one step along ray */
typedef struct RayStepStruct {
	float t;		/* time */
	float x,z;		/* x,z coordinates */
	float q1,p1,q2,p2;	/* Cerveny's dynamic ray tracing solution */
	int kmah;		/* KMAH index */
	float c,s;		/* cos(angle) and sin(angle) */
	float v,dvdx,dvdz;	/* velocity and its derivatives */
} RayStep;

/* one ray */
typedef struct RayStruct {
	int nrs;		/* number of ray steps */
	RayStep *rs;		/* array[nrs] of ray steps */
	int nc;			/* number of circles */
	int ic;			/* index of circle containing nearest step */
	void *c;		/* array[nc] of circles */
} Ray;

/* functions (raytracer)*/
Ray *makeRay (int SV, float x0, float z0, float a0, int nt, float dt, float ft,
	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);
void freeRay (Ray *ray);

/* functions (velocity interpolation)*/
void* vel2Alloc (int nx, float dx, float fx,
	int nz, float dz, float fz, float **v);
void vel2Free (void *vel2);
void vel2Interp (void *vel2, float x, float z,
	float *v, float *vx, float *vz, float *vxx, float *vxz, float *vzz);

int main(int argc, char **argv)
{
	int nx,nz,nt,ix,iz,ia1111=1,ia1133=1,iVS0=1,ia1113=1,ia3313=1,it,
		nxf,nxe,nzf,nze,ixo,izo,i,NRS,SV,
		idelta=1,iepsilon=1,cntt,isx,ia,nxo,nzo,nsx,na;
	float sx,dx,sz,dz,dt,a,fxo,dxo,fzo,dzo,fx,fz,ft,v0,vd,vd1,px,pz,exo,r2,
		r1,r1x,t1x,t2x,t2xz,xoc,zoc,xc,zc,terr,tc,norm2,curv,x1,c1,s1,
		z1,Det,fac,r2x,v1,t1,t2,a1,aperx,*dvdx1,*dvdz1,
		ezo,xo,zo,sd,fsx,dsx,amin,amax,da,fa,**VP0,**a1111,**a3333,
		**VS0,**a1133,**a1313,**a1113,**a3313,**delta,**epsilon,
		*tt,*xx,*zz,*t,*s,*vv,*cc,*ss,**X,**Y,**N;
	char *VP0file="",*VS0file="",*a1111file="",*a1133file="",
		*a1113file="",*a3313file="",*deltafile="",*epsilonfile="",
		*tfile="stdout";
	FILE *a1111fp,*a1133fp,*VP0fp,*VS0fp,
		*a1113fp,*a3313fp,*deltafp,*epsilonfp,*tfp;
	Ray *ray;

	/* hook up getpar */
	initargs(argc,argv);
	requestdoc(0);

	/* get required parameters */
	if (!getparstring("VP0file",&VP0file)) 
			err("Must specify VP0file!");
	if (!getparint("nt",&nt)) err("Must specify nt!");
	if (!getparint("nx",&nx)) err("Must specify nx!");
	if (!getparint("nz",&nz)) err("Must specify nz!");

	/* get optional parameters */
	/* get parameters for the parameter field*/
	if (!getparint("SV",&SV)) SV=0;
	if (!getparfloat("fx",&fx)) fx = 0.0;
	if (!getparfloat("fz",&fz)) fz = 0.0;
	if (!getparfloat("dx",&dx)) dx = 100.0;
	if (!getparfloat("dz",&dz)) dz = 100.0;
	if (!getparfloat("dt",&dt)) dt = 0.008;
	if (!getparfloat("fac",&fac)) fac = 0.01;

	/* get parameters for the take-off angle of the ray at the source*/
	if (!getparfloat("amin",&amin)) amin = 0;
	if (!getparfloat("amax",&amax)) amax = 90;
	if (!getparfloat("fa",&fa)) fa=-60;
	if (!getparint("na",&na)) na=61;
	if (!getparfloat("da",&da)) da=2;

	/* get parameters defining the output traveltime table*/
	if (!getparfloat("fxo",&fxo)) fxo = fx;
	if (!getparint("nxo",&nxo))   nxo = nx;
	if (!getparfloat("dxo",&dxo)) dxo = dx;
	if (!getparfloat("fzo",&fzo)) fzo = fz;
	if (!getparint("nzo",&nzo))   nzo = nz;
	if (!getparfloat("dzo",&dzo)) dzo = dz;

        /* get parameters defining source positions */
        if (!getparfloat("fsx",&fsx))  fsx = fx;
        if (!getparint("nsx",&nsx)) nsx = 1;
        if (!getparfloat("dsx",&dsx)) dsx = 2*dxo;
        if (!getparfloat("aperx",&aperx)) aperx = 0.5*nx*dx;

	/* output traveltime file */
	if (!getparstring("tfile",&tfile)) {
		tfp = stdout;
	} else {
		tfp = fopen(tfile,"w");
	}	
	
	/* ensure takeoff point is within model */
        if (fsx<fx || fsx>(fx+dx*(nx-1)) || (fsx+dsx*(nsx-1))<fx 
			|| (fsx+dsx*(nsx-1))>(fx+dx*(nx-1))) 
		err("Shot point locations exceed model dimensions");	
		
	/* allocate workspace */
	VP0 = alloc2float(nz,nx);
	VS0 = alloc2float(nz,nx);
	a1111 = alloc2float(nz,nx);
	a3333 = alloc2float(nz,nx);
	a1133 = alloc2float(nz,nx);
	a1313 = alloc2float(nz,nx);
	a1113 = alloc2float(nz,nx);
	a3313 = alloc2float(nz,nx);
	delta = alloc2float(nz,nx);
	epsilon = alloc2float(nz,nx);
	tt = alloc1float(nt);
	xx = alloc1float(nt);
	zz = alloc1float(nt);
	t  = alloc1float(nxo*nzo);
	s  = alloc1float(nxo*nzo);
	cc = alloc1float(nt);
	ss = alloc1float(nt);
	vv = alloc1float(nt);
	dvdx1 = alloc1float(nt);
	dvdz1 = alloc1float(nt);
	N  = alloc2float(2,2);
	X  = alloc2float(2,2);
	Y  = alloc2float(2,2);

	/* anisotropic parameters */
	if (!getparstring("a1111file",&a1111file)) 
			ia1111=0;
	if (!getparstring("a1133file",&a1133file)) 
			ia1133=0;
	if (!getparstring("VS0file",&VS0file)) 
			iVS0=0;
	if (!getparstring("a1113file",&a1113file)) 
			ia1113=0;
	if (!getparstring("a3313file",&a3313file)) 
			ia3313=0;
	if (!getparstring("deltafile",&deltafile)) 
			idelta=0;
	if (!getparstring("epsilonfile",&epsilonfile)) 
			iepsilon=0;


	/* read required velocities*/
	if ((VP0fp=fopen(VP0file,"r"))==NULL)
		err("error opening VP0file=%s!\n",VP0file);
	if (fread(VP0[0],sizeof(float),nz*nx,VP0fp)!=nz*nx)
		err("error reading VP0file=%s!\n",VP0file);

	/* compute a3333 */
	for (ix=0; ix<nx; ++ix)
		for (iz=0; iz<nz; ++iz)
			a3333[ix][iz] = VP0[ix][iz]*VP0[ix][iz];
	

	/* read velocities*/
	if(ia1111){
		if ((a1111fp=fopen(a1111file,"r"))==NULL)
			err("error opening a1111file=%s!\n",a1111file);
		if (fread(a1111[0],sizeof(float),nz*nx,a1111fp)!=nz*nx)
			err("error reading a1111file=%s!\n",a1111file);
	}else{
		for (ix=0; ix<nx; ++ix)
			for (iz=0; iz<nz; ++iz)
				a1111[ix][iz] = a3333[ix][iz];
	}		
	
	if(ia1133){
		if ((a1133fp=fopen(a1133file,"r"))==NULL)
			err("error opening a1133file=%s!\n",a1133file);
		if (fread(a1133[0],sizeof(float),nz*nx,a1133fp)!=nz*nx)
			err("error reading a1133file=%s!\n",a1133file);
	}else{
		for (ix=0; ix<nx; ++ix)
			for (iz=0; iz<nz; ++iz)
				a1133[ix][iz] = .4*a3333[ix][iz];
	}

	if(iVS0){
		if ((VS0fp=fopen(VS0file,"r"))==NULL)
			err("error opening VS0file=%s!\n",VS0file);
		if (fread(VS0[0],sizeof(float),nz*nx,VS0fp)!=nz*nx)
			err("error reading VS0file=%s!\n",VS0file);
	}else{
		for (ix=0; ix<nx; ++ix)
			for (iz=0; iz<nz; ++iz)
				VS0[ix][iz] = .5477*VP0[ix][iz];
	}

	/* compute a1313 */
	for (ix=0; ix<nx; ++ix)
		for (iz=0; iz<nz; ++iz)
			a1313[ix][iz] = VS0[ix][iz]*VS0[ix][iz];

	if(ia1113){
		if ((a1113fp=fopen(a1113file,"r"))==NULL)
			err("error opening a1113file=%s!\n",a1113file);
		if (fread(a1113[0],sizeof(float),nz*nx,a1113fp)!=nz*nx)
			err("error reading a1113file=%s!\n",a1113file);
	}else{
		for (ix=0; ix<nx; ++ix)
			for (iz=0; iz<nz; ++iz)
				a1113[ix][iz] = 0;
	}

	if(ia3313){
		if ((a3313fp=fopen(a3313file,"r"))==NULL)
			err("error opening a3313file=%s!\n",a3313file);
		if (fread(a3313[0],sizeof(float),nz*nx,a3313fp)!=nz*nx)
			err("error reading a3313file=%s!\n",a3313file);
	}else{
		for (ix=0; ix<nx; ++ix)
			for (iz=0; iz<nz; ++iz)
				a3313[ix][iz] = 0;
	}

	/* if specified read Thomsen's parameters*/
	if(idelta==1 || iepsilon==1) {
		if(idelta){
			if ((deltafp=fopen(deltafile,"r"))==NULL)
				err("error opening deltafile=%s!\n",deltafile);
			if (fread(delta[0],sizeof(float),nz*nx,deltafp)!=nz*nx)
				err("error reading deltafile=%s!\n",deltafile);
		}else{
			for (ix=0; ix<nx; ++ix)
				for (iz=0; iz<nz; ++iz)
					delta[ix][iz] = 0;
		}
		if(iepsilon) {
			if ((epsilonfp=fopen(epsilonfile,"r"))==NULL)
			err("error opening epsilonfile=%s!\n",epsilonfile);
			if (fread(epsilon[0],sizeof(float),nz*nx,epsilonfp)!=nz*nx)
			err("error reading epsilonfile=%s!\n",epsilonfile);
		}else{
			for (ix=0; ix<nx; ++ix)
				for (iz=0; iz<nz; ++iz)
					epsilon[ix][iz] = 0;
		}

		for (ix=0; ix<nx; ++ix)
			for (iz=0; iz<nz; ++iz){
			a1111[ix][iz] = a3333[ix][iz]+2*epsilon[ix][iz]
					*a3333[ix][iz];
			a1133[ix][iz] = sqrt(2*delta[ix][iz]*a3333[ix][iz]*
				(a3333[ix][iz]-a1313[ix][iz])+(a3333[ix][iz]-
				a1313[ix][iz])*(a3333[ix][iz]-a1313[ix][iz]))
					-a1313[ix][iz];
			}
	}

	/* iterate over sources	*/
        for(isx=0;isx<nsx;isx++) {
	sx=fsx+isx*dsx;

	/* initialize time and raypath length*/
	for (ixo=0;ixo<nxo;ixo++)
		for(izo=0;izo<nzo;izo++) {
		     i = izo+ixo*nzo;
		     t[i] = INITIAL;
		     s[i] = INITIAL;
		}

		/* iterate over takeoff angle */
		for(ia=0;ia<na;ia++) {
				a=fa+ia*da;	
				printf("sx = %f angle = %f\n",sx,a);
				/* Assume 0s as initial time */
        			ft=0; sz=0.0;
				ray = makeRay(SV,sx,sz,a,nt,dt,ft,a1111,a3333,
					a1133,a1313,a1113,a3313,nx,
					dx,fx,nz,dz,fz,amax,amin);
					/* store the ray */
					for(cntt=0;cntt<nt;cntt++){
						tt[cntt]=ray->rs[cntt].t;
						xx[cntt]=ray->rs[cntt].x;
						zz[cntt]=ray->rs[cntt].z;
						cc[cntt]=ray->rs[cntt].c;
						ss[cntt]=ray->rs[cntt].s;
						vv[cntt]=ray->rs[cntt].v;
						dvdx1[cntt]=ray->rs[cntt].dvdx;
						dvdz1[cntt]=ray->rs[cntt].dvdz;
					      }
		NRS=ray->nrs;

		/* free ray */
		freeRay(ray);

		/* Extrapolate traveltimes around the central ray 
		(Similar to Zhenyue Liu's code, rayt2d)*/
		v0 = vv[0];
		r2 = 0.0;
		xc = sx;
		zc = sz;
		sd = 0;
		vd1= v0;
		
		/* trace the auxillary ray */		
		a1=a+0.1;
		ray=makeRay(SV,sx,sz,a1,NRS,dt,ft,a1111,a3333,a1133,a1313,
			a1113,a3313,nx,dx,fx,nz,dz,fz,amax,amin);

		for (it=1;it<NRS;++it) {
			xo = xx[it];
			zo = zz[it];
			vd = vv[it];
			sd = sd+0.5*dt*(vd+vd1);
			vd1= vd;
			curv = fabs(-(dvdx1[it]*cc[it]-dvdz1[it]*ss[it]))/vv[it];