Implemented (2,2) mode of potential tidal love number for stars with …

…a radiative envelope
TeamCOMPAS · Apr 4, 2024 · 8341cf0 · 8341cf0
1 parent 14ef044
commit 8341cf0
Show file tree

Hide file tree

Showing 4 changed files with 17 additions and 9 deletions.
diff --git a/src/BaseStar.h b/src/BaseStar.h
@@ -175,6 +175,8 @@ class BaseStar {
     virtual void            CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams) { }                                                                                      // Default is NO-OP
     virtual void            CalculateGBParams()                                                                 { CalculateGBParams(m_Mass0, m_GBParams); }                         // Use class member variables
 
+    virtual double          CalculateImK22Tidal(const double p_Omega)                                           { return 0.0; }                                                     // Default is 0.0
+
             void            CalculateLambdas()                                                                  { CalculateLambdas(m_Mass - m_CoreMass); }                          // Use class member variables
             void            CalculateLambdas(const double p_EnvMass);
 

diff --git a/src/GiantBranch.cpp b/src/GiantBranch.cpp
@@ -1135,22 +1135,28 @@ double GiantBranch::CalculateMomentOfInertia() const {
  *
  * Zahn, 1977, Eq. (5.5) , with the value of E_2 coming from Kushnir et al., 2017, by comparing Eq. (8) to Eq. (1)
  *
- *
- * double CalculateImKlmTidal(const double p_Omega)
+ * double CalculateImK22Tidal(const double p_Omega)
  *
  * @param   [IN]    p_Omega                     Orbital angular frequency (1/yr)
  * @return                                      Imaginary component of pontential tidal love number (unitless)
  */
-double GiantBranch::CalculateImKlmTidal(const double p_Omega) {
-    double beta2Dynamical = 1;
+double GiantBranch::CalculateImK22Tidal(const double p_Omega) {
+    double beta2Dynamical = 1.0;
     double rhoFactorDynamcial = 0.1;
     double radiusAU = m_Radius * RSOL_TO_AU;
+    double coreRadiusAU = CalculateCoreRadius() * RSOL_TO_AU;
+    double coreRadius_over_radius = coreRadiusAU / radiusAU;
+    double coreRadius_over_radius_3 = coreRadius_over_radius * coreRadius_over_radius * coreRadius_over_radius;
+    double coreRadius_over_radius_9 = coreRadius_over_radius_3 * coreRadius_over_radius_3 * coreRadius_over_radius_3;
+    double mass_over_coreMass = m_Mass / m_CoreMass;
 
-    double E2Dynamical = (2.0 / 3.0) * PPOW((CalculateCoreRadius() / radiusAU), 9) * PPOW((m_Mass / m_CoreMass), (4.0 / 3.0)) * beta2Dynamical * rhoFactorDynamcial;
 
-    double s22 = 2.0 * (p_Omega - m_Omega) * std::sqrt(radiusAU * radiusAU * radiusAU / G_AU_Msol_yr / m_Mass);
+    double E2Dynamical = (2.0 / 3.0) * coreRadius_over_radius_9 * mass_over_coreMass * std::cbrt(mass_over_coreMass) * beta2Dynamical * rhoFactorDynamcial;
 
-    double k22Dynamical = E2Dynamical * PPOW(s22, 8.0/3.0);
+    double s22 = 2.0 * (p_Omega - Omega()) * std::sqrt(radiusAU * radiusAU * radiusAU / G_AU_Msol_yr / m_Mass);
+    double s22_4_3 = s22 * std::cbrt(s22);
+    double s22_8_3 = s22_4_3 * s22_4_3;
+    double k22Dynamical = E2Dynamical * s22_8_3;
 
     return k22Dynamical;
 }

diff --git a/src/GiantBranch.h b/src/GiantBranch.h
@@ -52,7 +52,7 @@ class GiantBranch: virtual public BaseStar, public MainSequence {
             void            CalculateGBParams()                                                             { CalculateGBParams(m_Mass0, m_GBParams); }                         // Use class member variables
 
     static  double          CalculateHRateConstant_Static(const double p_Mass);
-            double          CalculateImKlmTidal(const double p_Omega);
+            double          CalculateImK22Tidal(const double p_Omega);
     virtual double          CalculateInitialSupernovaMass() const                                           { return m_Mass; }                                                  // Use class member variables
 
             double          CalculateLifetimeToHeIgnition(const double p_Mass, const double p_Tinf1_FGB, const double p_Tinf2_FGB);

diff --git a/src/Star.h b/src/Star.h
@@ -162,7 +162,7 @@ class Star {
     double          CalculateConvectiveEnvelopeMass()                                                               { return m_Star->CalculateConvectiveEnvelopeMass(); }
 
     double          CalculateEddyTurnoverTimescale()                                                                { return m_Star->CalculateEddyTurnoverTimescale(); }
-    double          CalculateImKlmTidal(const double p_Omega)                                                       { return m_Star->CalculateImKlmTidal(p_Omega); }
+    double          CalculateImK22Tidal(const double p_Omega)                                                       { return m_Star->CalculateImK22Tidal(p_Omega); }
     void            CalculateLambdas()                                                                              { m_Star->CalculateLambdas(); }
     void            CalculateLambdas(const double p_EnvMass)                                                        { m_Star->CalculateLambdas(p_EnvMass); }