lincoeff.c

su 的源代码库

	SndKK_kappa = SndKK_0 + ((SndKK_90-SndKK_0)/D)*(FstKK_kappa-FstKK_0);
      }
      
      (*rp_1st).global[0]=Fst0_kappa;
      (*rp_2nd).global[0]=Snd0_kappa;
      (*rp_1st).global[1]=Fst90_kappa;
      (*rp_2nd).global[1]=Snd90_kappa;
      (*rp_1st).global[2]=FstK_kappa;
      (*rp_2nd).global[2]=SndK_kappa;
      (*rp_1st).global[3]=FstKK_kappa;
      (*rp_2nd).global[3]=SndKK_kappa;      

    }
  
			     
  return (0);
  
}

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


float P_term(float Da_a, float Dz_z, float DG_G, 
	     float Delta1_1, float Delta1_2, float Delta1_3,
	     float Delta2_1, float Delta2_2, float Delta2_3,  
	     float Eps1_1, float Eps1_2, float Eps2_1, float Eps2_2, 
	     float Gamma1_44, float Gamma2_44, float alpha, float beta, 
	     float ang, float azim, float a_k)

     /* This subroutine evaluates individual angular terms (inc. angle)  
	for Rpp 1st-order approximation. The input represents the medium parameters 
	(Thomsen's type aniso parameters (GammaX_44 represents GammaX_s), plus 
	S-wave impedance contrast ratio DG_G, density contrast ratio Dr_r and 
	background P- and S-wave velocities), and incidence angle "ang",
	azimuth "azim" and rotation angle "kappa", all being phase angles in radians.
     */

{
  float absolute, sn_2, sn_2_tg_2;
  
  absolute = Dz_z;

  sn_2 = 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_44*pow(sin(a_k),2) - Gamma1_44*pow(sin(azim),2)));

  sn_2_tg_2 = 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));
  
  /*fprintf(stderr,"dela1_1=%f delta1_2=%f gamma1s=%f delta2_1=%f delta2_2=%f gamma2s=%f \n",
    Delta1_1,Delta1_2,Gamma1_44,Delta2_1,Delta2_2,Gamma2_44);*/
  

  return(absolute + sn_2*pow(sin(ang),2) + sn_2_tg_2*pow(sin(ang),2)*pow(tan(ang),2));

}

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



  float Iso_exact(int type, float vp1, float vs1, float rho1, 
		  float vp2, float vs2, float rho2, float ang)  

     /* The subroutine returns the exact Rpp (type=0) or Rps (type=1) 
	real-valued reflection coefficient evaluated for a horizontal 
	interface separating two isotropic halfspaces. If the coefficient 
	is to be complex, the value -999.9 is returned. The subroutine 
	returns the value 1000 if the wave type "type" is chosen incorrectly.

	The other quantities are:
	vp1 ... P-wave velocity in the incidence halfspace
	vs1 ... S-wave velocity in the incidence halfspace
	rho1 ... density of the incidence halfspace
	vp2 ... P-wave velocity in the reflecting halfspace
	vs2 ... S-wave velocity in the reflecting halfspace
	rho2 ... density of the reflecting halfspace
	ang ... incidence angle in radians
     */

{
  float p, P1, P2, P3, P4;
  float q, X, Y, Z, D;
  float COEF;
  
  
  /***** some auxiliary quantities ***/

  p=sin(ang)/vp1;
  if((P1=1-pow(vp1*p,2)) < 0.) return(-999.9);
  if((P2=1-pow(vs1*p,2)) < 0.) return(-999.9);
  if((P3=1-pow(vp2*p,2)) < 0.) return(-999.9);
  if((P4=1-pow(vs2*p,2)) < 0.) return(-999.9);
  P1=sqrt(P1);
  P2=sqrt(P2);
  P3=sqrt(P3);
  P4=sqrt(P4);

  q=2*(rho2*pow(vs2,2) - rho1*pow(vs1,2));
  X=rho2 - q*pow(p,2);
  Y=rho1 + q*pow(p,2);
  Z=rho2 - rho1 - q*pow(p,2);
  
  D=pow(q*p,2)*P1*P2*P3*P4 +
    rho1*rho2*(vp1*vs2*P2*P3 + vs1*vp2*P1*P4) +
    vp1*vs1*P3*P4*pow(Y,2) +
    vp2*vs2*P1*P2*pow(X,2) +
    vp1*vp2*vs1*vs2*pow(p,2)*pow(Z,2);
  
  /*** the coefficient ***/

  if(type == 0)
    {
      COEF = (pow(q*p,2)*P1*P2*P3*P4 +
	      rho1*rho2*(vs1*vp2*P1*P4 - vp1*vs2*P2*P3) -
	      vp1*vs1*P3*P4*pow(Y,2) +
	      vp2*vs2*P1*P2*pow(X,2) -
	      vp1*vp2*vs1*vs2*pow(p,2)*pow(Z,2))/D;    
    }
  else if(type == 1) 
    {
      COEF = (-2*vp1*p*P1*(q*P3*P4*Y + vp2*vs2*X*Z))/D;
    }
  else 
    {
      COEF = 1000.;
    }
  
  return (COEF);
  
}


int polar(float ang, float azim, float *cphi, float *sphi, int i_hsp)

     /* cos(PHI) and sin(PHI) of rotation angle PHI describing
	a rotation of the polarization vector of S1 (S2)  wave
	in the isotropic slowness plane:

	ang = INCIDENCE phase angle (from the vertical) in rad
	azim = azimuth in rad
	cphi=cos(PHI)
	sphi=sin(PHI)
	i_hsp=0,1,2,3 for ISO,VTI,HTI,ORT respectively
	Phi si counted in right-hand Cartesian c.s. SV-SH,
	SV pointing downward. Positive Phi is from SV->SH.

	the subroutine requires the generalized medium parameters
	as the function LinRps, plus the parameter gamma1_1, which 
	enter the function as extern 

        perr=0 O.K.
	=-1 a singular point in VTI or HTI, computation completed
	=-2 a singular point in ORT, computation terminated 
	=-100 if a wrong choice of the inc. halfspace */


{
  extern float vp1,vs1;
  extern float eps1_1,eps1_2,delta1_1,delta1_2,delta1_3,gamma1_1,gamma1_44,gamma1_55;

  int c,perr=0;
  float B12, B11_B22, D;
  float A55Gs, A55G1, A33;
  float VS,sin_j, cos_j, sin_k, cos_k;
  
  
  /* do not do anything, if ISO inc. halfspace */

  if (i_hsp == 0)
    {
      *cphi=1.;
      *sphi=0.;
      return(perr=0);
    }  

  /* auxiliary quantities */

  A55Gs = pow(vs1,2)*(gamma1_44 - gamma1_55);
  A55G1 = pow(vs1,2)*(1 + 2*gamma1_55)*gamma1_1;
  A33 = pow(vp1,2);

  /* get the reflection angle of the iso S-wave (in rad), and azimuth angles */

  VS = sqrt( (2*pow(vs1,2)*(1+2*gamma1_44) + pow(vs1,2)*(1+2*gamma1_55)*(1+gamma1_1))/3 );
  
  sin_j = (VS/vp1)*sin(ang);
  cos_j = 1 - 0.5*pow(sin_j,2);
  cos_k = cos(azim);
  sin_k = sin(azim);

  /*fprintf(stderr,"vp=%f vs=%f ang=%f azim=%f a3=%f a55gs=%f a55g1=%f Sj=%f Cj=%f Saz=%f Caz=%f \n",vp1,vs1,ang,azim,A33,
	  A55Gs,A55G1,sin_j,cos_j,sin_k,cos_k);
  */
  
  /* get the polarization coefficients */

  B11_B22 = 2*A33*(eps1_2*pow(cos_k,2)-delta1_2)*pow(cos_k,2)*pow(cos_j,2)*pow(sin_j,2)
          + 2*A33*(eps1_1*pow(sin_k,2)-delta1_1)*pow(sin_k,2)*pow(cos_j,2)*pow(sin_j,2)
    + 2*A33*(1+2*eps1_2)*delta1_3*pow(cos_k,2)*pow(sin_k,2)*pow(sin_j,2)*(pow(cos_j,2)+1)
    + 4*A33*eps1_2*pow(cos_k,2)*pow(sin_k,2)*pow(cos_j,2)*pow(sin_j,2) 
    + 2*A33*(eps1_2-eps1_1)*pow(cos_k,2)*pow(sin_k,2)*pow(sin_j,2)
    - 2*A55G1*pow(sin_j,2)
    - 2*A55Gs*(pow(cos_k,2)-pow(sin_k,2))
    + 2*A55Gs*pow(cos_k,2)*pow(sin_j,2);
  
  B12 = A33*(1+2*eps1_2)*delta1_3*cos_k*sin_k*(pow(cos_k,2)-pow(sin_k,2))*cos_j*pow(sin_j,2)
    + A33*(delta1_2-delta1_1)*cos_k*sin_k*cos_j*pow(sin_j,2)
    + 2*A33*(eps1_1-eps1_2)*cos_k*pow(sin_k,3)*cos_j*pow(sin_j,2)
    + 2*A55Gs*cos_k*sin_k*cos_j;
  
  D = sqrt(pow(B11_B22,2) + 4*pow(B12,2));

  /*  fprintf(stderr,"B1122=%f B12=%f \n",B11_B22,B12);
   */

  if (D <= LC_TINY) perr=-1;
  
  /* do the logic for different halfspaces */

  if (i_hsp == 1){
    *cphi=1.;
    *sphi=0.;
    return(perr);     
  }
  if (i_hsp == 2){
    *cphi=(cos_j*cos_k)/sqrt(1-pow(sin_j*cos_k,2));
    *sphi=-sin_k/sqrt(1-pow(sin_j*cos_k,2));
    return(perr);
  }
  if (i_hsp == 3){
    if (perr == -1) return(perr=-2);
    if (fabs(B12) <= LC_TINY)
      {
      if (fabs(B11_B22) <= LC_TINY) return(perr=-2);
      switch (c=floor(B11_B22/fabs(B11_B22))) 
	{
	case 1: case 2:
	  *cphi=1;
	  *sphi=0;
	  break;
	case 0: case -1: case -2:
	  *cphi=0;
	  *sphi=1;
	  break;
	default : return(perr=-2);
	}
      }
    else {
      *cphi=sqrt( 0.5*(1+B11_B22/D)  );
      *sphi=(B12/fabs(B12))*sqrt( 0.5*(1-B11_B22/D) );
    }
    
    return(0);
   
  }
  
  return(-100);

}

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


int S_err_2nd_order(ErrorFlag *rsv_1st, ErrorFlag *rsv_2nd, ErrorFlag *rsh_1st, ErrorFlag *rsh_2nd, 
		    float true_kappa, int index)


     /* The subroutine evaluates first order (exactly) and second order (approximately)
	Rpsv and Rpsh coefficients for a given incidence angle and azimuths. The second order evaluation
	is based on the numerical modeling and comparing of the exact Rpsv (Rpsh) and its first order 
	approximation (so called semi-2nd order). The comparison has been done for the incidence 
	angle of 30 deg and it should be most accurate for this angle. However, such an estimated 
	semi-2nd order Rpsv (Rpsh) usually resulted in a global increase of accuracy of the first-order Rpsv 
	(Rpsh) so the subroutine can be used for other angles, as a rough error estimate, as well.

	INPUT:
	true_kappa ... mutual rotation of the upper and lower halfspaces in rad. This angle affects
	               the global error evaluation only.
	index ... =1, =2: if the true_kappa=0, then the subroutine should be called just once
		          with index=1. If the true_kappa is not zero, then the
		          subroutine needs to be called the second time (index=2) in
		          order to properly evaluate the GLOBAL error. No other values of
			  the "index" parameter are allowed. If some other value is on the
			  entry, the subroutine is terminated immediately with the value -1.
	err_ang ... as #define. The incidence angle for which the error analysis is
	                        performed. Not recommended to change the current value of
				30deg by a large ammount.
	lkappa, kkappa ... in declaration. Two azimuths for which the global error for Rpsv is evaluated,
	                                   besides 0deg and 90deg; their values can be changed arbitrary, 
					   but do not forget to change your output information as well 
					   (i.e. give proper azimuths for the output file corresponding to
					   *rp_1st(2nd).global[2],[3]). In the case of Rpsh, the angles are
					   given as 20deg, 45deg and 70deg, and can be changed in the
					   paragraph GLOBAL -> SH -> err_kappa=90deg section.

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

        	struct ErrorFlag
	         {