Got rid of more global variables

config.cpp

@@ -121,7 +121,7 @@ void Config::error_checking()
         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++)
+    for (int i=0; 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");
     if (ext_units_physical && ((norm_mass == 0) || (norm_length == 0)))

phigrape.cpp

@@ -325,26 +325,6 @@ double  r2, rv_ij,
         // NOTE this is a mess. Looks like eps_bh is never used, and the m_bh from outside the block is never used.
 #endif
 
-double pot_ext[N_MAX], pot_act_ext[N_MAX];    // for all EXTPOT
-
-
-int    i_bh, i_bh1, i_bh2,
-    num_bh = 0, num_bh1 = 0, num_bh2 = 0;
-
-double m_bh, m_bh1, m_bh2, b_bh,
-    r, r2,
-    x_ij, y_ij, z_ij,
-    vx_ij, vy_ij, vz_ij, rv_ij;
-
-// #ifdef BBH_INF
-int    inf_event[N_MAX];
-double DR2, tmp_r2;
-double3 x_bbhc, v_bbhc;
-double DV2, EB, SEMI_a, SEMI_a2;
-// #endif // BBH_INF
-
-// #ifdef ADD_N_BH
-
 // double x_bh1[3], x_bh2[3], v_bh1[3], v_bh2[3];
 double3 x_bh1, x_bh2, v_bh1, v_bh2;
 
@@ -358,99 +338,11 @@ double3 a_bh1, adot_bh1, a_bh2, adot_bh2;
 
 #include "n_bh.c"
 
-/*
-int calc_force_n_BH(double m1, double xx1[], double vv1[],
-                    double m2, double xx2[], double vv2[],
-                    double eps_BH,
-                    double pot_n1, double a_n1[], double adot_n1[],
-                    double pot_n2, double a_n2[], double adot_n2[])
-*/
-
-/*
-                INPUT
-
-m1         - mass of the 1 BH
-xx1[0,1,2] - coordinate of the 1 BH
-vv1[0,1,2] - velocity of the 1 BH
-
-m2         - mass of the 2 BH
-xx2[0,1,2] - coordinate of the 2 BH
-vv2[0,1,2] - velocity of the 2 BH
-
-eps_BH     - force softening, can be even exactly 0.0 !
-
-                OUTPUT
-
-pot_n1          for the 1 BH
-a_n1    [0,1,2] for the 1 BH
-adot_n1 [0,1,2] for the 1 BH
-
-pot_n2          for the 2 BH
-a_n2    [0,1,2] for the 2 BH
-adot_n2 [0,1,2] for the 2 BH
-
-return -   0 if everything OK
-*/
-
-// #endif // ADD_N_BH
-
-// #ifdef ADD_PN_BH
-
-// 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];
 double3 a_pn1[7], adot_pn1[7], a_pn2[7], adot_pn2[7];
 
-double s_bh1[3] = {0.0, 0.0, 1.0};
-double s_bh2[3] = {0.0, 0.0, 1.0};
 
 #include "pn_bh_spin.c"
 
-/*
-int calc_force_pn_BH(double m1, double xx1[], double vv1[], double ss1[],
-                    double m2, double xx2[], double vv2[], double ss2[],
-                    double CCC_NB, double dt_bh,
-                    int usedOrNot[],
-                    double a_pn1[][3], double adot_pn1[][3],
-                    double a_pn2[][3], double adot_pn2[][3])
-*/
-
-/*
-                INPUT
-
-m1         - mass of the 1 BH
-xx1[0,1,2] - coordinate of the 1 BH
-vv1[0,1,2] - velocity of the 1 BH
-spin1[0,1,2] - normalized spin of the 1 BH
-
-m2         - mass of the 2 BH
-xx2[0,1,2] - coordinate of the 2 BH
-vv2[0,1,2] - velocity of the 2 BH
-spin2[0,1,2] - normalized spin of the 2 BH
-
-CCC_NB           - Speed of light "c" in internal units
-dt_BH         - timestep of the BH's, needed for the SPIN integration
-
-usedOrNot[PN0, PN1, PN2, PN2.5, PN3, PN3.5, SPIN] - different PN term usage: PN1, PN2, PN2.5, PN3, PN3.5, SPIN
-        0    1    2    3      4    5      6
-
-                OUTPUT
-
-a_pn1   [0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5, 6 - SPIN] [3]    for the 1 BH
-adot_pn1[0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5, 6 - SPIN] [3]    for the 1 BH
-
-a_pn2   [0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5, 6 - SPIN] [3]    for the 2 BH
-adot_pn2[0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5, 6 - SPIN] [3]    for the 2 BH
-
-return -   0 if everything OK
-    - 505 if BH's separation < 4 x (RSwarch1 + RSwarch2)
-*/
-
-// #endif // ADD_N_BH
-
 #ifdef ETICS
 double t_exp, dt_exp=ETICS_DTSCF; // t_exp is just the initial value
 #ifdef ETICS_CEP
@@ -460,28 +352,25 @@ int grapite_cep_index;
 
 void get_CPU_time(double *time_real, double *time_user, double *time_syst)
 {
-struct rusage xxx;
+    struct rusage xxx;
-double sec_u, microsec_u, sec_s, microsec_s;
+    double sec_u, microsec_u, sec_s, microsec_s;
-struct timeval tv;
+    struct timeval tv;
 
-getrusage(RUSAGE_SELF,&xxx);
+    getrusage(RUSAGE_SELF,&xxx);
 
-sec_u = xxx.ru_utime.tv_sec;
+    sec_u = xxx.ru_utime.tv_sec;
-sec_s = xxx.ru_stime.tv_sec;
+    sec_s = xxx.ru_stime.tv_sec;
 
-microsec_u = xxx.ru_utime.tv_usec;
+    microsec_u = xxx.ru_utime.tv_usec;
-microsec_s = xxx.ru_stime.tv_usec;
+    microsec_s = xxx.ru_stime.tv_usec;
 
-*time_user = sec_u + microsec_u * 1.0E-06;
+    *time_user = sec_u + microsec_u * 1.0E-06;
-*time_syst = sec_s + microsec_s * 1.0E-06;
+    *time_syst = sec_s + microsec_s * 1.0E-06;
 
-//  *time_real = time(NULL);
+    gettimeofday(&tv, NULL);
-
+    *time_real = tv.tv_sec + 1.0E-06 * tv.tv_usec;
-gettimeofday(&tv, NULL);
-*time_real = tv.tv_sec + 1.0E-06 * tv.tv_usec;
-
-*time_user = *time_real;
 
+    *time_user = *time_real;
 }
 
 void read_data(char inp_fname[], int *diskstep, int *N, double *time_cur, int ind[], double m[], double3 x[], double3 v[])
@@ -885,7 +774,7 @@ DT_EXT_GRAV += (CPU_tmp_user - CPU_tmp_user0);
 
 }
 
-void energy_contr(const double time_cur, const double timesteps, const double n_act_sum, const double g6_calls, double& rcm_sum, double& vcm_sum, double& E_tot_0, double& E_tot_corr_0, double& E_tot_corr_sd_0, int &skip_con, int N, double m[], double3 x[], double3 v[], double pot[])
+void energy_contr(const double time_cur, const double timesteps, const double n_act_sum, const double g6_calls, double& rcm_sum, double& vcm_sum, double& E_tot_0, double& E_tot_corr_0, double& E_tot_corr_sd_0, int &skip_con, int N, double m[], double3 x[], double3 v[], double pot[], double pot_ext[])
 {
     // TODO maybe use static variables here for the previous step energy?
     double E_pot = 0.0;
@@ -1058,6 +947,21 @@ int main(int argc, char *argv[])
     double C_NB;
     int    usedOrNot[7];
 
+    double pot_ext[N_MAX], pot_act_ext[N_MAX];    // for all EXTPOT
+
+    int    i_bh1, i_bh2;
+
+    double m_bh1, m_bh2;
+
+    int    inf_event[N_MAX];
+    double DR2, tmp_r2;
+    double3 x_bbhc, v_bbhc;
+    double DV2, EB, SEMI_a, SEMI_a2;
+
+    double s_bh1[3] = {0.0, 0.0, 1.0};
+    double s_bh2[3] = {0.0, 0.0, 1.0};
+
+
     /* INIT the rand() !!! */
     srand(19640916);                 /* it is just my birthday :-) */
 
@@ -1503,7 +1407,7 @@ int main(int argc, char *argv[])
 
     /* Energy control... */
     if (myRank == rootRank) {
-        energy_contr(time_cur, timesteps, n_act_sum, g6_calls, rcm_sum, vcm_sum, E_tot_0, E_tot_corr_0, E_tot_corr_sd_0, skip_con, N, m, x, v, pot);
+        energy_contr(time_cur, timesteps, n_act_sum, g6_calls, rcm_sum, vcm_sum, E_tot_0, E_tot_corr_0, E_tot_corr_sd_0, skip_con, N, m, x, v, pot, pot_ext);
     } /* if (myRank == rootRank) */
 
 #ifdef ETICS_DUMP
@@ -2130,7 +2034,7 @@ int main(int argc, char *argv[])
         if (time_cur >= t_contr) {
             if (myRank == rootRank) {
 
-                energy_contr(time_cur, timesteps, n_act_sum, g6_calls, rcm_sum, vcm_sum, E_tot_0, E_tot_corr_0, E_tot_corr_sd_0, skip_con, N, m, x, v, pot);
+                energy_contr(time_cur, timesteps, n_act_sum, g6_calls, rcm_sum, vcm_sum, E_tot_0, E_tot_corr_0, E_tot_corr_sd_0, skip_con, N, m, x, v, pot, pot_ext);
 
                 /* write cont data */