elaray.c

su 的源代码库

		err("\n Error in findnewmode/transmission .\n");

	    /* isotropic case P-wave */
 	    if(modet == 0){

	      if(findPiso(a3333t,pl,gx,gz,&px,&pz,&vgx,&vgz,&g11,&g13,&g33,0)
			  !=1)
		return NULL;

	    /* isotropic case SV/SH-wave */
	    } else if(modet == 1 || modet == 2){

	      if(findSiso(a1313t,pl,gx,gz,&px,&pz,&vgx,&vgz,&g11,&g13,&g33,0)
			!=1)
			return NULL;
	      else if( modet == 2) {
		g11=g13=g33=0;
	      }

	    /* anisotropic case qP/qSV */
	    } else if(modet == 3 || modet == 4){
		rt=1;

		if(findqPqSV(gz,gx,pl,a1111t,a3333t,a1133t,a1313t,a1113t,
			a3313t,modet,&px,&pz,&vgx,&vgz,&g11,&g13,&g33,rt,ifp)
			!=1)
			return NULL;

	    /* anisotropic case SH */
	    } else if(modet == 5){

		if(findqSH(gz,gx,pl,a1212t,a1223t,a2323t,&px,&pz,&vgx,
		&vgz,0) != 1)
			return NULL;
		else {
			g11=g13=g33=0;
		}


	    } else
		err(" ERROR. Wrong mode in RayAcrossEdge");

	    /* This is the final check */
	    pl = gz*pz + px*gx;
	    if((pl-plold)*(pl-plold) > 0.000000001)
		err(" ERROR; Snell's law not satisfied in transmission \n");

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

	    /* check if we need to compute transmission coeff */
            if(reftrans !=0 )
            {
		
	    /* initialize */
	    coeff = 1;
	    phase = 0;

	    /* real/complex isotropic ref/transm coefficient */
	    if((modei==0 || modei==1 || modei==2) &&
	       (modet==0 || modet==1 || modet==2) ){

		temp = rt_iso_real(sqrt(a1111i),sqrt(a1111t),sqrt(a1313i),
		sqrt(a1313t),rhoi,rhot,pl,modei,modet,0,&coeff);

		/* complex coeff; needs to be tested */
		if(temp==0) {
			if(rt_iso_cmplx(sqrt(a1111i),sqrt(a1313t),
				sqrt(a1111t), sqrt(a1313t),rhoi,rhot,pl,
				modei,modet,0,&coeff,&phase) ==-1){

				warn("problems in rt_iso_cmplx/trans");
				return NULL;
			}

		/* problems with rt_iso_real/transm */
		}else if (temp==-1){

			warn(" Problems in rt_iso_real/transm");
			return NULL;
		}

	    /* real/complex coeff for SH-ani propagation */
	    } else if(modei==5 && modet==5){

		/* get the reflected wave root */
		if(findqSH(gz,gx,plold,a1212i,a1223i,a2323i,&pxr,&pzr,&vgxr,
		&vgzr,1) != 1){

			warn(" No real SH-reflection ");
			return NULL;

		}

                if(rt_SHa_real(a1212i*rhoi,a2323i*rhoi,a1223i*rhoi,
			a1212t*rhot,a2323t*rhot,a1223t*rhot,pxi,pzi,
			px,pz,pxr,pzr,0,&coeff,gz,gx) !=1 ){

				warn("problems in SHa_real");
				return NULL;
			}


	    /* real/complex anisotropic coeff */
	    } else {

		if(rt_ani_real(gz,gx,a1111i,a3333i,a1133i,a1313i,a1113i, 	
			a3313i,rhoi,a1111t,a3333t,a1133t,a1313t,a1113t,a3313t,
		        rhot,pxi,pzi,g11i,g13i,g33i,px,pz,g11,g13,g33, 	
	                modei,modet,0,&coeff,ifp) !=1 ){

				warn("problems in rt_ani_real");
				return NULL;
		}

	   
            }

	    ampli *= coeff;
	    ampliphase += phase;

	    /* optional output transmitted values */
            if (ifp!=NULL){
		 fprintf(ifp," ampli=%g \t ampliphase=%g \n",
			     ampli,ampliphase);
	    }

	    /* fprintf(junkfp,"%f %f \n", 
		   (atan(pxi/pzi)+asin(gz))*180/PI,coeff); */

	    if(ifp!=NULL)
		fprintf(ifp," transmission coeff=%g \n",coeff);


	    } /* close if reftrans */

	} /* close if boundary conditions */


	/* 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 = eum;
	rs->f = f;
	rs->vgx = vgx;
	rs->vgz = vgz;
	rs->rsNext = NULL; 
        rs->ampli = ampli;
        rs->ampliphase = ampliphase;
	rs->mode = modet;
        rs->g11 = g11;
        rs->g13 = g13;
        rs->g33 = g33;
        rs->rho = rhot;


	return rs;
}

static RayStep* reflectRay (RayStep *rs,FILE *ifp, int conv,FILE *junkfp,
	int reftrans)
/* Reflect ray from edge and return pointer to new RayStep. */
{
	int kmah,nref,modei,modet,temp,mindexi,mindext,rort;
	float x,z,px,pz,t,q1,p1,q2,p2,rt,
	      scale,gx,gz,frac,dx,dz,plold;
	float g11,g13,g33,g11i,g13i,g33i;
	float rhoi,rhot,coeff,phase;
	float pxi,pzi,vgxi,vgzi,vgxt,vgzt,pxt,pzt ;
        float ampli,ampliphase,sigma,pl,vgx,vgz;
	float a1111i,a3333i,a1133i,a1313i,a1113i,a3313i;
	float a1111t,a3333t,a1133t,a1313t,a1113t,a3313t;
	float a1212i,a1223i,a2323i,a1212t,a1223t,a2323t;

	EdgeUse *eu,*eum;
	EdgeUseAttributes *eua,*euma;
	Face *f,*fn;
	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;
	modei = rs->mode;
	g11 = g11i = rs->g11;
	g13 = g13i = rs->g13;
	g33 = g33i = rs->g33;


	/* reflection or free surface */
	/* check for boundary */
	if (eu->euMate->f==NULL) 
		rort=2;
	else 
		rort=1;

	/* determine stiffness on this side of edge */
	fa = f->fa;

	a1111t=a1111i= fa->a1111;
	a3333t=a3333i= fa->a3333;
	a1133t=a1133i= fa->a1133;
	a1313t=a1313i= fa->a1313;
	a1113t=a1113i= fa->a1113;
	a3313t=a3313i= fa->a3313;
        a1212t=a1212i = fa->a1212;            
        a2323t=a2323i = fa->a2323;             
        a1223t=a1223i = fa->a1223;
        rhot = rhoi   = fa->rho;
	mindext = mindexi= fa->mindex;


	if(rort != 2){
		/* determine stiffness on other side of edge */
		eum = eu->euMate;
		fn = eum->f;
		fa = fn->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;

		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;
		gx = frac*euma->tx-(1.0-frac)*eua->tx;
		gz = frac*euma->tz-(1.0-frac)*eua->tz;
		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  */
		pl = plold = 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);
		}
	} else {
		gx = 1.0;
		gz=0.0;
		pl = plold = px;
	}	

	/* find the ray mode in the new triangle */
 	if(findnewmode(modei,&modet,conv,mindexi) != 1)
		err("\n Error in findnewmode/transmission .\n");

	/* isotropic case P-wave */
 	if(modet == 0){

	if(findPiso(a3333i,pl,gx,gz,&px,&pz,&vgx,&vgz,&g11,&g13,&g33,1)
			  !=1)
		return NULL;

	/* isotropic case SV/SH-wave */
	} else if(modet == 1 || modet == 2){

		if(findSiso(a1313i,pl,gx,gz,&px,&pz,&vgx,&vgz,&g11,&g13,&g33,1)
			!=1)
			return NULL;
		else if( modet == 2) {
			g11=g13=g33=0;
		}

	/* anisotropic case qP/qSV */
	} else if(modet == 3 || modet == 4){
		rt = -1;

		if(findqPqSV(gz,gx,pl,a1111i,a3333i,a1133i,a1313i,a1113i,
			a3313i,modet,&px,&pz,&vgx,&vgz,&g11,&g13,&g33,rt,ifp)
			!=1)
			return NULL;

	/* anisotropic case SH */
	} else if(modet == 5){

		if(findqSH(gz,gx,pl,a1212i,a1223i,a2323i,&px,&pz,&vgx,
		&vgz,1) != 1)
			return NULL;
		else {
			g11=g13=g33=0;
		}


	} else
		err(" ERROR. Wrong mode in ReflectRay");

	/* This is the final check */
	pl = gz*pz + px*gx;
	if((pl-plold)*(pl-plold) > 0.000000001)
		err(" ERROR; Snell's law not satisfied in reflection \n");

	/* reflection coefficients */
	if(reftrans != 0){

	    /* initialize */
	    coeff = 1;
	    phase = 0;

	    /* real/complex isotropic ref/transm coefficient */
	    if((modei==0 || modei==1 || modei==2) &&
	       (modet==0 || modet==1 || modet==2) ){

		temp = rt_iso_real(sqrt(a1111i),sqrt(a1111t),sqrt(a1313i),
		sqrt(a1313t),rhoi,rhot,pl,modei,modet,rort,&coeff);

		/* real isotropic coeff */
		if(temp==1)
			ampli*=coeff;

		/* complex iso coeff NOT FINISHED YET*/
		else if(temp==0) {

			if(rt_iso_cmplx(sqrt(a1111i),sqrt(a1313t),
				sqrt(a1111t), sqrt(a1313t),rhoi,rhot,pl,
				modei,modet,rort,&coeff,&phase) ==-1){

				       warn("problems in rt_iso_cmplx/trans");
				       return NULL;
			}

		/* problems with rt_iso_real/transm */
		} else if (temp==-1){
			warn(" Problems in rt_iso_real/transm");
			return NULL;
		}
	
	    /* real/complex coeff for SH-ani propagation */
	    } else if(modei==5 && modet==5 && rort !=2 ){

		/* get the transmitted wave root */
		if(findqSH(gz,gx,plold,a1212t,a1223t,a2323t,&pxt,&pzt,&vgxt,
		&vgzt,0) != 1){
			warn(" No real SH-transmission ");
			return NULL;
		}

		/* get SH reflection coeff */
                if(rt_SHa_real(a1212i*rhoi,a2323i*rhoi,a1223i*rhoi,
			a1212t*rhot,a2323t*rhot,a1223t*rhot,pxi,pzi,
			pxt,pzt,px,pz,rort,&coeff,gz,gx)!= 1){
				warn(" problems in rt_SHa_real ");
				return NULL;
		}

	    /* get free surface SH reflection coeff */
	    } else if(modei==5 && modet==5 && rort == 2 ){
                if(rt_SHa_real(a1212i*rhoi,a2323i*rhoi,a1223i*rhoi,
			0,0,0,pxi,pzi,
			0,0,px,pz,rort,&coeff,gz,gx)!= 1){
				warn(" problems in rt_SHa_real ");
				return NULL;
		}

	    /* real/complex anisotropic coeff */
	    } else {

		temp=rt_ani_real(gz,gx,a1111i,a3333i,a1133i,a1313i,a1113i, 	
			a3313i,rhoi,a1111t,a3333t,a1133t,a1313t,a1113t,a3313t,
		        rhot,pxi,pzi,g11i,g13i,g33i,px,pz,g11,g13,g33, 	
	                modei,modet,1,&coeff,ifp);
		if(temp!=1) 
			return NULL;

	    }

	    ampli *= coeff;
	    ampliphase += phase;

	    /* fprintf(junkfp,"%f %f\n", 
		   (atan(pxi/pzi)+asin(gz))*180/PI,coeff); */

	}

	/* optional output reflected values */
        if (ifp!=NULL){
		 fprintf(ifp,"\n reflected values :\n");
		 fprintf(ifp," mode=%i \t vgx=%g vgz=%g \n",
			     modet,vgx,vgz);
		 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);
		 fprintf(ifp," reflection coeff=%g\n",coeff);

	}


	/* 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+1;
	rs->eu = eu;
	rs->f = f;
	rs->rsNext = NULL;
        rs->ampli = ampli;
        rs->ampliphase = ampliphase;
        rs->vgx = vgx;
        rs->vgz = vgz;
	rs->mode =modet;
	rs->g11 =g11;
	rs->g33 =g33;
	rs->g13 =g13;
        rs->rho = rhoi;

	return rs;
}


int testIfSourceOnEdge(Face *tris, float *z, float *x)
/*****************************************************************************
 check if source is on an edge or on a vertex. If so, move the source
 not the most elegant program, feel free to change it
******************************************************************************
Input:
*xs		source coordinate
*zs		source coordinate
*tris		face containing source coordinate
******************************************************************************
Output:
0		source is not on edge/vertex
1		source is on edge/vertex
This is a test version
******************************************************************************
Author:   Andreas Rueger, Colorado School of Mines, 02/01/94
******************************************************************************/
{
    float x1,x2,x3,z1,z2,z3,eps,zmax,xs,zs;
    float xmax;
    float m12,m13,m23,b12,b13,b23;
    EdgeUse *eu;

    eu=tris->eu;
    eps=tris->m->eps;
    zmax=tris->m->xmax;
    xmax=tris->m->ymax;
    xs=*x;
    zs=*z;
    m12=m13=m23=0.0;


    /* get vertices */
    x1=eu->vu->v->y;
    z1=eu->vu->v->x;
    x2=eu->euCW->vu->v->y;
    z2=eu->euCW->vu->v->x;
    x3=eu->euCCW->vu->v->y;
    z3=eu->euCCW->vu->v->x;

    /* source is sitting on vertex */
    if( (xs-x1)*(xs-x1) + (zs-z1)*(zs-z1) < eps*eps ||
	(xs-x2)*(xs-x2) + (zs-z2)*(zs-z2) < eps*eps ||
	(xs-x3)*(xs-x3) + (zs-z3)*(zs-z3) < eps*eps ){
             *z = (zs < zmax-eps ? zs + 2*eps : zs - 2*eps);
             *x = (xs < xmax-eps ? xs + 2*eps : xs - 2*eps);
	     return 1;
    /* check the most typical cases */
    } else if((x3-x2)*(x3-x2) < eps*eps && (x3-xs)*(x3-xs) < eps*eps){
             *x = (xs < xmax-eps ? xs + eps : xs - eps);
	     return 1;