Refactored giant branch helpers into a single class

src/SSE/Landmarks/BaseOfGiantBranch.cpp

@@ -13,7 +13,7 @@
 #include "BaseOfGiantBranch.h"
 
 #include "Constants.h"
-#include "GiantBranchHelpers.h"
+#include "GiantBranchHelper.h"
 #include "Utilities.h"
 
 #include <Exceptions/PreconditionError.h>
@@ -136,7 +136,7 @@ void BaseOfGiantBranch::ComputeMetallicityDependents( Herd::Generic::Metallicity
   rB[ 7 ] = 9.301992e+00;
 
   // RBGB
-  m_ZDependents.m_pRGBComputer = std::make_unique< Herd::SSE::GiantBranchRadius >( i_Z );
+  m_ZDependents.m_pGBHelper = std::make_unique< Herd::SSE::GiantBranchHelper >( i_Z );
 }
 
 /**
@@ -188,7 +188,7 @@ Herd::Generic::Luminosity BaseOfGiantBranch::ComputeLuminosity( Herd::Generic::M
  */
 Herd::Generic::Radius BaseOfGiantBranch::ComputeRadius( Herd::Generic::Mass i_Mass ) const
 {
-  return m_ZDependents.m_pRGBComputer->Compute( i_Mass, ComputeLuminosity( i_Mass ) );
+  return m_ZDependents.m_pGBHelper->Radius( i_Mass, ComputeLuminosity( i_Mass ) );
 }
 
 }

src/SSE/Landmarks/BaseOfGiantBranch.h

@@ -25,7 +25,7 @@
 namespace Herd::SSE
 {
 
-class GiantBranchRadius;
+class GiantBranchHelper;
 
 /**
  * @brief Computations for characteristic properties at BGB
@@ -58,7 +58,7 @@ class BaseOfGiantBranch : public Herd::SSE::ILandmark
     std::array< double, 5 > m_TBGB; ///< \f$ t_{BGB} \f$ calculations
     std::array< double, 8 > m_LBGB;  ///< \f$ L_{BGB} \f$ calculations
 
-    std::unique_ptr< Herd::SSE::GiantBranchRadius > m_pRGBComputer; ///< Computes the giant branch radius
+    std::unique_ptr< Herd::SSE::GiantBranchHelper > m_pGBHelper; ///< Helper functions for the GB computations
   };
 
   MetallicityDependents m_ZDependents;  ///< Metallicity-dependent quantities

src/SSE/Landmarks/CMakeLists.txt

@@ -3,7 +3,7 @@ set(HEADER_LIST BaseOfGiantBranch.h
 								BaseOfGiantBranch.h
 								Constants.h
 								CriticalMassValues.h
-								GiantBranchHelpers.h
+								GiantBranchHelper.h
 								HeliumIgnition.h						
 								TerminalMainSequence.h
 								Utilities.h
@@ -13,7 +13,7 @@ set(HEADER_LIST BaseOfGiantBranch.h
 
 set(SOURCE_LIST BaseOfGiantBranch.cpp
 								CriticalMassValues.cpp
-								GiantBranchHelpers.cpp
+								GiantBranchHelper.cpp
 								HeliumIgnition.cpp
 								TerminalMainSequence.cpp
 								ZeroAgeMainSequence.cpp

src/SSE/Landmarks/GiantBranchHelper.cpp

@@ -0,0 +1,212 @@
+/**
+ * @file GiantBranchHelper.cpp
+ * @author Evren Imre
+ * @date 1 Sep 2021	
+ */
+/* This file is a part of HeRD, a stellar evolution library
+ * Copyright © 2021 Evren Imre 
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at https://mozilla.org/MPL/2.0/. */
+
+#include "GiantBranchHelper.h"
+
+#include <Generic/MathHelpers.h>
+#include <SSE/Landmarks/Constants.h>
+
+#include <algorithm>
+#include <cmath>
+#include <iterator>
+
+#include <range/v3/algorithm.hpp>
+
+namespace
+{
+
+// @formatter:off
+
+std::array< double, 24 > s_ZRGB { 9.960283e-01, 8.164393e-01, 2.383830e+00, 2.223436e+00, 8.638115e-01, 1.231572e-01,
+  2.561062e-01, 7.072646e-02, -5.444596e-02, -5.798167e-02, -1.349129e-02, 0.,
+  1.157338e+00, 1.467883e+00, 4.299661e+00, 3.130500e+00, 6.992080e-01, 1.640687e-02,
+  4.022765e-01, 3.050010e-01, 9.962137e-01, 7.914079e-01, 1.728098e-01, 0.
+};
+// @formatter:on
+}
+
+namespace Herd::SSE
+{
+
+/**
+ * @param i_Z Initial metallicity
+ * @pre \c i_Z>0
+ * @throws PreconditionError If preconditions are violated
+ */
+GiantBranchHelper::GiantBranchHelper( Herd::Generic::Metallicity i_Z )
+{
+  Herd::Generic::ThrowIfNotPositive( i_Z, "i_Z" );
+
+  ComputeMetallicityDependents( i_Z );
+}
+
+/**
+ * @remarks Without a user-defined destructor forward declaration and unique_ptr do not work together
+ */
+GiantBranchHelper::~GiantBranchHelper() = default;
+
+/**
+ * @param i_Mass Mass
+ * @param i_Luminosity Luminosity
+ * @return Giant branch radius
+ * @pre \c i_Mass>0
+ * @pre \c i_Luminosity>0
+ * @throws PreconditionError If preconditions are violated
+ */
+Herd::Generic::Radius GiantBranchHelper::Radius( Herd::Generic::Mass i_Mass, Herd::Generic::Luminosity i_Luminosity )
+{
+  Herd::Generic::ThrowIfNotPositive( i_Mass, "i_Mass" );
+  Herd::Generic::ThrowIfNotPositive( i_Luminosity, "i_Luminosity" );
+
+  ComputeA( i_Mass );
+
+  const auto& rB = m_ZDependents.m_RGB;
+  double a = m_MDependents.m_A.second;
+
+  double term1 = std::pow( i_Luminosity, rB[ 1 ] );
+  double term2 = Herd::Generic::BXhC( i_Luminosity, rB[ 0 ], rB[ 2 ] );
+  return Herd::Generic::Radius( a * ( term1 + term2 ) );
+}
+
+/**
+ * @param i_Mass Mass
+ * @param i_Luminosity Luminosity
+ * @param i_Metallicity Metallicity
+ * @return Core mass
+ * @pre \c i_Mass>0
+ * @pre \c i_Luminosity>0
+ * @throws PreconditionError If preconditions are violated
+ */
+Herd::Generic::Mass GiantBranchHelper::CoreMass( Herd::Generic::Mass i_Mass, Herd::Generic::Luminosity i_Luminosity )
+{
+  Herd::Generic::ThrowIfNotPositive( i_Mass, "i_Mass" );
+  Herd::Generic::ThrowIfNotPositive( i_Luminosity, "i_Luminosity" );
+
+  ComputeGB( i_Mass );
+
+  // Derived from Eq. 37 and Eq. 38
+  const auto& rMcGB = m_MDependents.m_McGB;
+  if( i_Luminosity < rMcGB.m_Lx )
+  {
+    return Herd::Generic::Mass( Herd::Generic::BXhC( i_Luminosity, 1.0 / rMcGB.m_D, 1.0 / rMcGB.m_P ) );
+  } else
+  {
+    return Herd::Generic::Mass( Herd::Generic::BXhC( i_Luminosity, 1.0 / rMcGB.m_B, 1.0 / rMcGB.m_Q ) );
+  }
+}
+
+/**
+ * @param i_Z Metallicity
+ */
+void GiantBranchHelper::ComputeMetallicityDependents( Herd::Generic::Metallicity i_Z )
+{
+  std::array< double, 6 > zetaPowers5;
+  double zeta = std::log10( i_Z / Herd::SSE::Constants::s_SolarMetallicityTout96 );
+  Herd::Generic::ComputePowers( zetaPowers5, zeta );
+
+  std::array< double, 3 > zetaPowers2;
+  ranges::cpp20::copy_n( zetaPowers5.begin(), 3, zetaPowers2.begin() );
+
+  double logZ = std::log10( i_Z );
+
+  // RGB calculations
+  auto& rB = m_ZDependents.m_RGB;
+  rB[ 0 ] = std::max( std::pow( 10., -4.6739 - 0.9394 * logZ ), -0.04167 + 55.67 * i_Z );
+  rB[ 0 ] = std::min( rB[ 0 ], 0.4771 - Herd::Generic::BXhC( i_Z, 9329.21, 2.94 ) );
+
+  rB[ 1 ] = std::min( 0.54, Herd::Generic::ComputeInnerProduct( { 0.397, 0.28826, 0.5293 }, zetaPowers2 ) );
+
+  rB[ 2 ] = std::pow( 10., std::max( -0.1451, -2.2794 - 1.5175 * logZ - 0.254 * logZ * logZ ) );
+  if( i_Z > 0.004 )
+  {
+    rB[ 2 ] = std::max( rB[ 2 ], Herd::Generic::ApBXhC( i_Z, 0.7307, 14265.1, 3.395 ) );
+  }
+
+  std::array< double, 4 > tempB;
+  Herd::Generic::MultiplyMatrixVector( tempB, s_ZRGB, zetaPowers5 );
+  ranges::cpp20::copy( tempB, std::next( rB.begin(), 3 ) ); // @suppress("Invalid arguments")
+
+  // D0
+  m_ZDependents.m_D0 = 5.37 + zeta * 0.135;
+}
+
+/**
+ * @param i_Mass Mass
+ */
+void GiantBranchHelper::ComputeA( Herd::Generic::Mass i_Mass )
+{
+  // Compute only if the mass value is different
+  auto& [ rKey, rA ] = m_MDependents.m_A;
+  if( rKey && *rKey == i_Mass )
+  {
+    return;
+  }
+
+  rKey = i_Mass;
+
+  const auto& rB = m_ZDependents.m_RGB;
+  double left = Herd::Generic::BXhC( i_Mass, rB[ 3 ], -rB[ 4 ] );
+  double right = Herd::Generic::BXhC( i_Mass, rB[ 5 ], -rB[ 6 ] );
+  rA = std::min( left, right );  // Eq. 46
+}
+
+/**
+ * @param i_Mass Mass
+ */
+void GiantBranchHelper::ComputeGB( Herd::Generic::Mass i_Mass )
+{
+  // Compute only if the mass value is different
+  auto& rEvaluatedAt = m_MDependents.m_McGB.m_EvaluatedAt;
+  if( rEvaluatedAt && *rEvaluatedAt == i_Mass )
+  {
+    return;
+  }
+
+  rEvaluatedAt = i_Mass;
+
+  auto& rMcGB = m_MDependents.m_McGB;
+  rMcGB.m_B = std::max( 3.0e+04, Herd::Generic::ApBXhC( i_Mass, 500, 17500, 0.6 ) );
+
+  // AMASS.SSE implementation for D, p and q replaces MHeF with 2.0
+  double blendWeight = std::clamp( Herd::Generic::ComputeBlendWeight( i_Mass, 2.0, 2.5 ), 0.0, 1.0 );
+  rMcGB.m_P = 6.0 - blendWeight;
+  rMcGB.m_Q = 3.0 - blendWeight;
+
+  // D
+  double d0 = m_ZDependents.m_D0;
+  double logD = 0;
+  if( i_Mass <= 2.0 )
+  {
+    logD = d0;
+  }
+
+  if( i_Mass > 2.0 && i_Mass < 2.5 )
+  {
+    logD = d0 - ( 0.025 * d0 + 0.45 ) * blendWeight; // This is linear interpolation between d0 and 0.975 * d0 - 0.18 * i_Mass at i_Mass = 2.5
+  }
+
+  if( i_Mass >= 2.5 )
+  {
+    logD = std::max( { -1.0, 0.5 * d0 - 0.06 * i_Mass, 0.975 * d0 - 0.18 * i_Mass } );
+  }
+
+  rMcGB.m_D = std::pow( 10.0, logD );
+
+  // Mx
+  rMcGB.m_Mx.Set( std::pow( rMcGB.m_B / rMcGB.m_D, 1.0 / ( rMcGB.m_P - rMcGB.m_Q ) ) );  // Eq. 38
+
+  // Lx
+  rMcGB.m_Lx.Set( std::pow( rMcGB.m_D * rMcGB.m_Mx, rMcGB.m_P ) );
+}
+
+}
+

src/SSE/Landmarks/GiantBranchHelper.h

@@ -0,0 +1,93 @@
+/**
+ * @file GiantBranchHelper.h
+ * @author Evren Imre
+ * @date 1 Sep 2021	
+ */
+/* This file is a part of HeRD, a stellar evolution library
+ * Copyright © 2021 Evren Imre 
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at https://mozilla.org/MPL/2.0/. */
+
+#ifndef HF72D459C_7D4D_4F08_84B3_4DE9A7E08995
+#define HF72D459C_7D4D_4F08_84B3_4DE9A7E08995
+
+#include <Generic/Quantity.h>
+
+#include <array>
+#include <optional>
+#include <utility>
+
+namespace Herd::SSE
+{
+/**
+ * @brief Giant branch radius computations
+ * @cite Hurley00
+ */
+class GiantBranchHelper
+{
+public:
+  GiantBranchHelper( Herd::Generic::Metallicity i_Z ); ///< Constructor
+  ~GiantBranchHelper(); ///< Destructor
+
+  Herd::Generic::Radius Radius( Herd::Generic::Mass i_Mass, Herd::Generic::Luminosity i_Luminosity ); ///< Computes \f$ R_{GB} \f$
+  Herd::Generic::Mass CoreMass( Herd::Generic::Mass i_Mass, Herd::Generic::Luminosity i_Luminosity ); ///<  Computes \f$ M_{c, GB}\f$ as per Eq. 37
+
+private:
+
+  void ComputeMetallicityDependents( Herd::Generic::Metallicity i_Z ); ///< Computes various metallicity-dependent quantities
+
+  // Mass dependents
+  void ComputeA( Herd::Generic::Mass i_Mass );  ///< Computes f$ A \f$ in Eq. 46
+  void ComputeGB( Herd::Generic::Mass i_Mass ); ///< Computes various giant branch parameters
+
+  /**
+   * @brief Various quantities and values that depend on metallicity only
+   */
+  struct MetallicityDependents
+  {
+    std::array< double, 7 > m_RGB;  ///< \f$ R_{GB}\f$ calculations
+
+    double m_D0;  ///< \f$ D_0 \f$ in \f$ M_{c, GB}\f$ calculations
+  };
+
+  MetallicityDependents m_ZDependents;  ///< Metallicity-dependent quantities
+
+  /**
+   * @brief Various quantities and values that depend on mass
+   * @remarks Note that these quantities also depend on metallicity via the equation coefficients
+   */
+  struct MassDependents
+  {
+    using TKey = std::optional< Herd::Generic::Mass >;
+
+    std::pair< TKey, double > m_A; ///< \f$ A \f$ in Eq. 46
+
+    /**
+     * @brief These parameters are used in Eq. 37 and Eq. 38
+     */
+    struct CoreMassParameters
+    {
+      TKey m_EvaluatedAt;  ///< Mass at which the vakyes are computed
+
+      double m_D = 0; ///< \f$ D \f$ in Eq. 37
+      double m_B = 0; ///< \f$ B \f$ in Eq. 37
+      double m_P = 0; ///< \f$ p \f$ in Eq. 37
+      double m_Q = 0; ///< \f$ q \f$ in Eq. 37
+
+      Herd::Generic::Mass m_Mx; ///< \f$ M_x \f$ in Eq. 38
+      Herd::Generic::Luminosity m_Lx; ///< Luminosity at \f$ M_x \f$
+    };
+    CoreMassParameters m_McGB;
+
+  };
+
+  MassDependents m_MDependents; ///< Mass-dependent quantities
+};
+
+}
+
+
+
+#endif /* HF72D459C_7D4D_4F08_84B3_4DE9A7E08995 */