elaray.c

su 的源代码库

			polar(px,pz,g11,g13,g33,&polx,&polz,mode);



		/* initialize linked list of ray steps */
		rs = rsp[iangle] = (RayStep*)emalloc(sizeof(RayStep));
		rs->sigma = 0.0;
		rs->x = xs;
		rs->z = zs;
		rs->px = px;
		rs->pz = pz;
		rs->q1 = 1.0;
		rs->p1 = 0.0;
		rs->q2 = 0.0;
		rs->p2 = 1.0;
		rs->t = 0.0;
		rs->kmah = 0;
		rs->nref = 0;
		rs->eu = NULL;
		rs->f = tris;
                rs->ampliphase = 0;
                rs->vgx = vgx;
                rs->vgz = vgz;
                rs->g11 = g11;
                rs->g13 = g13;
                rs->g33 = g33;
                rs->rho = rhob;

		rs->rsNext = NULL;
		rs->mode = mode;
		
		/******* initialize ray amplitude **********/
		if(mode==2 || mode==5){
			rs->ampli = fy;
		} else {
                	rs->ampli = fx*polx + fz*polz;
		}

                /* initialize iresets for new ray */
                for(i=0;i<kk;i++)
		       ireset[i] = 0;


		if (ifp!=NULL){
			fprintf(ifp,"\n"
				"****************************************\n"
				"RAY WITH TAKEOFF ANGLE: %g (degrees)\n"
				"****************************************\n",
				angle*180.0/PI);
			fprintf(ifp,"\n initial values : \n");
			fprintf(ifp,"\n"
				"px=%g \t pz=%g \t vgx=%g \t vgz=%g \n",
				px,pz,vgx,vgz);
			fprintf(ifp,"g11=%g \t g33=%g \t g13=%g \t ampli=%g\n",
				g11,g33,g13,rs->ampli);
			fprintf(ifp,"xs=%g \t zs=%g \t mode=%i \t rho=%g\n ",
				xs,zs,mode,rhob);
		}
			
		/* trace ray */
		do {
			/* trace ray in triangle */
			rs = RayInAnTri(rs);
			
			/* remember last ray step */
			rslast = rs;
			
			/* edge attributes */
			ea = rs->eu->e->ea;
                       

			/* null edges are transmitting */
			refl = stop = temp = refmod = transmod = free= 0;

			if (ea!=NULL)
			{

			        /* check if edge is reflecting or stopping */
				for (i=0,temp=ea->k;i<kk; ++i) {

					if (temp==rc[i].k) {

						if (rc[i].hitseq[0] == 1)
							refl = 1;
						else if (rc[i].hitseq[0] == -1)
							stop = 1;
						else if (rc[i].hitseq[0] == 2)
							transmod = 1;
						else if (rc[i].hitseq[0] == 3)
							refmod = 1;						
						else if (rc[i].hitseq[0] == 4)
							free = 1;						
						else if (rc[i].hitseq[0] != 0)
						warn("wrong mode in refseq\n");

						++(rc[i].nhits);

						if(rc[i].nhits<nrefseq[i]){
							++(rc[i].hitseq);
							++(ireset[i]);
						}

						break;
				        }
				}
			}

			if (ifp!=NULL){
				fprintf(ifp,"^^^\n"
					"Interaction with interface %d at "
					"(x=%g,z=%g).\n",
					temp,rslast->x,rslast->z);

                         	if(refl){
				 fprintf(ifp," Ray is reflected ");
				} else if(refmod){
				 fprintf(ifp," Ray is reflected and "
				 "converted.\n");
			 	} else if(transmod){
				 fprintf(ifp," Ray is transmitted and " 	
				 "converted.\n");
				} else if(stop){
				 fprintf(ifp,"Hit stopping edge.  "
					     "Ray stopped.\n");
				} else if(free){
				 fprintf(ifp,"Free surface.  "
					     "Ray is reflected.\n");
				} else {
				 fprintf(ifp," Ray is transmitted ");
			        }
			}



			/* if ray at stopping edge, stop */
			if (stop) {
			 break;

			/* else if ray at reflecting edge, reflect */
			} else if (ea!=NULL && (refl)) {
				rs = reflectRay(rs,ifp,0,junkfp,reftrans);
				nameref=temp;

			/* else reflecting edge, mode conv */
			} else if (ea!=NULL && (refmod)) {
				rs = reflectRay(rs,ifp,1,junkfp,reftrans);
				nameref=temp;

			/* else transm edge, mode conv */
			} else if (ea!=NULL && (transmod)) {
				rs = RayAcrossEdge(rs,ifp,1,junkfp,reftrans);

			/* else free surface reflection */
			} else if (ea!=NULL && (free)) {
				rs = reflectRay(rs,ifp,2,junkfp,reftrans);

			/* else trace ray across edge */
			} else {
				rs = RayAcrossEdge(rs,ifp,0,junkfp,reftrans);
			}
			
		} while(rs!=NULL);
		

                /* reinitialize RayChecks for new ray */
                for(i=0;i<kk;i++){
		       rc[i].hitseq=rc[i].hitseq -ireset[i];
		       rc[i].nhits=0;

                }

		/* save ray end parameters 
		fa = rslast->f->fa;*/

		re[iangle].sigma = rslast->sigma;
		re[iangle].x = rslast->x;
		re[iangle].z = rslast->z;
		re[iangle].px = rslast->px;
		re[iangle].pz = rslast->pz;
		re[iangle].t = rslast->t;
		re[iangle].q1 = rslast->q1;
		re[iangle].p1 = rslast->p1;
		re[iangle].q2 = rslast->q2;
		re[iangle].p2 = rslast->p2;
		re[iangle].kmah = rslast->kmah;
		re[iangle].nref = rslast->nref;
		re[iangle].vgx = rslast->vgx;
		re[iangle].vgz = rslast->vgz;
		re[iangle].ab = angle;
                re[iangle].kend = temp;
                re[iangle].ampli = rslast->ampli;
                re[iangle].ampliphase = rslast->ampliphase;
                re[iangle].dangle = dangle;
                re[iangle].mode = rslast->mode;
                re[iangle].g11 = rslast->g11;
                re[iangle].g33 = rslast->g33;
                re[iangle].g13 = rslast->g13;
                re[iangle].num = iangle;
                re[iangle].nameref = nameref;
                re[iangle].rhob = rhob;
                re[iangle].rhoe = rslast->rho;

		
		/* optional output of rayendinformation */
                if (ifp!=NULL){
		 fprintf(ifp,"\n---------------------------"
				"------------------------ \n");
		 fprintf(ifp,"\n The following information is "
				"stored at the rayend \n");
		 fprintf(ifp,"\n Ray stops at "
			     "(%g,%g) at interface %i\n", 	
			     re[iangle].x,re[iangle].z,re[iangle].kend);
	         fprintf(ifp," Slowness components px= "
			    "%g pz=%g\n",re[iangle].px,re[iangle].pz);
		 fprintf(ifp," Traveltime = "
			     "%g  \t Sigma=%g\n",re[iangle].t,re[iangle].sigma);
		 fprintf(ifp," kmah index=%i \t Number of reflections=%i\n", 	
			     re[iangle].kmah,re[iangle].nref);
		 fprintf(ifp," g11=%g \t \t g33=%g \t g13=%g\n",
			     re[iangle].g11, re[iangle].g33,re[iangle].g13);
		 fprintf(ifp," Refl/Transm Amplitude =%g Phase=%g\n",
			    re[iangle].ampli,re[iangle].ampliphase); 
		 fprintf(ifp," ray number =%i reflected from=%i\n",
			    re[iangle].num,re[iangle].nameref); 
		 fprintf(ifp," density at begin =%g dens at end=%g\n",
			    re[iangle].rhob,re[iangle].rhoe); 
		 fprintf(ifp,"\n ------------------------" 	
			      "--------------------------- \n");
               } 
	}
}

    



static RayStep* RayInAnTri (RayStep *rs)
/*****************************************************************************
Trace rays in anisotropic triangles
******************************************************************************
Input:
* rs		pointer to last ray step
Output:
*rs		new ray step pointer
******************************************************************************
Notes:
Rays are traced from triangle to triangle. Intersection is found with
next triangle edge. Kinematic and dynamic parameters are updated.
******************************************************************************
Credits:  Andreas Rueger, Colorado School of Mines, 02/06/94

  The program is based on :
 	(gbray.c)traceRayInTri.c, AUTHOR: Andreas Rueger, 08/12/93
 	(sdray.c)traceRayInTri.c, AUTHOR: Dave Hale, CSM, 02/26/91
******************************************************************************/

{
	int kmah,nref,mode;
	float x,z,px,pz,t,q1,p1,q2,p2,ampli,
	      xa,za,xb,zb,dx,dz,b,c,
	      sigma,dt,dt1,small,frac,dsigma,
              ampliphase,rho;
	float vgx,vgz,g11,g33,g13;
	EdgeUse *eu,*eut,*eusmall;
	Face *f;
	
	/* get input parameters */
	sigma =	rs->sigma;
	x = rs->x;
	z = rs->z;
	px = rs->px;
	pz = rs->pz;
	t = rs->t;
	q1 = rs->q1;
	p1 = rs->p1;
	q2 = rs->q2;
	p2 = rs->p2;
	kmah = rs->kmah;
	nref = rs->nref;
	eu = rs->eu;
	f = rs->f;
	ampli = rs->ampli;
        ampliphase = rs->ampliphase;
	mode = rs->mode;
	vgx = rs->vgx;
	vgz = rs->vgz;
	g11 = rs->g11;
	g33 = rs->g33;
	g13 = rs->g13;
	rho = rs->rho;

	/* determine edge intersection with smallest positive dt */
	eusmall = NULL;
	small = FLT_MAX;
	eut = f->eu;
	do {
		/* edge endpoints */
		xa = eut->vu->v->y;  za = eut->vu->v->x;
		xb = eut->euCW->vu->v->y;  zb = eut->euCW->vu->v->x;
		
		/* coefficients b and c of equation to be solved for dt */
		dx = xb-xa;  dz = zb-za;
		b = dx*vgz - dz*vgx;
		c = (eut==eu ? 0.0 : dx*(z-za)-dz*(x-xa)); 
		
		if (b!=0.0)
			dt1 = -c/b;
		else
			dt1 = FLT_MAX;

		
		/* remember edge with smallest positive dt */
		if (0.0 <dt1 && dt1<small) {
			small = dt1;
			eusmall = eut;
		}
		
		/* next edge use */
		eut = eut->euCW;
		
	} while (eut!=f->eu);
	
	/* ray exits at edge with smallest positive dt */
	dt = small;
	eu = eusmall;
		
	/* change this if you include 2.5D spreading */
	dsigma=0;

	/* update ray parameters */
	sigma += dsigma;
	t +=dt;
	x += dt*vgx;
	z += dt*vgz;
	
	/* buggy place. Need to find a save fix. 
	   origionally, 0.0001 0.9999 */

	/* don't let ray exit too close to a vertex */
	xa = eu->vu->v->y;  dx = eu->euCW->vu->v->y-xa;
	za = eu->vu->v->x;  dz = eu->euCW->vu->v->x-za;
	frac = (ABS(dx)>ABS(dz)?(x-xa)/dx:(z-za)/dz);
	if (frac<0.001) {
		x = xa+0.001*dx;
		z = za+0.001*dz;
	} else if (frac>0.999) {
		x = xa+0.999*dx;
		z = za+0.999*dz;
	}

	/* return new raystep */
	rs->rsNext = (RayStep*)emalloc(sizeof(RayStep));
	rs = rs->rsNext;
	rs->sigma = sigma;
	rs->x = x;
	rs->z = z;
	rs->px = px;
	rs->pz = pz;
	rs->t = t;
	rs->q1 = q1;
	rs->p1 = p1;
	rs->q2 = q2;
	rs->p2 = p2;
	rs->kmah = kmah;
	rs->nref = nref;
	rs->eu = eu;
	rs->f = f;
	rs->rsNext = NULL;
        rs->vgx = vgx;
        rs->vgz = vgz;
        rs->ampli = ampli;
        rs->ampliphase = ampliphase;
	rs->mode = mode;
        rs->g11 = g11;
        rs->g13 = g13;
        rs->g33 = g33;
        rs->rho = rho;

	return rs;
}


static RayStep* RayAcrossEdge (RayStep *rs,FILE *ifp, int conv,FILE *junkfp,
	int reftrans)
/*****************************************************************************
rays across edges in anisotropic models
******************************************************************************
Input:
* rs		pointer to last ray step
conv		mode conversion=1
reftrans	=1 include transmission coeff

Output:
*rs		new ray step pointer
******************************************************************************
Notes:
Rays are traced across triangle edges. Kinematic and dynamic
boundary conditions are applied.If successful, return pointer 
to new RayStep.
If unsuccessful, return NULL.
Failure to trace across an edge is because:
(1) the ray was incident with angle greater than the critical angle, or
(2) the ray is incident at a boundary edge 
(3) grazing incidence
******************************************************************************
Credits:  Andreas Rueger, Colorado School of Mines, 02/06/94

  The program is based on :
 	(gbray.c)traceRayAcrossEdge.c, AUTHOR: Andreas Rueger, 08/12/93
 	(sdray.c)traceRayAcrossEdge.c, AUTHOR: Dave Hale, CSM, 02/26/91
******************************************************************************/
{
	int kmah,nref,modei,temp,mindexi,mindext,modet;
	float a1111i,a3333i,a1133i,a1313i,a1113i,a3313i;
	float a1111t,a3333t,a1133t,a1313t,a1113t,a3313t;
	float a1212i,a1223i,a2323i,a1212t,a1223t,a2323t;
	float rhoi,rhot,scale,gx,gz,frac,dx,dz,vgx,vgz;
	float vgxr,vgzr,pxr,pzr,pxi,pzi,coeff,phase;
	float sigma,x,z,px,pz,t,q1,p1,q2,p2,pl,rt,vgxi,vgzi;
        float ampli,ampliphase,g11,g13,g33,g11i,g13i,g33i,plold;
	EdgeUse *eu,*eum;
	EdgeUseAttributes *eua,*euma;
	Face *f;
	FaceAttributes *fa;

	/* get input parameters */
	sigma = rs->sigma;
	x = rs->x;
	z = rs->z;
	pxi = px = rs->px;
	pzi = pz = rs->pz;
	t = rs->t;
	q1 = rs->q1;
	p1 = rs->p1;
	q2 = rs->q2;
	p2 = rs->p2;
	kmah = rs->kmah;
	nref = rs->nref;
	eu = rs->eu;
	f = rs->f;
        ampli = rs->ampli;
	ampliphase = rs->ampliphase;
	vgxi = vgx = rs->vgx;
	vgzi = vgz = rs->vgz;
	modet = modei = rs->mode;
	g11 = g11i = rs->g11;
	g13 = g13i = rs->g13;
	g33 = g33i = rs->g33;

	/* check for boundary */
	if (eu->euMate->f==NULL) return NULL;
	
	/* determine stiffness on this side of edge */
	fa = f->fa;

	a1111i = fa->a1111;
	a3333i = fa->a3333;
	a1133i = fa->a1133;
	a1313i = fa->a1313;
	a1113i = fa->a1113;
	a3313i = fa->a3313;
        a1212i = fa->a1212;            
        a2323i = fa->a2323;             
        a1223i = fa->a1223;
        rhoi   = fa->rho;
	mindexi= fa->mindex;

	/* determine stiffness on other side of edge */
	eum = eu->euMate;
	f = eum->f;
	fa = f->fa;
	a1111t = fa->a1111;
	a3333t = fa->a3333;
	a1133t = fa->a1133;
	a1313t = fa->a1313;
	a1113t = fa->a1113;
	a3313t = fa->a3313;
        a1212t = fa->a1212;            
        a2323t = fa->a2323;             
        a1223t = fa->a1223;
        rhot   = fa->rho;
	mindext= fa->mindex;

	/* check if we need to apply boundary conditions */
	if(a1111i != a1111t || a3333i != a3333t ||
	   a1133i != a1133t || a1313i != a1313t ||
	   a1113i != a1113t || a3313i != a3313t ||
	   a1212i != a1212t || a1223i != a1223t ||
	   a2323i != a2323t || rhoi != rhot || conv)
        {
	    /* edge vector */
	    dx = eum->vu->v->y-eu->vu->v->y;
	    dz = eum->vu->v->x-eu->vu->v->x;

	    /* fractional distance along edge */
	    frac = (ABS(dx)>ABS(dz)?
		(x-eu->vu->v->y)/dx:
		(z-eu->vu->v->x)/dz);

	    /* linearly interpolate unit vector g tangent to edge */
	    eua = eu->eua;
	    euma = eum->eua;
	    if (eua!=NULL && euma!=NULL) {
		gx = frac*euma->tx-(1.0-frac)*eua->tx;
		gz = frac*euma->tz-(1.0-frac)*eua->tz;
	    } else {
		gx = -dx;
		gz = -dz;
	    }
	    scale = 1.0/sqrt(gx*gx+gz*gz);
	    gx *= scale;
	    gz *= scale;
	
	    /* gx = cos(horizontal/interface). The angle is
	       measured with respect to a right handed coordinate
	       system on the interface. gz = sin accordingly.
	       pl is measured relative to the same coordinate
	       system. Note that the sign differs for incidence
	       rays from both sides  */
	
	    /* tangent slowness component  */
	    plold = pl = gx*px+pz*gz;

	    /* optional output incidence values */
            if (ifp!=NULL){
		 fprintf(ifp,"\n incidence values :\n");
		 fprintf(ifp," mode=%i \t vgx=%g vgz=%g \n",
			     modei,vgxi,vgzi);
		 fprintf(ifp," g11=%g \t g13=%g g33=%g \n",
			     g11,g13,g33);
		 fprintf(ifp," ampli=%g \t ampliphase=%g \n",
			     ampli,ampliphase);
		 fprintf(ifp," tangential slowness component=%g \n",
			     pl);
	    }

	    /* find the ray mode in the new triangle */
 	    if(findnewmode(modei,&modet,conv,mindext) != 1)