From 22f983cf17964338f132659d3fbf6591dffad8fa Mon Sep 17 00:00:00 2001
From: Yohai Meiron <sahmes@gmail.com>
Date: Wed, 18 Mar 2020 11:55:24 -0400
Subject: [PATCH] Moved some variables from global to main and removed some
 unused variables

---
 .gitignore   |   2 +
 phigrape.cpp | 265 ++++++++++++++++++++++++---------------------------
 pn_bh_spin.c |  15 +--
 3 files changed, 130 insertions(+), 152 deletions(-)

diff --git a/.gitignore b/.gitignore
index 9f184f3..773bfd7 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,3 +6,5 @@
 *.cfg
 *.inp
 *.mask
+grapite-dev-exec-threshold
+phigrape
diff --git a/phigrape.cpp b/phigrape.cpp
index b194e0b..b545836 100644
--- a/phigrape.cpp
+++ b/phigrape.cpp
@@ -306,39 +306,12 @@ double3 operator-(const double3& a, const double3& b)
 }
 
 
-int    name_proc, n_proc=1, myRank=0, rootRank=0, cur_rank,
-    i, j, k, ni, nj, diskstep=0, power, jjj, iii,
-    skip_con=0, tmp_i;
-
-double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
-    dt_bh, t_bh=0.0, dt_bh_tmp,
-    t_end, time_cur, dt_min, dt_max, min_t, min_t_loc, dt_new, min_dt,
-    eta_s, eta, eta_bh,
-    E_pot, E_pot_ext, E_kin, E_tot, E_tot_0, DE_tot,
-    E_tot_corr, E_tot_corr_0, DE_tot_corr,
-    E_tot_corr_sd, E_tot_corr_sd_0, DE_tot_corr_sd,
-    E_corr = 0.0, E_sd = 0.0, t_diss_on = 0.125,
-    mcm, rcm_mod, vcm_mod,
-    rcm_sum=0.0, vcm_sum=0.0,
-    eps=0.0, eps2,
-    a2_mod, adot2_mod,
-    dt_tmp, dt2half, dt3over6, dt4over24, dt5over120,
-    dtinv, dt2inv, dt3inv,
-    a0mia1, ad04plad12, ad0plad1,
-    a1abs, adot1abs, a2dot1abs, a3dot1abs,
-    timesteps=0.0, n_act_sum=0.0, n_act_distr[N_MAX], g6_calls=0.0, g6_calls_sum=0.0,
-    tmp, tmp_r, tmp_v, tmp_rv, tmp_cpu,
-    tmp_pot, tmp_a, tmp_adot,
-    tmp_a_bh, tmp_adot_bh,
-    tmp_a_bh_pn0, tmp_a_bh_pn1, tmp_a_bh_pn2, tmp_a_bh_pn2_5, tmp_a_bh_pn3, tmp_a_bh_pn3_5, tmp_a_bh_spin,
-    tmp_adot_bh_pn0, tmp_adot_bh_pn1, tmp_adot_bh_pn2, tmp_adot_bh_pn2_5, tmp_adot_bh_pn3, tmp_adot_bh_pn3_5, tmp_adot_bh_spin;
-
 
 double L[3]; // needed in pn_bh_spin.c
 // Needed for things related to BHs
 #include "debug.h"
-int    ind_sort[N_MAX];
-double var_sort[N_MAX];
+// int    ind_sort[N_MAX];
+// double var_sort[N_MAX];
 
 
 double3 xcm, vcm, mom,
@@ -611,28 +584,24 @@ gettimeofday(&tv, NULL);
 
 }
 
-void read_data()
+void read_data(int *diskstep, int *N, double *time_cur, int ind[], double m[], double3 x[], double3 v[])
 {
-inp = fopen(inp_fname,"r");
-fscanf(inp,"%d \n", &diskstep);
-
-fscanf(inp,"%d \n", &N);
-
-fscanf(inp,"%lE \n", &time_cur);
-
-for (i=0; i<N; i++) fscanf(inp,"%d %lE %lE %lE %lE %lE %lE %lE \n", &ind[i], &m[i], &x[i][0], &x[i][1], &x[i][2], &v[i][0], &v[i][1], &v[i][2]);
-
-fclose(inp);
+    auto inp = fopen(inp_fname, "r");
+    fscanf(inp, "%d \n", diskstep);
+    fscanf(inp, "%d \n", N);
+    printf("zzzzzz reading *%s* %d\n", inp_fname, *N);
+    fscanf(inp, "%lE \n", time_cur);
+    for (int i=0; i<*N; i++) fscanf(inp,"%d %lE %lE %lE %lE %lE %lE %lE \n", &ind[i], &m[i], &x[i][0], &x[i][1], &x[i][2], &v[i][0], &v[i][1], &v[i][2]);
+    fclose(inp);
 }
 
-void write_snap_data()
+void write_snap_data(char out_fname[], int diskstep, int N, double time_cur, int ind[], double m[], double3 x[], double3 v[])
 {
-    sprintf(out_fname, "%06d.dat", diskstep);
     out = fopen(out_fname, "w");
     fprintf(out,"%06d \n", diskstep);
     fprintf(out,"%07d \n", N);
     fprintf(out,"%.10E \n", time_cur);
-    for (i=0; i<N; i++) {
+    for (int i=0; i<N; i++) {
         fprintf(out,"%07d \t %.10E \t % .10E % .10E % .10E \t % .10E % .10E % .10E \n",
             ind[i],
             m[i],
@@ -642,35 +611,34 @@ void write_snap_data()
     fclose(out);
 }
 
-void write_cont_data() // virtually identical to write_snap_data()
-{
-out = fopen("data.con","w");
-fprintf(out,"%06d \n", diskstep);
+// // // // // // // // // // // void write_cont_data() // virtually identical to write_snap_data()
+// // // // // // // // // // // {
+// // // // // // // // // // // out = fopen("data.con","w");
+// // // // // // // // // // // fprintf(out,"%06d \n", diskstep);
+// // // // // // // // // // // 
+// // // // // // // // // // // fprintf(out,"%07d \n", N);
+// // // // // // // // // // // 
+// // // // // // // // // // // fprintf(out,"%.16E \n", time_cur);
+// // // // // // // // // // // 
+// // // // // // // // // // // for (i=0; i<N; i++)
+// // // // // // // // // // //     {
+// // // // // // // // // // //     fprintf(out,"%07d \t %.10E \t % .10E % .10E % .10E \t % .10E % .10E % .10E \n",
+// // // // // // // // // // //         ind[i],
+// // // // // // // // // // //         m[i],
+// // // // // // // // // // //         x[i][0], x[i][1], x[i][2],
+// // // // // // // // // // //         v[i][0], v[i][1], v[i][2]);
+// // // // // // // // // // //     }
+// // // // // // // // // // // 
+// // // // // // // // // // // fclose(out);
+// // // // // // // // // // // }
 
-fprintf(out,"%07d \n", N);
-
-fprintf(out,"%.16E \n", time_cur);
-
-for (i=0; i<N; i++)
-    {
-    fprintf(out,"%07d \t %.10E \t % .10E % .10E % .10E \t % .10E % .10E % .10E \n",
-        ind[i],
-        m[i],
-        x[i][0], x[i][1], x[i][2],
-        v[i][0], v[i][1], v[i][2]);
-    }
-
-fclose(out);
-}
-
-void write_bh_data()
+void write_bh_data(double time_cur, double m[], double3 x[], double3 v[])
 {
     if (config->live_smbh_count == 2) {
 
-        out = fopen("bh.dat", "a");
+        auto out = fopen("bh.dat", "a");
 
         for (int i=0; i < 2; i++) {
-//             double (*a_pn)[3], (*adot_pn)[3], pot_bh, *a_bh, *adot_bh;
             double3 *a_pn, *adot_pn, a_bh, adot_bh;
             double pot_bh;
             if (i==0) {
@@ -687,33 +655,33 @@ void write_bh_data()
                 adot_bh = adot_bh2;
             }
 
-            tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) );
-            tmp_v = sqrt( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
+            double tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) );
+            double tmp_v = sqrt( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
 
-            tmp_a = sqrt( SQR(a[i][0]) + SQR(a[i][1]) + SQR(a[i][2]) );
-            tmp_adot = sqrt( SQR(adot[i][0]) + SQR(adot[i][1]) + SQR(adot[i][2]) );
+            double tmp_a = sqrt( SQR(a[i][0]) + SQR(a[i][1]) + SQR(a[i][2]) );
+            double tmp_adot = sqrt( SQR(adot[i][0]) + SQR(adot[i][1]) + SQR(adot[i][2]) );
 
             if (config->live_smbh_custom_eps >= 0) {
 
-                tmp_a_bh = sqrt( SQR(a_bh[0]) + SQR(a_bh[1]) + SQR(a_bh[2]) );
-                tmp_adot_bh = sqrt( SQR(adot_bh[0]) + SQR(adot_bh[1]) + SQR(adot_bh[2]) );
+                double tmp_a_bh = sqrt( SQR(a_bh[0]) + SQR(a_bh[1]) + SQR(a_bh[2]) );
+                double tmp_adot_bh = sqrt( SQR(adot_bh[0]) + SQR(adot_bh[1]) + SQR(adot_bh[2]) );
 
                 if (config->binary_smbh_pn) {
-                    tmp_a_bh_pn0   = sqrt( SQR(a_pn[0][0]) + SQR(a_pn[0][1]) + SQR(a_pn[0][2]) );
-                    tmp_a_bh_pn1   = sqrt( SQR(a_pn[1][0]) + SQR(a_pn[1][1]) + SQR(a_pn[1][2]) );
-                    tmp_a_bh_pn2   = sqrt( SQR(a_pn[2][0]) + SQR(a_pn[2][1]) + SQR(a_pn[2][2]) );
-                    tmp_a_bh_pn2_5 = sqrt( SQR(a_pn[3][0]) + SQR(a_pn[3][1]) + SQR(a_pn[3][2]) );
-                    tmp_a_bh_pn3   = sqrt( SQR(a_pn[4][0]) + SQR(a_pn[4][1]) + SQR(a_pn[4][2]) );
-                    tmp_a_bh_pn3_5 = sqrt( SQR(a_pn[5][0]) + SQR(a_pn[5][1]) + SQR(a_pn[5][2]) );
-                    tmp_a_bh_spin  = sqrt( SQR(a_pn[6][0]) + SQR(a_pn[6][1]) + SQR(a_pn[6][2]) );
+                    double tmp_a_bh_pn0   = sqrt( SQR(a_pn[0][0]) + SQR(a_pn[0][1]) + SQR(a_pn[0][2]) );
+                    double tmp_a_bh_pn1   = sqrt( SQR(a_pn[1][0]) + SQR(a_pn[1][1]) + SQR(a_pn[1][2]) );
+                    double tmp_a_bh_pn2   = sqrt( SQR(a_pn[2][0]) + SQR(a_pn[2][1]) + SQR(a_pn[2][2]) );
+                    double tmp_a_bh_pn2_5 = sqrt( SQR(a_pn[3][0]) + SQR(a_pn[3][1]) + SQR(a_pn[3][2]) );
+                    double tmp_a_bh_pn3   = sqrt( SQR(a_pn[4][0]) + SQR(a_pn[4][1]) + SQR(a_pn[4][2]) );
+                    double tmp_a_bh_pn3_5 = sqrt( SQR(a_pn[5][0]) + SQR(a_pn[5][1]) + SQR(a_pn[5][2]) );
+                    double tmp_a_bh_spin  = sqrt( SQR(a_pn[6][0]) + SQR(a_pn[6][1]) + SQR(a_pn[6][2]) );
 
-                    tmp_adot_bh_pn0   = sqrt( SQR(adot_pn[0][0]) + SQR(adot_pn[0][1]) + SQR(adot_pn[0][2]) );
-                    tmp_adot_bh_pn1   = sqrt( SQR(adot_pn[1][0]) + SQR(adot_pn[1][1]) + SQR(adot_pn[1][2]) );
-                    tmp_adot_bh_pn2   = sqrt( SQR(adot_pn[2][0]) + SQR(adot_pn[2][1]) + SQR(adot_pn[2][2]) );
-                    tmp_adot_bh_pn2_5 = sqrt( SQR(adot_pn[3][0]) + SQR(adot_pn[3][1]) + SQR(adot_pn[3][2]) );
-                    tmp_adot_bh_pn3   = sqrt( SQR(adot_pn[4][0]) + SQR(adot_pn[4][1]) + SQR(adot_pn[4][2]) );
-                    tmp_adot_bh_pn3_5 = sqrt( SQR(adot_pn[5][0]) + SQR(adot_pn[5][1]) + SQR(adot_pn[5][2]) );
-                    tmp_adot_bh_spin  = sqrt( SQR(adot_pn[6][0]) + SQR(adot_pn[6][1]) + SQR(adot_pn[6][2]) );
+                    double tmp_adot_bh_pn0   = sqrt( SQR(adot_pn[0][0]) + SQR(adot_pn[0][1]) + SQR(adot_pn[0][2]) );
+                    double tmp_adot_bh_pn1   = sqrt( SQR(adot_pn[1][0]) + SQR(adot_pn[1][1]) + SQR(adot_pn[1][2]) );
+                    double tmp_adot_bh_pn2   = sqrt( SQR(adot_pn[2][0]) + SQR(adot_pn[2][1]) + SQR(adot_pn[2][2]) );
+                    double tmp_adot_bh_pn2_5 = sqrt( SQR(adot_pn[3][0]) + SQR(adot_pn[3][1]) + SQR(adot_pn[3][2]) );
+                    double tmp_adot_bh_pn3   = sqrt( SQR(adot_pn[4][0]) + SQR(adot_pn[4][1]) + SQR(adot_pn[4][2]) );
+                    double tmp_adot_bh_pn3_5 = sqrt( SQR(adot_pn[5][0]) + SQR(adot_pn[5][1]) + SQR(adot_pn[5][2]) );
+                    double tmp_adot_bh_spin  = sqrt( SQR(adot_pn[6][0]) + SQR(adot_pn[6][1]) + SQR(adot_pn[6][2]) );
 
                     fprintf(out,"%.16E \t %.8E \t % .16E % .16E % .16E \t %.16E \t % .16E % .16E % .16E \t %.16E \t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t %.8E \t\t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t\t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \n",
                             time_cur, m[i],
@@ -768,10 +736,10 @@ void write_bh_data()
         fclose(out);
     } else if (config->live_smbh_count == 1) {
         out = fopen("bh.dat", "a");
-        tmp_r = sqrt( SQR(x[0][0]) + SQR(x[0][1]) + SQR(x[0][2]) );
-        tmp_v = sqrt( SQR(v[0][0]) + SQR(v[0][1]) + SQR(v[0][2]) );
-        tmp_a = sqrt( SQR(a[0][0]) + SQR(a[0][1]) + SQR(a[0][2]) );
-        tmp_adot = sqrt( SQR(adot[0][0]) + SQR(adot[0][1]) + SQR(adot[0][2]) );
+        double tmp_r = sqrt( SQR(x[0][0]) + SQR(x[0][1]) + SQR(x[0][2]) );
+        double tmp_v = sqrt( SQR(v[0][0]) + SQR(v[0][1]) + SQR(v[0][2]) );
+        double tmp_a = sqrt( SQR(a[0][0]) + SQR(a[0][1]) + SQR(a[0][2]) );
+        double tmp_adot = sqrt( SQR(adot[0][0]) + SQR(adot[0][1]) + SQR(adot[0][2]) );
         fprintf(out,"%.16E \t %.8E \t % .16E % .16E % .16E \t %.16E \t % .16E % .16E % .16E \t %.16E \t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t %.8E \n",
                 time_cur, m[0],
                 x[0][0], x[0][1], x[0][2], tmp_r,
@@ -785,19 +753,18 @@ void write_bh_data()
     }
 }
 
-void write_bh_nb_data()
+void write_bh_nb_data(double time_cur, double m[], double3 x[], double3 v[])
 {
-    int i_bh, nb = config->live_smbh_neighbor_number;
-    double tmp, tmp_r, tmp_v;
+    int nb = config->live_smbh_neighbor_number;
+    int    ind_sort[N_MAX];
+    double var_sort[N_MAX];
 
     out = fopen("bh_neighbors.dat", "a");
 
     /* 1st BH */
 
-    i_bh = 0;
-
     for (int i_bh=0; i_bh < config->live_smbh_count; i_bh++) {
-        for (i=0; i<N; i++) var_sort[i] = (x[i]-x[i_bh]).norm();
+        for (int i=0; i<N; i++) var_sort[i] = (x[i]-x[i_bh]).norm();
         std::iota(ind_sort, ind_sort+N, 0);
         std::partial_sort(ind_sort, ind_sort + nb, ind_sort + N, [&](int i, int j) {return var_sort[i] < var_sort[j];});
 
@@ -809,7 +776,7 @@ void write_bh_nb_data()
                 v[i_bh][0], v[i_bh][1], v[i_bh][2]);
 
         for (int j=1; j < nb; j++) {
-            i = ind_sort[j];
+            int i = ind_sort[j];
             fprintf(out,"%02d %07d   %.8E   % .8E % .8E % .8E   %.8E   % .8E % .8E % .8E   %.8E \t",
                     j, i,
                     m[i],
@@ -1072,7 +1039,7 @@ for (i=0; i<ni; i++)
 
 }
 
-void calc_self_grav(double t)
+void calc_self_grav(double t, double eps2, double &g6_calls)
 {
     /* calc the new grav for the active particles */
 
@@ -1082,11 +1049,11 @@ void calc_self_grav(double t)
 
     g6_set_ti(clusterid, t);
 
-    ni = n_act;
+    int ni = n_act; // TODO why is this needed?
 
     /* define the local phi, a, adot for these active particles */
 
-    for (i=0; i<ni; i+=npipe) {
+    for (int i=0; i<ni; i+=npipe) {
         nn = npipe;
         if (ni-i < npipe) nn = ni - i;
         for (ii=0; ii<nn; ii++) {
@@ -1100,8 +1067,8 @@ void calc_self_grav(double t)
 
     /* Store the new value of the local partial force etc... */
 
-    for (i=0; i<n_act; i++) {
-        iii = ind_act[i];
+    for (int i=0; i<n_act; i++) {
+        int iii = ind_act[i];
         pot_act_tmp_loc[iii] = pot_act_tmp[i];
         a_act_tmp_loc[iii] = a_act_tmp[i];
         adot_act_tmp_loc[iii] = adot_act_tmp[i];
@@ -1370,34 +1337,35 @@ DT_EXT_GRAV += (CPU_tmp_user - CPU_tmp_user0);
 #endif        // EXTPOT
 }
 
-void energy_contr()
+void energy_contr(const double time_cur, const double timesteps, int n_act_sum, int 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) // TODO how do we know N?
 {
-    E_pot = 0.0;
-    for (i=0; i<N; i++) E_pot += m[i]*pot[i];
+    // TODO maybe use static variables here for the previous step energy?
+    double E_pot = 0.0;
+    for (int i=0; i<N; i++) E_pot += m[i]*pot[i];
     E_pot *= 0.5;
 
-    E_kin = 0.0;
-    for (i=0; i<N; i++) E_kin += m[i]*( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
+    double E_kin = 0.0;
+    for (int i=0; i<N; i++) E_kin += m[i]*( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
     E_kin *= 0.5;
 
-    E_pot_ext = 0.0;
+    double E_pot_ext = 0.0;
 #ifdef EXTPOT
     for (i=0; i<N; i++) E_pot_ext += m[i]*pot_ext[i];
 #endif
 
-    if (timesteps == 0.0) {
+    if (timesteps == 0.0) { //TODO what is this thing?
         rcm_sum = 0.0;
         vcm_sum = 0.0;
     }
 
-    mcm     = 0.0;
-    rcm_mod = 0.0;
-    vcm_mod = 0.0;
+    double mcm     = 0.0;
+    double rcm_mod = 0.0;
+    double vcm_mod = 0.0;
 
-    xcm = {0, 0, 0};
-    vcm = {0, 0, 0};
+    double3 xcm = {0, 0, 0};
+    double3 vcm = {0, 0, 0};
 
-    for (i=0; i<N; i++) {
+    for (int i=0; i<N; i++) {
         mcm += m[i];
         xcm += m[i] * x[i];
         vcm += m[i] * v[i];
@@ -1414,7 +1382,7 @@ void energy_contr()
 
     mom = {0, 0, 0};
 
-    for (i=0; i<N; i++) {
+    for (int i=0; i<N; i++) {
         mom[0] += m[i] * ( x[i][1] * v[i][2] - x[i][2] * v[i][1] );
         mom[1] += m[i] * ( x[i][2] * v[i][0] - x[i][0] * v[i][2] );
         mom[2] += m[i] * ( x[i][0] * v[i][1] - x[i][1] * v[i][0] );
@@ -1422,11 +1390,10 @@ void energy_contr()
 
     get_CPU_time(&CPU_time_real, &CPU_time_user, &CPU_time_syst);
 
-    E_tot = E_pot + E_kin + E_pot_ext;
-
-    E_tot_corr    = E_tot + E_corr;
-    E_tot_corr_sd = E_tot + E_corr + E_sd;
-
+    double E_tot = E_pot + E_kin + E_pot_ext;
+    double E_tot_corr    = E_tot /* + E_corr */; // TODO what is this?
+    
+    
     if (timesteps == 0.0) {
         /* initialize the E_tot_0 + etc... */
 
@@ -1439,6 +1406,7 @@ void energy_contr()
         /* read the E_tot_0 + etc... */
         if (skip_con == 0) {
             inp = fopen("contr.con","r");
+            double tmp;
             fscanf(inp,"%lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE",
                    &tmp, &tmp, &tmp, &tmp,
                    &tmp, &tmp, &tmp, &tmp, &tmp, &tmp, &tmp,
@@ -1450,9 +1418,7 @@ void energy_contr()
             skip_con = 1;
         }
     }
-    DE_tot         = E_tot         - E_tot_0;
-    DE_tot_corr    = E_tot_corr    - E_tot_corr_0;
-    DE_tot_corr_sd = E_tot_corr_sd - E_tot_corr_sd_0;
+    double DE_tot         = E_tot         - E_tot_0;
 
     printf("%.3E %.3E  % .4E %.4E % .4E  % .4E % .4E  %.2E\n",
             time_cur, timesteps,
@@ -1483,11 +1449,31 @@ void energy_contr()
 
     E_tot_0 = E_tot;
     E_tot_corr_0 = E_tot_corr;
-    E_tot_corr_sd_0 = E_tot_corr_sd;
+    E_tot_corr_sd_0 = 0;
 }
 
 int main(int argc, char *argv[])
 {
+
+// These variables are moved from global
+int    name_proc, n_proc=1, myRank=0, rootRank=0, cur_rank,
+    i, j, k, nj, diskstep=0, power, jjj, iii,
+    skip_con=0, tmp_i;
+
+double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
+    dt_bh, t_bh=0.0, dt_bh_tmp,
+    t_end, time_cur, dt_min, dt_max, min_t, min_t_loc, dt_new,
+    eta_s, eta, eta_bh,
+    E_tot_0, E_tot_corr_0, E_tot_corr_sd_0,
+    rcm_sum=0.0, vcm_sum=0.0,
+    eps=0.0, eps2,
+    a2_mod, adot2_mod,
+    dt_tmp, dt2half, dt3over6, dt4over24, dt5over120,
+    dtinv, dt2inv, dt3inv,
+    a1abs, adot1abs, a2dot1abs, a3dot1abs,
+    timesteps=0.0, n_act_sum=0.0, n_act_distr[N_MAX], g6_calls=0.0, g6_calls_sum=0.0,
+    tmp;
+
     /* INIT the rand() !!! */
     srand(19640916);                 /* it is just my birthday :-) */
 
@@ -1584,7 +1570,7 @@ int main(int argc, char *argv[])
 
         /* read the global data for particles to the rootRank */
 
-        read_data();
+        read_data(&diskstep, &N, &time_cur, ind, m, x, v);
 
         /* possible coordinate & velocity limits for ALL particles !!! */
 
@@ -1867,7 +1853,7 @@ int main(int argc, char *argv[])
     // NOTE this is where calc_self_grav_zero() used to be. Hopefully the copying from *_act to *_act_new is just a temporary measure.
     std::copy(&x_act[0][0], &x_act[0][0]+3*N, &x_act_new[0][0]);
     std::copy(&v_act[0][0], &v_act[0][0]+3*N, &v_act_new[0][0]);
-    calc_self_grav(time_cur);
+    calc_self_grav(time_cur, eps2, g6_calls);
 
     /* Wait to all processors to finish his works... */
     MPI_Barrier(MPI_COMM_WORLD);
@@ -1946,7 +1932,7 @@ int main(int argc, char *argv[])
                                     m_bh2, x_bh2, v_bh2, s_bh2,
                                     C_NB,  dt_bh_tmp, usedOrNot,
                                     (double(*)[3])a_pn1, (double(*)[3])adot_pn1,
-                                    (double(*)[3])a_pn2, (double(*)[3])adot_pn2);
+                                    (double(*)[3])a_pn2, (double(*)[3])adot_pn2, myRank, rootRank);
 
             for (k=0;k<3;k++) {
                 a[i_bh1][k] += a_pn1[1][k] + a_pn1[2][k] + a_pn1[3][k] + a_pn1[4][k] + a_pn1[5][k] + a_pn1[6][k];
@@ -1978,7 +1964,7 @@ int main(int argc, char *argv[])
 
     /* Energy control... */
     if (myRank == rootRank) {
-        energy_contr();
+        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);
     } /* if (myRank == rootRank) */
 
 #ifdef ETICS_DUMP
@@ -2075,10 +2061,10 @@ int main(int argc, char *argv[])
     if (myRank == rootRank) {
 
         /* Write BH data... */
-        if (config->live_smbh_output) write_bh_data();
+        if (config->live_smbh_output) write_bh_data(time_cur, m, x, v);
 
         /* Write BH NB data... */
-        if (config->live_smbh_neighbor_output) write_bh_nb_data();
+        if (config->live_smbh_neighbor_output) write_bh_nb_data(time_cur, m, x, v);
 
     } /* if (myRank == rootRank) */
 
@@ -2087,7 +2073,6 @@ int main(int argc, char *argv[])
     if (myRank == rootRank) {
         get_CPU_time(&CPU_time_real0, &CPU_time_user0, &CPU_time_syst0);
         get_CPU_time(&CPU_time_real,  &CPU_time_user,  &CPU_time_syst);
-        tmp_cpu = CPU_time_real-CPU_time_real0;
     } /* if (myRank == rootRank) */
 
     timesteps = 0.0;
@@ -2275,7 +2260,7 @@ int main(int argc, char *argv[])
         DT_ACT_PRED += (CPU_tmp_user - CPU_tmp_user0);
 #endif
 
-        calc_self_grav(min_t);
+        calc_self_grav(min_t, eps2, g6_calls);
 
         /* Reduce the "global" vectors from "local" on all the nodes */
 
@@ -2351,7 +2336,7 @@ int main(int argc, char *argv[])
                                         m_bh2, x_bh2, v_bh2, s_bh2,
                                         C_NB,  dt_bh_tmp, usedOrNot,
                                         (double(*)[3])a_pn1, (double(*)[3])adot_pn1,
-                                        (double(*)[3])a_pn2, (double(*)[3])adot_pn2);
+                                        (double(*)[3])a_pn2, (double(*)[3])adot_pn2, myRank, rootRank);
 
                 a_act_new[i_bh1] += a_pn1[1] + a_pn1[2] + a_pn1[3] + a_pn1[4] + a_pn1[5] + a_pn1[6];
                 a_act_new[i_bh2] += a_pn2[1] + a_pn2[2] + a_pn2[3] + a_pn2[4] + a_pn2[5] + a_pn2[6];
@@ -2594,10 +2579,10 @@ int main(int argc, char *argv[])
         if (time_cur >= t_bh) {
             if (myRank == rootRank) {
                 /* Write BH data... */
-                if (config->live_smbh_output) write_bh_data();
+                if (config->live_smbh_output) write_bh_data(time_cur, m, x, v);
 
                 /* Write BH NB data... */
-                if (config->live_smbh_neighbor_output) write_bh_nb_data();
+                if (config->live_smbh_neighbor_output) write_bh_nb_data(time_cur, m, x, v);
 
             } /* if (myRank == rootRank) */
 
@@ -2607,11 +2592,11 @@ int main(int argc, char *argv[])
         if (time_cur >= t_contr) {
             if (myRank == rootRank) {
 
-                energy_contr();
+                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);
 
                 /* write cont data */
 
-                write_cont_data();
+                write_snap_data((char*)"data.con", diskstep, N, time_cur, ind, m, x, v);
 
                 /* possible OUT for timing !!! */
 
@@ -2665,7 +2650,9 @@ int main(int argc, char *argv[])
         if (time_cur >= t_disk) {
             if (myRank == rootRank) {
                 diskstep++;
-                write_snap_data();
+                char out_fname[256];
+                sprintf(out_fname, "%06d.dat", diskstep);
+                write_snap_data(out_fname, diskstep, N, time_cur, ind, m, x, v);
             } /* if (myRank == rootRank) */
 
 #ifdef ETICS_DUMP
diff --git a/pn_bh_spin.c b/pn_bh_spin.c
index 54822f3..c0e54da 100644
--- a/pn_bh_spin.c
+++ b/pn_bh_spin.c
@@ -10,7 +10,7 @@ int calc_force_pn_BH(double m1, double xx1[], double vv1[], double spin1[],
                      double CCC_NB, double dt_bh, 
                      int usedOrNot[], 
                      double a_pn1[][3], double adot_pn1[][3],
-                     double a_pn2[][3], double adot_pn2[][3])
+                     double a_pn2[][3], double adot_pn2[][3], int myRank, int rootRank)
 {
 
 /*
@@ -52,7 +52,7 @@ double PI2 = 9.86960440108935;
 double c_1, c_2, c_4, c_5, c_6, c_7, RS_DIST;
 
 double M, eta, r, r2, r3, MOR;
-double V1_V22,VWHOLE, RP, RPP, VA;
+double V1_V22, RP, RPP, VA;
 double N[3], x[3], v[3], A[3];
 double A1, B1, A2, B2, A2_5, B2_5, AK2, BK2, AK4, BK4, AK5, BK5;
 double A1D, A2D, A2_5D, B1D, B2D, B2_5D, ADK2, BDK2, ADK4, BDK4, ADK5, BDK5;
@@ -72,7 +72,6 @@ double SS1aux[3],SS2aux[3],SU[3],SV[3],XAD[3],XSD[3];
 double NDOTCS[3], NCSU[3], NDOTCSIG[3], NCSV[3], ACS[3], VCSU[3], ACSIG[3], VCSV[3], SNVDOT,  
        SIGNVDOT, NSDOT, NSIGDOT, VSDOT, VSIGDOT, NXSDOT,  NXADOT;
 double C1_5D[3],C2D[3], C2_5D[3];
-double ADK, BDK, AD[3], KSAK, KSBK;
 double nu, Spin1Abs2, Spin2Abs2, rS1, rS2, S1Dir[3], S2Dir[3], QM[3];
 double Spin1Abs, Spin2Abs, QMAux2_1[3], QMAux2_2[3], QMAux1[3] , QMD[3], SPINPrev[3][2], SpinPrev2_1, 
        SpinPrev2_2, SPSPP1, SPSPP2, Spin1AbsNew2, Spin2AbsNew2, Spin1AbsNew, Spin2AbsNew, S1DirNew[3], 
@@ -129,7 +128,6 @@ r3 = r2*r;
 
 MOR    = M/r;
 V1_V22 = v[0]*v[0]+v[1]*v[1]+v[2]*v[2];
-VWHOLE = sqrt(V1_V22);
 RP     = (x[0]*v[0]+x[1]*v[1]+x[2]*v[2])/r;
 
 
@@ -667,15 +665,6 @@ if(usedOrNot[5] == 1)		// PN3.5	~1/c^7
     } /* if(Van_Spin==1) */
 
 
-   ADK = ADK2+ADK4+ADK5+ADK6+ADK7;
-   BDK = BDK2+BDK4+BDK5+BDK6+BDK7;
-
-   KSAK = AK2+AK4+AK5+AK6+AK7;
-   KSBK = BK2+BK4+BK5+BK6+BK7;
-
-   for(k=0;k<3;k++) AD[k] = -2.0*MOR*RP*(KSAK*N[k]+KSBK*v[k])/r2 + MOR*(ADK*N[k]+BDK*v[k])/r + MOR*(KSAK*(v[k]-N[k]*RP)/r+KSBK*AT[k])/r + C1_5D[k]/c_2 + C2D[k]/c_4 +C2_5D[k]/c_4 + QMD[k]/c_4;
-
-
 for(k=0;k<3;k++)	// new values of the BH's spins, returned back to the main program... 
   {
   spin1[k] = SPIN[k][0];