diff --git a/phigrape.cpp b/phigrape.cpp
index 681089e..23e7bf4 100644
--- a/phigrape.cpp
+++ b/phigrape.cpp
@@ -562,6 +562,110 @@ public:
     double *m, *t, *dt, *pot;
 };
 
+inline void calc_high_derivatives(const double dt_tmp, const double3 a_act, const double3 a_act_new, const double3 adot_act, const double3 adot_act_new, double3 *a2, double3 *a3)
+{
+    double dtinv = 1.0/dt_tmp;
+    double dt2inv = dtinv*dtinv;
+    double dt3inv = dt2inv*dtinv;
+
+    double3 a0mia1 = a_act-a_act_new;
+    double3 ad04plad12 = 4.0*adot_act + 2.0*adot_act_new;
+    double3 ad0plad1 = adot_act + adot_act_new;
+
+    *a2 = -6.0*a0mia1*dt2inv - ad04plad12*dtinv;
+    *a3 = 12.0*a0mia1*dt3inv + 6.0*ad0plad1*dt2inv;
+}
+
+inline void corrector(const double dt_tmp, const double3 a2, const double3 a3, double3* x_act_new, double3* v_act_new)
+{
+    double dt3over6 = dt_tmp*dt_tmp*dt_tmp/6.0;
+    double dt4over24 = dt3over6*dt_tmp/4.0;
+    double dt5over120 = dt4over24*dt_tmp/5.0;
+
+    *x_act_new += dt4over24*a2 + dt5over120*a3;
+    *v_act_new += dt3over6*a2 + dt4over24*a3;
+}
+
+inline double aarseth_step(const double eta, const double dt_tmp, const double3 a_act_new, const double3 adot_act_new, const double3 a2, const double3 a3)
+{
+    double a1abs = a_act_new.norm();
+    double adot1abs = adot_act_new.norm();
+    double3 a2dot1 = a2 + dt_tmp*a3;
+    double a2dot1abs = a2dot1.norm();
+    double a3dot1abs = a3.norm();
+
+    return sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
+}
+
+void binary_smbh_influence_sphere_output(int i_bh1, int i_bh2, int ind_act[], double m_act[], double3 x_act[], double3 v_act[], double pot_act[], double dt_act[], int n_act, double timesteps, double time_cur,  double factor, int inf_event[])
+{
+    //TODO !!IMPORTANT!! only open the file IF THERE IS ANYTHING TO WRITE!!!
+    //TODO inf_event to be static or something?
+    auto out = fopen("bbh_inf.dat", "a");
+
+    double m_bh1 = m_act[i_bh1];
+    double m_bh2 = m_act[i_bh2];
+
+    double3 x_bh1 = x_act[i_bh1];
+    double3 x_bh2 = x_act[i_bh2];
+
+    double3 v_bh1 = v_act[i_bh1];
+    double3 v_bh2 = v_act[i_bh2];
+
+    double3 x_bbhc = (m_bh1*x_bh1 + m_bh2*x_bh2)/(m_bh1 + m_bh2);
+    double3 v_bbhc = (m_bh1*v_bh1 + m_bh2*v_bh2)/(m_bh1 + m_bh2);
+
+    double DR2 = SQR(x_bh1[0] - x_bh2[0]) + SQR(x_bh1[1] - x_bh2[1]) + SQR(x_bh1[2] - x_bh2[2]);
+    double DV2 = SQR(v_bh1[0] - v_bh2[0]) + SQR(v_bh1[1] - v_bh2[1]) + SQR(v_bh1[2] - v_bh2[2]);
+
+    double EB = -(m_bh1 + m_bh2) / sqrt(DR2) + 0.5 * DV2;
+
+    double SEMI_a = -0.5 * (m_bh1 + m_bh2)/EB;
+    double SEMI_a2 = SQR(SEMI_a);
+
+    for (int i=2; i<n_act; i++) {
+
+        double tmp_r2 = SQR(x_act[i][0] - x_bbhc[0]) + SQR(x_act[i][1] - x_bbhc[1]) + SQR(x_act[i][2] - x_bbhc[2]);
+
+        int iii = ind_act[i];
+
+        if (tmp_r2 < SEMI_a2*SQR(factor)) {
+
+            if (inf_event[iii] == 0) {
+
+                fprintf(out,"INF1 %.6E %.16E %07d %07d   %.6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E \n",
+                    timesteps, time_cur, i, ind_act[i],
+                    sqrt(DR2), x_bbhc[0], x_bbhc[1], x_bbhc[2], v_bbhc[0], v_bbhc[1], v_bbhc[2],
+                    m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_act[0],
+                    m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_act[1],
+                    sqrt(tmp_r2),
+                    m_act[i], x_act[i][0], x_act[i][1], x_act[i][2], v_act[i][0], v_act[i][1], v_act[i][2], pot_act[i],
+                    dt_act[i]);
+
+            inf_event[iii] = 1;
+            }
+
+        } else {
+            if (inf_event[iii] == 1) {
+                fprintf(out,"INF2 %.6E %.16E %07d %07d   %.6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E \n",
+                    timesteps, time_cur, i, ind_act[i],
+                    sqrt(DR2), x_bbhc[0], x_bbhc[1], x_bbhc[2], v_bbhc[0], v_bbhc[1], v_bbhc[2],
+                    m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_act[0],
+                    m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_act[1],
+                    sqrt(tmp_r2),
+                    m_act[i], x_act[i][0], x_act[i][1], x_act[i][2], v_act[i][0], v_act[i][1], v_act[i][2], pot_act[i],
+                    dt_act[i]);
+            }
+
+            inf_event[iii] = 0;
+
+        } /* if (tmp_r2 < DR2*R_INF2) */
+
+    } /* i */
+    fclose(out);
+
+}
+
 int main(int argc, char *argv[])
 {
     int name_proc, n_proc=1, myRank=0, rootRank=0, cur_rank,
@@ -1104,14 +1208,14 @@ int main(int argc, char *argv[])
         get_CPU_time(&CPU_time_real,  &CPU_time_user,  &CPU_time_syst);
     } /* if (myRank == rootRank) */
 
-    timesteps = 0.0;
+    timesteps = 0.0; // Why won't those two be long long instead of double + should include the zeroth step
     n_act_sum = 0.0;
 
     for (int i=1; i<N+1; i++) { // TODO why like this?
         n_act_distr[i-1] = 0.0;
     }
 
-    g6_calls = 0.0;
+    g6_calls = 0.0; //TODO this should include the calls at the zeroth step, so move it further up.
 
 #ifdef TIMING
     DT_TOT      = 0.0;
@@ -1351,51 +1455,41 @@ int main(int argc, char *argv[])
             // TODO split this loop into the three tasks it is doing, and remove the redundancy.
             dt_tmp = min_t - t_act[i];
 
-            dt3over6 = dt_tmp*dt_tmp*dt_tmp/6.0;
-            dt4over24 = dt3over6*dt_tmp/4.0;
-            dt5over120 = dt4over24*dt_tmp/5.0;
+            double3 a2, a3;
+            calc_high_derivatives(dt_tmp, a_act[i], a_act_new[i], adot_act[i], adot_act_new[i], &a2, &a3);
 
-            dtinv = 1.0/dt_tmp;
-            dt2inv = dtinv*dtinv;
-            dt3inv = dt2inv*dtinv;
+            corrector(dt_tmp, a2, a3, &x_act_new[i], &v_act_new[i]);
 
-            double3 a0mia1 = a_act[i]-a_act_new[i];
-            double3 ad04plad12 = 4.0*adot_act[i] + 2.0*adot_act_new[i];
-            double3 ad0plad1 = adot_act[i] + adot_act_new[i];
 
-            a2 = -6.0*a0mia1*dt2inv - ad04plad12*dtinv;
-            a3 = 12.0*a0mia1*dt3inv + 6.0*ad0plad1*dt2inv;
+            //TODO make beautiful
+            double eta_curr;
+            if ((config->live_smbh_count > 0) && (ind_act[i] < config->live_smbh_count)) eta_curr = eta_bh;
+            else eta_curr = eta;
 
-            x_act_new[i] += dt4over24*a2 + dt5over120*a3;
-            v_act_new[i] += dt3over6*a2 + dt4over24*a3;
-
-            a1abs = sqrt(a_act_new[i][0]*a_act_new[i][0]+a_act_new[i][1]*a_act_new[i][1]+a_act_new[i][2]*a_act_new[i][2]);
-            adot1abs = sqrt(adot_act_new[i][0]*adot_act_new[i][0]+adot_act_new[i][1]*adot_act_new[i][1]+adot_act_new[i][2]*adot_act_new[i][2]);
-
-            a2dot1 = a2 + dt_tmp*a3;
-
-            a2dot1abs = a2dot1.norm();
-            a3dot1abs = a3.norm();
-
-            if ((config->live_smbh_count > 0) && (ind_act[i] < config->live_smbh_count))
-                dt_new = sqrt(eta_bh*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
-            else
-                dt_new = sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
+            double dt_new = aarseth_step(eta_curr, dt_tmp, a_act_new[i], adot_act_new[i], a2, a3);
 
+            //TODO the below should be moved to a function
             if (dt_new < dt_min) dt_tmp = dt_min;
 
             if ((dt_new < dt_tmp) && (dt_new > dt_min)) {
                 power = log(dt_new)/log(2.0) - 1;
-
-                dt_tmp = pow(2.0, (double)power);
+                dt_tmp = pow(2.0, (double)power); // TODO why is this casting needed here?
             }
 
-            if ((dt_new > 2.0*dt_tmp) && (fmod(min_t, 2.0*dt_tmp) == 0.0) && (2.0*dt_tmp <= dt_max)) {
-                dt_tmp *= 2.0;
+            if ((dt_new > 2*dt_tmp) && (fmod(min_t, 2*dt_tmp) == 0) && (2*dt_tmp <= dt_max)) {
+                dt_tmp *= 2;
             }
 
+            if (config->dt_min_warning && (myRank == 0)) {
+                if (dt_act[i] == dt_min) {
+                    printf("!!! Warning1: dt_act = dt_min = %.6E \t ind_act = %07d \n", dt[i], ind_act[i]);
+                    fflush(stdout);
+                }
+            }
+
+            /* BEGIN copy of everything */
+
             dt_act[i] = dt_tmp;
-
             t_act[i] = min_t;
 
             pot_act[i] = pot_act_new[i];
@@ -1404,13 +1498,10 @@ int main(int argc, char *argv[])
             v_act[i]    =    v_act_new[i];
             a_act[i]    =    a_act_new[i];
             adot_act[i] = adot_act_new[i];
+            
+            /* END copy of everything */
+
 
-            if (config->dt_min_warning && (myRank == 0)) {
-                if (dt_act[i] == dt_min) {
-                    printf("!!! Warning1: dt_act = dt_min = %.6E \t ind_act = %07d \n", dt[i], ind_act[i]);
-                    fflush(stdout);
-                }
-            }
 
         } /* i */
 
@@ -1422,74 +1513,9 @@ int main(int argc, char *argv[])
             if (config->live_smbh_count==2) dt_act[i_bh2] = min_dt;
         }
 
-        if (config->binary_smbh_influence_sphere_output) {
-            if (myRank == rootRank) {
-
-                out = fopen("bbh_inf.dat", "a");
-
-                m_bh1 = m_act[i_bh1];
-                m_bh2 = m_act[i_bh2];
-
-                x_bh1 = x_act[i_bh1];
-                x_bh2 = x_act[i_bh2];
-
-                v_bh1 = v_act[i_bh1];
-                v_bh2 = v_act[i_bh2];
-
-                x_bbhc = (m_bh1*x_bh1 + m_bh2*x_bh2)/(m_bh1 + m_bh2);
-                v_bbhc = (m_bh1*v_bh1 + m_bh2*v_bh2)/(m_bh1 + m_bh2);
-
-                DR2 = SQR(x_bh1[0] - x_bh2[0]) + SQR(x_bh1[1] - x_bh2[1]) + SQR(x_bh1[2] - x_bh2[2]);
-                DV2 = SQR(v_bh1[0] - v_bh2[0]) + SQR(v_bh1[1] - v_bh2[1]) + SQR(v_bh1[2] - v_bh2[2]);
-
-                EB = -(m_bh1 + m_bh2) / sqrt(DR2) + 0.5 * DV2;
-
-                SEMI_a = -0.5 * (m_bh1 + m_bh2)/EB;
-                SEMI_a2 = SQR(SEMI_a);
-
-                for (int i=2; i<n_act; i++) {
-
-                    tmp_r2 = SQR(x_act[i][0] - x_bbhc[0]) + SQR(x_act[i][1] - x_bbhc[1]) + SQR(x_act[i][2] - x_bbhc[2]);
-
-                    iii = ind_act[i];
-
-                    if (tmp_r2 < SEMI_a2*SQR(config->binary_smbh_influence_radius_factor)) {
-
-                        if (inf_event[iii] == 0) {
-
-                            fprintf(out,"INF1 %.6E %.16E %07d %07d   %.6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E \n",
-                                timesteps, time_cur, i, ind_act[i],
-                                sqrt(DR2), x_bbhc[0], x_bbhc[1], x_bbhc[2], v_bbhc[0], v_bbhc[1], v_bbhc[2],
-                                m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_act[0],
-                                m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_act[1],
-                                sqrt(tmp_r2),
-                                m_act[i], x_act[i][0], x_act[i][1], x_act[i][2], v_act[i][0], v_act[i][1], v_act[i][2], pot_act[i],
-                                dt_act[i]);
-
-                        inf_event[iii] = 1;
-                        }
-
-                    } else {
-                        if (inf_event[iii] == 1) {
-                            fprintf(out,"INF2 %.6E %.16E %07d %07d   %.6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E \n",
-                                timesteps, time_cur, i, ind_act[i],
-                                sqrt(DR2), x_bbhc[0], x_bbhc[1], x_bbhc[2], v_bbhc[0], v_bbhc[1], v_bbhc[2],
-                                m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_act[0],
-                                m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_act[1],
-                                sqrt(tmp_r2),
-                                m_act[i], x_act[i][0], x_act[i][1], x_act[i][2], v_act[i][0], v_act[i][1], v_act[i][2], pot_act[i],
-                                dt_act[i]);
-                        }
-
-                        inf_event[iii] = 0;
-
-                    } /* if (tmp_r2 < DR2*R_INF2) */
-
-                } /* i */
-
-                fclose(out);
-
-            } /* if (myRank == rootRank) */
+        if (config->binary_smbh_influence_sphere_output && (myRank == rootRank)) {
+            //TODO clean up this mass. We don't want to have all these _act arrays; they are only needed for this lame function.
+            binary_smbh_influence_sphere_output(i_bh1, i_bh2, ind_act, m_act, x_act, v_act, pot_act, dt_act, n_act, timesteps, time_cur, config->binary_smbh_influence_radius_factor, inf_event);
         }
 
         /* Return back the new coordinates + etc... of active part. to the global data... */