lincoeff.c

su 的源代码库

						 delta2_1*sin2a_k -
						 delta1_2*cos2a -
						 delta1_1*sin2a);

  CS3 = -alpha*beta/(pow(alpha,2)-pow(beta,2))*(delta1_1*sin2a + 
						delta1_2*cos2a -
						delta1_3*cos2a*sin2a -
						eps1_2*(cos2a*cos2a+2*cos2a*sin2a) -
						eps1_1*sin2a*sin2a +
						eps2_2*(cos2a_k*cos2a_k+2*sin2a_k*cos2a_k) +
						eps2_1*sin2a_k*sin2a_k -
						delta2_1*sin2a_k -
						delta2_2*cos2a_k +
						delta2_3*cos2a_k*sin2a_k);
   
  /*
    S3_C = pow(beta,2)/pow(alpha,2)*Dr_r +
    pow(beta,2)/(2*(pow(alpha,2)-pow(beta,2)))*(delta1_1*sin2a +
						delta1_2*cos2a -
						delta2_1*sin2a_k -
						delta2_2*cos2a_k) +
    pow(alpha,2)/(pow(alpha,2)-pow(beta,2))*(delta1_1*sin2a +
					     delta1_2*cos2a -
					     delta1_3*cos2a*sin2a -
					     eps1_1*sin2a*sin2a -
					     eps1_2*(cos2a*cos2a+2*cos2a*sin2a) +
					     eps2_1*sin2a_k*sin2a_k +
					     eps2_2*(cos2a_k*cos2a_k+2*sin2a_k*cos2a_k) -
					     delta2_1*sin2a_k -
					     delta2_2*cos2a_k +
					     delta2_3*cos2a_k*sin2a_k) +
    2*pow(beta/alpha,2)*(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:  */
 
  
  S3_C = pow(beta,2)/pow(alpha,2)*DG_G +
    pow(beta,2)/(2*(pow(alpha,2)-pow(beta,2)))*(delta1_1*sin2a +
						delta1_2*cos2a -
						delta2_1*sin2a_k -
						delta2_2*cos2a_k) +
    pow(alpha,2)/(pow(alpha,2)-pow(beta,2))*(delta1_1*sin2a +
					     delta1_2*cos2a -
					     delta1_3*cos2a*sin2a -
					     eps1_1*sin2a*sin2a -
					     eps1_2*(cos2a*cos2a+2*cos2a*sin2a) +
					     eps2_1*sin2a_k*sin2a_k +
					     eps2_2*(cos2a_k*cos2a_k+2*sin2a_k*cos2a_k) -
					     delta2_1*sin2a_k -
					     delta2_2*cos2a_k +
					     delta2_3*cos2a_k*sin2a_k) +
    2*pow(beta/alpha,2)*(gamma2_55*cos2a_k +
			 gamma2_44*sin2a_k -
			 gamma1_55*cos2a -
			 gamma1_44*sin2a);
  
  

  S5_C = -pow(beta,2)/(pow(alpha,2)-pow(beta,2))*(delta1_1*sin2a +
						  delta1_2*cos2a -
						  delta1_3*cos2a*sin2a -
						  eps1_1*sin2a*sin2a -
						  eps1_2*(cos2a*cos2a+2*cos2a*sin2a) +
						  eps2_1*sin2a_k*sin2a_k +
						  eps2_2*(cos2a_k*cos2a_k+2*sin2a_k*cos2a_k) -
						  delta2_1*sin2a_k -
						  delta2_2*cos2a_k +
						  delta2_3*cos2a_k*sin2a_k);
  
  S = pow(alpha,2)/(4*(pow(alpha,2)-pow(beta,2)))*((delta2_2-delta2_1)*sin(2*(azim-kappa)) +
						   (delta1_1-delta1_2)*sin(2*azim));

  CS_C = alpha*beta/(4*(pow(alpha,2)-pow(beta,2)))*((delta2_1-delta2_2)*sin(2*(azim-kappa)) +
						    (delta1_2-delta1_1)*sin(2*azim)) +
    beta/alpha*((gamma2_44-gamma2_55)*sin(2*(azim-kappa)) -
		(gamma1_44-gamma1_55)*sin(2*azim));
  
  CS3_C = alpha*beta/(2*(pow(alpha,2)-pow(beta,2)))*((delta1_1-delta1_2-delta1_3*(cos2a-sin2a) +
						      2*(eps1_2-eps1_1)*sin2a)*sin(azim)*cos(azim) -
						     (delta2_1-delta2_2-delta2_3*(cos2a_k-sin2a_k) +
						      2*(eps2_2-eps2_1)*sin2a_k)*sin(azim-kappa)*cos(azim-kappa));
  
  S3 = pow(alpha,2)/(2*(pow(alpha,2)-pow(beta,2)))*((-delta1_1+delta1_2+delta1_3*(cos2a-sin2a) -
					       2*(eps1_2-eps1_1)*sin2a)*sin(azim)*cos(azim) +
					      (delta2_1-delta2_2-delta2_3*(cos2a_k-sin2a_k) +
					       2*(eps2_2-eps2_1)*sin2a_k)*sin(azim-kappa)*cos(azim-kappa));


  *sv= CS * cos(ang) * sin(ang)
    + S_C * sin(ang)
    + CS3 * cos(ang)*pow(sin(ang),3)
    + (S3_C + S_C*0.5*pow(beta/alpha,2)) * pow(sin(ang),3)
    + (S5_C + S3_C*0.5*pow(beta/alpha,2)) * pow(sin(ang),5)
    + S5_C*0.5*pow(beta/alpha,2) * pow(sin(ang),7);
  

  *sh= S * sin(ang)
    + CS_C * cos(ang)*sin(ang)
    + (CS3_C + CS_C*0.5*pow(beta/alpha,2)) * cos(ang)*pow(sin(ang),3)
    + S3 * pow(sin(ang),3)
    + CS3_C*0.5*pow(beta/alpha,2) * cos(ang)*pow(sin(ang),5);
    
  /*  fprintf(stderr,"vp=%f vs=%f \n",alpha,beta); 
  fprintf(stderr,"ang=%f azim=%f k=%f sv=%f S=%f CSCj=%f S3=%f S3CCj=%f \n",ang*180/PI,azim*180/PI,kappa*180/PI,*sv,S,CS_C,S3,CS3_C);
  
  
    fprintf(stderr,"ang=%f azim=%f k=%f sv=%f CS=%f SCj=%f CS3=%f S3Cj=%f S5Cj=%f \n",ang*180/PI,azim*180/PI,kappa*180/PI,*sv,CS,S_C,CS3,S3_C,S5_C);  
  */

  /*    fprintf(stderr,"IN_ANG=%f AZIM=%f Da=%f \t Dr=%f \t Db=%f \t a/b=%f \n",ang*180./PI,azim*180./PI,Da_a,Dr_r,Db_b,beta/alpha);
   */
  

  /* finnaly, complete the Rps1, Rps2 */
  
  /*  fprintf(stderr,"ang=%f azim=%f S3/S1=%f \n",ang*180/PI,azim*180/PI,(CS3+S3_C+S_C*0.5*pow(beta/alpha,2)-0.5*CS)/(CS+S_C));  
   */

  *rps1 = (*sv)*(*cphi) + (*sh)*(*sphi);
  *rps2 = -(*sv)*(*sphi) + (*sh)*(*cphi);
  
  return(err);
}

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



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


     /* The subroutine evaluates first order (exactly) and second order (approximately)
	Rpp coefficients for a given incidence angle and azimuths. The second order evaluation
	is based on the numerical modeling and comparing of the exact Rpp 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 Rpp usually resulted in a global increase of accuracy of the first-order Rpp, 
	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 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]).

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

        	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 azimuths 0deg and 90 deg - global model
	            float angle[4];    ... inc. angle, two azimuths and kappa for which the quantities
		                           above are evaluated
	          }, 
 	 where  
         iso[0]-iso[4] ... isotropic part of Rpp for the following inc. angles (due to velocity
	                   and density contrasts): 15deg, 20deg, 25deg, 30deg and 35deg. 
		 	   If a change is required, seek the paragraph "ISO part" in this subroutine.
        upper[0],upper[1] ... Rpp coefficients for azimuths 0deg and 90deg, respectively, computed
                              for the ANISO/ISO configuration of the halfspaces.
        lower[0],lower[1] ... Rpp coefficients for azimuths 0deg and 90deg, respectively, computed
                              for the ISO/ANISO configuration of the halfspaces.
        global[0]-[3] ... Rpp coefficients for ANISO/ANISO configuration evaluated in the following
                          azimuths, respectively: 0deg, 90deg, lkappa, kkappa.
        angle[0]-[3] ... inc. angle (err_ang), 0deg, 90deg, kappa.

      P_err_2nd_order = 0 ... usual termination
                      =-1 ... unusual termination (wrong input)
       
     */
{

  extern float vp1,vp2,vs1,vs2,rho1,rho2;
  extern float eps1_1,eps1_2,delta1_1,delta1_2,delta1_3,gamma1_44;
  extern float eps2_1,eps2_2,delta2_1,delta2_2,delta2_3,gamma2_44;
  
  float Da_a, Dz_z, DG_G, Db_b, Dr_r, alpha, beta, azim;
  float Delta1_1, Delta1_2, Delta2_1, Delta2_2, Gamma1_44, Gamma2_44;  
  float D1_1, D1_2, D2_1, D2_2, GS_1, GS_2;
  float ISO_2nd;
  float d11_plus=0., d11_min=0., d12_plus=0., d12_min=0., gs1_plus=0., gs1_min=0.;
  float d21_plus=0., d21_min=0., d22_plus=0., d22_min=0., gs2_plus=0., gs2_min=0.;
  float a_plus=0., a_min=0, b_plus=0., b_min=0., r_plus=0., r_min=0. ;

  float Fst0_0, Fst0_90, Fst90_0, Fst90_90, FstK_0, FstK_90, FstKK_0, FstKK_90;
  float Fst0_kappa, Fst90_kappa, FstK_kappa, FstKK_kappa;
  float Snd0_0, Snd0_90, Snd90_0, Snd90_90, SndK_0, SndK_90, SndKK_0, SndKK_90;
  float Snd0_kappa, Snd90_kappa, SndK_kappa, SndKK_kappa;  
  float err_kappa=0.0, lkappa= 30.*PI/180, kkappa= 60.*PI/180, D;
  
  int N=0;  
    
  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);

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

  /* some useful quantities in advance */

  Da_a=(vp2-vp1)/(0.5*(vp2+vp1));
  Db_b=(vs2-vs1)/(0.5*(vs2+vs1));
  Dz_z=(vp2*rho2 - vp1*rho1)/(vp2*rho2 + vp1*rho1);
  DG_G=(pow(vs2,2)*rho2-pow(vs1,2)*rho1)/(0.5*(pow(vs2,2)*rho2+pow(vs1,2)*rho1));
  Dr_r=(rho2-rho1)/(0.5*(rho2+rho1));
  alpha=0.5*(vp2+vp1);
  beta=0.5*(vs2+vs1);
    
  if(index == 1) err_kappa=0.;
  if(index == 2) err_kappa=PI/2.;
  if((index != 1) && (index != 2)) return(-1);

  if(delta1_1) d11_plus=(fabs(delta1_1)+delta1_1)/(2*fabs(delta1_1));  
  if(delta1_1) d11_min=(fabs(delta1_1)-delta1_1)/(2*fabs(delta1_1));
  if(delta1_2) d12_plus=(fabs(delta1_2)+delta1_2)/(2*fabs(delta1_2));
  if(delta1_2) d12_min=(fabs(delta1_2)-delta1_2)/(2*fabs(delta1_2));
  if(gamma1_44) gs1_plus=(fabs(gamma1_44)+gamma1_44)/(2*fabs(gamma1_44));
  if(gamma1_44) gs1_min=(fabs(gamma1_44)-gamma1_44)/(2*fabs(gamma1_44));

  if(delta2_1) d21_plus=(fabs(delta2_1)+delta2_1)/(2*fabs(delta2_1));  
  if(delta2_1) d21_min=(fabs(delta2_1)-delta2_1)/(2*fabs(delta2_1));
  if(delta2_2) d22_plus=(fabs(delta2_2)+delta2_2)/(2*fabs(delta2_2));
  if(delta2_2) d22_min=(fabs(delta2_2)-delta2_2)/(2*fabs(delta2_2));
  if(gamma2_44) gs2_plus=(fabs(gamma2_44)+gamma2_44)/(2*fabs(gamma2_44));
  if(gamma2_44) gs2_min=(fabs(gamma2_44)-gamma2_44)/(2*fabs(gamma2_44));

  if(Da_a) a_plus=(fabs(Da_a)+Da_a)/(2*fabs(Da_a));
  if(Da_a) a_min=(fabs(Da_a)-Da_a)/(2*fabs(Da_a));
  if(Db_b) b_plus=(fabs(Db_b)+Db_b)/(2*fabs(Db_b));
  if(Db_b) b_min=(fabs(Db_b)-Db_b)/(2*fabs(Db_b));
  if(Dr_r) r_plus=(fabs(Dr_r)+Dr_r)/(2*fabs(Dr_r));
  if(Dr_r) r_min=(fabs(Dr_r)-Dr_r)/(2*fabs(Dr_r));


  /* the angles used for the anisotropic accuracy evaluation */

  (*rp_1st).angle[0]=(*rp_2nd).angle[0] = err_ang;        /* incidence angle */
  (*rp_1st).angle[1]=(*rp_2nd).angle[1] = 0.*PI/180;      /* first azimuth */
  (*rp_1st).angle[2]=(*rp_2nd).angle[2] = 90.*PI/180;     /* second azimuth */ 
  (*rp_1st).angle[3]=(*rp_2nd).angle[3] = true_kappa;      /* kappa */



  /* ISO part first */
  
  for(N = 0; N < 5; N++)
    {
      (*rp_1st).iso[N]=P_term(Da_a, Dz_z, DG_G, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 
			      0., 0., alpha, beta, (15+5*N)*PI/180, 0., 0.);
      (*rp_2nd).iso[N]=Iso_exact(0, vp1, vs1, rho1, vp2, vs2, rho2, (15+5*N)*PI/180);
    }

  /* get the ISO semi-2nd order for the inc. angle err_ang */
  ISO_2nd = Iso_exact(0, vp1, vs1, rho1, vp2, vs2, rho2, err_ang) - 
            P_term(Da_a, Dz_z, DG_G, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 
		   0., 0., alpha, beta, err_ang, 0., 0.);


  /* UPPER halfspace contribution */

  Delta1_1=(-0.45*d11_plus  -0.50*d11_min)*pow(delta1_1,2) + 
    (-0.22*d11_plus*gs1_plus+0*d11_plus*gs1_min-0.93*d11_min*gs1_plus-0.33*d11_min*gs1_min)*delta1_1*gamma1_44 + 
    (0.74*d11_plus*a_plus + 0.59*d11_plus*a_min + 2.56*d11_min*a_plus + 0.31*d11_min*a_min)*delta1_1*Da_a + 
    (0.19*d11_plus*b_plus + 0*d11_plus*b_min + 1.27*d11_min*b_plus + 0.31*d11_min*b_min)*delta1_1*Db_b + 
    (0*d11_plus*r_plus + 0.18*d11_plus*r_min + 1.27*d11_min*r_plus + 0.95*d11_min*r_min)*delta1_1*Dr_r;

  D1_1=Delta1_1;
  
  Delta1_2=(-0.45*d12_plus -0.50*d12_min)*pow(delta1_2,2) +
    (0.74*d12_plus*a_plus + 0.59*d12_plus*a_min + 2.56*d12_min*a_plus + 0.31*d12_min*a_min)*delta1_2*Da_a + 
    (0.19*d12_plus*b_plus + 0*d12_plus*b_min + 1.27*d12_min*b_plus + 0.31*d12_min*b_min)*delta1_2*Db_b + 
    (0*d12_plus*r_plus + 0.18*d12_plus*r_min + 1.27*d12_min*r_plus + 0.95*d12_min*r_min)*delta1_2*Dr_r;
							  

  D1_2=Delta1_2;  

  Delta2_1=0.;
  Delta2_2=0.;

  Gamma1_44=((1./3.)*gs1_plus + 0.55*gs1_min)*pow(gamma1_44,2) +
    (0.92*gs1_plus*a_plus + 0.9*gs1_plus*a_min + 1.27*gs1_min*a_plus + 0.98*gs1_min*a_min)*gamma1_44*Da_a +
    (-1.74*gs1_plus*b_plus -2.1*gs1_plus*b_min -1.84*gs1_min*b_plus -1.91*gs1_min*b_min)*gamma1_44*Db_b +
    (-1.02*gs1_plus*r_plus -1.04*gs1_plus*r_min -0.9*gs1_min*r_plus -0.9*gs1_min*r_min)*gamma1_44*Dr_r;
								
  GS_1=Gamma1_44;