refaniso.c

su 的源代码库

Notes:
	routine returns (-1) if evanescent energy present
	it is assumed that a1111,a3333,a1313 >0 
	Currently no SH scattering supported

******************************************************************************
Technical reference: Graebner, M.; Geophysics, Vol 57, No 11:
		Plane-wave reflection and transmission coefficients
		for a transversely isotropic solid.

******************************************************************************
Author: CWP: Andreas Rueger, Colorado School of Mines, 01/26/96
******************************************************************************/
{
	double K1,K2,K3;
	double qa1,qb1,sqr,p2;
	double qasq,qbsq;
	double a1111,a3333,a1133;
	double a1313,oa3333,oa1313;

	p2=pl*pl;

	/* note: we need only density normalized
	coefficients to get vertical slowness and
	eigenvectors */


	/**********************  Medium 1  ***********************************/
	a1111=spar1->a1111;
	a3333=spar1->a3333;
	a1133=spar1->a1133;
 	a1313=spar1->a1313;

	oa3333=1./a3333;
	oa1313=1./a1313;

	sqr=a1111*oa1313+a1313*oa3333-(a1313+a1133)*(a1313+a1133)*oa3333*oa1313;
	K1=oa3333+oa1313 - sqr*p2;

	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," evanescent energy generated  \n");
		return (-1);
	}

	/* vertical slowness in medium 1 */
	qasq=0.5*((K1 - sqrt(sqr)));
        if(qasq < 0){
		fprintf(stderr,"evanescent energy generated  \n");
		return (-1);
	}

	qa1=sqrt(qasq);
	qbsq=0.5*(K1 + sqrt(sqr));
	qb1=sqrt(qbsq);

	if(modei==0)
	   *p_vert =qa1;
	else
           *p_vert =qb1;
	
	/* ddprint(sqr); */
	return 1;
}

void rottens2D (Stiff2D *spar, double phi)
/*****************************************************************************
rottens2D - rotate 2-D elastic stiffness tensor
******************************************************************************
Input:
aijkl		input/output stiffness elements
phi 		rotation angle (counterclock wise)

******************************************************************************
Notes:
Warning!: double precision required for this routine!

******************************************************************************
Author:  Andreas Rueger, Colorado School of Mines, Jan 07, 1995
******************************************************************************/
{
		double ss,cc,sc,a,c,f,l,m,n,o,p,q;
		double si=sin(phi),co=cos(phi);

		ss=si*si; cc=co*co;
		sc=si*co;
		a=spar->a1111;
		c=spar->a3333;
		f=spar->a1133;
		l=spar->a1313;
		m=spar->a1113;
		n=spar->a3313;
		o=spar->a1212;
		p=spar->a2323;
		q=spar->a1223;


		spar->a1111=a*cc*cc+c*ss*ss+2*f*cc*ss+4*l*cc*ss
			-4*m*cc*sc-4*n*ss*sc;
		spar->a3333=a*ss*ss+c*cc*cc+2*f*cc*ss+4*l*cc*ss+
			4*m*ss*sc+4*n*cc*sc;
		spar->a1133=a*ss*cc+c*ss*cc+f*(ss*ss+cc*cc)-4*l*ss*cc-2*m*
			(ss*sc-cc*sc)-2*n*(cc*sc- ss*sc);
		spar->a1313=a*ss*cc+c*ss*cc-2*f*ss*cc+l*(cc*cc+ss*ss-2*ss*cc)
			-2*m*(ss*sc-cc*sc)-2*n*(cc*sc- ss*sc);
		spar->a1113=a*cc*sc-c*ss*sc+f*(ss*sc-cc*sc)+2*l*(ss*sc-cc*sc)+
			m*(cc*cc -3*cc*ss) +(3*cc*ss - ss*ss )*n;
		spar->a3313=a*ss*sc-c*cc*sc+f*(cc*sc-ss*sc)+2*l*(cc*sc-ss*sc)+
			m*( 3*cc*ss - ss*ss ) + n*( cc*cc - 3*ss*cc);
		spar->a1212=cc*o-2*sc*q+ss*p;
		spar->a2323=ss*o+2*sc*q+cc*p;
		spar->a1223=sc*o+(cc-ss)*q-sc*p;
}

int thom2stiffTI (double vp, double vs, double eps, double delta, double gamma,
	double phi, Stiff2D *spar, int sign)
/*****************************************************************************
thom2stiffTI - convert Thomsen's parameters into density normalized 
		stiffness components for transversely isotropic material
		with in-plane-tilted axes of symmetry
******************************************************************************
Input:
vp,vs	P and S wave velocity along symmetry axes
eps	Thomsen's parameter with respect to symmetry axes
delta	Thomsen's parameter with respect to symmetry axes
gamma	Thomsen's parameter with respect to symmetry axes
phi	angle(rad) vertical --> symmetry axes (clockwise) 
sign	sign of c11+c44 ( for most materials sign=1)

Output:
*aijkl	density normalized stiffness components
*******************************************************************************
Notes:
subroutine returns (-1) if unphysical parameter detected, otherwise
returns (1).

(1)-axes horizontal (2)-axes out-of-plane (3)-axis vertical

For vertical axes of symmetry see thom2stiffVTI.c
*******************************************************************************
Author: Andreas Rueger, Colorado School of Mines, Jan 07, 1995
******************************************************************************/
{
	double temp1,temp2;

	if(vp<=0. || vs<0. )
		return (-1);
	
	spar->a3333 = temp1 = vp*vp;
	spar->a1313 = temp2 = vs*vs;
	spar->a1111 = 2.* temp1*eps + temp1;

	temp1 = 2.* delta *temp1 *(temp1-temp2)
	        +(temp1-temp2)*(temp1-temp2);

	if(temp1 < 0.)
		return (-1);
	 
	if (sign == 1)
		spar->a1133 = sqrt(temp1)-temp2;
	else if(sign == -1)
		spar->a1133 = - sqrt(temp1)-temp2;
	else
		return (-1);

	spar->a1212 = 2.*temp2*gamma + temp2;

	spar->a1113 = spar->a3313 = spar->a1223 = 0.0;

	spar->a2323 = temp2; /* a1313=a2323 */

	rottens2D(spar,phi);


	return (1);

}



int stiff2thomVTI (double a1111, double a3333, double a1133, double a1313, 	
	double a1212, double *vp, double *vs, double *eps,
        double *delta, double *gamma)
/*****************************************************************************
stiff2thomVTI - convert density normalized stiffness components 
		parameters into Thomsen's for transversely isotropic 
		material with vertical axes of symmetry
******************************************************************************
Input:
*aijkl	density normalized stiffness components

Output:
vp,vs	vertical P and S wave velocity
eps	Thomsen's parameter 
delta	Thomsen's parameter 
gamma	Thomsen's parameter 

*******************************************************************************
Notes:
subroutine returns (-1) if unphysical parameter detected, otherwise
returns (1).

(1)-axes horizontal (2)-axes out-of-plane (3) axis vertical

For tilted in-plane axes of symmetry see stiff2thomTI.c
or libTest/stiff2thoma.ma
*******************************************************************************
Author: Andreas Rueger, Colorado School of Mines, Jan 08, 1995
******************************************************************************/
{
	double temp;

	if(a1111<=0. || a1313<0.)
		return (-1);
	 
	
	*vp  = sqrt(a3333);
	*vs  = sqrt(a1313);
	*eps = (a1111-a3333)*0.5/a3333;

	temp   = (a1133+a1313)*(a1133+a1313)-(a3333-a1313)*(a3333-a1313);
	*delta = temp/(a3333*(a3333-a1313)) * 0.5;
	*gamma = (a1212-a1313)/a1313*0.5;

	return (1);

}

int stiff2tv(Stiff2D *spar,double *alpha,double *beta,double *ev,
double *dv,double *gv)
/*****************************************************************************
stiff2tv - Convert stiffnesses into Thomsen's parameter of the equivalent
VTI model, P-wave reflections and azimuthal AVO in TI-media
******************************************************************************
Input:
spar            stiffnesses (normalized or non-normal.)

Output:
alpha           fracture plane compressional velocity
beta            S_parallel vertical velocity
ev              eps of equiv. VTI
dv              delta of ..
gv              gamma of ..

******************************************************************************
Author: CWP: Andreas Rueger, Colorado School of Mines, 01/27/96
******************************************************************************/
{



   	   *alpha=sqrt(spar->a3333);
	   *beta=sqrt(spar->a2323);
	   *ev=(spar->a1111 - spar-> a3333)/(2.*spar->a3333);
	   *dv=((spar->a1133+spar->a1313)*(spar->a1133+spar->a1313)-
		   (spar->a3333-spar->a1313)*(spar->a3333-spar->a1313))/
		  (2.*spar->a3333*(spar->a3333-spar->a1313));
	   *gv=(spar->a1212 - spar->a2323)/(2.*spar->a2323);
	   
	   return 1;
	   
	   
}

int thom2tv(double vp,double vs,double eps,double delta,double gamma,
	    double *alpha,double *beta,double *ev,double *dv,double *gv)
/*****************************************************************************
thom2tv - Convert generic Thomsen parameters into Thomsen's parameter 
          of the equivalent VTI model
******************************************************************************
Input:

vp              symm. axis compressional velocity
vs              symm. axis shear velocity
eps             Thomsen's generic parameter as defined with resp. to sym. axis
delta           Thomsen's ..
gamma           Thomsen's ..

Output:
alpha           fracture plane compressional velocity
beta            S_parallel vertical velocity
ev              eps of equiv. VTI
dv              delta of ..
gv              gamma of ..

******************************************************************************
Notes:
Technical Reference:
Andreas Rueger, P-wave reflections and azimuthal AVO in TI-media
******************************************************************************
Author: Andreas Rueger, Colorado School of Mines, 01/27/96
******************************************************************************/
{
   double f;
   
   *ev    = -eps/(1.+2.*eps);

   *gv    = -gamma/(1.+2.*gamma);

   *alpha = vp*sqrt(1.+2.*eps);
   
   *beta  = vs*sqrt(1.+2.*gamma);
    
   f = vs/vp;

   f = 1 - f*f;
   
   *dv    = (delta -2.*eps*(1. + eps/f))/((1.+2.*eps)*(1.+2.*eps/f));

   return 1;
   
}

int v_phase2DVTI(Stiff2D *spar1, double sangle, int mode, double *v_phase)
/*****************************************************************************
v_phase2DVTI - Given phase angle, compute phase-velocity for TI media
******************************************************************************
Input:
spar1            (density normalized) stiffnesses
sangle           sin(phase_angle)
mode             mode (0=P 1=SV)

Output:
v_phase          phase-velocity for angle
******************************************************************************
Notes: Currently no SH mode supported
******************************************************************************
Technical reference: White, J. E..; Underground Sound.

******************************************************************************
Author: Andreas Rueger, Colorado School of Mines, 01/26/96
******************************************************************************/
{
	double a1111,a3333,a1133,a1313;
	double s2,c2,sqr,dummy;
	
	s2=sangle*sangle;
        c2=1.-s2;
	

	/***density normalized stiffnesses *************/
	a1111=spar1->a1111;
	a3333=spar1->a3333;
	a1133=spar1->a1133;
 	a1313=spar1->a1313;

	sqr=(a1111-a1313)*s2 - (a3333-a1313)*c2;
	sqr *=sqr;
	sqr +=4.*(a1133+a1313)*(a1133+a1313)*s2*c2;
	
	if(ABS(sqr)<FLT_EPSILON)
		sqr=0.;

	else if(sqr < 0){
		fprintf(stderr," error in v_phase \n");
		return (-1);
	}

        sqr=sqrt(sqr);
	
        dummy=(a1111+a1313)*s2+(a3333+a1313)*c2;
	
	if(mode==0)
	   *v_phase = sqrt(0.5*(dummy+sqr));
	else
           *v_phase = sqrt(0.5*(dummy-sqr));
	
	return 1;
}