lincoeff.c

su 的源代码库

		    float iso[5];     
	            float upper[2];   
	            float lower[2];   
	            float global[4]; 
	            float angle[4]; 
		                    
	          }, 
 	 where  
         iso[0]-iso[4] ... Rpsv: isotropic part 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.
			   Rpsh: iso[n]=0 for all n.
        upper[0],upper[1] ... Rpsv: coefficient for azimuths 0deg and 90deg, respectively, computed
                              for the ANISO/ISO configuration of the halfspaces. If required,
			      changes can be made in UPPER -> SV paragraph.
			      Rpsh: the same as above. However, only upper[0] is activated for 45deg.
			      If required, changes can be made in UPPER -> SH paragraph.
        lower[0],lower[1] ... Rpsv: coefficients for azimuths 0deg and 90deg, respectively, computed
                              for the ISO/ANISO configuration of the halfspaces. If required,
			      changes can be made in LOWER -> SV paragraph. 
			      Rpsh: the same as above. However, only lower[0] is activated for 45deg.
			      If required, changes can be made in LOWER -> SH paragraph.
        global[0]-[3] ... Rpsv: coefficients for ANISO/ANISO configuration evaluated in the following
                          azimuths, respectively: 0deg, 90deg, lkappa, kkappa. If required,
			  changes can be made using lkappa and kkappa variables.  
			  Rpsh: the same meaning as above. However, only the global[0]-global[2] are
			  activated for the angles 20deg, 45deg and 70deg. If required,
			  changes can be made in GLOBAL -> SH -> err_kappa=PI/2. paragraph. 
        angle[0]-[3] ... inc. angle (err_ang), 0deg, 90deg, kappa.

      S_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, Db_b, Dr_r, DG_G, alpha, beta;   
  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, SHD1_1, SHD1_2, SHD2_1, SHD2_2, SHGS_1, SHGS_2;  
  float ISO_2nd, azim;               
  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 SHFst20_0, SHFst20_90, SHFst45_0, SHFst45_90, SHFst70_0, SHFst70_90;
  float SHFst20_kappa, SHFst45_kappa, SHFst70_kappa;
  float SHSnd20_0, SHSnd20_90, SHSnd45_0, SHSnd45_90, SHSnd70_0, SHSnd70_90;
  float SHSnd20_kappa, SHSnd45_kappa, SHSnd70_kappa; 
  float err_kappa=0.0, lkappa=30.*PI/180, kkappa=60.*PI/180, D;

  int N=0;  

  float SV_term(float DG_G, float Dr_r, 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 kappa);
  
  float SH_term(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 kappa);

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

  /* some quantities in advance */

  Da_a=(vp2-vp1)/(0.5*(vp2+vp1));
  Db_b=(vs2-vs1)/(0.5*(vs2+vs1));
  Dr_r=(rho2-rho1)/(0.5*(rho2+rho1));  
  DG_G=(pow(vs2,2)*rho2-pow(vs1,2)*rho1)/(0.5*(pow(vs2,2)*rho2+pow(vs1,2)*rho1));
  alpha=0.5*(vp2+vp1);
  beta=0.5*(vs2+vs1);
  /*kkappa=PI/2.-true_kappa;*/

  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 */

  (*rsv_1st).angle[0]=(*rsv_2nd).angle[0] = err_ang;        /* incidence angle */
  (*rsv_1st).angle[1]=(*rsv_2nd).angle[1] = 0.;             /* first azimuth */
  (*rsv_1st).angle[2]=(*rsv_2nd).angle[2] = 90.;            /* second azimuth */ 
  (*rsv_1st).angle[3]=(*rsv_2nd).angle[3] = true_kappa;      /* kappa */

  (*rsh_1st).angle[0]=(*rsh_2nd).angle[0] = err_ang;        /* incidence angle */
  (*rsh_1st).angle[1]=(*rsh_2nd).angle[1] = 0.;             /* first azimuth */
  (*rsh_1st).angle[2]=(*rsh_2nd).angle[2] = 90.;            /* second azimuth */ 
  (*rsh_1st).angle[3]=(*rsh_2nd).angle[3] = true_kappa;      /* kappa */


  /* ISO part first */

  /* SV */

  for(N = 0; N < 5; N++)
    {
      (*rsv_1st).iso[N]=SV_term(DG_G, Dr_r, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
				0., 0., alpha, beta, (15+5*N)*PI/180, 0., 0.);
      (*rsv_2nd).iso[N]=Iso_exact(1, vp1, vs1, rho1, vp2, vs2, rho2, (15+5*N)*PI/180);
    }

  /* evaluate the ISO 2nd order for the inc. angle err_ang */

  ISO_2nd = Iso_exact(1, vp1, vs1, rho1, vp2, vs2, rho2, err_ang) -
            SV_term(DG_G, Dr_r, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
		    0., 0., alpha, beta, err_ang, 0., 0.);
  
  /* SH */
  (*rsh_1st).iso[0]=(*rsh_2nd).iso[0]=0.;
  (*rsh_1st).iso[1]=(*rsh_2nd).iso[1]=0.;
  (*rsh_1st).iso[2]=(*rsh_2nd).iso[2]=0.;
  (*rsh_1st).iso[3]=(*rsh_2nd).iso[3]=0.;
  (*rsh_1st).iso[4]=(*rsh_2nd).iso[4]=0.;


  /* UPPER halfspace contribution */

  /* SV */
  Delta1_1=(-1*d11_plus -2.5*d11_min)*pow(delta1_1,2) + 
    (-0.22*d11_plus*gs1_plus+0*d11_plus*gs1_min-0.5*d11_min*gs1_plus-0.39*d11_min*gs1_min)*delta1_1*gamma1_44 + 
    (1.28*d11_plus*a_plus +0.75*d11_plus*a_min +1.93*d11_min*a_plus +2.93*d11_min*a_min)*delta1_1*Da_a + 
    (0.64*d11_plus*b_plus +0.23*d11_plus*b_min +0.18*d11_min*b_plus +1.87*d11_min*b_min)*delta1_1*Db_b + 
    (1.15*d11_plus*r_plus +0.66*d11_plus*r_min +3.51*d11_min*r_plus +2.69*d11_min*r_min)*delta1_1*Dr_r;
  

  D1_1=Delta1_1;
  
  Delta1_2=(-1*d12_plus -2.5*d12_min)*pow(delta1_2,2) +
    (1.28*d12_plus*a_plus +0.75*d12_plus*a_min +1.93*d12_min*a_plus +2.93*d12_min*a_min)*delta1_2*Da_a + 
    (0.64*d12_plus*b_plus +0.23*d12_plus*b_min +0.18*d12_min*b_plus +1.87*d12_min*b_min)*delta1_2*Db_b + 
    (1.15*d12_plus*r_plus +0.66*d12_plus*r_min +3.51*d12_min*r_plus +2.69*d12_min*r_min)*delta1_2*Dr_r;
						       

  D1_2=Delta1_2;  

  Delta2_1=0.;
  Delta2_2=0.;

  Gamma1_44=(-0.6*gs1_plus -1.1*gs1_min)*pow(gamma1_44,2) +
    (-0.385*gs1_plus*a_plus -0.36*gs1_plus*a_min -0.55*gs1_min*a_plus -0.225*gs1_min*a_min)*gamma1_44*Da_a +
    (0*gs1_plus*b_plus +0*gs1_plus*b_min +0*gs1_min*b_plus -0.3*gs1_min*b_min)*gamma1_44*Db_b +
    (-0.33*gs1_plus*r_plus -0.32*gs1_plus*r_min -0.64*gs1_min*r_plus -0.57*gs1_min*r_min)*gamma1_44*Dr_r;
							   
  GS_1=Gamma1_44;
  
  Gamma2_44=0.;

  if(index == 1)
    {     
      azim=0.;
      (*rsv_1st).upper[0]=SV_term(DG_G, Dr_r, delta1_1, delta1_2, delta1_3, 0., 0., 0.,
				  eps1_1, eps1_2, 0., 0.,
				  gamma1_44, 0., alpha, beta, err_ang, azim, 0.); 
      (*rsv_2nd).upper[0]=SV_term(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;  
      (*rsv_1st).upper[1]=SV_term(DG_G, Dr_r, delta1_1, delta1_2, delta1_3, 0., 0., 0.,
				  eps1_1, eps1_2, 0., 0.,
				  gamma1_44, 0., alpha, beta, err_ang, azim, 0.);
      (*rsv_2nd).upper[1]=SV_term(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;
      
    }
  
  /* SH */
  Delta1_1=(-0.8*d11_plus -2.0*d11_min)*pow(delta1_1,2) + 
    (-0.22*d11_plus*gs1_plus+0*d11_plus*gs1_min-0.4*d11_min*gs1_plus-0.22*d11_min*gs1_min)*delta1_1*gamma1_44 + 
    (0*d11_plus*d12_plus + 0*d11_plus*d12_min + 0*d11_min*d12_plus + 0*d11_min*d12_min)*delta1_1*delta1_2 +
    (1.32*d11_plus*a_plus +1.09*d11_plus*a_min +2.93*d11_min*a_plus +2.4*d11_min*a_min)*delta1_1*Da_a + 
    (0.95*d11_plus*b_plus +0.42*d11_plus*b_min +0.73*d11_min*b_plus +0.75*d11_min*b_min)*delta1_1*Db_b + 
    (1.67*d11_plus*r_plus +1.68*d11_plus*r_min +3.71*d11_min*r_plus +3.12*d11_min*r_min)*delta1_1*Dr_r;
  

  SHD1_1=Delta1_1;
  
  Delta1_2=(-0.8*d12_plus -2.0*d12_min)*pow(delta1_2,2) +
    (0.22*d12_plus*gs1_plus-0.33*d12_plus*gs1_min+0.22*d12_min*gs1_plus+0.32*d12_min*gs1_min)*delta1_2*gamma1_44 +
    (1.32*d12_plus*a_plus +1.09*d12_plus*a_min +2.93*d12_min*a_plus +2.4*d12_min*a_min)*delta1_2*Da_a + 
    (0.95*d12_plus*b_plus +0.42*d12_plus*b_min +0.73*d12_min*b_plus +0.75*d12_min*b_min)*delta1_2*Db_b + 
    (1.67*d12_plus*r_plus +1.68*d12_plus*r_min +3.71*d12_min*r_plus +3.12*d12_min*r_min)*delta1_2*Dr_r;
  

  SHD1_2=Delta1_2;  

  Delta2_1=0.;
  Delta2_2=0.;

  Gamma1_44=(-0.45*gs1_plus -0.78*gs1_min)*pow(gamma1_44,2) +
    (-0.17*gs1_plus*a_plus -0.16*gs1_plus*a_min -0.28*gs1_min*a_plus -0.09*gs1_min*a_min)*gamma1_44*Da_a +
    (-0.27*gs1_plus*b_plus -0.45*gs1_plus*b_min -0.55*gs1_min*b_plus -0.52*gs1_min*b_min)*gamma1_44*Db_b +
    (-0.38*gs1_plus*r_plus -0.48*gs1_plus*r_min -0.57*gs1_min*r_plus -0.61*gs1_min*r_min)*gamma1_44*Dr_r;
							     

  SHGS_1=Gamma1_44;
  
  Gamma2_44=0.;


  if(index == 1 )
    {
      azim=45.*PI/180;  
      (*rsh_1st).upper[0]=SH_term(delta1_1, delta1_2, delta1_3, 0., 0., 0.,
				  eps1_1, eps1_2, 0., 0.,
				  gamma1_44, 0., alpha, beta, err_ang, azim, 0.);
      (*rsh_2nd).upper[0]=SH_term(Delta1_1, Delta1_2, 0., Delta2_1, Delta2_2, 0.,
				  0., 0., 0., 0.,
				  Gamma1_44, Gamma2_44, alpha, beta, err_ang, azim, 0.);  
      
      (*rsh_2nd).upper[1]=0.0;
      (*rsh_1st).upper[1]=0.0;
    }
  
 
  /* LOWER halfspace contribution */

  /* SV */
  Delta1_1=0.;
  Delta1_2=0.;

  Delta2_1=(-0.72*d21_plus -1.89*d21_min)*pow(delta2_1,2) + 
    (0.61*d21_plus*gs2_plus+1.26*d21_plus*gs2_min+0.83*d21_min*gs2_plus+2.28*d21_min*gs2_min)*delta2_1*gamma2_44+
    (-0.28*d21_plus*a_plus +0*d21_plus*a_min -1.1*d21_min*a_plus -1.05*d21_min*a_min)*delta2_1*Da_a + 
    (0.82*d21_plus*b_plus +0.83*d21_plus*b_min +2.93*d21_min*b_plus +1.73*d21_min*b_min)*delta2_1*Db_b + 
    (-0.07*d21_plus*r_plus +0.1*d21_plus*r_min +0*d21_min*r_plus +0*d21_min*r_min)*delta2_1*Dr_r;
  

  D2_1=Delta2_1;
  
  Delta2_2=(-0.72*d22_plus -1.89*d22_min)*pow(delta2_2,2) +
    (-0.28*d22_plus*a_plus +0*d22_plus*a_min -1.1*d22_min*a_plus -1.05*d22_min*a_min)*delta2_2*Da_a + 
    (0.82*d22_plus*b_plus +0.83*d22_plus*b_min +2.93*d22_min*b_plus +1.73*d22_min*b_min)*delta2_2*Db_b + 
    (-0.07*d22_plus*r_plus +0.1*d22_plus*r_min +0*d22_min*r_plus + 0*d22_min*r_min)*delta2_2*Dr_r;
							   

  D2_2=Delta2_2;
  
  Gamma1_44=0.;
  Gamma2_44=(-0.28*gs2_plus -(2./3.)*gs2_min)*pow(gamma2_44,2) +
    (-0.37*gs2_plus*a_plus -0.3*gs2_plus*a_min -0.55*gs2_min*a_plus -0.45*gs2_min*a_min)*gamma2_44*Da_a +
    (1.1*gs2_plus*b_plus +0.69*gs2_plus*b_min +0.92*gs2_min*b_plus +0.6*gs2_min*b_min)*gamma2_44*Db_b +
    (0*gs2_plus*r_plus +0*gs2_plus*r_min +0.06*gs2_min*r_plus +0.20*gs2_min*r_min)*gamma2_44*Dr_r;
								

  GS_2=Gamma2_44;

  if(index == 1)
    {
      azim=0.;
      (*rsv_1st).lower[0]=SV_term(DG_G, Dr_r, 0., 0., 0., delta2_1, delta2_2, delta2_3,
				  0., 0., eps2_1, eps2_2,
				  0., gamma2_44, alpha, beta, err_ang, azim, 0.);
      (*rsv_2nd).lower[0]=SV_term(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;  
      (*rsv_1st).lower[1]=SV_term(DG_G, Dr_r, 0., 0., 0., delta2_1, delta2_2, delta2_3,
				  0., 0., eps2_1, eps2_2,
				  0., gamma2_44, alpha, beta, err_ang, azim, 0.);
      (*rsv_2nd).lower[1]=SV_term(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;
      
    }
  
  /* SH */
  Delta1_1=0.;
  Delta1_2=0.;

  Delta2_1=(-0.61*d21_plus -1.56*d21_min)*pow(delta2_1,2) + 
    (0.28*d21_plus*gs2_plus+1.1*d21_plus*gs2_min+0.44*d21_min*gs2_plus+1.56*d21_min*gs2_min)*delta2_1*gamma2_44 + 
    (0*d21_plus*d22_plus + 0*d21_plus*d22_min + 0*d21_min*d22_plus + 0*d