lincoeff.c

su 的源代码库

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


#include "par.h"

#define err_ang 30.*PI/180
#define LC_TINY FLT_EPSILON

int P_err_2nd_order(ErrorFlag *rp_1st, ErrorFlag *rp_2nd, 
				float true_kappa, int index);

float lincoef_Rp(float ang, float azim, float kappa, float *rpp, ErrorFlag *rp_1st, ErrorFlag *rp_2nd, int count)

     /* The subroutine for computation of weak-contrast-weak-anisotropy
	Rpp reflection coefficient for an arbitrary combination
	of ISO, VTI, HTI or ORTHORHOMBIC halfspaces, see adopted
	convention described in the main LinRORT.

	INPUT:
	ang = phase incidence angle (in rad)
	azim = phase azimuthal angle (in rad)
	kappa = mutual rotation of the reflecting halfspace with respect
	        to the incidence halfspace (in rad)
	count ... calling number of the subroutine Rp. If, and only if, count=1 
	          (first call) then the error quantities *rp_1st and *rp_2nd are
		  evaluated.

	External medium parameters of the both halfspaces: they must be
	pre-processed, i.e. correctly assigned for particular halfspaces. 
	There is no an appropriate check guard in this subroutine.

	OUTPUT:
	*rpp = Rpp(ang,azim,kappa, el.param.) reflection coefficient
	*rp_1st ... ErrorFlag structure: contains information on first-order
	            Rpp reflection coefficient evaluated for the purpose of
		    error analysis.
	*rp_2nd ... ErrorFlag structure: contains information on semi 
	            second-order Rpp reflection coefficient evaluated for 
		    the purpose of error analysis. 

	The error estimate is based on the comparison of the 1st and 2nd order 
	quantities from the structures *rp_1st and *rp_2nd.

	struct ErrorFlag
	        {
		   float iso[5];      ... isotropic part of Rpp for 5 different inc. angles
	           float upper[2];    ... ANISO/ISO Rpp for azimuths 0deg and 90 deg - upper hsp contribution
	           float lower[2];    ... ISO/ANISO Rpp for azimuths 0deg and 90 deg - lower hsp contribution
	           float global[4];   ... ANISO/ANISO Rpp for four azimuths - global model
	           float angle[4];    ... inc. angle, two azimuths and kappa for which the quantities
		                          above are evaluated
	        }; 


	err = 0...usual output
            = ANG_CR...incidence angle approaches the crytical angle ANG_CR
	      (approximately ang>=0.7*ANG_CR)
	      */
{

  extern float vp1,vp2,vs1,vs2,rho1,rho2;
  extern float eps1_1,eps1_2,delta1_1,delta1_2,delta1_3,gamma1_44,gamma1_55;
  extern float eps2_1,eps2_2,delta2_1,delta2_2,delta2_3,gamma2_44,gamma2_55;
  
  float err=0., ANG_CR;                               /* to control thinks */
  float Da_a, DG_G, alpha, beta, a_k;     /* to simplify thinks */
  float abs_lincoeff, S2, S2T2;                                /* to compute thinks */
 

  
  /* first, an approximate check for the crytical angle */

  if (vp1 <= vp2)
    {
      if (ang >= 0.7*(ANG_CR=asin(vp1/vp2)))
	err=ANG_CR;
    }

  /* some quantities in advance */

  Da_a=(vp2-vp1)/(0.5*(vp2+vp1));
  DG_G=(pow(vs2,2)*rho2-pow(vs1,2)*rho1)/(0.5*(pow(vs2,2)*rho2+pow(vs1,2)*rho1));
  /*-- DG_G=(rho2-rho1)/(0.5*(rho2+rho1)) + 2*(vs2-vs1)/2.85; --*/
  alpha=0.5*(vp2+vp1);
  beta=0.5*(vs2+vs1);
  

  /* first call => error quantities determination */

  if(count==1) 
    {
    P_err_2nd_order(rp_1st, rp_2nd, kappa, 1);
    if(kappa != 0.) P_err_2nd_order(rp_1st, rp_2nd, kappa, 2);
    }


  /*----- tests for different ratios alpha/beta -------
    fprintf(stderr,"alpha_orig=%f \t beta_orig=%f \t a/b_orig=%f \n",alpha,beta,beta/alpha);
    alpha=3.5;
    beta=1.9214265;
    vole
    fprintf(stderr,"alpha=%f \t beta=%f \t a/b=%f \n",alpha,beta,beta/alpha);
    ----------------------------------------------------*/


  /* now, the individual parts of the coefficient */

  a_k=azim-kappa;

  abs_lincoeff = (vp2*rho2 - vp1*rho1)/(vp2*rho2 + vp1*rho1);
  
  S2 = 0.5*(
	    Da_a - pow(2*beta/alpha,2)*DG_G + delta2_1*pow(sin(a_k),2) + delta2_2*pow(cos(a_k),2) -
	    delta1_1*pow(sin(azim),2) - delta1_2*pow(cos(azim),2) -  8*pow(beta/alpha,2)*(
				 gamma2_55*pow(cos(a_k),2) + gamma2_44*pow(sin(a_k),2) -
				 gamma1_55*pow(cos(azim),2) - gamma1_44*pow(sin(azim),2))
				 );
  
  /* -----------  no iso terms here 
     S2 = 0.5*(
	    delta2_1*pow(sin(a_k),2) + delta2_2*pow(cos(a_k),2) -
	    delta1_1*pow(sin(azim),2) - delta1_2*pow(cos(azim),2) -  8*pow(beta/alpha,2)*(
				 gamma2_55*pow(cos(a_k),2) + gamma2_44*pow(sin(a_k),2) -
				 gamma1_55*pow(cos(azim),2) - gamma1_44*pow(sin(azim),2))
				 );
  ----------------------------------*/

  S2T2 = 0.5*(
	      Da_a + 
	      eps2_1*pow(sin(a_k),4) + eps2_2*(pow(cos(a_k),4)+2*pow(cos(a_k),2)*pow(sin(a_k),2)) +
	      delta2_3*pow(cos(a_k),2)*pow(sin(a_k),2) -
	      eps1_1*pow(sin(azim),4) - eps1_2*(pow(cos(azim),4)+2*pow(sin(azim),2)*pow(cos(azim),2)) -
	      delta1_3*pow(cos(azim),2)*pow(sin(azim),2)
	      );
  

  /* finally, make the coefficient */

  /*  fprintf(stderr,"IN_ANG=%f AZIM=%f abs_lincoeff=%f \t S2=%f \t S2T2=%f \n",ang*180./PI,azim*180./PI,abs_lincoeff,S2,S2T2);
   */

  /*fprintf(stderr,"ang=%f azim=%f S3/S1=%f \n",ang*180/PI,azim*180/PI,S2T2/S2);*/

  *rpp = abs_lincoeff + S2*pow(sin(ang),2) + S2T2*pow(sin(ang),2)*pow(tan(ang),2);
  

  return(err*180/PI);
}

/********************************************/
/****************** the end *****************/
/********************************************/


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)

     /* The subroutine for computation of weak-contrast-weak-anisotropy
	Rps reflection coefficients for an arbitrary combination
	of ISO, VTI, HTI or ORTHORHOMBIC halfspaces, see adopted
	convention described in the main LinRORT.

	INPUT:
	ang = phase incidence angle (in rad)
	azim = phase azimuthal angle (in rad)
	kappa = mutual rotation of the reflecting halfspace with respect
	        to the incidence halfspace (in rad)
	i_hsp=0...iso, =1...VTI, =2...HTI, =3...ORT; the type of upper (incidence)
	      halfspace (strictly, this parameter would not be necesary, but knowing
	      the type of the halfspace in advance simplifies this subroutines)
	count ... calling number of the subroutine Rs. If, and only if, count=1 
	          (first call) then the error quantities *rsv_1st, *rsv_2nd, *rsh_1st  
		  and *rsh_2nd are evaluated.

	external medium parameters of the both halfspaces: they must be
	pre-processed, i.e. correctly assigned for all possible halfspaces. 
	There is no a check for that in this subroutine.

	OUTPUT:
	*rps1 = Rps1(ang,azim,kappa,el.parameters) reflection coefficient
	*rps2 = Rps2(ang,azim,kappa,el.parameters) reflection coefficient
	*sv = Rpsv(ang,azim,kappa,el.parameters) reflection coefficient
	*sh = Rpsh(ang,azim,kappa,el.parameters) reflection coefficient
	*cphi = cosine of rotation angle PHI, characterizing an approximate
	        rotation of *rps1 from *sv direction 
	*sphi = sine of rotation angle PHI, characterizing an approximate
	        rotation of *rps1 from *sv direction

	           the following holds:
	           rps1 = sv*cos(phi) + sh*sin(phi)
	           rps2 = - sv*sin(phi) + sh*cos(phi)


	*rsv_1st, *rsh_1st ... ErrorFlag structure: contains information on first-order
	                       Rpsv, Rpsh reflection coefficients evaluated for the purpose of
		               error analysis.
	*rsv_2nd, *rsh_2nd ... ErrorFlag structure: contains information on semi 
	                       second-order Rpsv, Rpsh reflection coefficient evaluated for 
		               the purpose of error analysis. 

	The error estimate is based on the comparison of the 1st and 2nd order 
	quantities from the structures *rsv_1st and *rsv_2nd, and *rsh_1st and *rsh_2nd.

	struct ErrorFlag
	        {
		   float iso[5];      ... isotropic part of Rpsv for 5 different inc. angles (for Rpsh, iso[n]=0)
	           float upper[2];    ... ANISO/ISO Rpsv (Rpsh) for azimuths 0deg and 90 deg (45deg) 
		                           - upper halfspace contribution
	           float lower[2];    ... ISO/ANISO Rpsv (Rpsh) for azimuths 0deg and 90 deg (45deg)
		                          - lower halfspace contribution
	           float global[4];   ... ANISO/ANISO Rpsv (Rpsh) for four (three) azimuths
		                          - global model
	           float angle[4];    ... inc. angle, azimuths and kappa for which the quantities
		                          above are evaluated
	        }; 


	err = 0...usual output
            = ANG_CR...inc. ang approaches the crytical angle ANG_CR 
	      (approximately ang>=0.7*ANG_CR)
	    =-1...inc.angle is a singular point (or very close to it: Vs1^2-Vs2^2<=1E-6 Km/s).
	          However, the incidence halfspace is ISO or VTI or HTI, so the 
		  computation of Rps is completed.
            =-2...inc.angle is a singular point (or very close to it: Vs1^2-Vs2^2<=1E-6 Km/s).
	          Incidence halfspace is ORT, so the computation of Rps is terminated
		  (-1 and -2 have priority over ANG_CR)
      */

{

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

  float err=0.,ANG_CR;  
  int perr;

  float alpha, beta, Da_a, Dr_r, Db_b, CS, S_C, CS3, S3_C, S5_C, S, CS_C, CS3_C, S3;
  float cos2a, sin2a, cos2a_k, sin2a_k, DG_G;
  
  
  int polar(float ang, float azim, float *cphi, float *sphi, int i_hsp);
  int S_err_2nd_order(ErrorFlag *rsv_1st, ErrorFlag *rsv_2nd, ErrorFlag *rsh_1st, ErrorFlag *rsh_2nd, 
		       float true_kappa, int index);

  /* first, an approximate check for the crytical angle */

  if (vp1 <= vp2)
    {
      if (ang >= 0.7*(ANG_CR=asin(vp1/vp2)))
	err=ANG_CR*180/PI;
    }

  /* compute cosine and sine of the rotation angle PHI */

  perr=polar(ang, azim, cphi, sphi, i_hsp);
  if (perr) err=perr;

  if (err == -2) return(err);
  
  /* now, determine the SV and SH terms in Rps coefficients */

  alpha = 0.5*(vp2+vp1);
  beta = 0.5*(vs2+vs1);
  Dr_r = (rho2-rho1)/(0.5*(rho2+rho1));
  Db_b = (vs2-vs1)/beta;
  Da_a = (vp2-vp1)/alpha;  
  cos2a_k = pow(cos(azim-kappa),2);
  sin2a_k = pow(sin(azim-kappa),2);
  cos2a = pow(cos(azim),2);
  sin2a = pow(sin(azim),2);

  /*** for DG_G insted of Dr_R + 2Db_b ****/
  DG_G=(pow(vs2,2)*rho2-pow(vs1,2)*rho1)/(0.5*(pow(vs2,2)*rho2+pow(vs1,2)*rho1)); 

  if(count==1) {
    S_err_2nd_order(rsv_1st, rsv_2nd, rsh_1st, rsh_2nd, kappa, 1);
    if(kappa != 0.) S_err_2nd_order(rsv_1st, rsv_2nd, rsh_1st, rsh_2nd, kappa, 2);
     ;
  }
  
  /*------ tests for different ratios alpha/beta -------
    fprintf(stderr,"alpha_orig=%f \t beta_orig=%f \t a/b_orig=%f \n",alpha,beta,beta/alpha);
  alpha=3.5;
  beta=1.9214265;
  vole
  delta2_3=-0.22;
  fprintf(stderr,"alpha=%f \t beta=%f \t a/b=%f \n",alpha,beta,beta/alpha);
      ------------------------------------------------------------------*/


  /*------------------      
  CS = -beta/alpha*Dr_r - 
    (alpha*beta)/(2*(pow(alpha,2)-pow(beta,2)))*(delta2_2*cos2a_k +
						 delta2_1*sin2a_k -
						 delta1_1*sin2a -
						 delta1_2*cos2a )
    - 2*beta/alpha*(gamma2_55*cos2a_k + gamma2_44*sin2a_k -
		    gamma1_55*cos2a - gamma1_44*sin2a +	Db_b);
  -----------------*/

  /*  for DG_G insted of Dr_R + 2Db_b: */
  CS = -beta/alpha*DG_G - 
      (alpha*beta)/(2*(pow(alpha,2)-pow(beta,2)))*(delta2_2*cos2a_k +
      delta2_1*sin2a_k -
      delta1_1*sin2a -
      delta1_2*cos2a )
      -2*beta/alpha*(gamma2_55*cos2a_k + gamma2_44*sin2a_k -
      gamma1_55*cos2a - gamma1_44*sin2a );
  

  
  S_C = -0.5*Dr_r +
    pow(alpha,2)/(2*(pow(alpha,2)-pow(beta,2)))*(delta2_2*cos2a_k +