elaray.c

su 的源代码库

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

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

/* ELARAY: $Test Release: 1.2 $ ; $Date: 1997/03/05 16:49:53 $	*/


#include "par.h"
#include "tri.h"
#include "elastic.h"

#define diprint(expr) printf(#expr " = %i\n",expr)
typedef struct RSStruct {
	float sigma,x,z,px,pz,t;
	float q1,p1,q2,p2;
	float vgx,vgz; 
	float g11,g13,g33;    
	int kmah,nref,mode;
	EdgeUse *eu;
        float atten;
        float ampli,ampliphase;
	Face *f;
	struct RSStruct *rsNext;
	float rho;
} RayStep;
/*********************** self documentation **********************/
char *sdoc[] = {
" ELARAY - ray tracing for elastic anisotropic models",
"                                              				",
" elaray <modelfile >rayends [optional parameters]			",
"									",
" Optional Parameters:							",
" xs=(max-min)/2 x coordinate of source (default is halfway across model)",
" zs=min         z coordinate of source (default is at top of model)	",
" nangle=101     number of takeoff angles				",
" fangle=-45     first takeoff angle (in degrees)			",
" langle=45      last takeoff angle (in degrees)			",
" nxz=101        number of (x,z) in optional rayfile (see notes below)	",
" mode=0         shoot P-rays						",
"	      =1 shoot SV-rays						",
"	      =2 shoot SH-rays						",
" prim        =1 only reflected rays are plotted 		",     
"             =0 only direct hits are displayed  			",
" refseq=1,0,0   index of reflector followed by sequence of:		",
"		 transmission(0)					",
"		 reflection (1)						",
"		 transmission with mode conversion (2)			",					
"		 reflection with mode conversion (3)			",					
"                ray stops(-1).						",
" krecord        if specified, only rays incident at interface with index",
"                krecord are displayed and stored			",
" f0=1	         force impact strenght					",
" fdip=0         force dip with respect to vertical			",
" fazi=0         force azimuth with respect to positive x-axis 		",
" reftrans=0	 =1 include reflec/transm coeff(currently only for P)	",
" rayfile        file of ray x,z coordinates of ray-edge intersections	",
" wavefile       file of ray x,z coordinates uniformly sampled in time	",
" nt=		 number of (x,z) in optional wavefile (see notes below)	",
" tw=		 traveltime associated with wavefront (alternative to nt)",	
" infofile       ASCII-file to store useful information 		",
" outparfile     contains parameters for the plotting software. 	",
"                default is <outpar> 					",
" NOTES:								",
" The rayends file contains ray parameters for the locations at which	",
" the rays terminate.  							",
"									",
" The rayfile is useful for making plots of ray paths.			",
" nxz should be larger than twice the number of triangles intersected	",
" by the rays.								",
"									",
" The wavefile is useful for making plots of wavefronts.		",
" The time sampling interval in the wavefile is tmax/(nt-1),		",
" where tmax is the maximum time for all rays. Alternatively, ",
" one wavefront at time tw is plotted.	",
"									",
" The infofile is useful for collecting information along the		",
" individual rays. 							",
" The outparfile stores information used for the plotting software	",
"									",
NULL};

/*
 * AUTHORS:  Andreas Rueger, Colorado School of Mines, 01/02/95
  The program is based on :
 	        gbray.c, AUTHOR: Andreas Rueger, 08/12/93
 	       	sdray.c, AUTHOR Dave Hale, CSM, 02/26/91
 */
/**************** end self doc ***********************************/


#define TAPER 0.9



/* prototypes for functions defined and used internally */
void shootRaysAni (Model *m, float xs, float zs, int nangle, int reftrans,
	float dangle, float fangle, RefCheck *rc, int kk, RayStep *rsp[], 
	RayEnd re[], int *nrefseq,FILE *ifp, int mode, float fdip, float f0,
	float fazi);
void writeRays (Model *m,FILE *rayfp, int nxz, int nray, RayStep *rs[], RayEnd  
        re[], int krecord, int prim, FILE *ifp, FILE *outparfp);
void writeWaves (FILE *wavefp, int nt, int nray, RayStep *rs[]);
void writeWaves2 (FILE *wavefp, float tw, int nray, RayStep *rs[]);
static RayStep* RayInAnTri (RayStep *rs);
static RayStep* RayAcrossEdge (RayStep *rs,FILE *ifp, int conv,FILE *junkfp,
	int reftrans);
static RayStep* reflectRay (RayStep *rs,FILE *ifp, int conv,FILE *junkfp,
	int reftrans);
int testIfSourceOnEdge(Face *tris, float *zs, float *xs);
static void gvelreal (float a1111, float a3333, float a1133, float a1313,
	float a1113, float a3313, float px, float pz, float *vgx,  float *vgz,
	float *g11n, float *g13n, float *g33n);
void polar(float px, float pz, float g11, float g13,
	float g33, float *polx, float *polz, int mode );
int findnewmode(int mode , int* newmode, int conv, int mindex);
int findqPqSV(float s, float c, float pl, float a1111, float a3333,
	float a1133,float a1313, float a1113, float a3313, int mode, float
	*pxnew, float *pznew, float *vgx, float *vgz, float *g11, float *g13,  
	float *g33, float rt, FILE *ifp);


/* the main program */
int main (int argc, char **argv)
{
	Model *m;
	int i,ii,nangle,nxz,nt,kk,nseq,krecord,prim;
	int **hitseqa,mode,reftrans,*nrefseq;
	float xs,zs,fangle,langle,dangle,f0,fdip,fazi;
	/* float t0,t1,t2,tw; */
        float tw;
	
	char *rayfile,*wavefile,*infofile,*outparfile;
	FILE *rayfp=NULL,*wavefp=NULL,*ifp=NULL;
        FILE *outparfp=NULL;

	RayStep **rsp;
	RayEnd *re;	
        RefCheck *rc;

        /* t0=cpusec(); */

	/* hook up getpar to handle the parameters */
	initargs(argc,argv);
	requestdoc(0);
	
	/* read model */
	m = readModel(stdin);

	/* check raymode */
	if (!getparint("mode",&mode)) mode = 0;
	
	if (mode != 0 && mode != 1 && mode != 2)
		err("ERROR Wrong raymode defined \n");


	/* get optional parameters */
	if (!getparfloat("xs",&xs)) xs = 0.5*(m->ymin+m->ymax);
	if (!getparfloat("zs",&zs)) zs = m->xmin;
	if (!getparint("nangle",&nangle)) nangle = 101;
	if (!getparint("reftrans",&reftrans)) reftrans = 0;
	if (!getparfloat("fangle",&fangle)) fangle = -45.0;
	if (!getparfloat("langle",&langle)) langle = 45.0;
	if (!getparfloat("f0",&f0)) f0 = 1.0;
	if (!getparfloat("fdip",&fdip)) fdip = 0;
	if (!getparfloat("fazi",&fazi)) fazi = 0;
	if (getparstring("rayfile",&rayfile)) rayfp = efopen(rayfile,"w");
	if (getparstring("infofile",&infofile))ifp = efopen(infofile,"w");
	if (!getparint("nxz",&nxz)) nxz = 101;
	if (getparstring("wavefile",&wavefile)) wavefp = efopen(wavefile,"w");
	if (!getparint("nt",&nt)) nt = INT_MAX;
	if (!getparfloat("tw",&tw)) tw = FLT_MAX;
	if (!getparint("krecord",&krecord)) krecord = INT_MAX;
	if (!getparint("prim",&prim)) prim = INT_MAX;
	if (getparstring("outparfile",&outparfile)) {
                 outparfp = efopen(outparfile,"w");
        } else {
                 outparfp = efopen("outpar","w");
        }

	/* how many interfaces are of interest */
        kk=countparname("refseq");

        /* errorchecks */
        if(prim == 1 && kk == 0)
           	 err("\n first define target reflector \n ");

        /* allocate space for raychecks */
	rc=(RefCheck*)emalloc(kk*sizeof(RefCheck));

	/* Warning */
	if(reftrans != 0){
	  warn("\n Refl/Transm Coeff. need more work and testing !!!!\n");
	}
	
        /* allocate space for hitseq */
	if(kk>0){
	  hitseqa=(int**)emalloc(kk*sizeof(int*));
	  nrefseq=alloc1int(kk);
	} else {
	  hitseqa=(int**)emalloc(sizeof(int*));
	  nrefseq=alloc1int(1);
  	}

        if (ifp!=NULL) {
	 fprintf(ifp,"\n THIS FILE CONTAINS RAY TRACING INFORMATION \n"); 		
	 fprintf(ifp,"\n *********************************************\n");
 	 fprintf(ifp," defined reflection/transmission sequences \n");
	 fprintf(ifp," *********************************************\n");
        }

	/* determine interface modes */
        for(i=0;i<kk;i++){

	  /* how many parameters are defined for <refseq> ? */
	  nseq=countnparval(i+1,"refseq");
	  nrefseq[i]=nseq-1;

	  hitseqa[i]=alloc1int(nseq);

	  /* get the address and strip off the interface index */
	  getnparint(i+1,"refseq",hitseqa[i]);
	  rc[i].k=hitseqa[i][0];
	  rc[i].nhits=0;

	  /* now get the defined sequence of trans/ref/stops */
          rc[i].hitseq=hitseqa[i]+1;


          if (ifp!=NULL){
		fprintf(ifp,"\n interface %i :",rc[i].k);
		for(ii=0;ii<nseq-1;ii++){
		  fprintf(ifp,"%i ",rc[i].hitseq[ii]);
		}
          }
        }
        if (ifp!=NULL)
	 fprintf(ifp,
		"\n interfaces without defined refseq are transmitting \n");

	
	/* convert angles to radians and determine angle increment */
	fangle *= PI/180.0;
	langle *= PI/180.0;
	fdip   *= PI/180.0;
	fazi   *= PI/180.0;

	if(nangle==1){
	   dangle=0;
	} else {
	   dangle = (langle-fangle)/(nangle-1);
	}

	/* allocate space for ray step lists and ray ends */
	rsp = (RayStep**)alloc1(nangle,sizeof(RayStep*));
	re = (RayEnd*)alloc1(nangle,sizeof(RayEnd));

	/* shoot rays */
	shootRaysAni(m,xs,zs,nangle,reftrans,dangle,fangle,rc,kk,rsp,
			re,nrefseq,ifp,mode,fdip,f0,fazi);
	
        /* t0 = cpusec() - t0;*/

	/* if requested, write rays to rayfile */
        if (rayfp!=NULL) writeRays(m,rayfp,nxz,nangle,rsp,re,krecord,
		prim,ifp,outparfp);

	/* if requested, write waves to wavefile */
	if (wavefp!=NULL) {
		if(nt==INT_MAX && tw!=FLT_MAX)
			writeWaves2(wavefp,tw,nangle,rsp);
		else if(nt!=INT_MAX && tw==FLT_MAX)
			writeWaves(wavefp,nt,nangle,rsp);
		else 
			fprintf(stderr,"\n define either nt or tw \n");
	}
			

	/* t2 = cpusec()- t0;*/

	
	/* write ray ends */
	if (efwrite(re,sizeof(RayEnd),nangle,stdout)!=nangle)
		err("Error writing ray ends to stdout!\n");

	/* t1 = cpusec()- t2 +t0; */

	/* fprintf(stderr,"\n ray tracing cpu time= %g sec \n",t1);
	fprintf(stderr,"\n plotting cpu time= %g sec \n\n",t2); */


	return 1;
}

void shootRaysAni (Model *m, float xs, float zs, int nangle, int reftrans, 
	float dangle, float fangle, RefCheck *rc, int kk, RayStep *rsp[],  
	RayEnd re[], int *nrefseq, FILE *ifp, int mode, float fdip,
        float f0, float fazi)
/*****************************************************************************
Shoot rays in anisotropic models
******************************************************************************
Input:
m		triangulated model
xs		horizontal coordinate of source
zs		vertical coordinate (depth) of source
nangle		number of ray takeoff angles
dangle		increment in ray takeoff angle
fangle		first ray takeoff angle
rc      	pointer to reflector structure
nrefseq		pointer to refseq-size
kk	        number of reflectors 
mode		ray mode at takeoff
fdip		force direction with vertical
fazi		force azimuth with positive x-axis
f0		magnitude of impact source
reftrans	=1 refl/transm coeff included

Output:
rsp		array[nangle] of pointers to linked list of RaySteps
re		array[nangle] of RayEnds
******************************************************************************
Notes:
Rays are traced from triangle to triangle.
Whenever a ray crosses a triangle edge, two RaySteps are added
to the linked list of RaySteps, one for each side of the edge.
The RaySteps are useful if ray parameters along the entire raypath
are required.  The RayEnds are useful if parameters at the ray
endpoint are required.
******************************************************************************
AUTHORS:  Andreas Rueger, Colorado School of Mines, 02/02/94
  The program is based on :
 	(gbray.c)shootRays.c, AUTHOR: Andreas Rueger, 08/12/93
 	(sdray.c)shootRays.c, AUTHOR Dave Hale, CSM, 02/26/91
******************************************************************************/
{
	int iangle,refl,stop,i,temp,nameref,free;
	int refmod,transmod,*ireset,mindex;
	float a1111,a3333,a1133,a1313,a1113,a3313;
	float a1212,a1223,a2323,polx,polz,rhob;
	float s,ss,c,cc,sc,px,pz,sa,soa,angle;
	float gamm11,gamm33,gamm13,sqr,vp2,vp,ovp;
	float g11,g33,g13,vgx,vgz,fx,fz,fy;
	FILE *junkfp;

	Tri *tris;
	FaceAttributes *fa;
	EdgeAttributes *ea;
	RayStep *rs,*rslast;
	
	ireset = alloc1int(kk);

	/* junkfp=efopen("trans.dat","w");*/
        junkfp=0;
	


	/* initializations */
	px=pz=0.0;
	nameref=0;

	/* determine triangle in model containing source location */
	tris = insideTriInModel(m,NULL,zs,xs);
	temp = 1;
	i = 0;

	while(temp == 1){

		temp = testIfSourceOnEdge(tris,&zs,&xs);
		if(temp){
		    tris = insideTriInModel(m,NULL,zs,xs);
		    ++i;
		    if(i==5)
			err("\n ERROR in testIfSourceOnEdge\n");
		}
	}	

        
	/* determine stiffness at source location */
	fa = tris->fa;

	/* get stiffness at source point  */

	a1111 = fa->a1111;
	a3333 = fa->a3333;
	a1133 = fa->a1133;
	a1313 = fa->a1313;
	a1113 = fa->a1113;
	a3313 = fa->a3313;
	a1212 = fa->a1212;
	a1223 = fa->a1223;
	a2323 = fa->a2323;
	rhob  = fa->rho;
	mindex  = fa->mindex;

	/* error check */
	if(a3333 < 0) 
		err("\n ERROR: a3333 cannot be negative !! \n");
	if(a1313 < 0) 
		err("\n ERROR: a1313 cannot be negative !! \n");
	if(a1212 < 0) 
		err("\n ERROR: a1212 cannot be negative !! \n");
	if(a2323 < 0) 
		err("\n ERROR: a1313 cannot be negative !! \n");

	/* isotropic case P-rays*/
	if(mindex == 0 && mode == 0) mode=0;

	/* isotropic case SV-rays*/
	else if(mindex == 0 && mode == 1) mode=1;

	/* isotropic case SH-rays*/
	else if(mindex == 0 && mode == 2) mode=2;

	/* anisotropic case qP-rays*/
	else if(mindex != 0 && mode == 0) mode=3;

	/* anisotropic case qSV-rays*/
	else if(mindex != 0 && mode == 1) mode=4;

	/* anisotropic case qSH-rays*/
	else if(mindex != 0 && mode == 2) mode=5;


	/* force components */
	fx=f0*sin(fdip)*cos(fazi);
	fz=f0*cos(fdip);
	fy=f0*sin(fdip)*sin(fazi);

	if (ifp!=NULL){
	  fprintf(ifp,"\n ***** Force components *****\n");
	  fprintf(ifp,"\n fx=%g \t fy=%g \t fz=%g \n",fx,fy,fz);
	}

	/* loop over takeoff angles */
	for (iangle=0,angle=fangle; iangle<nangle; ++iangle,angle+=dangle) {

		/************** initialize ray tracing *************/

		c = cos(angle);
		s = sin(angle);

		/* isotropic case P-ray*/
		if(mode == 0 ){
			sa = sqrt(a3333);
			soa=1/sa;

			px  = s*soa;
			pz  = c*soa;
			vgx = px*a3333;
			vgz = pz*a3333;
			g11 = vgx*px;
			g33 = vgz*pz;
			g13 = vgx*pz;

		/* isotropic case SV-ray*/
		} else if(mode == 1){
			sa = sqrt(a1313);
			soa=1/sa;

			px  = s*soa;
			pz  = c*soa;
			vgx = px*a1313;
			vgz = pz*a1313;
			g11 = vgz*pz;
			g33 = vgx*px;
			g13 = -vgx*pz;

		/* isotropic case SH-ray */
		} else if(mode == 2){
			sa = sqrt(a1313);
			soa=1/sa;

			px  = s*soa;
			pz  = c*soa;
			vgx = px*a1313;
			vgz = pz*a1313;
			g11 = 0;
			g33 = 0;
			g13 = 0;


		/* anisotropic case qP/qSV-ray*/
		} else if(mode == 3 || mode == 4){
			cc = c*c;
			ss = s*s;
			sc = s*c;

			/*computing phase velocity for angle */
			gamm11 = a1111*ss+ a1313*cc +2*a1113*sc;
			gamm33 = a3333*cc + a1313*ss+2*a3313*sc;
			gamm13 = (a1133+a1313)*sc+ a1113*ss+ a3313*cc;
			sqr    = sqrt((gamm11+gamm33)*(gamm11+gamm33)-
				4*(gamm11*gamm33-gamm13*gamm13));
			vp2    = gamm11+gamm33;
			vp2    =((mode == 3)? vp2+sqr : vp2-sqr);

			/* error check */
			if(vp2 < 0) 
			   err("\n ERROR: ray initialization \n");

			vp     = sqrt(vp2*.5);
			ovp    = 1/vp;
			px     = s*ovp;
			pz     = c*ovp;

			gvelreal (a1111,a3333,a1133,a1313,a1113,a3313,px,
				pz,&vgx,&vgz,&g11,&g13,&g33);

		/* anisotropic case SH-ray*/
		} else if(mode == 5){
			cc = c*c;
			ss = s*s;
			sc = s*c;

			/*computing phase velocity for angle */
			sqr    = a1212*ss+ 2*a1223*sc + a2323*cc;

			/* error check */
			if(sqr < 0) 
				err("\n ERROR: ray initialization \n");

			vp     = sqrt(sqr);
			ovp    = 1/vp;
			px     = s*ovp;
			pz     = c*ovp;
				
			vgx = a1212*px + a1223*pz;
			vgz = a1223*px + a2323*pz;
			
			g11=g13=g33=0;

		} else 
			err("\n ERROR: wrong ray mode in initialization \n");


		/* get polariation with respect to ray direction */
		if( mode != 2 && mode !=5)