refaniso.c

su 的源代码库

	K2=a1111*oa3333*p2 -oa3333;

	K3=p2-oa1313;

	sqr= K1*K1 -4.*K2*K3;
	if(ABS(sqr)<FLT_EPSILON)
		sqr=0.;

	if(sqr < 0){
		fprintf(stderr,"ERROR slowness medium 2 \n");
		return (-1);
	}

	/* vertical slowness in medium 2 */
	qasq=0.5*((K1 - sqrt(sqr)));

	if(ABS(qasq)<FLT_EPSILON)
		qasq=0.;

	if(qasq < 0){
		fprintf(stderr,"\n COMPLEX vertical P- slowness medium 2 \n");
		return (-1);
	}

	qa2=sqrt(qasq);
	qbsq=0.5*(K1 + sqrt(sqr));
	qb2=sqrt(qbsq);

	/* ddprint(qa2);ddprint(qb2); */

	sqr=a1111*p2 + a1313*qasq + a3333*qasq + a1313*p2 -2.;
	if (ABS(sqr) < FLT_EPSILON){
		fprintf(stderr,"ERROR P eigenvector medium 2 \n");
		return (-1);
	}

	sqr=1./sqr;

	la2=a3333*qasq + a1313*p2 -1;
	
	la2=la2*sqr;
	if(ABS(la2)<FLT_EPSILON)
		la2=0.;

	if(la2 <0){
		fprintf(stderr,"ERROR la2<0\n");
		return (-1);
	}

	/* P eigenvector component medium 2*/
	la2=sqrt(la2);

	ma2=a1313*qasq + a1111*p2 -1;
	ma2=ma2*sqr;

	if(ABS(ma2)<FLT_EPSILON)
		ma2=0.;

	if(ma2 <0){
		fprintf(stderr,"ERROR ma2<0\n");
		return (-1);
	}

	/* P eigenvector component medium 2*/
	ma2=sqrt(ma2);

	sqr=a1111*p2 + a1313*qbsq + a3333*qbsq + a1313*p2 -2.;
	if (ABS(sqr) < FLT_EPSILON){
		fprintf(stderr,"ERROR SV eigenvector medium 2 \n");
		return (-1);
	}

	sqr=1./sqr;

	lb2=a1313*qbsq + a1111*p2 -1.;
	lb2=lb2*sqr;
	if(ABS(lb2)<FLT_EPSILON)
		lb2=0.;

	if(lb2 < 0){
		fprintf(stderr,"ERROR lb2<0\n");
		return (-1);
	}

	/* SV-eigenvector component medium 2*/
	lb2=sqrt(lb2);

	mb2=a3333*qbsq + a1313*p2 -1.;
	mb2=mb2*sqr;
	if(ABS(mb2)<FLT_EPSILON)
		mb2=0.;

	if(mb2 <0){
		fprintf(stderr,"ERROR mb2<0\n");
		return (-1);
	}

	/* SV-eigenvector component medium 1*/
	mb2=sqrt(mb2);

	/* ddprint(la2);ddprint(lb2);
	   ddprint(ma2);ddprint(mb2);*/

	/* convert into stiffness */
	a1313=a1313*rho2;
	a3333=a3333*rho2;
	a1133=a1133*rho2;

	a2=a1313*(qa2*la2 + pl*ma2);
	b2=a1313*(qb2*mb2 - pl*lb2);
	c2=pl*la2*a1133 + qa2*ma2*a3333;
	d2=pl*mb2*a1133 - qb2*lb2*a3333;

	/* ddprint(a2);ddprint(b2);
	   ddprint(c2);ddprint(d2); */

	m11=la1; m12=mb1; m13=-la2; m14=-mb2;
	m21=c1;  m22=d1;  m23=-c2;  m24=-d2;
	m31=ma1; m32=-lb1;m33=ma2;  m34=-lb2;
	m41=a1;  m42=b1;  m43=a2;   m44=b2; 

	/*ddprint(m11); ddprint(m12); ddprint(m13); ddprint(m14);
	ddprint(m21); ddprint(m22); ddprint(m23); ddprint(m24);
	ddprint(m31); ddprint(m32); ddprint(m33); ddprint(m34);
	ddprint(m41); ddprint(m42); ddprint(m43); ddprint(m44);*/

	sqr=(m14*m23*m32*m41 - m13*m24*m32*m41 - m14*m22*m33*m41 + 
      	     m12*m24*m33*m41 + m13*m22*m34*m41 - m12*m23*m34*m41 - 
     	     m14*m23*m31*m42 + m13*m24*m31*m42 + m14*m21*m33*m42 - 
     	     m11*m24*m33*m42 - m13*m21*m34*m42 + m11*m23*m34*m42 + 
     	     m14*m22*m31*m43 - m12*m24*m31*m43 - m14*m21*m32*m43 + 
      	     m11*m24*m32*m43 + m12*m21*m34*m43 - m11*
      	     m22*m34*m43 -m13*m22*m31*m44 + m12*m23*m31*m44 + 
	     m13*m21*m32*m44 - m11*m23*m32*m44 - m12*m21*m33*m44 + 	
	     m11*m22*m33*m44);
	/* ddprint(sqr); */

	if(ABS(sqr)<FLT_EPSILON){
		fprintf(stderr,"ERROR Denominator == 0 \n");
		return (-1);
	}

	*coeff=1./sqr;


	/* PP transmission*/
	if(modei==0 && modet==0 && rort==0)
		*coeff *=2*(-(m14*m21*m32*m41)+m11*m24*m32*m41+
			m12*m21*m34*m41-m11*m22*m34*m41+m14*m21*m31*m42-
			m11*m24*m31*m42-m12*m21*m31*m44+m11*m22*m31*m44);

	/* PP reflection*/
	else if(modei==0 && modet==0 && rort==1)
		*coeff *=(m14*m23*m32-m13*m24*m32-m14*m22*m33+m12*m24*m33+
       			  m13*m22*m34-m12*m23*m34)*m41+m31*(-(m14*m23*m42)+
		   	  m13*m24*m42+m14*m22*m43-m12*m24*m43-m13*m22*m44+
 			  m12*m23*m44)+m21*(-(m14*m33*m42)+m13*m34*m42+
			  m14*m32*m43-m12*m34*m43-m13*m32*m44+m12*m33*m44)+
    			  m11*(m24*m33*m42-m23*m34*m42-m24*m32*m43+
			  m22*m34*m43+m23*m32*m44-m22*m33*m44);

	/* PS transmission*/
	else if(modei==0 && modet==1 && rort==0)
		*coeff *=2*(m13*m21*m32*m41-m11*m23*m32*m41-m12*m21*m33*m41+
      			 m11*m22*m33*m41-m13*m21*m31*m42+m11*m23*m31*m42+
			 m12*m21*m31*m43-m11*m22*m31*m43);

	/* PS reflection*/
	else if(modei==0 && modet==1 && rort==1)
		*coeff *=2*(m14*m21*m33*m41-m11*m24*m33*m41-m13*m21*m34*m41+
		         m11*m23*m34*m41-m14*m21*m31*m43+m11*m24*m31*m43+
			 m13*m21*m31*m44-m11*m23*m31*m44);

	/* SP transmission*/
	else if(modei==1 && modet==0 && rort==0)
		*coeff *=2*(-(m14*m22*m32*m41)+m12*m24*m32*m41+
			 m14*m22*m31*m42-m12*m24*m31*m42+m12*m21*m34*m42-
			 m11*m22*m34*m42-m12*m21*m32*m44+m11*m22*m32*m44);

	/* SP reflection*/
	else if(modei==1 && modet==0 && rort==1)
		*coeff *=2*(-(m14*m22*m33*m42)+m12*m24*m33*m42+
			 m13*m22*m34*m42-m12*m23*m34*m42+m14*m22*m32*m43-
			 m12*m24*m32*m43-m13*m22*m32*m44+m12*m23*m32*m44);

	/* SS transmission*/
	else if(modei==1 && modet==1 && rort==0)
		*coeff *=2*(m13*m22*m32*m41-m12*m23*m32*m41-m13*m22*m31*m42+		
			 m12*m23*m31*m42-m12*m21*m33*m42+m11*m22*m33*m42+
			 m12*m21*m32*m43-m11*m22*m32*m43);

	/* SS reflection*/
	else if(modei==1 && modet==1 && rort==1)
		*coeff *=(-(m14*m23*m31)+m13*m24*m31+m14*m21*m33-
		         m11*m24*m33-m13*m21*m34+m11*m23*m34)*m42+
			 m32*(m14*m23*m41-m13*m24*m41-m14*m21*m43+		
			 m11*m24*m43+m13*m21*m44-m11*m23*m44)+
			 m22*(m14*m33*m41-m13*m34*m41-m14*m31*m43+
			 m11*m34*m43+m13*m31*m44-m11*m33*m44)+	
			 m12*(-(m24*m33*m41)+m23*m34*m41+m24*m31*m43-		
			 m21*m34*m43 - m23*m31*m44 + m21*m33*m44);
	else 
		return (-1);

	sqr=*coeff;

	/* ddprint(sqr); */
	return 1;
}


void gvelpolSH(double a1212, double a2323, double a1223, double px, 
		double pz, double *vgx, double *vgz, double *g11, 
		double *g13, double *g33)
/*****************************************************************************
gvelpolSH - compute SH group velocity and polarization for TI medium
******************************************************************************
Input:
aijkl		stiffness elements
px,pz		slowness elements

Output:
*vgx, *vgz	group velocities
*g11,*g13,*g33  polarizations

******************************************************************************
Author:  Andreas Rueger, Colorado School of Mines, 01/09/95
******************************************************************************/
{

	*vgx = a1212*px + a1223*pz;
	*vgz = a1223*px + a2323*pz;

	*g11=0.;
	*g13=0.;
	*g33=0.;
}

int gvelpolTI (double a1111, double a3333, double a1133, double a1313,
	double a1113, double a3313, double px, double pz, double *vgx,
	double *vgz, double *g11n, double *g13n, double *g33n)
/*****************************************************************************
gvelpolTI - compute P/SV group velocity and polarization for TI medium
******************************************************************************
Input:
aijkl		stiffness elements
px,pz		slowness elements

Output:
*vgx, *vgz	group velocities
*g11,*g13,*g33  polarizations

******************************************************************************
Notes:
Input px/pz values determine ray mode.
Subroutine returns:
(-1) if error in g11, g33 detected,
(1) otherwise
******************************************************************************
Author:  Andreas Rueger, Colorado School of Mines, 01/09/95
******************************************************************************/
{
	double gamm11,gamm13,gamm33;
	double px2,pz2,pxz,den,g11,g13,g33;

	px2   = px*px;
	pz2   = pz*pz;
	pxz   = px*pz;


	/*anisotropy parameters*/
	gamm11 = a1111*px2+ a1313*pz2 +2*a1113*pxz;
	gamm33 = a3333*pz2 + a1313*px2+2*a3313*pxz;
	gamm13 = (a1133+a1313)*pxz+ a1113*px2+ a3313*pz2;
	den     = 1/(gamm11+gamm33-2);
	g11     = (gamm33-1)*den;
	g33     = (gamm11-1)*den;
	g13     = -gamm13*den;

	/* computing ray (group) velocities */
	*vgx =  (a1111*px*g11+(a1133+a1313)*pz*g13+a3313*pz*g33+
		a1113*(pz*g11+2*px*g13)+a1313*g33*px);
	*vgz =  (a3333*pz*g33+(a1133+a1313)*px*g13+a1113*px*g11+
		+a3313*(px*g33+2*pz*g13)+a1313*g11*pz);

	/* kill round-off */
	if( g11 < - FLT_EPSILON || g33 < - FLT_EPSILON) 
		return (-1);

	else if( g11 < 0.0 )
		g11 = 0;

	else if( g33 < 0.0 )
		g33 = 0;


	*g11n=g11;
	*g13n=g13;
	*g33n=g33;

	return (1);

}


int p_hor2DTI (Stiff2D *spar, double s, int mode, double *p)
/*****************************************************************************
p_hor2DTI - compute horizontal slowness for in-plane propagation in TI medium
******************************************************************************
Input:
spar		stiffness elements
s		sin(incidence angle)
mode		0=qP-Wave, 1=qSV-wave

Output:
*p	horizontal slowness component p_x

******************************************************************************
Notes:
(-1) if error detected,
(1) otherwise
******************************************************************************
Author:  Andreas Rueger, Colorado School of Mines, 01/26/95
******************************************************************************/
{

        double cc,c,ss,sc,gamm11,gamm33,gamm13,sqr,vp2;

	if(mode != 0 && mode !=1)
		return (-1);	/* wrong mode */

	cc = 1.-s*s;

	/* ddprint(s); */

	if(ABS(cc)<FLT_EPSILON)
		cc=0.;
	if(cc <0)
		return (-1);	/* evanescent energy generated */


 	c  = sqrt(cc);
	ss = s*s;
	sc = s*c;

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

	/* error check */
	if(vp2 < 0) 
		err("\n ERROR: ray initialization \n");
	/* p_hor = sin(pangle)/ABS(phasevel) */

	*p     = s/sqrt(vp2*.5);

	return (1);
	
}

int p_vert2DVTI(Stiff2D *spar1, double pl, int modei, double *p_vert)
/*****************************************************************************
p_vert2DVTI - Given the horizontal slowness, compute vertical slowness component
******************************************************************************
Input:
spar1            (density normalized) stiffnesses
pl               horizontal slowness
modei            mode (0=P 1=SV)

Output:
p_vert           vertical slowness component

******************************************************************************