lincoeff.c

su 的源代码库

  
  Gamma2_44=0.;

  if(index == 1)
    {      
      azim=0.;
      (*rp_1st).upper[0]=P_term(Da_a, Dz_z, DG_G, delta1_1, delta1_2, delta1_3, 0., 0., 0.,
				eps1_1, eps1_2, 0., 0.,
				gamma1_44, 0., alpha, beta, err_ang, azim, 0.);
      (*rp_2nd).upper[0]=P_term(0., 0., 0., Delta1_1, Delta1_2, 0., Delta2_1, Delta2_2, 0., 0., 0., 0., 0., 
				Gamma1_44, Gamma2_44, alpha, beta, err_ang, azim, 0.) +
	                 ISO_2nd;
      
      azim=90.*PI/180;  
      (*rp_1st).upper[1]=P_term(Da_a, Dz_z, DG_G, delta1_1, delta1_2, delta1_3, 0., 0., 0.,
				eps1_1, eps1_2, 0., 0.,
				gamma1_44, 0., alpha, beta, err_ang, azim, 0.);
      (*rp_2nd).upper[1]=P_term(0., 0., 0., Delta1_1, Delta1_2, 0., Delta2_1, Delta2_2, 0., 0., 0., 0., 0.,
				Gamma1_44, Gamma2_44, alpha, beta, err_ang, azim, 0.) +
                         ISO_2nd;
      
    }
  


  /* LOWER halfspace contribution */

  Delta1_1=0.;
  Delta1_2=0.;

  Delta2_1=(-0.44*d21_plus -0.89*d21_min)*pow(delta2_1,2) + 
    (-(1./3.)*d21_plus*gs2_plus+0*d21_plus*gs2_min-0.65*d21_min*gs2_plus+0*d21_min*gs2_min)*delta2_1*gamma2_44+ 
    (1.83*d21_plus*a_plus +1.43*d21_plus*a_min +2.56*d21_min*a_plus +1.35*d21_min*a_min)*delta2_1*Da_a + 
    (-0.37*d21_plus*b_plus -0.59*d21_plus*b_min -0.55*d21_min*b_plus -0.07*d21_min*b_min)*delta2_1*Db_b + 
    (0*d21_plus*r_plus +0*d21_plus*r_min -0.13*d21_min*r_plus +0.10*d21_min*r_min)*delta2_1*Dr_r;


  D2_1=Delta2_1;
  
  Delta2_2=(-0.44*d22_plus -0.89*d22_min)*pow(delta2_2,2) +
    (1.83*d22_plus*a_plus +1.43*d22_plus*a_min +2.56*d22_min*a_plus +1.35*d22_min*a_min)*delta2_2*Da_a + 
    (-0.37*d22_plus*b_plus -0.59*d22_plus*b_min -0.55*d22_min*b_plus -0.07*d22_min*b_min)*delta2_2*Db_b + 
    (0*d22_plus*r_plus + 0*d22_plus*r_min -0.13*d22_min*r_plus + 0.10*d22_min*r_min)*delta2_2*Dr_r;
							   
  D2_2=Delta2_2;
  
  Gamma1_44=0.;
  Gamma2_44=(-0.37*gs2_plus -0.5*gs2_min)*pow(gamma2_44,2) +
    (0*gs2_plus*a_plus +1.05*gs2_plus*a_min +0.83*gs2_min*a_plus +0.83*gs2_min*a_min)*gamma2_44*Da_a +
    (-0.46*gs2_plus*b_plus +0.23*gs2_plus*b_min +0.24*gs2_min*b_plus +0.38*gs2_min*b_min)*gamma2_44*Db_b +
    (-0.13*gs2_plus*r_plus +0.07*gs2_plus*r_min +0*gs2_min*r_plus +0.24*gs2_min*r_min)*gamma2_44*Dr_r;
							    
  GS_2=Gamma2_44;
  
  if(index == 1)
    {     
      azim=0.;
      (*rp_1st).lower[0]=P_term(Da_a, Dz_z, DG_G, 0., 0., 0., delta2_1, delta2_2, delta2_3,
				0., 0., eps2_1, eps2_2,
				0., gamma2_44, alpha, beta, err_ang, 0., azim);
      (*rp_2nd).lower[0]=P_term(0., 0., 0., Delta1_1, Delta1_2, 0., Delta2_1, Delta2_2, 0.,
				0., 0., 0., 0.,
				Gamma1_44, Gamma2_44, alpha, beta, err_ang, 0., azim) +
	                 ISO_2nd;
      
      
      azim=90.*PI/180;  
      (*rp_1st).lower[1]=P_term(Da_a, Dz_z, DG_G, 0., 0., 0., delta2_1, delta2_2, delta2_3,
				0., 0., eps2_1, eps2_2,
				0., gamma2_44, alpha, beta, err_ang, 0., azim);
      (*rp_2nd).lower[1]=P_term(0., 0., 0., Delta1_1, Delta1_2, 0., Delta2_1, Delta2_2, 0.,
				0., 0., 0., 0.,
				Gamma1_44, Gamma2_44, alpha, beta, err_ang, 0., azim) +
	                 ISO_2nd;
      
      
    }
  
  /* GLOBAL contribution*/
  
  /* The following terms are derived on the assumption that Rpp coefficients for
     two identicaly anisotropic halfspaces are zero (with 0 velocity contrats)
  */
  
  Delta1_1=D1_1 +
           (0.45*d11_plus*d21_plus -0.45*d11_plus*d21_min -0.50*d11_min*d21_plus +0.50*d11_min*d21_min)*
                                                             delta1_1*delta2_1*pow(cos(err_kappa),2) + 
           (0.45*d11_plus*d22_plus -0.45*d11_plus*d22_min -0.50*d11_min*d22_plus +0.50*d11_min*d22_min)* 
                                                             delta1_1*delta2_2*pow(sin(err_kappa),2) +
           (0.22*d11_plus*gs2_plus +0.*d11_plus*gs2_min +0.93*d11_min*gs2_plus +0.33*d11_min*gs2_min)*
	                                                     delta1_1*gamma2_44*pow(cos(err_kappa),2);
  
  Delta1_2=D1_2 + 
           (0.45*d12_plus*d21_plus -0.45*d12_plus*d21_min -0.50*d12_min*d21_plus +0.50*d12_min*d21_min)*
                                                             delta1_2*delta2_1*pow(sin(err_kappa),2) + 
           (0.45*d12_plus*d22_plus -0.45*d12_plus*d22_min -0.50*d12_min*d22_plus +0.50*d12_min*d22_min)* 
                                                             delta1_2*delta2_2*pow(cos(err_kappa),2) +
           (0.22*d12_plus*gs2_plus +0.*d12_plus*gs2_min +0.93*d12_min*gs2_plus +0.33*d12_min*gs2_min)*
                                                             delta1_2*gamma2_44*pow(sin(err_kappa),2); 

  Delta2_1=D2_1 + 
           (0.44*d21_plus*d11_plus -0.44*d21_plus*d11_min -0.89*d21_min*d11_plus +0.89*d21_min*d11_min)*
                                                             delta2_1*delta1_1*pow(cos(err_kappa),2) + 
           (0.44*d21_plus*d12_plus -0.44*d21_plus*d12_min -0.89*d21_min*d12_plus +0.89*d21_min*d12_min)* 
                                                             delta2_1*delta1_2*pow(sin(err_kappa),2) +
           ((1./3.)*d21_plus*gs1_plus + 0.*d21_plus*gs1_min +0.65*d21_min*gs1_plus +0.*d21_min*gs1_min)*
                                                             delta2_1*gamma1_44*pow(cos(err_kappa),2);
  

  Delta2_2=D2_2 + 
           (0.44*d22_plus*d11_plus -0.44*d22_plus*d11_min -0.89*d22_min*d11_plus +0.89*d22_min*d11_min)*
                                                             delta2_2*delta1_1*pow(sin(err_kappa),2) + 
           (0.44*d22_plus*d12_plus -0.44*d22_plus*d12_min -0.89*d22_min*d12_plus +0.89*d22_min*d12_min)* 
                                                             delta2_2*delta1_2*pow(cos(err_kappa),2) +
           ((1./3.)*d22_plus*gs1_plus + 0.*d22_plus*gs1_min +0.65*d22_min*gs1_plus +0.*d22_min*gs1_min)*
                                                             delta2_2*gamma1_44*pow(sin(err_kappa),2);


  Gamma1_44=GS_1 + 
            (0.0*gs1_plus*d21_plus +0.0*gs1_plus*d21_min +0.0*gs1_min*d21_plus +0.0*gs1_min*d21_min)*
                                                             gamma1_44*delta2_1*pow(cos(err_kappa),2) +
            (0.0*gs1_plus*d22_plus +0.0*gs1_plus*d22_min +0.0*gs1_min*d22_plus +0.0*gs1_min*d22_min)*
                                                             gamma1_44*delta2_2*pow(sin(err_kappa),2) +
            (-(1./3.)*gs1_plus*gs2_plus -0.75*gs1_plus*gs2_min -0.95*gs1_min*gs2_plus -0.55*gs1_min*gs2_min)*
                                                             gamma1_44*gamma2_44*pow(cos(err_kappa),2);
  
  Gamma2_44=GS_2 +
            (0.37*gs2_plus*gs1_plus +0.*gs2_plus*gs1_min +0.*gs2_min*gs1_plus +0.5*gs2_min*gs1_min)*
                                                            gamma2_44*gamma1_44*pow(cos(err_kappa),2);

  /* The following terms are based on error minimization for 2-parameter model only.
     Notice that such 2nd order terms do not result in Rpp=0 for two identically
     anisotropic halfspaces (probably not convenient to use).
  
  
   Delta1_1=D1_1 +
           (0.45*d11_plus*d21_plus -0.45*d11_plus*d21_min -0.50*d11_min*d21_plus +0.50*d11_min*d21_min)*
                                                             delta1_1*delta2_1*pow(cos(err_kappa),2) + 
           (0.45*d11_plus*d22_plus -0.45*d11_plus*d22_min -0.50*d11_min*d22_plus +0.50*d11_min*d22_min)* 
                                                             delta1_1*delta2_2*pow(sin(err_kappa),2) +
           (0*d11_plus*gs2_plus +0.1*d11_plus*gs2_min +0.8*d11_min*gs2_plus +1.1*d11_min*gs2_min)*
	                                                     delta1_1*gamma2_44*pow(cos(err_kappa),2);
  
  Delta1_2=D1_2 + 
           (0.45*d12_plus*d21_plus -0.45*d12_plus*d21_min -0.50*d12_min*d21_plus +0.50*d12_min*d21_min)*
                                                             delta1_2*delta2_1*pow(sin(err_kappa),2) + 
           (0.45*d12_plus*d22_plus -0.45*d12_plus*d22_min -0.50*d12_min*d22_plus +0.50*d12_min*d22_min)* 
                                                             delta1_2*delta2_2*pow(cos(err_kappa),2) +
           (0*d12_plus*gs2_plus +0.1*d12_plus*gs2_min +0.8*d12_min*gs2_plus +1.1*d12_min*gs2_min)*
                                                             delta1_2*gamma2_44*pow(sin(err_kappa),2); 

  Delta2_1=D2_1 + 
           (0.44*d21_plus*d11_plus -0.44*d21_plus*d11_min -0.89*d21_min*d11_plus +0.89*d21_min*d11_min)*
                                                             delta2_1*delta1_1*pow(cos(err_kappa),2) + 
           (0.44*d21_plus*d12_plus -0.44*d21_plus*d12_min -0.89*d21_min*d12_plus +0.89*d21_min*d12_min)* 
                                                             delta2_1*delta1_2*pow(sin(err_kappa),2);
  

  Delta2_2=D2_2 + 
           (0.44*d22_plus*d11_plus -0.44*d22_plus*d11_min -0.89*d22_min*d11_plus +0.89*d22_min*d11_min)*
                                                             delta2_2*delta1_1*pow(sin(err_kappa),2) + 
           (0.44*d22_plus*d12_plus -0.44*d22_plus*d12_min -0.89*d22_min*d12_plus +0.89*d22_min*d12_min)* 
                                                             delta2_2*delta1_2*pow(cos(err_kappa),2);


  Gamma1_44=GS_1 + 
            (0.1*gs1_plus*d21_plus +0.22*gs1_plus*d21_min +0.2*gs1_min*d21_plus +0.55*gs1_min*d21_min)*
                                                             gamma1_44*delta2_1*pow(cos(err_kappa),2) +
            (0.1*gs1_plus*d22_plus +0.22*gs1_plus*d22_min +0.2*gs1_min*d22_plus +0.55*gs1_min*d22_min)*
                                                             gamma1_44*delta2_2*pow(sin(err_kappa),2) +
            (-(1./3.)*gs1_plus*gs2_plus -0.75*gs1_plus*gs2_min -0.95*gs1_min*gs2_plus -0.55*gs1_min*gs2_min)*
                                                             gamma1_44*gamma2_44*pow(cos(err_kappa),2);
  
  Gamma2_44=GS_2 +
            (0.37*gs2_plus*gs1_plus +0*gs2_plus*gs1_min +0*gs2_min*gs1_plus +0.5*gs2_min*gs1_min)*
                                                            gamma2_44*gamma1_44*pow(cos(err_kappa),2);
  
  */
  
  if(index == 1)

    /* here, err_kappa=0: evaluate the error for alligned halfspaces */

    {
      azim=0.;
      (*rp_1st).global[0]=P_term(Da_a, Dz_z, DG_G, delta1_1, delta1_2, delta1_3, delta2_1, delta2_2, delta2_3,
				 eps1_1, eps1_2, eps2_1, eps2_2,
				 gamma1_44, gamma2_44, alpha, beta, err_ang, azim, azim-err_kappa);
      (*rp_2nd).global[0]=P_term(0., 0., 0., Delta1_1, Delta1_2, 0., Delta2_1, Delta2_2, 0.,
				 0., 0., 0., 0.,
				 Gamma1_44, Gamma2_44, alpha, beta, err_ang, azim, azim-err_kappa) +
	                  ISO_2nd;
      

      azim=90.*PI/180;
      (*rp_1st).global[1]=P_term(Da_a, Dz_z, DG_G, delta1_1, delta1_2, delta1_3, delta2_1, delta2_2, delta2_3,
				 eps1_1, eps1_2, eps2_1, eps2_2,
				 gamma1_44, gamma2_44, alpha, beta, err_ang, azim, azim-err_kappa);
      (*rp_2nd).global[1]=P_term(0., 0., 0., Delta1_1, Delta1_2, 0., Delta2_1, Delta2_2, 0.,
				 0., 0., 0., 0.,
				 Gamma1_44, Gamma2_44, alpha, beta, err_ang, azim, azim-err_kappa) +
	                  ISO_2nd;
      
      azim=lkappa;
      (*rp_1st).global[2]=P_term(Da_a, Dz_z, DG_G, delta1_1, delta1_2, delta1_3, delta2_1, delta2_2, delta2_3,
				 eps1_1, eps1_2, eps2_1, eps2_2,
				 gamma1_44, gamma2_44, alpha, beta, err_ang, azim, azim-err_kappa);
      (*rp_2nd).global[2]=P_term(0., 0., 0., Delta1_1, Delta1_2, 0., Delta2_1, Delta2_2, 0.,
				 0., 0., 0., 0.,
				 Gamma1_44, Gamma2_44, alpha, beta, err_ang, azim, azim-err_kappa) +
	                  ISO_2nd;
      
      
      azim=kkappa;
      (*rp_1st).global[3]=P_term(Da_a, Dz_z, DG_G, delta1_1, delta1_2, delta1_3, delta2_1, delta2_2, delta2_3,
				 eps1_1, eps1_2, eps2_1, eps2_2,
				 gamma1_44, gamma2_44, alpha, beta, err_ang, azim, azim-err_kappa);
      (*rp_2nd).global[3]=P_term(0., 0., 0., Delta1_1, Delta1_2, 0., Delta2_1, Delta2_2, 0.,
				 0., 0., 0., 0.,
				 Gamma1_44, Gamma2_44, alpha, beta, err_ang, azim, azim-err_kappa) +
	                  ISO_2nd;
      
    }
  
  Fst0_0=(*rp_1st).global[0];
  Snd0_0=(*rp_2nd).global[0];
  Fst90_0=(*rp_1st).global[1];
  Snd90_0=(*rp_2nd).global[1];
  FstK_0=(*rp_1st).global[2];
  SndK_0=(*rp_2nd).global[2];
  FstKK_0=(*rp_1st).global[3];
  SndKK_0=(*rp_2nd).global[3];

  if(index == 2)

    /* here, err_kappa = 90.: evaluate the error for the halfspaces rotated by 90 deg */

    {    
      azim=0.;  
      Fst0_90=P_term(Da_a, Dz_z, DG_G, delta1_1, delta1_2, delta1_3, delta2_1, delta2_2, delta2_3,
				 eps1_1, eps1_2, eps2_1, eps2_2,
				 gamma1_44, gamma2_44, alpha, beta, err_ang, azim, azim-err_kappa);
      Snd0_90=P_term(0., 0., 0., Delta1_1, Delta1_2, 0., Delta2_1, Delta2_2, 0.,
				 0., 0., 0., 0.,
				 Gamma1_44, Gamma2_44, alpha, beta, err_ang, azim, azim-err_kappa) +
	     ISO_2nd;
      Fst0_kappa=P_term(Da_a, Dz_z, DG_G, delta1_1, delta1_2, delta1_3, delta2_1, delta2_2, delta2_3,
				 eps1_1, eps1_2, eps2_1, eps2_2,
				 gamma1_44, gamma2_44, alpha, beta, err_ang, azim, azim-true_kappa);
      
      
      azim=90.*PI/180;  
      Fst90_90=P_term(Da_a, Dz_z, DG_G, delta1_1, delta1_2, delta1_3, delta2_1, delta2_2, delta2_3,
				 eps1_1, eps1_2, eps2_1, eps2_2,
				 gamma1_44, gamma2_44, alpha, beta, err_ang, azim, azim-err_kappa);  
      Snd90_90=P_term(0., 0., 0., Delta1_1, Delta1_2, 0., Delta2_1, Delta2_2, 0.,
				 0., 0., 0., 0.,
				 Gamma1_44, Gamma2_44, alpha, beta, err_ang, azim, azim-err_kappa) +
	       ISO_2nd;
      Fst90_kappa=P_term(Da_a, Dz_z, DG_G, delta1_1, delta1_2, delta1_3, delta2_1, delta2_2, delta2_3,
				 eps1_1, eps1_2, eps2_1, eps2_2,
				 gamma1_44, gamma2_44, alpha, beta, err_ang, azim, azim-true_kappa);


      azim=lkappa;
      FstK_90=P_term(Da_a, Dz_z, DG_G, delta1_1, delta1_2, delta1_3, delta2_1, delta2_2, delta2_3,
				 eps1_1, eps1_2, eps2_1, eps2_2,
				 gamma1_44, gamma2_44, alpha, beta, err_ang, azim, azim-err_kappa);
      SndK_90=P_term(0., 0., 0., Delta1_1, Delta1_2, 0., Delta2_1, Delta2_2, 0.,
				 0., 0., 0., 0.,
				 Gamma1_44, Gamma2_44, alpha, beta, err_ang, azim, azim-err_kappa) +
	      ISO_2nd;
      FstK_kappa=P_term(Da_a, Dz_z, DG_G, delta1_1, delta1_2, delta1_3, delta2_1, delta2_2, delta2_3,
				 eps1_1, eps1_2, eps2_1, eps2_2,
				 gamma1_44, gamma2_44, alpha, beta, err_ang, azim, azim-true_kappa);
       

      azim=kkappa;
      FstKK_90=P_term(Da_a, Dz_z, DG_G, delta1_1, delta1_2, delta1_3, delta2_1, delta2_2, delta2_3,
				 eps1_1, eps1_2, eps2_1, eps2_2,
				 gamma1_44, gamma2_44, alpha, beta, err_ang, azim, azim-err_kappa);
      SndKK_90=P_term(0., 0., 0., Delta1_1, Delta1_2, 0., Delta2_1, Delta2_2, 0.,
				 0., 0., 0., 0.,
				 Gamma1_44, Gamma2_44, alpha, beta, err_ang, azim, azim-err_kappa) +
	       ISO_2nd;
      FstKK_kappa=P_term(Da_a, Dz_z, DG_G, delta1_1, delta1_2, delta1_3, delta2_1, delta2_2, delta2_3,
				 eps1_1, eps1_2, eps2_1, eps2_2,
				 gamma1_44, gamma2_44, alpha, beta, err_ang, azim, azim-true_kappa);
       
      /* interpolate the error for the halfspaces rotated by the angle true_kappa */

      if(fabs(D=Fst0_90-Fst0_0) < LC_TINY){
	Snd0_kappa = 0.5*(Snd0_0 + Snd0_90);
      }
      else {
	Snd0_kappa = Snd0_0 + ((Snd0_90-Snd0_0)/D)*(Fst0_kappa-Fst0_0);
      }

      if(fabs(D=Fst90_90-Fst90_0) < LC_TINY) {
	Snd90_kappa = 0.5*(Snd90_0 + Snd90_90);
      }
      else {
	Snd90_kappa = Snd90_0 + ((Snd90_90-Snd90_0)/D)*(Fst90_kappa-Fst90_0);
      }

      if(fabs(D=FstK_90-FstK_0) < LC_TINY) {
	SndK_kappa =0.5*(SndK_90 + SndK_0);
      }
      else {
	SndK_kappa = SndK_0 + ((SndK_90-SndK_0)/D)*(FstK_kappa-FstK_0);
      }
      
      if(fabs(D=FstKK_90-FstKK_0) < LC_TINY) {
	SndKK_kappa = 0.5*(SndKK_90 + SndKK_0);
      }
      else {