better utils/attenuation/conversion_from_Qkappa_Qmu_to_Qp_Qs_from_Dah…

…len_Tromp_959_960.f90

utils/attenuation/conversion_from_Qkappa_Qmu_to_Qp_Qs_from_Dahlen_Tromp_959_960.f90

@@ -1,13 +1,20 @@
 
   program conversion
 
-! Dimitri Komatitsch, CNRS Marseille, France, October 2015
+! Dimitri Komatitsch, CNRS Marseille, France, October 2015 and March 2017
 
 ! see formulas 9.59 and 9.60 in the book of Dahlen and Tromp, 1998
-! (in that book, P is called alpha and S is called beta)
+! (in that book, P is called alpha and S is called beta).
+! The formulas in Dahlen and Tromp, 1998 are for the 3D case,
+! in the 2D plane strain case the 4/3 coefficient must be changed to 1,
+! see more details in file doc/Qkappa_Qmu_versus_Qp_Qs_relationship_in_2D_plane_strain.pdf
 
   implicit none
 
+! coefficient for the 3D case is 4/3, for the 2D plane strain case it is 1,
+! see more details in file doc/Qkappa_Qmu_versus_Qp_Qs_relationship_in_2D_plane_strain.pdf
+  double precision, parameter :: coefficient = 4.d0/3.d0   !! 1.d0
+
   double precision :: Qkappa,Qmu,Qp,Qs,inverse_of_Qp,cp,cs
 
 !!! this is for the Carcione et al. 1988 example
@@ -34,7 +41,7 @@ program conversion
   Qs = Qmu
 
 ! for Qp the formula is more complex
-  inverse_of_Qp = (1.d0 - (4.d0/3.d0)*(cs**2)/(cp**2))/Qkappa + (4.d0/3.d0)*(cs**2)/(cp**2)/Qmu
+  inverse_of_Qp = (1.d0 - coefficient*(cs**2)/(cp**2))/Qkappa + coefficient*(cs**2)/(cp**2)/Qmu
   Qp = 1.d0/inverse_of_Qp
 
 ! print the result