refaniso.c

su 的源代码库

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


#include "par.h"
#include "anisotropy.h"

/*********************** self documentation **********************/
/***********************************************************************
REFANISO - Reflection coefficients for Anisotropic media

graebner - Reflection/Transmission coeff. (for VTI and TIH media)
           Coefficients are based on Graebner's paper.
gvelpolSH - compute SH group velocity and polarization for TI medium
gvelpolTI - compute P/SV group velocity and polarization for TI medium
p_hor2DTI - compute horizontal slowness for in-plane propagation in TI medium
p_vert2DVTI - Given the horizontal slowness, compute vertical slowness component
rottens2D - rotate 2-D elastic stiffness tensor
stiff2thomVTI - convert density normalized stiffness components 
		parameters into Thomsen's for transversely isotropic 
		material with vertical axes of symmetry
stiff2tv - Convert stiffnesses into Thomsen's parameter of the equivalent
          VTI model, P-wave reflections and azimuthal AVO in TI-media
thom2stiffTI - convert Thomsen's parameters into density normalized 
		stiffness components for transversely isotropic material
		with in-plane-tilted axes of symmetry
thom2tv - Convert generic Thomsen parameters into Thomsen's parameter 
          of the equivalent VTI model
v_phase2DVTI - Given phase angle, compute phase-velocity for TI media

******************************************************************************
graebner:
Input:
spar1,spar2      stiffnesses for both media
rho1,rho2        densities
pl               horizontal slowness
modei            incident mode (0=P 1=SV)
modet            scattered mode
rort             =1 reflection =0 transmission

Output:
coeff            ref/trans coeff

******************************************************************************
gvelpolSH:
Input:
aijkl		stiffness elements
px,pz		slowness elements

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

******************************************************************************
gvelpolTI:
Input:
aijkl		stiffness elements
px,pz		slowness elements

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

******************************************************************************
p_hor2DTI: 
Input:
spar		stiffness elements
s		sin(incidence angle)
mode		0=qP-Wave, 1=qSV-wave

Output:
*p	horizontal slowness component p_x

******************************************************************************
p_vert2DVTI:
Input:
spar1            (density normalized) stiffnesses
pl               horizontal slowness
modei            mode (0=P 1=SV)

Output:
p_vert           vertical slowness component

******************************************************************************
rottens2D:
Input:
aijkl		input/output stiffness elements
phi 		rotation angle (counterclock wise)

******************************************************************************
stiff2thomVTI :
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 

******************************************************************************
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 ..

******************************************************************************
thom2tv:
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 ..

******************************************************************************
v_phase2DVTI:
Input:
spar1            (density normalized) stiffnesses
sangle           sin(phase_angle)
mode             mode (0=P 1=SV)

Output:
v_phase          phase-velocity for angle
******************************************************************************


************************************************************************
Function prototypes:
int graebner2D(Stiff2D *spar1, double rho1, Stiff2D *spar2, double rho2,
	 double pl, int modei, int modet, int rort, double *coeff);
void gvelpolSH(double a1212, double a2323, double a1223, double px, 
		double pz, double *vgx, double *vgz, double *g11, 
		double *g13, double *g33)
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);
int p_hor2DTI (Stiff2D *spar, double s, int mode, double *p);
int p_vert2DVTI(Stiff2D *spar1, double pl, int modei, double *p_vert);
void rottens2D (Stiff2D *spar, double phi);
int stiff2thomVTI (double a1111, double a3333, double a1133, double a1313, 	
	double a1212, double *vp, double *vs, double *eps,
        double *delta, double *gamma);
int stiff2tv(Stiff2D *spar,double *alpha,double *beta,double *ev,
			double *dv,double *gv);
int thom2stiffTI (double vp, double vs, double eps, double delta, double gamma,
	double phi, Stiff2D *spar, int sign)
int thom2tv(double vp,double vs,double eps,double delta,double gamma,
	    double *alpha,double *beta,double *ev,double *dv,double *gv);
int v_phase2DVTI(Stiff2D *spar1, double sangle, int mode, double *v_phase);

************************************************************************
Author: CWP: Andreas Rueger, 1994-1996, Colorado School of Mines
***********************************************************************/
/**************** end self doc ********************************/



#define diprint(expr) printf(#expr " = %i\n",expr)
#define dfprint(expr) printf(#expr " = %f\n",expr)
#define ddprint(expr) printf(#expr " = %g\n",expr)


int graebner2D(Stiff2D *spar1, double rho1, Stiff2D *spar2, double rho2,
	 double pl, int modei, int modet, int rort, double *coeff)
/*****************************************************************************
graebner - Reflection/Transmission coeff. (for VTI and TIH media)
           Coefficients are based on Graebner's paper.
******************************************************************************
Input:
spar1,spar2      stiffnesses for both media
rho1,rho2        densities
pl               horizontal slowness
modei            incident mode (0=P 1=SV)
modet            scattered mode
rort             =1 reflection =0 transmission

Output:
coeff            ref/trans coeff

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

 Note the mistype in the equation for K1. The algorithm can be 
 used for VTI and TIH media on the incidence and scattering side

******************************************************************************
 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, 02/01/94
******************************************************************************/
{
	double K1,K2,K3;
	double qa1,qa2,qb1,qb2,sqr,p2;
	double la1,la2,lb1,lb2;
	double ma1,ma2,mb1,mb2;
	double a1,a2,b1,b2;
	double c1,c2,d1,d2;
	double qasq,qbsq;
	double a1111,a3333,a1133;
	double a1313,oa3333,oa1313;
	double m11,m12,m13,m14,m21,m22,m23,m24;
	double m31,m32,m33,m34,m41,m42,m43,m44;

	p2=pl*pl;

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

	/* diprint(modei);
	diprint(modet);
	diprint(rort);
	ddprint(spar1->a1111);
	ddprint(spar1->a3333);
	ddprint(spar1->a1133);
	ddprint(spar1->a1313);
	ddprint(spar2->a1111);
	ddprint(spar2->a3333);
	ddprint(spar2->a1133);
	ddprint(spar2->a1313); */


	/**********************  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,"ERROR in slowness medium1 \n");
		return (-1);
	}

	/* vertical slowness in medium 1 */
	qasq=0.5*((K1 - sqrt(sqr)));
	qa1=sqrt(qasq);
	qbsq=0.5*(K1 + sqrt(sqr));
	qb1=sqrt(qbsq);

	/* ddprint(qa1);ddprint(qb1); */

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

	sqr=1./sqr;

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

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

	/* P-eigenvector component medium 1*/
	la1=sqrt(la1);

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

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

	/* P-eigenvector component medium 1*/
	ma1=sqrt(ma1);

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

	sqr=1./sqr;

	lb1=a1313*qbsq + a1111*p2 -1.;
	
	lb1=lb1*sqr;

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

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

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

	mb1=a3333*qbsq + a1313*p2 -1.;
	mb1=mb1*sqr;

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

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

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

	/* ddprint(la1);ddprint(lb1);
	   ddprint(ma1);ddprint(mb1);*/

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

	a1=a1313*(qa1*la1 + pl*ma1);
	b1=a1313*(qb1*mb1 - pl*lb1);
	c1=pl*la1*a1133 + qa1*ma1*a3333;
	d1=pl*mb1*a1133 - qb1*lb1*a3333;

	/* ddprint(a1);ddprint(b1);
	   ddprint(c1);ddprint(d1);*/


	/**********************  Medium 2  ***********************************/

	a1111=spar2->a1111;
	a3333=spar2->a3333;
	a1133=spar2->a1133;
	a1313=spar2->a1313;

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

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