linrort.c

su 的源代码库

	err(" incidence angle from <0;90) only; error in fangle or langle or dangle.  Program terminated!");
	return (-1);
	}
  if (fazim > lazim)
	{
	err(" an error in fazim or lazim.  Program terminated!");
	return (-1);
	}
	

  /* dumping some info into the file info.out  */

  fprintf(out_inf_p," ************************ \n *** INFORMATION FILE *** \n ************************ \n \n INCIDENCE HALFSPACE 1: %s \n",hspace1);
  if (HSP1==0) fprintf(out_inf_p," Vp1=%6.3f \t Vs1=%6.3f \t rho1=%6.3f",vp1,vs1,rho1);
  if (HSP1==1) fprintf(out_inf_p," Vp1(Vp33)=%6.3f \t Vs1(Vs44)=%6.3f \t rho1=%6.3f eps1=%6.3f \t\t delta1=%6.3f ",vp1,vs1,rho1,eps1,delta1);

  if (HSP1==2) fprintf(out_inf_p," Vp1(Vp33)=%6.3f \t Vs1(Vs44)=%6.3f \t rho1=%6.3f eps1_v=%6.3f \t\t delta1_v=%6.3f \t gamma1_v=%6.3f",vp1,vs1,rho1,eps1_v,delta1_v,gamma1_v);
  if (HSP1==3) fprintf(out_inf_p," Vp1(Vp33)=%6.3f \t Vs1(Vs44)=%6.3f \t rho1=%6.3f eps1_1=%6.3f \t\t delta1_1=%6.3f \t gamma1_1=%6.3f eps1_2=%6.3f \t\t delta1_2=%6.3f \t gamma1_2=%6.3f delta1_3=%6.3f",vp1,vs1,rho1,eps1_1,delta1_1,gamma1_1,eps1_2,delta1_2,gamma1_2,delta1_3);

  fprintf(out_inf_p,"\n\n REFLECTING HALFSPACE 2: %s \n", hspace2);
  if (HSP2==0) fprintf(out_inf_p," Vp2=%6.3f \t Vs2=%6.3f \t rho2=%6.3f",vp2,vs2,rho2);
  if (HSP2==1) fprintf(out_inf_p," Vp2(Vp33)=%6.3f \t Vs2(Vs44)=%6.3f \t rho2=%6.3f eps2=%6.3f \t\t delta2=%6.3f ",vp2,vs2,rho2,eps2,delta2);
  if (HSP2==2) fprintf(out_inf_p," Vp2(Vp33)=%6.3f \t Vs2(Vs44)=%6.3f \t rho2=%6.3f eps2_v=%6.3f \t\t delta2_v=%6.3f \t gamma2_v=%6.3f",vp2,vs2,rho2,eps2_v,delta2_v,gamma2_v);
  if (HSP2==3) fprintf(out_inf_p," Vp2(Vp33)=%6.3f \t Vs2(Vs44)=%6.3f \t rho2=%6.3f eps2_1=%6.3f \t\t delta2_1=%6.3f \t gamma2_1=%6.3f eps2_2=%6.3f \t\t delta2_2=%6.3f \t gamma2_2=%6.3f delta2_3=%6.3f",vp2,vs2,rho2,eps2_1,delta2_1,gamma2_1,eps2_2,delta2_2,gamma2_2,delta2_3); 

  fprintf(out_inf_p,"\n \n INCIDENCE ANGLE AND AZIMUTH RANGE:\n");
  if (strcmp(a_file,"-1")) 
	fprintf(out_inf_p," Incidence angles and azimuths are specified in the file \"%s\", \n kappa=%7.3f \n",a_file,kappa);
 else
	fprintf(out_inf_p," fangle=%7.3f \t langle=%7.3f \t dangle=%7.3f \t fazim=%7.3f	\t lazim=%7.3f	\t dazim=%7.3f \t rotation of the halfspace 2:  kappa=%7.3f \n", fangle,langle,dangle,fazim,lazim,dazim,kappa); 


  /* initialization of anisotropic parameters */

  /* ISO media */
  if (HSP1==0)
	{
	eps1_1=0.;
	eps1_2=0.;
	delta1_1=0.;
	delta1_2=0.;
	delta1_3=0.;
	gamma1_44=0.;
	gamma1_55=0.;
	gamma1_1=0.;	
	}
  if (HSP2==0)
	{
	eps2_1=0.;
	eps2_2=0.;
	delta2_1=0.;
	delta2_2=0.;
	delta2_3=0.;
	gamma2_44=0.;
	gamma2_55=0.;
	}


  /* VTI media */
  if (HSP1==1)
	{
	eps1_1=eps1;
	eps1_2=eps1;
	delta1_1=delta1;
	delta1_2=delta1;
	delta1_3=0.;
	gamma1_44=0.;
	gamma1_55=0.;
	gamma1_1=gamma1;
	}
  if (HSP2==1)
	{
	eps2_1=eps2;
	eps2_2=eps2;
	delta2_1=delta2;
	delta2_2=delta2;
	delta2_3=0.;
	gamma2_44=0.;
	gamma2_55=0.;
	}
  

  /* HTI media */
  if (HSP1==2)
	{
	eps1_1=0.;
	eps1_2=eps1_v;
	delta1_1=0.;
	delta1_2=delta1_v;
	delta1_3=delta1_v-2*eps1_v;
	/*  exact delta_3 gives usually worse results: 
	delta1_3=(delta1_v-2*eps1_v)/(1+2*eps1_v); */
	gamma1_1=0.;
	if (HTIback1==44)
	{
	  gamma1_44=0.;
	  gamma1_55=gamma1_v;
	  count+=1; 
	  back=1;	  
	}
	if (HTIback1==55)
	{
	  gamma1_44=-gamma1_v/(1+2*gamma1_v);
	  gamma1_55=0.;
	  vs1=vs1*sqrt(1+2*gamma1_v);	  
	}
	}
  if (HSP2==2)
	{
	eps2_1=0.;
	eps2_2=eps2_v;
	delta2_1=0.;
	delta2_2=delta2_v;
	delta2_3=delta2_v-2*eps2_v;
	/*exact delta_3 gives usually worse results: 
	delta2_3=(delta2_v-2*eps2_v)/(1+2*eps2_v); */
	if (HTIback2==44)
	{
	  gamma2_44=0.;
	  gamma2_55=gamma2_v;
	  count+=1;	  
	  back=1;
	}
	if (HTIback2==55)
	{
	  gamma2_44=-gamma2_v/(1+2*gamma2_v);
	  gamma2_55=0.;
	  vs2=vs2*sqrt(1+2*gamma2_v);
	}
	}
  
  /* ORT media */
  if (HSP1==3)
	{
	if (ORTback1==44)
	{
	  gamma1_44=0.;
	  gamma1_55=gamma1_2-(gamma1_1*(1+2*gamma1_2)/(1+2*gamma1_1));
	  count+=1;	
	  back=1;	  
	}
	if (ORTback1==55)
	{
	  gamma1_44=gamma1_1-(gamma1_2*(1+2*gamma1_1)/(1+2*gamma1_2));
	  gamma1_55=0.;
	  vs1=vs1*sqrt((1+2*gamma1_2)/(1+2*gamma1_1));
	}
	}
  if(HSP2==3)
	{
	if (ORTback2==44)
	{
	  gamma2_44=0.;
	  gamma2_55=gamma2_2-(gamma2_1*(1+2*gamma2_2)/(1+2*gamma2_1));
	  count+=1;
	  back=1;
	}
	if (ORTback2==55)
	{
	  gamma2_44=gamma2_1-(gamma2_2*(1+2*gamma2_1)/(1+2*gamma2_2));
	  gamma2_55=0.;
	  vs2=vs2*sqrt((1+2*gamma2_2)/(1+2*gamma2_1));
	}
	}
  

  /* some more dumping into the file info.out  */
  
  fprintf(out_inf_p,"\n INPUT PARAMETERS FOR THE GENERAL Rpp AND Rps SUBROUTINES: \n");
  fprintf(out_inf_p," alpha1=%6.3f \t\t beta1=%6.5f \t\t rho1=%6.3f eps1_1=%6.3f \t\t delta1_1=%6.3f eps1_2=%6.3f \t\t delta1_2=%6.3f delta1_3=%6.3f gamma1_44=%6.5f \t gamma1_55=%6.5f \t gamma1_1=%6.5f \n",vp1,vs1,rho1,eps1_1, delta1_1,eps1_2,delta1_2,delta1_3,gamma1_44,gamma1_55,gamma1_1);

  fprintf(out_inf_p,"\n alpha2=%6.3f \t\t beta2=%6.5f \t\t rho2=%6.3f eps2_1=%6.3f \t\t delta2_1=%6.3f eps2_2=%6.3f \t\t delta2_2=%6.3f delta2_3=%6.3f gamma2_44=%6.5f \t gamma2_55=%6.5f \n",vp2,vs2,rho2,eps2_1, delta2_1,eps2_2,delta2_2,delta2_3,gamma2_44,gamma2_55);

  fprintf(out_inf_p,"\n ERROR MESSAGES: \n");
  
  /* some useful combinations for inversion */
  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));
  DZ_Z=(vp2*rho2-vp1*rho1)/(0.5*(vp2*rho2+vp1*rho1));
  alpha=0.5*(vp2+vp1);
  beta=0.5*(vs2+vs1);
  Db_b=(vs2-vs1)/(0.5*(vs2+vs1));

  PS2=0.5*((delta1_2-delta1_1+8*pow(beta/alpha,2)*gamma1_44) -
	(delta2_2-delta2_1+8*pow(beta/alpha,2)*gamma2_44)*cos(2*kappa*PI/180));	
  
		PSC=0.5*(delta2_2-delta2_1+8*pow(beta/alpha,2)*gamma2_44)*sin(2*kappa*PI/180);
  PABS=0.5*(-(delta2_2-delta2_1+8*pow(beta/alpha,2)*gamma2_44)*pow(sin(kappa*PI/180),2) +
		delta2_2-delta1_2+ Da_a - pow(2*beta/alpha,2)*DG_G);

  SS2=alpha/(2*(alpha+beta))*(delta1_2-delta1_1) + 2*beta/alpha*gamma1_44 - 
	(alpha/(2*(alpha+beta))*(delta2_2-delta2_1)+ 2*beta/alpha*gamma2_44)*cos(2*kappa*PI/180);
  SSC=(alpha/(2*(alpha+beta))*(delta2_2-delta2_1)+ 2*beta/alpha*gamma2_44)*sin(2*kappa*PI/180);
  /*  SABS=-(alpha/(2*(alpha+beta))*(delta2_2-delta2_1)+ 2*beta/alpha*gamma2_44)*pow(sin(kappa*PI/180),2) +
	alpha/(2*(alpha+beta))*(delta2_2-delta1_2)-(alpha+2*beta)/(2*alpha)*(rho2-rho1)/(0.5*(rho2+rho1))-
	2*beta/alpha*Db_b;
  */
  SABS=-(alpha/(2*(alpha+beta))*(delta2_2-delta2_1)+ 2*beta/alpha*gamma2_44)*pow(sin(kappa*PI/180),2) +
	alpha/(2*(alpha+beta))*(delta2_2-delta1_2)-0.5*(rho2-rho1)/(0.5*(rho2+rho1))-
	beta/alpha*DG_G;

  /*  some auxiliary outputs, may be useful

fprintf(stderr,"\n b_a=%f \t Da_a=%f \t Dr_r=%f \t Db_b=%f \n DG_G=%f \t 0.5*DZ_Z=%f \n PS2=%f \t PSC=%f \t PABS=%f \n SS2=%f \t SSC=%f \t SABS=%f \n", beta/alpha,Da_a,(rho2-rho1)/(0.5*(rho2+rho1)),Db_b,
DG_G,0.5*DZ_Z,PS2,PSC,PABS,SS2,SSC,SABS);
  */

  /* finally, the evaluation of the coefficients */
  
  if (!strcmp(a_file,"-1"))
	{
	for( azim = fazim; azim <= lazim; azim += dazim){
	 for( ang = fangle; ang <= langle; ang += dangle)
		{	  
		if((errp=lincoef_Rp(ang*PI/180, azim*PI/180, kappa*PI/180, &RPP, &Rp_1st, &Rp_2nd, count))!=0)
		fprintf(out_inf_p,"\n !!! WARNING: incidence ANGLE=%4.2fdeg probably approaches, or may have exceeded, the crytical angle ANGL_CR=%4.2fdeg (this is an approximate value only).\n", ang, errp);
		errs=lincoef_Rs(ang*PI/180, azim*PI/180, kappa*PI/180, &RPS1, &RPS2, &SV, &SH, &cosPHI, &sinPHI, HSP1,
			&Rsv_1st, &Rsv_2nd, &Rsh_1st, &Rsh_2nd, count);		
		fprintf(out_P_p," %f \t %f \t %f \n",ang,azim,RPP);
 count+=1;
		if (errs >= 0)
		{
		 fprintf(out_S_p," %f \t %f \t %f \t %f \t %f \t %f \n",ang,azim,RPS1,RPS2,cosPHI,sinPHI);
		 fprintf(out_SVSH_p," %f \t %f \t %f \t %f \t %f \t %f \n",ang,azim,SV,SH,cosPHI,sinPHI);	
		}
		else if (errs==-1)
		{
		fprintf(out_S_p," %f \t %f \t %f \t %f \t %f \t %f \n",ang,azim,RPS1,RPS2,cosPHI,sinPHI);
		fprintf(out_SVSH_p," %f \t %f \t %f \t %f \t %f \t %f \n",ang,azim,SV,SH,cosPHI,sinPHI);		
		fprintf(out_inf_p,"\n !!! WARNING: incidence ANGLE=%4.2fdeg is essentially at a singular point of S-wave propagation. However, the Rps computation is completed. \n", ang);
		}				
		else if (errs==-2)
		{
		 fprintf(out_inf_p,"\n !!! WARNING: incidence ANGLE=%4.2f is essentially at a singular point of S-wave propagation, the computation of Rps is terminated for this angle (ORT inc. halfspace).\n", ang);
		}
		}	 
	}
	}  
	else
	{
	while(fscanf(a_file_p,"%f %f",&ang,&azim)!=EOF)
	 {
		if (ang >= 90.) continue;
		if ((errp=lincoef_Rp(ang*PI/180, azim*PI/180, kappa*PI/180, &RPP, &Rp_1st, &Rp_2nd, count))!=0.)
		fprintf(out_inf_p,"\n !!! WARNING: incidence ANGLE=%4.2fdeg probably approaches, or may have exceeded, the crytical angle ANGL_CR=%4.2fdeg (this is an approximate value only).\n", ang, errp);
		errs=lincoef_Rs(ang*PI/180, azim*PI/180, kappa*PI/180, &RPS1, &RPS2, &SV, &SH, &cosPHI, &sinPHI, HSP1, &Rsv_1st, &Rsv_2nd, &Rsh_1st, &Rsh_2nd, count);
		fprintf(out_P_p," %f \t %f \t %f \n",ang,azim,RPP);
		count+=1;
		if (errs >= 0.)
		{
		fprintf(out_S_p," %f \t %f \t %f \t %f \t %f \t %f \n",ang,azim,RPS1,RPS2,cosPHI,sinPHI);
		fprintf(out_SVSH_p," %f \t %f \t %f \t %f \t %f \t %f \n",ang,azim,SV,SH,cosPHI,sinPHI);
		}
		else if (errs==-1)
		{
		fprintf(out_S_p," %f \t %f \t %f \t %f \t %f \t %f \n",ang,azim,RPS1,RPS2,cosPHI,sinPHI);
		fprintf(out_SVSH_p," %f \t %f \t %f \t %f \t %f \t %f \n",ang,azim,SV,SH,cosPHI,sinPHI);		
		fprintf(out_inf_p,"\n !!! WARNING: incidence ANGLE=%4.2fdeg is essentially at a singular point of S-wave propagation. However, the Rps computation is completed. \n", ang);
		}				
		else if (errs==-2)
		{
		fprintf(out_inf_p,"\n !!! WARNING: incidence ANGLE=%4.2fdeg azimuth AZIM=%4.2fdeg is essentially at a singular point of S-wave propagation, the computation of Rps is terminated for this angle (ORT inc. halfspace). \n", ang,azim);
		}		
	 }
	}
  

  /* Almost done - filling up the out_Error file */

  fprintf(out_Error_p,"************************************************************* \n");
  fprintf(out_Error_p,"******************* FAST TEST OF ACCURACY ******************* \n");
  fprintf(out_Error_p,"************************************************************* \n");
  fprintf(out_Error_p,"(The test is based on the comparison of 1st and 2nd order terms.) \n");
  fprintf(out_Error_p,"(The test applies for the S-wave projection \"back=0\" only.) \n");


  if(back == 0){	


	/* P COEFFICIENT FIRST */

	fprintf(out_Error_p,"\n Rpp COEFFICIENT:  \n");
	fprintf(out_Error_p,"****************** \n");

	/* P isotropic */
	fprintf(out_Error_p,"\n A) Isotropic component (due to velocity and density contrasts):  \n");
	fprintf(out_Error_p,"\n Inc. Angle: \t \t	15deg	20deg	25deg	30deg	35deg \n");
	fprintf(out_Error_p," ----------------------------------------------------------------------------- \n");
	fprintf(out_Error_p," Exact: \t \t %7.4f	%7.4f	%7.4f	%7.4f	%7.4f \n",
		Rp_2nd.iso[0],Rp_2nd.iso[1],Rp_2nd.iso[2],Rp_2nd.iso[3], Rp_2nd.iso[4]);
	fprintf(out_Error_p," Approximation: \t %7.4f	%7.4f	%7.4f	%7.4f	%7.4f	\n",
		Rp_1st.iso[0],Rp_1st.iso[1],Rp_1st.iso[2],Rp_1st.iso[3], Rp_1st.iso[4]);
	fprintf(out_Error_p," ------------------------------------------------------------------------------ \n");

	Max=0.;  
	for(n=0; n < 5 ; n++){
	if(fabs(Rp_2nd.iso[n]) > 0.0001){
	Iso_e[n]=(Rp_1st.iso[n]/Rp_2nd.iso[n]-1)*100.;
	Max = fabs(Max) > fabs(Iso_e[n]) ? Max : Iso_e[n];
	}
	else{
	Iso_e[n]=999.9;
	}
	}