Service release 1.5.1

This fixes issues #26, issue #27, and issue #28

Author: Wolfgang Kastaun

ET_interface/thorns/RePrimAnd/dist/RePrimAnd.tar

@@ -917,6 +917,10 @@ All quantities are in the same geometric units used by the EOS.
              py::vectorize(&etk::spherical_star::lambda_from_rc),
              "Metric function lambda at circumf. radius",
              py::arg("rc"))     
+        .def("gm1_from_rc", 
+             py::vectorize(&etk::spherical_star::gm1_from_rc),
+             "Pseudo enthalpy g - 1 at circumf. radius",
+             py::arg("rc"))
         .def("mbary_from_rc", 
              py::vectorize(&etk::spherical_star::mbary_from_rc),
              "Baryonic mass within circumf. radius",

bindings/python/pyreprimand.cpp

@@ -107,7 +107,7 @@ class units {
   static units geom_udensity(double udensity, 
                                        double g_si=G_SI)
   {
-    return geom_ulength( c_SI / std::sqrt(g_si * udensity));
+    return geom_ulength( c_SI / std::sqrt(g_si * udensity), g_si);
   }
 
   /**\brief Compute units with G=c=1 and given mass unit
@@ -118,7 +118,7 @@ class units {
   **/  
   static constexpr units geom_umass(double umass, double g_si=G_SI)
   {
-    return geom_ulength(umass * g_si /(c_SI*c_SI));
+    return geom_ulength(umass * g_si /(c_SI*c_SI), g_si);
   }
 
   /**\brief Compute units with G=c=1 and the mass unit is the solar mass

library/BasicStuff/include/unitconv.h

@@ -68,12 +68,11 @@ real_t eos_barotr_gpoly::hm1(real_t gm1) const
 }
 
 /**
-\return Soundspeed squared \f$ c_s^2 = \frac{h-1}{nh} \f$
+\return Soundspeed squared \f$ c_s^2 = \frac{g-1}{ng} \f$
 */
 real_t eos_barotr_gpoly::csnd(real_t gm1) const
 {
-  real_t h1 = hm1(gm1);
-  real_t cs2 = h1/(n*(h1+1));
+  real_t cs2 = gm1/(n*(gm1+1));
   return sqrt(cs2);
 }
 

library/EOS_Barotropic/eos_barotr_gpoly.cc

@@ -64,26 +64,47 @@ tidal_ode::tidal_ode(eos_barotr eos_, real_t gm1_center_,
                              "computed for isentropic EOS"));
   }
 
-  std::vector<real_t> revrho;
-  for (auto i = dnu_.rbegin(); i != dnu_.rend(); ++i) 
+  std::vector<real_t> revrho, revlambda, revdnu, revrsqr, revmbr3;
+  
+  auto ilambda = lambda_.rbegin();
+  auto irsqr = rsqr_.rbegin();
+  for (auto idnu = dnu_.rbegin(); idnu != dnu_.rend(); ++idnu) 
   {
-    real_t gm1{ gm1_from_dnu(*i) };
+    assert(ilambda != lambda_.rend());
+    assert(irsqr != rsqr_.rend());
+    
+    real_t dnu{ *idnu };
+    real_t lambda{ *ilambda++ };
+    real_t rsqr{ *irsqr++ };
+
+    real_t gm1{ gm1_from_dnu(dnu) };
     auto s{ eos.at_gm1(eos.range_gm1().limit_to(gm1)) };
+    
     assert(s);
-    revrho.push_back(s.rho());
+    assert(lambda >= 0.0);
+    assert(rsqr >= 0.0);
+    
+    real_t rho{ s.rho() };
+    real_t rho_e{ rho * (1.0 + s.eps()) };
+    real_t mbr3{ m_by_r3(rsqr, lambda, rho_e) };
+    
+    revrho.push_back(rho);
+    revlambda.push_back(lambda);
+    revdnu.push_back(dnu);
+    revrsqr.push_back(rsqr);
+    revmbr3.push_back(mbr3);
   }
 
 
-  const std::vector<real_t> revdnu(dnu_.rbegin(), dnu_.rend());
-
-  const std::vector<real_t> revlambda(lambda_.rbegin(), lambda_.rend());
-  const std::vector<real_t> revrsqr(rsqr_.rbegin(), rsqr_.rend());
-  
   dnu_rho = make_interpol_pchip_spline(revrho, revdnu);
 
   lambda_rho = make_interpol_pchip_spline(revrho, revlambda);
 
   rsqr_rho = make_interpol_pchip_spline(revrho, revrsqr);
+
+  mbr3_rho = make_interpol_pchip_spline(revrho, revmbr3);
+
+
   assert(x_start()>x_end());
 }
 
@@ -111,21 +132,20 @@ auto tidal_ode::drho_y(real_t rho_, real_t ym2) const -> real_t
   real_t cs2{ std::pow(s.csnd(), 2) }; 
 
   real_t rho_e{ rho * (1.0 + eps) };
-  
-  real_t lambda{ lambda_rho(rho) };
+  assert(lambda_rho.range_x().contains(rho));
+  real_t lambda{ std::max(0.0, lambda_rho(rho)) };
   real_t e2l{ std::exp(2.0 * lambda) };
   real_t rsqr{ rsqr_rho(rho) };
   real_t wtfac{ cs2 / rho };
- 
-      
-  real_t mbyr3{ m_by_r3(rsqr, lambda, rho_e) };
+  real_t mbyr3{ mbr3_rho(rho) };
 
   real_t a{ 4.0*PI * p + mbyr3 };
   real_t b{ 2.0 * rsqr * (mbyr3 + 2.0*PI * (p - rho_e)) };
   real_t c{ 4.0*PI * (3.0 * rho_e + 11.0 * p) - 8.0 * mbyr3 };
   real_t d{ 4.0 * rsqr * wtfac * e2l * a };
   real_t g{ (ym2+2.0 + 3.0) / e2l + b };
 
+  
   real_t f{0.0};
   if (rsqr>0) {
     f = wtfac * ym2 / rsqr;
@@ -134,9 +154,11 @@ auto tidal_ode::drho_y(real_t rho_, real_t ym2) const -> real_t
     f = (-4.0*PI/7.0) * (
                 (11.0 * h - (32.0/3.0) * (1.0 + eps) ) * cs2 + h);
   }
-  
-  
-  return (4.0*PI * h +  (f * g + wtfac * c)) / a - d;
+
+
+  real_t res{ (4.0*PI * h +  (f * g + wtfac * c)) / a - d };
+  assert(std::isfinite(res));
+  return res;
 }
 
 auto tidal_ode::initial_data() const -> state_t

library/NeutronStar/TOV/tidal_deform_ode.cc

@@ -17,6 +17,7 @@ class tidal_ode {
   interpolator dnu_rho;
   interpolator lambda_rho;
   interpolator rsqr_rho;
+  interpolator mbr3_rho;
 
   auto gm1_from_dnu(real_t dnu) const -> real_t;
 

library/NeutronStar/TOV/tidal_deform_ode.h

@@ -1,6 +1,6 @@
 
 project('RePrimAnd', ['cpp','c'], default_options : ['cpp_std=c++11'],
-        version : '1.5', license : 'CC BY-NC-SA 4.0')
+        version : '1.5.1', license : 'CC BY-NC-SA 4.0')
 
 project_headers_dest = 'reprimand'