refrealazihti.c

su 的源代码库

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

/* REFRELAAZIHTI: $Revision: 1.1 $ ; $Date: 2005/01/20 22:32:34 $	*/

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

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

#define EPS 0.0001

/* prototype of subroutine used internally */

int phasevel3DTIH(Stiff2D *spar, int mode, double d1, double d2, 	
	       double d3, double *v);
int p_hor3DTIH(Stiff2D *spar1, int mode, double sangle, double cangle, 	
	       double sazi, double cazi, double *p);
void polconv(Vector3D *d, int m, int rort);
int eigVect3DTIH(Stiff2D *spar, int mode, double p1, double p2, 
		double p3, Vector3D *d);
int p_vert3DTIH(Stiff2D *spar, int mode, double p, double s, double c, 
		int rort, double *q, Vector3D *d);
int graebner3D(Stiff2D *spar1, Stiff2D *spar2, double rho1, double rho2, 
	      int modei, int modet, int rort, double sazi, double cazi,
	      double p, double *b, double **a, int *ipvt, double *z,
	      double *rcond);
int testEikonal(Stiff2D *spar, int mode, double p1, double p2, 
		double p3, Vector3D *d);

/*********************** self documentation **********************/
char *sdoc[] = {
"									",
" REFREALAZIHTI -  REAL AZImuthal REFL/transm coeff for HTI media 	",
"									",
" refRealAziHTI  [optional parameters]	>coeff.data 			",
" 									",
" Optional parameters: 							",
" vp1=2         p-wave velocity medium 1 (with respect to symm.axes)	",
" vs1=1         s-wave velocity medium 1 (with respect to symm.axes)	",
" eps1=0        epsilon medium1						",
" delta1=0	delta medium 1						",
" gamma1=0	gamma medium 1						",
" rho1=2.7	density medium 1 					",
" vp2=2         p-wave velocity medium 2 (with respect to symm.axes)	",
" vs2=1         s-wave velocity medium 2 (with respect to symm.axes)	",
" eps2=0        epsilon medium 2					",
" delta2=0	delta medium 2						",
" gamma2=0	gamma medium 2						",
" rho2=2.7	density medium 2 					",
" modei=0 	incident mode is qP					",
"		=1 incident mode is qSV					",
"		=2 incident mode is SP					",
" modet=0 	scattered mode						",
" rort=1 	reflection(1)	or transmission (0)			",
" azimuth=0	azimuth with respect to x1-axis (clockwise)		",
" fangle=0	first incidence angle					",
" langle=45	last incidence angle					",
" dangle=1	angle increment						",
" iscale=0      default: angle in degrees				",
"		=1 angle-axis in rad                                    ",
"               =2 axis  horizontal slowness                            ",
"               =3 sin^2 of incidence angle                             ",
" ibin=1 	binary output 						",
"		=0 Ascci output						",
" outparfile	=outpar parameter file for plotting			",
" coeffile	=coeff.data coefficient-output file			",
" test=1 	activate testing routines in code			",
" info=0 	output intermediate results				",
"									",
" Notes:								",
" Axes of symmetry have to coincide in both media.  This code computes	",
" all 6 REAL reflection/transmissions coefficients on the fly. However,	",
" the set-up is such Real reflection/transmission coefficients in 	",
" HTI-media with coinciding symmetry axes.				",
" However, the set-up is such that currently only one coefficient is	",
" dumped into the output. This is easily changed.  The solution of the	",
" scattering problem is obtained numerically and involves the Gaussian	",
" elimination of a 6X6 matrix.						", 
"									",
NULL};

/*
 * AUTHOR:: Andreas Rueger, Colorado School of Mines, 02/10/95
 *                original name of code <graebnerTIH.c>
 *           modified, extended version of this code <refTIH3D>
 * 
 *  Technical references:
 *
 * 	Sebastian Geoltrain: Asymptotic solutions to direct
 *		and inverse scattering in anisotropic elastic media;
 *		CWP 082.
 *	Graebner, M.; Geophysics, Vol 57, No 11:
 *		Plane-wave reflection and transmission coefficients
 *		for a transversely isotropic solid.
 *	Cerveny, V., 1972, Seismic rays and ray intensities in inhomogeneous 	
 *		anisotropic media: Geophys. J. R. astr. Soc., 29, 1-13.
 *
 *	.. and some derivations by Andreas Rueger.
 *
 * If propagation is perpendicular or 
 * parallel to the symmetry axis, the solution is analytic (see 		",
 * graebner2D.c and rtRealIso.c		
 */
/**************** end self doc ********************************/

int test;
int info;

/* the main program */
int main (int argc, char **argv)
{
	double vp1,vp2,vs1,vs2,rho1,rho2;
	double eps1,eps2,delta1,delta2;
	double gamma1,gamma2,azimuth;
	float fangle,langle,dangle,angle;
	double *coeff,p=0;
	double sangle,cangle,sazi,cazi;
	float anglef,dummy;
	FILE *outparfp=NULL, *coeffp=NULL;
	int ibin,modei,modet,rort,iangle,iscale,index;
	char *outparfile=NULL,*coeffile=NULL;
	Stiff2D *spar1, *spar2;
	double **a,*rcond,*z;
	int *ipvt;

	/* allocate space for stiffness elements */
	spar1=(Stiff2D*)emalloc(sizeof(Stiff2D));
	spar2=(Stiff2D*)emalloc(sizeof(Stiff2D));

	/* allocate space for matrix system */
	a = alloc2double(6,6);
	coeff = alloc1double(6);
	ipvt=alloc1int(6);
	z = alloc1double(6);
	rcond=alloc1double(6);

	/* hook up getpar to handle the parameters */
	initargs(argc,argv);
	requestdoc(0);

	if (!getparint("ibin",&ibin)) ibin = 1;
	if (!getparint("modei",&modei)) modei = 0;
	if (!getparint("modet",&modet)) modet = 0;
	if (!getparint("rort",&rort)) rort = 1;
        if (!getparint("iscale",&iscale)) iscale = 0;
	if (!getparint("test",&test)) test = 1;
	if (!getparint("info",&info)) info = 0;

	if(modei != 0 && modei !=1 && modei !=2){
		fprintf(stderr," \n ERROR wrong incidence mode \n");
		return (-1);	/* wrong mode */
	}

	if(modet != 0 && modet !=1 && modet !=2){
		fprintf(stderr," \n ERROR wrong scattering mode \n");
		return (-1);	/* wrong mode */
	}

        if(rort != 0 && rort !=1){
                fprintf(stderr," ERROR wrong rort parameter \n");
                return (-1);    /* wrong mode */
        }

        if(iscale != 0 && iscale !=1 && iscale !=2 && iscale!=3 ){
                fprintf(stderr," ERROR wrong iscale parameter \n");
                return (-1);    /* wrong mode */
        }


	if (!getparfloat("fangle",&fangle)) fangle = 0.0;
	if (!getparfloat("langle",&langle)) langle = 45.0;
	if (!getparfloat("dangle",&dangle)) dangle = 1.0;

	if (!getpardouble("azimuth",&azimuth)) azimuth = 0.;
	if (!getpardouble("vp1",&vp1)) vp1 = 2.0;
	if (!getpardouble("vp2",&vp2)) vp2 = 2.0;
	if (!getpardouble("vs1",&vs1)) vs1 = 1.0;
	if (!getpardouble("vs2",&vs2)) vs2 = 1.0;
	if (!getpardouble("rho1",&rho1)) rho1 = 2.7;
	if (!getpardouble("rho2",&rho2)) rho2 = 2.7;

	if (!getpardouble("eps1",&eps1)) eps1 = 0.;
	if (!getpardouble("eps2",&eps2)) eps2 = 0.;

	if (!getpardouble("delta1",&delta1)) delta1 = 0.;
	if (!getpardouble("delta2",&delta2)) delta2 = 0.;
	if (!getpardouble("gamma1",&gamma1)) gamma1 = 0.;
	if (!getpardouble("gamma2",&gamma2)) gamma2 = 0.;

	if (getparstring("outparfile",&outparfile)) {
                 outparfp = efopen(outparfile,"w");
        } else {
                 outparfp = efopen("outpar","w");
        }

	if (getparstring("coeffile",&coeffile)) {
                 coeffp = efopen(coeffile,"w");
        } else {
                 coeffp = efopen("coeff.data","w");
        }


	/******   some debugging information ******************/
	if(info){
		ddprint(azimuth);
		ddprint(vp1); ddprint(vs1); ddprint(rho1);
		ddprint(eps1); ddprint(delta1); ddprint(gamma1);

		ddprint(vp2); ddprint(vs2); ddprint(rho2);
		ddprint(eps2); ddprint(delta2); ddprint(gamma2);
	}


	/* convert into rad */
	azimuth=azimuth*PI /180.;
	sazi=sin(azimuth);
	cazi=cos(azimuth);

	/******   convertion into cij's ************************/
	if (!thom2stiffTI(vp1,vs1,eps1,delta1,gamma1,PI/2.,spar1,1) ){
		fprintf(stderr," \n ERROR in thom2stiffTI (1) \n");
		return (-1);
	}

	if (!thom2stiffTI(vp2,vs2,eps2,delta2,gamma2,PI/2.,spar2,1) ){
		fprintf(stderr,"\n ERROR in thom2stiffTI (2) \n");
		return (-1);
	}

	/*****    more debugging output ************************/
	
	if(info){
		diprint(modei);
		diprint(modet);
		diprint(rort);
		ddprint(spar1->a1111);
		ddprint(spar1->a3333);
		ddprint(spar1->a1133);
		ddprint(spar1->a1313);
		ddprint(spar1->a2323);
		ddprint(spar1->a1212);

		ddprint(spar2->a1111);
		ddprint(spar2->a3333);
		ddprint(spar2->a1133);
		ddprint(spar2->a1313);
		ddprint(spar2->a2323);
		ddprint(spar2->a1212);
	}

	/********  find generated wave type-index     ************/
	/* reflect_P (0) reflect_S (1) transm_P (2) transm_S (3) */

	if(modet == 0 && rort==1)
		index = 0;
	else if(modet == 1 && rort==1)
		index = 1;
	else if(modet == 2 && rort==1)
		index = 2;
	else if(modet == 0 && rort==0)
		index = 3;
	else if(modet == 1 && rort==0)
		index = 4;
	else if(modet == 2 && rort==0)
		index = 5;
	else {
		fprintf(stderr,"\n ERROR wrong (index) \n ");
		return (-1);
	}


	/***************** LOOP OVER ANGLES ************************/
	for(angle=fangle,iangle=0;angle<=langle;angle+=dangle){
		
		if(info) ddprint(angle);
	
		sangle=(double) angle*PI/180;
		cangle=cos(sangle);
		sangle=sin(sangle);


		/* get horizontal slowness */
		if(p_hor3DTIH(spar1,modei,sangle,cangle,sazi,cazi,&p)!=1){
			fprintf(stderr,"\n ERROR in p_hor3DTIH \n ");
			return (-1);
		}

		/* compute reflection/transmission coefficient */
		if(graebner3D(spar1,spar2,rho1,rho2,modei,modet,rort,
		   sazi,cazi,p,coeff,a,ipvt,z,rcond)!=1){
			fprintf(stderr,"\n ERROR in p_hor3DTIH \n ");
			return (-1);
		}

		++iangle;

                if(iscale==0)
                     anglef=(float) angle;
                else if(iscale==1)
                     anglef=(float) angle*PI/180.;
                else if(iscale==2)
                     anglef=(float) p;
                else if(iscale==3) 
                     anglef=(float) sangle*sangle;
                  
                dummy= (float)coeff[index];
                  

                /* Binary output for x_t */
                if(ibin==1){

                        fwrite(&anglef,sizeof(float),1,coeffp);
                        fwrite(&dummy,sizeof(float),1,coeffp);

                /* ASCII output  */
                } else if(ibin==0){

                        fprintf(coeffp,"%f      %f\n",anglef,dummy);
		}
	}

	/*********  No of output pairs for plotting ********/
	if(ibin) fprintf(outparfp,"%i\n",iangle);

	return 1;
}


int graebner3D(Stiff2D *spar1, Stiff2D *spar2, double rho1, double rho2, 
	      int modei, int modet, int rort, double sazi, double cazi,
	      double p, double *b, double **a, int *ipvt, double *z,
	      double *rcond)
/*****************************************************************************
	Real reflection/transmission coefficients in
	TIH-media  with coinciding symmetry axes.
Input:
	spar1	density normalized stiffness components medium 1
	spar2	density normalized stiffness components medium 2
	modei	incident wave mode (=0 P; =1 SV; =2 SP)
	modet	scattered wave mode (=0 P; =1 SV; =2 SP)
	rort    reflection or transmission
	sazi	sin(azimuth from symmetry axis)
	cazi	cos(azimuth from symmetry axis)
	p	horizontal slowness component
Output:
	coeff	reflection/transmission coefficient

Technical references:
	Sebastian Geoltrain: Asymptotic solutions to direct
		and inverse scattering in anisotropic elastic media;
		CWP 082.
 	Graebner, M.; Geophysics, Vol 57, No 11:
		Plane-wave reflection and transmission coefficients
		for a transversely isotropic solid.
	Cerveny, V., 1972, Seismic rays and ray intensities in inhomogeneous 	
		anisotropic media: Geophys. J. R. astr. Soc., 29, 1-13.

	.. and some own derivations.

***************************************************************************
AUTHOR: Andreas Rueger, Colorado School of Mines, 01/29/95
***************************************************************************/
{
	double q_in,q_Pr,q_Pt,q_SVr,q_SVt,q_SPr,q_SPt;
	double c55i,c55t,c44i,c44t,c13i,c13t,c33i,c33t;
	Vector3D d_in,d_Pr,d_Pt,d_SVr,d_SVt,d_SPr,d_SPt;

	double p1=cazi*p;
	double p2=sazi*p;

	c55i=spar1->a1313 * rho1;
	c55t=spar2->a1313 * rho2;
	c44i=spar1->a2323 * rho1;
	c44t=spar2->a2323 * rho2;
	c13i=spar1->a1133 * rho1;
	c13t=spar2->a1133 * rho2;
	c33i=spar1->a3333 * rho1;
	c33t=spar2->a3333 * rho2;




	/* compute vertical slowness components and polarization vectors */

	/* incident wave */
	if(p_vert3DTIH(spar1,modei,p,sazi,cazi,0,&q_in,&d_in) !=1){
		fprintf(stderr,"\n ERROR in p_vert3DTIH (incident wave) \n ");
		return (-1);
	}
	if(test && testEikonal(spar1,modei,p1,p2,q_in,&d_in)!=1){
		fprintf(stderr,"\n ERROR in testEikonal (incident wave)");
		return (-1);
	}

	/* reflected qP */
	if(p_vert3DTIH(spar1,0,p,sazi,cazi,1,&q_Pr,&d_Pr) !=1){
		fprintf(stderr,"\n ERROR in p_vert3DTIH (reflected qP) \n ");
		return (-1);
	}
	if(test && testEikonal(spar1,0,p1,p2,q_Pr,&d_Pr)!=1){
		fprintf(stderr,"\n ERROR in testEikonal (reflected qP)");
		return (-1);
	}

	/* reflected SV */
	if(p_vert3DTIH(spar1,1,p,sazi,cazi,1,&q_SVr,&d_SVr) !=1){
		fprintf(stderr,"\n ERROR in p_vert3DTIH (reflected SV) \n ");
		return (-1);
	}
	if(test && testEikonal(spar1,1,p1,p2,q_SVr,&d_SVr)!=1){
		fprintf(stderr,"\n ERROR in testEikonal (reflected SV) \n");
		return (-1);
	}

	/* reflected SP */
	if(p_vert3DTIH(spar1,2,p,sazi,cazi,1,&q_SPr,&d_SPr) !=1){
		fprintf(stderr,"\n ERROR in p_vert3DTIH (reflected SP) \n ");
		return (-1);
	}
	if(test && testEikonal(spar1,2,p1,p2,q_SPr,&d_SPr)!=1){
		fprintf(stderr,"\n ERROR in testEikonal (reflected SP) \n");
		return (-1);
	}

	/* transmitted qP */
	if(p_vert3DTIH(spar2,0,p,sazi,cazi,0,&q_Pt,&d_Pt) !=1){
		fprintf(stderr,"\n ERROR in p_vert3DTIH (transmitted qP) \n ");
		return (-1);
	}
	if(test && testEikonal(spar2,0,p1,p2,q_Pt,&d_Pt)!=1){
		fprintf(stderr,"\n ERROR in testEikonal (transmitted qP) \n");
		return (-1);
	}

	/* transmitted SV */
	if(p_vert3DTIH(spar2,1,p,sazi,cazi,0,&q_SVt,&d_SVt) !=1){
		fprintf(stderr,"\n ERROR in p_vert3DTIH (transmitted SV) \n ");
		return (-1);
	}
	if(test && testEikonal(spar2,1,p1,p2,q_SVt,&d_SVt)!=1){
		fprintf(stderr,"\n ERROR in testEikonal (transmitted SV) \n");
		return (-1);
	}

	/* transmitted SP */
	if(p_vert3DTIH(spar2,2,p,sazi,cazi,0,&q_SPt,&d_SPt) !=1){
		fprintf(stderr,"\n ERROR in p_vert3DTIH (transmitted SP) \n ");
		return (-1);
	}
	if(test && testEikonal(spar2,2,p1,p2,q_SPt,&d_SPt)!=1){
		fprintf(stderr,"\n ERROR in testEikonal (transmitted SP) \n");
		return (-1);
	}


	/* compute matrix elements */
	a[0][0] = d_Pr.x; 	a[1][0] = d_SVr.x;	a[2][0] = d_SPr.x;
	a[3][0] = -d_Pt.x;	a[4][0] = -d_SVt.x;	a[5][0] = -d_SPt.x;
	
	a[0][1] = d_Pr.y; 	a[1][1] = d_SVr.y;	a[2][1] = d_SPr.y;
	a[3][1] = -d_Pt.y;	a[4][1] = -d_SVt.y;	a[5][1] = -d_SPt.y;

	a[0][2] = d_Pr.z; 	a[1][2] = d_SVr.z;	a[2][2] = d_SPr.z;
	a[3][2] = -d_Pt.z;	a[4][2] = -d_SVt.z;	a[5][2] = -d_SPt.z;

	a[0][3] =  c55i*(d_Pr.z*p1 +  d_Pr.x*q_Pr);
	a[1][3] =  c55i*(d_SVr.z*p1+d_SVr.x*q_SVr);
	a[2][3] =  c55i*(d_SPr.z*p1+d_SPr.x*q_SPr);
	a[3][3] = -c55t*(d_Pt.z*p1 +  d_Pt.x*q_Pt);
	a[4][3] = -c55t*(d_SVt.z*p1+d_SVt.x*q_SVt);
	a[5][3] = -c55t*(d_SPt.z*p1+d_SPt.x*q_SPt);

	a[0][4] =  c44i*(d_Pr.z*p2 +  d_Pr.y*q_Pr);