move all fundamentals to constexpr

include/fundamental.hpp

@@ -27,9 +27,47 @@
 #include "precision.hpp"
 
 namespace fargOCA {
-  real constexpr PI = 3.14159265358979323844;
-  //real const PI4 = 4*PI;
+  namespace details {
+    real constexpr
+    cstSqrtImpl(real x, real l, real r) {
+      real const epsilon = x/1e19;
+      while ( l <= r ) {
+        real const mid = (l + r) / 2;
+        if ( mid * mid > x ) {
+          r = mid;
+        } else {
+          if ( x - mid * mid < epsilon ) {
+            return mid;
+          }
+          l = mid;
+        }
+      }
+      return -1;
+    }
+    
+    real constexpr
+    cstSqrt(real x) {
+      if ( x > 0 && x < 1 ) {
+        return 1/cstSqrt(1/x);
+      } else {
+        return cstSqrtImpl(x, 0, x);
+      }
+    }
     
+    real constexpr
+    cstPow(real r, int e) {
+      if (e == 1) {
+        return r;
+      } else {
+        real const s = cstPow(r, e/2);
+        return e%2 == 0 ? s*s : r*s*s;
+      }
+    }
+  }
+  using details::cstSqrt;
+  using details::cstPow;
+  
+  real constexpr PI       = 3.14159265358979323844;
   real constexpr YEAR     = 365.25*24*3600; ///< one year in s
   real constexpr RSUN     = 6.96e+10;       ///< Sun radius in cm 
   real constexpr XMH      = 1.67e-24;  ///< Hydrogen mass
@@ -53,14 +91,14 @@ namespace fargOCA {
   /// paper : use C2^2 instead of C2
   real constexpr CVNR     = 1.41; 
 
-  real const TIME0      = std::sqrt(VOL0/GRAVC/XM0); // Unit of time in s
-  real const V0         = R0/TIME0;
-  real const XNU0       = R0*V0;
-  real const TEMP0      = V0*V0/RGAS;
-  real const RHO0       = XM0/VOL0;
-  real const C0         = CLIGHT/V0;
-  real const A0         = AR*std::pow(TEMP0,4)/(RHO0*V0*V0);
-  real const SIGMARAD   = A0*C0/4;   // Boltzmann constant
-  real const E0         = XM0*R0*R0/(TIME0*TIME0);
+  real constexpr TIME0    = cstSqrt(VOL0/GRAVC/XM0); // Unit of time in s
+  real constexpr V0       = R0/TIME0;
+  real constexpr XNU0     = R0*V0;
+  real constexpr TEMP0    = cstPow(V0,2)/RGAS;
+  real constexpr RHO0     = XM0/VOL0;
+  real constexpr C0       = CLIGHT/V0;
+  real constexpr A0       = AR*cstPow(TEMP0,4)/(RHO0*cstPow(V0,2));
+  real constexpr SIGMARAD = A0*C0/4;   // Boltzmann constant
+  real constexpr E0       = XM0*cstPow(R0,2)/cstPow(TIME0,2);
 }
 #endif