linrort.c

su 的源代码库

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

/* LINRORT: $Revision: 1.1 $ ; $Date: 2003/08/19 21:24:44 $	*/

#include <par.h>		/* definuje getchar, sdoc a error funkci */


/*********************** self documentation **********************/
char *sdoc[] = {
"									",
" LINRORT - linearized P-P, P-S1 and P-S2 reflection coefficients 	",
"		for a horizontal interface separating two of any of the	",
"		following halfspaces: ISOTROPIC, VTI, HTI and ORTHORHOMBIC. ",
"									",
" linrort [optional parameters]						",
"									",
" hspace1=ISO	medium type of the incidence halfspace:		 	",
"		=ISO ... isotropic					",
"		=VTI ... VTI anisotropy				 	",
"		=HTI ... HTI anisotropy				 	",
"		=ORT ... ORTHORHOMBIC anisotropy			",
" for ISO:								",
" vp1=2	 P-wave velocity, halfspace1					",
" vs1=1	 S-wave velocity, halfspace1					",
" rho1=2.7	density, halfspace1					",
"									",
" for VTI:								",
" vp1=2	 P-wave vertical velocity (V33), halfspace1			",
" vs1=1	 S-wave vertical velocity (V44=V55), halfspace1			",
" rho1=2.7	density, halfspace1					",
" eps1=0	Thomsen's generic epsilon, halfspace1			",
" delta1=0	Thomsen's generic delta, halfspace1			",
" gamma1=0	Thomsen's generic gamma, halfspace1			", 
"									",
" for HTI:								",
" vp1=2	 P-wave vertical velocity (V33), halfspace1			",
" vs1=1	 \"fast\" S-wave vertical velocity (V44), halfspace1		",
" rho1=2.7	density, halfspace1					",
" eps1_v=0	Tsvankin's \"vertical\" epsilon, halfspace1		",
" delta1_v=0	Tsvankin's \"vertical\" delta, halfspace1		",
" gamma1_v=0	Tsvankin's \"vertical\" gamma, halfspace1		",	
"									",
" for ORT:								",
" vp1=2	 P-wave vertical velocity (V33), halfspace1			",
" vs1=1	 x2-polarized S-wave vertical velocity (V44), halfspace1 	",
" rho1=2.7	density, halfspace1					",
" eps1_1=0	Tsvankin's epsilon in [x2,x3] plane, halfspace1	 	",
" delta1_1=0	Tsvankin's delta in [x2,x3] plane, halfspace1		",
" gamma1_1=0	Tsvankin's gamma in [x2,x3] plane, halfspace1	  	",
" eps1_2=0	Tsvankin's epsilon in [x1,x3] plane, halfspace1		",
" delta1_2=0	Tsvankin's delta in [x1,x3] plane, halfspace1		",
" gamma1_2=0	Tsvankin's gamma in [x1,x3] plane, halfspace1	  	",
" delta1_3=0	Tsvankin's delta in [x1,x2] plane, halfspace1		",
"									",
" hspace2=ISO	medium type of the reflecting halfspace (the same	",
"		convention as above)					",
"									",
" medium parameters of the 2nd halfspace follow the same convention	",
" as above:								",
"									",
" vp2=2.5		 vs2=1.2		rho2=3.0		",
" eps2=0		  delta2=0					",
" eps2_v=0		delta2_v=0		gamma2_v=0		",
" eps2_1=0		delta2_1=0		gamma2_1=0		",
" eps2_2=0		delta2_2=0		gamma2_2=0		",
" delta2_3=0								",
"									",
"	(note you do not need \"gamma2\" parameter for evaluation	",
"	of weak-anisotropy reflection coefficients)			",
"									",
" a_file=-1	the string '-1' ... incidence and azimuth angles are	",
"		generated automatically using the setup values below	",
"		a_file=file_name ... incidence and azimuth angles are	",
"		read from a file \"file_name\"; the program expects a	",
"		file of two columns [inc. angle, azimuth]		",
"									",
" in the case of a_file=-1:						",
" fangle=0	first incidence phase angle				",
" langle=30	last incidence angle					",
" dangle=1	incidence angle increment				",
" fazim=0	first azimuth (in deg)				  	",
" lazim=0	last azimuth  (in deg)				  	",
" dazim=1	azimuth increment (in deg)				",
"									",
" kappa=0.	azimuthal rotation of the lower halfspace2 (e.t. a	",
"		symmetry axis plane for HTI, or a symmetry plane for	",
"		ORTHORHOMBIC) with respect to the x1-axis		",
"									",
" out_inf=info.out	information output file				",
" out_P=Rpp.out	file with Rpp reflection coefficients			",
" out_S=Rps.out	file with Rps reflection coefficients			",
" out_SVSH=Rsvsh.out  file with SV and SH projections of reflection	",
"			coefficients					",
" out_Error=error.out file containing error estimates evaluated during  ",
"			the computation of the reflection coefficients;	",
"									",
"									",
" Output:								",
" out_P:								",
" inc. phase angle, azimuth, reflection coefficient; for a_file=-1, the ",
" inc. angle is the fast dimension					",
" out_S:								",
" inc. phase angle, azimuth, Rps1, Rps2, cos(PHI), sin(PHI); for	",
" a_file=-1, the inc. angle is the fast dimension			", 
" out_SVSH:								",
" inc. phase angle, azimuth, Rsv, Rsh, cos(PHI), sin(PHI); for	  ",
" a_file=-1, the inc. angle is the fast dimension			",
" out_Error:								",
" error estimates of Rpp, Rpsv and Rpsh approximations; global error is ",
" analysed as well as partial contributions to the error due to the	",
" isotropic velocity contrasts, and due to anisotropic  upper and lower ",
" halfspaces. The error file is self-explanatory, see also descriptions ",
" of subroutines P_err_2nd_order, SV_err_2nd_order and SH_err_2nd_order.",
"									",
"									",
" Adopted Convention:							",
"									",
" The right-hand Cartesian coordinate system with the x3-axis pointing  ",
" upward has been chosen. The upper halfspace (halfspace1)		",
" contains the incident P-wave. Incidence angles can vary from <0,PI/2),",
" azimuths are unlimited, +azimuth sense counted from x1->x2 axes	",
" (azimuth=0 corresponds to the direction of x1-axis). In the current	",
" version, the coordinate system is attached to the halfspace1 (e.t.	",
" the symmetry axis plane of HTI halfspace1, or one of symmetry planes  ",
" of ORTHORHOMBIC halfspace1, is aligned with the x1-axis), however, the",
" halfspace2 can be arbitrarily rotated along the x3-axis with respect  ",
" to the halfspace1. The positive weak-anisotropy polarization of the	",
" reflected P-P wave (e.t. positive P-P reflection coefficient) is close",
" to the direction of isotropic slowness vector of the wave (pointing	",
" outward the interface). Similarly, weak-anisotropy S-wave reflection  ",
" coefficients are described in terms of \"SV\" and \"SH\" isotropic	",
" polarizations, \"SV\" and \"SH\" being unit vectors in the plane	",
" perpendicular to the isotropic slowness vector. Then, the positive	",
" \"SV\" polarization vector lies in the incidence plane and points	",
" towards the interface, and positive \"SH\" polarization vector is	",
" perpendicular to the incidence plane, aligned with the positive	",
" x2-axis, if azimuth=0. Rotation angle \"PHI\", characterizing a	",
" rotation of \"the best projection\" of the S1-wave polarization	",
" vector in the isotropic SV-SH plane in the incidence halfspace1, is	",
" counted in the positive sense from \"SV\" axis (PHI=0) towards the	",
" \"SH\" axis (PHI=PI/2). Of course, S2 is perpendicular to S1, and	",
" the projection of S1 and S2 polarizations onto the SV-SH plane	",
" coincides with SV and SH directions, respectively, for PHI=0.		",
"									",
" The units for velocities are km/s, angles I/O are in degrees		",
"									",
" Additional Notes:							",
"	The coefficients are computed as functions of phase incidence	",
"	angle and azimuth (determined by the incidence slowness vector).",
"	Vertical symmetry planes of the HTI and				",
"	ORTHORHOMBIC halfspaces can be arbitrarily rotated along the	",
"	x3-axis. The linearization is based on the assumption of weak	", 
"	contrast in elastic medium parameters across the interface,	",
"	and the assumption of weak anisotropy in both halfspaces.	",
"	See the \"Adopted Convention\" paragraph below for a proper	",
"	input.								",
"									",
NULL};

/* 
 * Author: Petr Jilek, CSM-CWP, December 1999.
 */

/**************** end self doc ********************************/

/* this sets up the S-background; optimum for HTI1/HTI2 is HTIback1=a55, */
/* HTIback2=a44, optimum for ORT unknown yet at this point*/

#define HTIback1 55	/* a55 or a44 in HTI1 (44 keeps correct values	*/
			/*  at 90 deg of azimuth)			*/
#define HTIback2 55	/* a55 or a44 in HTI2				*/
#define ORTback1 55	/* a55 or a44 in ORT1 , a55 is standart		*/
#define ORTback2 55	/* a55 or a44 in ORT2, a55 is standart		*/

float vp1,vp2,vs1,vs2,rho1,rho2;
float eps1_1,eps1_2,delta1_1,delta1_2,delta1_3,gamma1_1,gamma1_44,gamma1_55;
float eps2_1,eps2_2,delta2_1,delta2_2,delta2_3,gamma2_44,gamma2_55;

int main(int argc, char **argv)
{

	extern float vp1,vp2,vs1,vs2,rho1,rho2;
	extern float eps1_1,eps1_2,delta1_1,delta1_2,delta1_3,gamma1_1,gamma1_44,gamma1_55;
  extern float eps2_1,eps2_2,delta2_1,delta2_2,delta2_3,gamma2_44,gamma2_55;

  float eps1,delta1,eps2,delta2,gamma1;
  float eps1_v,delta1_v,gamma1_v,eps2_v,delta2_v,gamma2_v;
  float gamma1_2,gamma2_1,gamma2_2;
  float fangle,langle,dangle,fazim,lazim,dazim,kappa,ang,azim;

  float RPP,RPS1,RPS2,SV,SH,cosPHI,sinPHI;
  float errp,errs, Max, Iso_e[4]={0.,0.,0.,0.};

  int HSP1=-1,HSP2=-1,count=1,n,back=0;
  
  char *hspace1=NULL, *hspace2=NULL;
  char *a_file=NULL, *out_inf=NULL, *out_P=NULL,*out_S=NULL, *out_SVSH=NULL, *out_Error=NULL;
  FILE *a_file_p=NULL,*out_inf_p=NULL,*out_P_p=NULL,*out_S_p=NULL,*out_SVSH_p=NULL,*out_Error_p=NULL;

  ErrorFlag Rp_1st, Rp_2nd, Rsv_1st, Rsv_2nd, Rsh_1st, Rsh_2nd;
  
  
  float lincoef_Rp(float ang, float azim, float kappa, float *rpp, ErrorFlag *rp_1st,
		ErrorFlag *rp_2nd, int count);
  float lincoef_Rs(float ang, float azim, float kappa, float *rps1, 
		float *rps2, float *sv, float *sh, float *cphi, float *sphi, int i_hsp,
		ErrorFlag *rsv_1st, ErrorFlag *rsv_2nd, ErrorFlag *rsh_1st, ErrorFlag *rsh_2nd, int count);


  /*auxiliar quantities for inversion info */
  float beta, alpha, Da_a, Db_b, DG_G, DZ_Z, PS2,PSC,PABS,SS2,SSC,SABS;


  

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

  if (!getparstring("hspace1",&hspace1)) hspace1 = "ISO";
  if (!getparfloat("vp1",&vp1)) vp1 = 2.0;
  if (!getparfloat("vs1",&vs1)) vs1 = 1.0;
  if (!getparfloat("rho1",&rho1)) rho1 = 2.7;
  if (!getparfloat("eps1",&eps1)) eps1 = 0;
  if (!getparfloat("delta1",&delta1)) delta1 = 0.;
  if (!getparfloat("gamma1",&gamma1)) gamma1 = 0.;
  if (!getparfloat("eps1_v",&eps1_v)) eps1_v = 0;
  if (!getparfloat("delta1_v",&delta1_v)) delta1_v = 0.;
  if (!getparfloat("gamma1_v",&gamma1_v)) gamma1_v = 0.;
  if (!getparfloat("eps1_1",&eps1_1)) eps1_1 = 0;
  if (!getparfloat("delta1_1",&delta1_1)) delta1_1 = 0.;
  if (!getparfloat("gamma1_1",&gamma1_1)) gamma1_1 = 0.;
  if (!getparfloat("eps1_2",&eps1_2)) eps1_2 = 0;
  if (!getparfloat("delta1_2",&delta1_2)) delta1_2 = 0.;
  if (!getparfloat("gamma1_2",&gamma1_2)) gamma1_2 = 0.;
  if (!getparfloat("delta1_3",&delta1_3)) delta1_3 = 0.;

  if (!getparstring("hspace2",&hspace2)) hspace2 = "ISO";
  if (!getparfloat("vp2",&vp2)) vp2 = 2.5;
  if (!getparfloat("vs2",&vs2)) vs2 = 1.2;
  if (!getparfloat("rho2",&rho2)) rho2 = 3.0;
  if (!getparfloat("eps2",&eps2)) eps2 = 0;
  if (!getparfloat("delta2",&delta2)) delta2 = 0.;
  if (!getparfloat("eps2_v",&eps2_v)) eps2_v = 0;
  if (!getparfloat("delta2_v",&delta2_v)) delta2_v = 0.;
  if (!getparfloat("gamma2_v",&gamma2_v)) gamma2_v = 0.;
  if (!getparfloat("eps2_1",&eps2_1)) eps2_1 = 0;
  if (!getparfloat("delta2_1",&delta2_1)) delta2_1 = 0.;
  if (!getparfloat("gamma2_1",&gamma2_1)) gamma2_1 = 0.;
  if (!getparfloat("eps2_2",&eps2_2)) eps2_2 = 0;
  if (!getparfloat("delta2_2",&delta2_2)) delta2_2 = 0.;
  if (!getparfloat("gamma2_2",&gamma2_2)) gamma2_2 = 0.;
  if (!getparfloat("delta2_3",&delta2_3)) delta2_3 = 0.;
  
  if (!getparstring("a_file",&a_file)) a_file = "-1";  
  if (!getparfloat("fangle",&fangle)) fangle = 0.;
  if (!getparfloat("langle",&langle)) langle = 30.;
  if (!getparfloat("dangle",&dangle)) dangle = 1.;
  if (!getparfloat("fazim",&fazim)) fazim = 0.;
  if (!getparfloat("lazim",&lazim)) lazim = 0.;
  if (!getparfloat("dazim",&dazim)) dazim = 1.;
  if (!getparfloat("kappa",&kappa)) kappa = 0.;
  
  if (!getparstring("out_inf", &out_inf))  out_inf = "info.out" ;
  if (!getparstring("out_P", &out_P))  out_P = "Rpp.out" ;
  if (!getparstring("out_S", &out_S))  out_S = "Rps.out" ; 
  if (!getparstring("out_SVSH", &out_SVSH))  out_SVSH = "Rsvsh.out" ; 
  if (!getparstring("out_Error", &out_Error))  out_Error = "error.out" ; 


  /* open the files and "ID" the halfspaces*/

  if (strcmp(a_file,"-1")) 
	a_file_p = efopen(a_file, "r");
  out_inf_p = efopen(out_inf, "w");
  out_P_p = efopen(out_P, "w");
  out_S_p = efopen(out_S, "w");
  out_SVSH_p = efopen(out_SVSH, "w");
  out_Error_p = efopen(out_Error, "w");


  if((!strcmp(hspace1,"ISO"))||(!strcmp(hspace1,"iso")))HSP1=0;
  if((!strcmp(hspace1,"VTI"))||(!strcmp(hspace1,"vti")))HSP1=1;
  if((!strcmp(hspace1,"HTI"))||(!strcmp(hspace1,"hti")))HSP1=2;
  if((!strcmp(hspace1,"ORT"))||(!strcmp(hspace1,"ort")))HSP1=3;
  if((!strcmp(hspace2,"ISO"))||(!strcmp(hspace2,"iso")))HSP2=0;
  if((!strcmp(hspace2,"VTI"))||(!strcmp(hspace2,"vti")))HSP2=1;
  if((!strcmp(hspace2,"HTI"))||(!strcmp(hspace2,"hti")))HSP2=2;
  if((!strcmp(hspace2,"ORT"))||(!strcmp(hspace2,"ort")))HSP2=3;


  fprintf(stdout,"LinRORT working ...\n");
  /* necessary checks */

  if((HSP1==-1)||(HSP2==-1))
	{
	err(" wrong halfspace denotation; error in: hspace1 or hspace2 parameter.  Program terminated!");
	return(-1);
	}  
  if ((fangle < 0.|| fangle >= 90.) || (langle <= 0.|| langle >= 90.)
	|| (fangle > langle) || (dangle > 1+(langle-fangle)) || (dangle == 0) ) 
	{