Added post-Newtonian parameters to the config file

config.cpp

@@ -2,6 +2,7 @@
 #include <stdexcept>
 #include <cstdio>
 #include <fstream>
+#include <algorithm>
 
 using Dictionary = std::unordered_map<std::string,std::string>;
 
@@ -54,7 +55,7 @@ template<> double Config::string_cast(const std::string str)
 template<> int Config::string_cast(const std::string str)
 {
     size_t idx;
-    auto value = std::stoi(str, &idx);
+    auto value = std::stoi(str, &idx); // WARNING stoi can throw
     if (idx == str.size()) return value;
     else throw std::runtime_error("Cannot convert \"" + str + "\" into an int");
 }
@@ -66,6 +67,24 @@ template<> bool Config::string_cast(const std::string str)
     throw std::runtime_error("Cannot convert \"" + str + "\" into a bool");
 }
 
+#include <iostream>
+template<> std::vector<int> Config::string_cast(const std::string str)
+{
+    auto error = std::runtime_error("Cannot convert \"" + str + "\" into an integer array");
+    if (!( (str.front()=='{') && (str.back()=='}'))) throw error;
+    std::string new_str = strip(str.substr(1, str.length()-2));
+    std::replace(new_str.begin(), new_str.end(), ',', ' ');
+
+    std::vector<int> result;
+    while (new_str.length() > 0) {
+        size_t idx;
+        auto value = std::stoi(new_str, &idx);
+        result.push_back(value);
+        new_str = new_str.substr(idx, new_str.length()-idx);
+    }
+    return result;
+}
+
 
 // For mandatory parameters
 template<typename T>
@@ -98,6 +117,11 @@ void Config::error_checking()
         throw std::runtime_error("Black hole neighbour output (live_smbh_neighbor_output=true) requires at least one live black hole (live_smbh_count)");
     if (binary_smbh_influence_sphere_output && (live_smbh_count != 2))
         throw std::runtime_error("Binary black hole influence sphere output (binary_smbh_influence_sphere_output=true) requires exactly two live black holes (live_smbh_count=2)");
+    if (pn_usage.size() != 7)
+        throw std::runtime_error("PN usage array (pn_usage) must have exactly seven components");
+    for (int i; i<7; i++)
+        if (!((pn_usage[i] == 0) || (pn_usage[i] == 1)))
+            throw std::runtime_error("PN usage array (pn_usage) must have ones and zeros only");
 }
 
 Config::Config(std::string file_name)
@@ -121,6 +145,8 @@ Config::Config(std::string file_name)
     binary_smbh_pn = get_parameter<bool>(dictionary, "binary_smbh_pn", false);
     binary_smbh_influence_sphere_output = get_parameter<bool>(dictionary, "binary_smbh_influence_sphere_output", false);
     binary_smbh_influence_radius_factor = get_parameter<double>(dictionary, "binary_smbh_influence_radius_factor", 10.);
+    pn_usage = get_parameter<std::vector<int>>(dictionary, "pn_usage", std::vector<int>({1,1,1,1,1,1,1}));
+    pn_c = get_parameter<double>(dictionary, "pn_c", 500);
 
     error_checking();
 }

config.h

@@ -1,27 +1,30 @@
 #pragma once
 #include <string>
 #include <unordered_map>
+#include <vector>
 
 class Config {
 public:
     Config(std::string file_name);
 
-    double      eps;
-    double      t_end;
-    double      dt_disk;
-    double      dt_contr;
-    double      dt_bh;
-    double      eta;
-    std::string input_file_name;
-    bool        dt_min_warning;
-    int         live_smbh_count;
-    double      live_smbh_custom_eps;
-    bool        live_smbh_output;
-    bool        live_smbh_neighbor_output;
-    int         live_smbh_neighbor_number;
-    bool        binary_smbh_pn;
-    bool        binary_smbh_influence_sphere_output;
-    double      binary_smbh_influence_radius_factor;
+    double           eps;
+    double           t_end;
+    double           dt_disk;
+    double           dt_contr;
+    double           dt_bh;
+    double           eta;
+    std::string      input_file_name;
+    bool             dt_min_warning;
+    int              live_smbh_count;
+    double           live_smbh_custom_eps;
+    bool             live_smbh_output;
+    bool             live_smbh_neighbor_output;
+    int              live_smbh_neighbor_number;
+    bool             binary_smbh_pn;
+    bool             binary_smbh_influence_sphere_output;
+    double           binary_smbh_influence_radius_factor;
+    std::vector<int> pn_usage;
+    double           pn_c;
 private:
     using Dictionary = std::unordered_map<std::string,std::string>;
     Dictionary read_config_file(const std::string file_name);

phigrape.conf

@@ -66,6 +66,13 @@ binary_smbh_influence_radius_factor =   3.162277660168379497918067e+03
 # Add post Newtonian terms to SMBH-SMBH gravity. [default: false]
 binary_smbh_pn =                        true
 
+# A mask array (zeros and ones) determining whether or not to use specific post-Newtonian terms.
+# The elements represent {0, 1, 2, 2.5, 3, 3.5, spin}
+pn_usage =                              {1, 1, 1, 1, 0, 0, 0}
+
+# The speed of light in N-body units [default: 500]
+pn_c =                                  477.12
+
 ####################################
 # Negative powers of two           #
 ####################################

phigrape.cpp

@@ -400,9 +400,9 @@ return -   0 if everything OK
 
 // #ifdef ADD_PN_BH
 
-double C_NB = 477.12;
-
-int    usedOrNot[7] = {1, 1, 1, 1, 0, 0, 0};
+// double C_NB = 477.12;
+// 
+// int    usedOrNot[7] = {1, 1, 1, 1, 0, 0, 0};
 
 // double a_pn1[7][3], adot_pn1[7][3],
 //     a_pn2[7][3], adot_pn2[7][3];
@@ -1339,6 +1339,8 @@ int main(int argc, char *argv[])
 
     FILE   *out;
 
+    double C_NB;
+    int    usedOrNot[7];
 
     /* INIT the rand() !!! */
     srand(19640916);                 /* it is just my birthday :-) */
@@ -1364,6 +1366,8 @@ int main(int argc, char *argv[])
     /* read the input parameters to the rootRank */
 
     config = new Config("phigrape.conf");
+    C_NB = config->pn_c;
+    std::copy(config->pn_usage.begin(), config->pn_usage.end(), usedOrNot);
 
     if (myRank == rootRank) {