Added #include <stdexcept>

Enabled triple (plus) SMBHs
Fixed PN force adjustment before PN parameters set
2022-05-30 10:28:12 -04:00 · 2022-05-19 20:35:24 -04:00 · 2022-01-14 17:54:38 -05:00 · 2021-12-26 21:40:50 -05:00 · 2021-12-26 17:20:20 -05:00 · 2021-08-30 09:31:53 -04:00
29 changed files with 2171 additions and 16819 deletions
--- a/.gitignore
+++ b/.gitignore
@ -6,3 +6,9 @@
 *.cfg
 *.inp
 *.mask
+grapite-dev-exec-threshold
+phigrape
+*.h5
+.*
+CUDA
+.gitignore
--- a/44
+++ b/44
@ -1,30 +1,34 @@
-CUDAHOME     ?= /usr/local/cuda
-CPPFLAGS     += -DYEBISU -DETICS
+CUDA_HOME    ?= /usr/local/cuda
+CPPFLAGS     += -DETICS
 OPTIMIZATION ?= 3

-ETICS_DTSCF  ?= 0.015625
+CUDAINC       = -I$(CUDA_HOME)/include -I$(CUDA_HOME)/samples/common/inc/
+CUDALIB       = -L$(CUDA_HOME)/lib64 -lcudart -lcudadevrt -lcuda

-CUDAINC       = -I$(CUDAHOME)/include -I$(CUDAHOME)/samples/common/inc/
-CUDALIB       = -L$(CUDAHOME)/lib64 -lcudart -lcudadevrt
+default: grapite
+grapite: GRAPE_HOME = ../grapite
+grapite: GRAPELIB   = -L$(GRAPE_HOME) -lgrapite
+yebisu:  GRAPE_HOME = ../yebisu
+yebisu:  GRAPELIB   = -L$(GRAPE_HOME) -lyebisug6
+sapporo: GRAPE_HOME = ../sapporo2/lib
+sapporo: GRAPELIB   = -L$(GRAPE_HOME) -lsapporo
+         GRAPEINC   = -I$(GRAPE_HOME)

-GRAPEHOME     = ../grapite
-GRAPELIB      = -L$(GRAPEHOME) -lgrapite
-yebisu: GRAPEHOME = ../yebisu
-yebisu: GRAPELIB  = -L$(GRAPEHOME) -lyebisug6
-GRAPEINC      = -I$(GRAPEHOME)
-
-CFLAGS       ?= -mcmodel=large
-CFLAGS       += -O$(OPTIMIZATION)
+CXXFLAGS     += -std=c++11 -O$(OPTIMIZATION)
 INC           = $(GRAPEINC) $(CUDAINC)
-LIB           = $(GRAPELIB) $(CUDALIB) -lm -lgcc -lgfortran -lstdc++
-MPICC        ?= mpicc
-EXECUTABLE   ?= phi-GRAPE.exe
+LIB           = $(GRAPELIB) $(CUDALIB) -lm
+MPICXX       ?= mpic++
+EXECUTABLE   ?= phigrape
+
+# HDF5
+#CPPFLAGS     += -DHAS_HDF5
+#LIB          += -lhdf5 -lz -ldl

 default:
-	$(MPICC) $(CPPFLAGS) $(CFLAGS) -DETICS_DTSCF=$(ETICS_DTSCF) $(INC) phi-GRAPE.c -o $(EXECUTABLE) $(LIB)
+	$(MPICXX) $(CPPFLAGS) $(CXXFLAGS) $(INC) black_holes.cpp external.cpp io.cpp config.cpp phigrape.cpp -o $(EXECUTABLE) $(LIB)

-yebisu: CPPFLAGS := $(filter-out -DETICS, $(CPPFLAGS))
-yebisu: default
+yebisu sapporo: CPPFLAGS := $(filter-out -DETICS, $(CPPFLAGS))
+yebisu sapporo grapite: default

 clean:
-	rm -f *.o phi-GRAPE.exe
+	rm -f CUDA *.o phigrape
--- a/TODO.md
+++ b/TODO.md
@ -0,0 +1,8 @@
+TODO
+====
+
+* Break main() into smaller chunks; operations that are timed should become independent functions.
+
+* Const everything
+
+* OpenMP
--- a/act_def_linklist.c
+++ b/act_def_linklist.c
@ -1,376 +0,0 @@
-/**************************************************************
-   File     : act_def_linklist.c
-   Func.    : provide linear linking list functions
-            : for active particle def.
-   CODED BY : Zhong Shiyan
-   START    : 2014-03-28, 12:30
-**************************************************************/
-
-//***********************************************************//
-/*             Definition of T/P-node                        */
-//***********************************************************//
-typedef struct PNODE
-{
-    int    Pid;                  // Particle's real ID
-    struct PNODE *NextPNODE;
-} PNODE;
-
-typedef struct TNODE
-{
-    double t_node;               // t_node = t + dt
-    int    n_node;               // number of P-nodes under this node
-
-    struct PNODE *PartList, *PartListEnd;
-    struct TNODE *NextTNODE;
-} TNODE;
-
-struct TNODE *CurrT=NULL;
-
-//***********************************************************//
-/*             Operations on T/P-node                        */
-//***********************************************************//
-struct TNODE *CreateTNODE( double t ){
-    
-    struct TNODE *ptr;
-    
-    ptr = (TNODE*)malloc(sizeof(*ptr));
-    if( ptr == NULL ){
-	printf("Fail to create a node.");
-	exit(-1);
-    }
-    
-    ptr->t_node = t;
-    ptr->n_node = 0;
-    
-    ptr->NextTNODE = NULL;
-    
-    ptr->PartList    = NULL;
-    ptr->PartListEnd = NULL;
-    
-    return ptr;
-}
-//***********************************************************//
-struct PNODE *DeletePNODE( struct PNODE *ptr ){
-    struct PNODE *next;
-    
-    next = ptr->NextPNODE;
-    free(ptr);
-    
-    return next;
-}
-//***********************************************************//
-void DeleteTNODE( struct TNODE *Tptr ){
-    
-    struct PNODE *Pptr;
-    
-    Pptr = Tptr->PartList;
-    
-    while( Pptr != NULL ){
-      Pptr = DeletePNODE( Pptr );
-    }
-    
-    free(Tptr);
-    return;
-}
-//***********************************************************//
-struct PNODE *CreatePNODE( int id ){
-    struct PNODE *ptr;
-    
-    ptr = (PNODE*)malloc(sizeof(*ptr));
-    if( ptr == NULL ){
-	printf("Fail to create a P-node.");
-	return NULL;
-    }
-    
-    ptr->Pid = id;
-    ptr->NextPNODE = NULL;
-    
-    return ptr;
-
-}
-//***********************************************************//
-void InsertTNODE( struct TNODE *front, 
-                  struct TNODE *rear, 
-                  struct TNODE *ptr )
-{    
-    front->NextTNODE = ptr;
-    ptr->NextTNODE = rear;
-
-}
-//***********************************************************//
-
-
-
-
-//***********************************************************//
-/*             Functions about act_def                       */
-//***********************************************************//
-/**************************************************************
-This function only called 1 time, at beginning of simulation.
-After all particle's dt are computed
-**************************************************************/
-
-void CreateLinkList(){
-    
-    struct TNODE *head,*tail,*Tp1,*Tp2,*newTNODE,*Tptr;
-    struct PNODE *Pptr,*Ptail;
-    
-    int i, iii;
-    double ttmp, t1, t2;
-    
-    head = CreateTNODE( -1.0 );    // will be discarded at the end of this routine
-    tail = CreateTNODE( 2.0*t_end );
-    
-    head->NextTNODE = tail;
-    tail->NextTNODE = NULL;
-    
-    // building link list
-    
-for(i=0;i<N;i++){
-      
-      ttmp = t[i] + dt[i];
-      iii = ind[i];
-      if( m[iii] == 0.0 ) ttmp = t_end + 0.125;
-      
-      Tp1 = head;
-      Tp2 = Tp1->NextTNODE;
-      t1 = Tp1->t_node;
-      t2 = Tp2->t_node;
-      
-      Pptr = CreatePNODE( iii );
-      
-    while( Tp1->NextTNODE != NULL ){
-      if( ttmp == t1 ){                 // if T-node exist, add a P-node
-        
-        if( Tp1->PartListEnd == NULL ){ // This is the first P-node under current T-node
-          Tp1->PartList    = Pptr;
-          Tp1->PartListEnd = Pptr;
-          Tp1->n_node = Tp1->n_node + 1;
-        }
-        else{                            // There are already many P-nodes under this T-node
-          Ptail = Tp1->PartListEnd;
-          Ptail->NextPNODE = Pptr;
-          Tp1->PartListEnd = Pptr;
-          Tp1->n_node = Tp1->n_node + 1;
-        }
-        break; // jump out of this "while" loop
-      }
-      
-      if( ttmp > t1 && ttmp < t2 ){ // Create a new T-node and insert between *Tp1 and *Tp2, then add P-node to it
-        
-        newTNODE = CreateTNODE(ttmp);
-        InsertTNODE( Tp1, Tp2, newTNODE);
-        
-        newTNODE->n_node = 1;
-        newTNODE->PartList = Pptr;
-        newTNODE->PartListEnd = Pptr;
-        break; // jump out of this "while" loop
-      }
-      
-      // move to next T-node
-      Tp1 = Tp1->NextTNODE;
-      t1 = Tp1->t_node;
-      
-      if(Tp2->NextTNODE != NULL){
-        Tp2 = Tp2->NextTNODE;
-        t2  = Tp2->t_node;
-      }
-      else{ break; }
-    
-    }// while( Tp1->NextTNODE != NULL )
-    
-}//for(i=0;i<N;i++)
-
-      CurrT = head->NextTNODE;
-      free(head);
-      
-}
-
-//End of CreateLinkList()
-/**************************************************************/
-
-
-/**************************************************************
-This Function is used to modify the link list,
-after get new dt for active particles. Then point *CurrT to next
-T-node.
-**************************************************************/
-
-void ModifyLinkList(){
-    
-    struct TNODE *Tp1,*Tp2,*newTNODE;
-    struct PNODE *Pptr,*Ptail;
-    
-    int i, iii;
-    double ttmp, t1, t2;
-    
-  for(i=0;i<n_act;i++){
-    
-    iii = ind_act[i];
-    ttmp = t[iii] + dt[iii];
-    
-    Tp1 = CurrT;
-    Tp2 = Tp1->NextTNODE;
-    t1 = Tp1->t_node;
-    t2 = Tp2->t_node;
-      
-    Pptr = CreatePNODE( iii );
-      
-    while( Tp1->NextTNODE != NULL ){
-      if( ttmp == t1 ){                 // if T-node exist, add a P-node
-        
-        if( Tp1->PartListEnd == NULL ){ // This is the first P-node under current T-node
-          Tp1->PartList    = Pptr;
-          Tp1->PartListEnd = Pptr;
-          Tp1->n_node = Tp1->n_node + 1;
-        }
-        else{                            // There are already many P-nodes under this T-node
-          Ptail = Tp1->PartListEnd;
-          Ptail->NextPNODE = Pptr;
-          Tp1->PartListEnd = Pptr;
-          Tp1->n_node = Tp1->n_node + 1;
-        }
-        break; // jump out of this "while" loop
-      }
-      
-      if( ttmp > t1 && ttmp < t2 ){ // Create a new T-node and insert between *Tp1 and *Tp2, then add P-node to it
-        
-        newTNODE = CreateTNODE(ttmp);
-        InsertTNODE( Tp1, Tp2, newTNODE);
-        
-        newTNODE->n_node = 1;
-        newTNODE->PartList = Pptr;
-        newTNODE->PartListEnd = Pptr;
-        
-        break; // jump out of this "while" loop
-      }
-      
-      // move to next T-node
-      Tp1 = Tp1->NextTNODE;
-      t1  = Tp1->t_node;
-    
-      if(Tp2->NextTNODE != NULL){
-        Tp2 = Tp2->NextTNODE;
-        t2  = Tp2->t_node;
-      }
-      else{ break; }
-      
-    }//while( Tp1->NextTNODE != NULL )
-
-  }//for(i=0;i<n_act;i++)
-
-  Tp1 = CurrT;
-  CurrT = CurrT->NextTNODE;
-  DeleteTNODE( Tp1 );
-}
-
-//End of ModifyLinkList()
-/***************************************************************/
-
-
-/***************************************************************/
-/*
-void i_swap(int *a, int *b)
-{
-register int tmp;
-tmp = *a; *a = *b; *b = tmp;
-}
-*/
-/***************************************************************/
-/***************************************************************/
-void ind_act_sort(int l, int r)
-{
-
-int i, j, cikl, tmp;
-
-i = l; j = r;
-tmp = ind_act[(l+r)/2];
-
-cikl = 1;
-
-while(cikl)
-  {
-    while (ind_act[i]<tmp) i++;
-    while (tmp<ind_act[j]) j--;
-
-    if (i<=j)
-     {
-     i_swap(&ind_act[i],&ind_act[j]);
-     i++; j--;
-     }
-    else
-     {
-     cikl = 0;
-     }
-  }
-  if (l<j) ind_act_sort(l, j);
-  if (i<r) ind_act_sort(i, r);
-}
-/**************************************************************/
-
-/**************************************************************
-Get active particle list
-**************************************************************/
-void get_act_plist()
-  {
-  struct PNODE *Pptr;
-  int i, iii, j,k, itmp, ipt, flag;
-
-  char idcFile[30];
-  FILE *idcomp;
-
-  i=0;
-  Pptr  = CurrT->PartList;
-  n_act = 0;
-
-  min_t = CurrT->t_node; // IMPORTANT !!
-
-  flag = 0;
-
-  while(Pptr != NULL)
-    {
-    iii = Pptr->Pid;
-    Pptr = Pptr->NextPNODE;
-    if( m[iii] != 0.0 )			// Do not put zero mass part. in active plist
-      { 
-      ind_act[i] = iii;
-//      if(ind_act[i]==N-1) flag=1;
-      i++; n_act++;
-      }
-    }
-
-  if( n_act > 2 ) ind_act_sort( 0, n_act-1 );
-
-//  printf("last pid: %06d\n", ind_act_ll[n_act-1]);
-//  if( flag == 0 ) printf("Warning: BH not in the ind_act array!\n");
-//  if( ind_act[n_act-1]!= N-1) printf("Warning: BH not in the last of ind_act array!\n");
-
-  }// End of get_act_plist()
-/**************************************************************/
-
-
-#ifdef DEBUG_extra
-/**************************************************************
-Check link list
-**************************************************************/
-void check_linklist(int Tstep)
-  {
-  FILE *listf;
-  struct TNODE *Tptr;
-
-  listf = fopen("Check_linklist.dat","a");
-  Tptr = CurrT;
-  fprintf(listf,"Timesteps = %04d\n", Tstep);
-  while(Tptr != NULL)
-    {
-    fprintf(listf,"% 8E   %04d\n", Tptr->t_node, Tptr->n_node);
-    Tptr = Tptr->NextTNODE;
-    }
-
-  fprintf(listf,"========================\n\n");
-
-  fclose(listf);
-  }
-/**************************************************************/
-#endif
-
--- a/black_holes.cpp
+++ b/black_holes.cpp
@ -0,0 +1,180 @@
+#include <cstdio>
+#include <numeric>
+#include "black_holes.h"
+
+/* BEGIN legacy inclusion */
+// I'm not going to touch this C file
+#define SQR(x) ((x)*(x))
+double L[3]; // needed in pn_bh_spin.c
+#include "pn_bh_spin.c"
+#undef SQR
+/* END legacy inclusion */
+
+void two_body_gravity(const Particle_ref& j, const Particle_ref& i, const double eps, double& pot, double3& acc, double3& jrk)
+{
+    double3 dx = i.x - j.x;
+    double3 dv = i.v - j.v;
+    double r2 = dx.norm2() + eps*eps;
+    double r  = sqrt(r2);
+    double r3 = r2*r;
+    double r4 = r2*r2;
+    double RP = 3*(dx*dv)/r;
+
+    pot = -j.m/r;
+    acc = -j.m*dx/r3;
+    jrk = -j.m*(dv/r3 - RP*dx/r4);
+}
+
+void Black_hole_physics::adjust_softening(const std::vector<Particle_ref>& particles)
+{
+    if (eps_new < 0) return;
+    double pot_old, pot_new;
+    double3 acc_old, acc_new, jrk_old, jrk_new;
+    for (int j = 0; j < particles.size(); j++) {
+        for (int i = 0; i < particles.size(); i++) {
+            if (i == j) continue;
+            two_body_gravity(particles[j], particles[i], eps_old, pot_old, acc_old, jrk_old);
+            two_body_gravity(particles[j], particles[i], eps_new, pot_new, acc_new, jrk_new);
+            particles[i].pot += pot_new - pot_old;
+            particles[i].acc += acc_new - acc_old;
+            particles[i].jrk += jrk_new - jrk_old;
+        }
+    }
+}
+
+#if 0
+void Black_hole_physics::adjust_post_newtonian(
+    const double dt_bh, // pn_usage should be const
+    double3& acc1, double3& acc2,
+    double3& jrk1, double3& jrk2)
+{
+    // calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk
+    // TODO maybe have the PN terms as local variables here?
+    int tmp;
+    tmp = calc_force_pn_BH(masses[0], x1, v1, bbh_grav.spin1,
+                           masses[1], x2, v2, bbh_grav.spin2,
+                           c,  dt_bh, pn_usage,
+                           (double(*)[3])bbh_grav.a_pn1, (double(*)[3])bbh_grav.adot_pn1,
+                           (double(*)[3])bbh_grav.a_pn2, (double(*)[3])bbh_grav.adot_pn2, myRank, rootRank);
+    if (tmp == 505) exit(-1); // Very ugly way to terminate
+
+    // NOTE we have these _corr variables accumulating the corrections before
+    // applying it. It's almost the same but different from applying each
+    // correction term in a loop.
+    double3 acc1_corr(0,0,0), acc2_corr(0,0,0), jrk1_corr(0,0,0), jrk2_corr(0,0,0);
+    for (int i=1; i<7; i++) {
+        acc1_corr += bbh_grav.a_pn1[i];
+        acc2_corr += bbh_grav.a_pn2[i];
+        jrk1_corr += bbh_grav.adot_pn1[i];
+        jrk2_corr += bbh_grav.adot_pn2[i];
+    }
+    acc1 += acc1_corr;
+    acc2 += acc2_corr;
+    jrk1 += jrk1_corr;
+    jrk2 += jrk2_corr;
+}
+#endif
+
+void Black_hole_physics::write_bh_data(const double time_cur, const int count, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, const std::vector<double3> &a, const std::vector<double3> &adot, const std::vector<double> &dt)
+{
+    auto out = fopen("bh.dat", "a");
+    for (int i = 0; i < count; i++) {
+        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[i],
+            x[i][0], x[i][1], x[i][2], x[i].norm(),
+            v[i][0], v[i][1], v[i][2], v[i].norm(),
+            pot[i],
+            a[i][0],    a[i][1],    a[i][2],    a[i].norm(),
+            adot[i][0], adot[i][1], adot[i][2], adot[i].norm(),
+            dt[i]);
+    }
+    fprintf(out, "\n");
+    fclose(out);
+}
+
+void Write_bh_nb_data::operator()(double time_cur)
+{
+    for (int i_bh=0; i_bh < smbh_count; i_bh++) {
+        for (int i=0; i<N; i++) var_sort[i] = (x[i]-x[i_bh]).norm();
+        std::iota(ind_sort.begin(), ind_sort.end(), 0);
+        std::partial_sort(ind_sort.begin(), ind_sort.begin() + nb, ind_sort.begin() + N, [&](int i, int j) {return var_sort[i] < var_sort[j];});
+
+        fprintf(out,"%.16E \t %07d \t %.8E \t % .8E % .8E % .8E \t % .8E % .8E % .8E \t",
+                time_cur,
+                i_bh,
+                m[i_bh],
+                x[i_bh][0], x[i_bh][1], x[i_bh][2],
+                v[i_bh][0], v[i_bh][1], v[i_bh][2]);
+
+        for (int j=1; j < nb; j++) {
+            int i = ind_sort[j];
+            fprintf(out,"%02d %07d   %.8E   % .8E % .8E % .8E   %.8E   % .8E % .8E % .8E   %.8E \t",
+                    j, i,
+                    m[i],
+                    x[i][0], x[i][1], x[i][2], (x[i]-x[i_bh]).norm(),
+                    v[i][0], v[i][1], v[i][2], (v[i]-v[i_bh]).norm());
+        }
+        fprintf(out,"\n");
+    }
+    fprintf(out, "\n"); // this is redundant
+    fflush(out);
+}
+
+void Binary_smbh_influence_sphere_output::operator()(const std::vector<int>& ind_act, int n_act, double timesteps, double time_cur)
+{
+    double m_bh1  = m[0];
+    double m_bh2  = m[1];
+    double3 x_bh1 = x[0];
+    double3 x_bh2 = x[1];
+    double3 v_bh1 = v[0];
+    double3 v_bh2 = v[1];
+
+    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 = (x_bh1 - x_bh2).norm2();
+    double DV2 = (v_bh1 - v_bh2).norm2();
+    double EB = -(m_bh1 + m_bh2) / sqrt(DR2) + 0.5 * DV2;
+    double SEMI_a = -0.5 * (m_bh1 + m_bh2)/EB;
+    double SEMI_a2 = SEMI_a*SEMI_a;
+
+    for (int i=0; i<n_act; i++) {
+        int j_act = ind_act[i];
+        if (j_act<2) continue;
+        const double& pot_bh1 = pot[0];
+        const double& pot_bh2 = pot[1];
+        const double& m_act   = m[j_act];
+        const double3& x_act  = x[j_act];
+        const double3& v_act  = v[j_act];
+        const double& dt_act  = dt[j_act];
+        const double& pot_act = pot[j_act];
+        double tmp_r2 = (x_act - x_bbhc).norm2();
+        if (tmp_r2 < SEMI_a2*factor*factor) {
+            if (inf_event[j_act] == 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, j_act,
+                    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_bh1,
+                    m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_bh2,
+                    sqrt(tmp_r2),
+                    m_act, x_act[0], x_act[1], x_act[2], v_act[0], v_act[1], v_act[2], pot_act,
+                    dt_act);
+
+            inf_event[j_act] = 1;
+            }
+        } else {
+            if (inf_event[j_act] == 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_bh1,
+                    m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_bh2,
+                    sqrt(tmp_r2),
+                    m_act, x_act[0], x_act[1], x_act[2], v_act[0], v_act[1], v_act[2], pot_act,
+                    dt_act);
+            }
+            inf_event[j_act] = 0;
+        } /* if (tmp_r2 < DR2*R_INF2) */
+    } /* i */
+    fflush(out);
+}
--- a/black_holes.h
+++ b/black_holes.h
@ -0,0 +1,104 @@
+#include <algorithm>
+#include <cstdio>
+#include <vector>
+#include "double3.h"
+
+struct Particle_ref {
+    Particle_ref(double& m, double3& x, double3& v, double& pot, double3& acc, double3& jrk) :
+        m(m), x(x), v(v), pot(pot), acc(acc), jrk(jrk) {}
+    const double&  m;
+    const double3& x;
+    const double3& v;
+    double&        pot;
+    double3&       acc;
+    double3&       jrk;
+};
+
+struct Bbh_gravity {
+    double pot1, pot2;
+    double3 a1, a2, adot1, adot2, a_pn1[7], a_pn2[7], adot_pn1[7], adot_pn2[7];
+    double spin1[3], spin2[3];
+};
+
+class Black_hole_physics {
+public:
+    Black_hole_physics()
+        : count(0), c(0) {}
+    Black_hole_physics(const int count, const int myRank, const int rootRank)
+        : count(count), c(0), myRank(myRank), rootRank(rootRank) {}
+    void set_post_newtonian(const double c, const int pn_usage[7])
+    {
+        this->c = c;
+        std::copy(pn_usage, pn_usage+7, this->pn_usage);
+    }
+    void set_spins(const double spin1[3], const double spin2[3])
+    {
+        std::copy(spin1, spin1+3, this->bbh_grav.spin1);
+        std::copy(spin2, spin2+3, this->bbh_grav.spin2);
+    }
+    void set_softening(const double eps_old, const double eps_new)
+    {
+        this->eps_old = eps_old;
+        this->eps_new = eps_new;
+    }
+    void adjust_softening(const std::vector<Particle_ref>& particles);
+
+    void adjust_post_newtonian(
+        const double dt_bh, // pn_usage should be const
+        double3& acc1, double3& acc2,
+        double3& jrk1, double3& jrk2);
+    void write_bh_data(const double time_cur, const int count, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, const std::vector<double3> &a, const std::vector<double3> &adot, const std::vector<double> &dt);
+public: //TODO make private
+    /////////////std::vector<double> masses;
+    int count;
+    int myRank, rootRank;
+    double eps_old, eps_new;
+    double c;
+    int pn_usage[7];
+    Bbh_gravity bbh_grav;
+};
+
+class Write_bh_nb_data {
+public:
+    Write_bh_nb_data(int nb, int smbh_count, int N, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v)
+        : nb(nb), smbh_count(smbh_count), N(N), m(m), x(x), v(v)
+    {
+        ind_sort.resize(N);
+        var_sort.resize(N);
+        out = fopen("bh_neighbors.dat", "w");
+    }
+    ~Write_bh_nb_data()
+    {
+        fclose(out);
+    }
+    void operator()(double time_cur);
+private:
+    int nb, smbh_count, N;
+    const std::vector<double> &m;
+    const std::vector<double3> &x, &v;
+    std::vector<int> ind_sort;
+    std::vector<double> var_sort;
+    FILE *out;
+};
+
+class Binary_smbh_influence_sphere_output {
+public:
+    Binary_smbh_influence_sphere_output(double factor, int N, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, const std::vector<double> &dt)
+        : factor(factor), m(m), x(x), v(v), pot(pot), dt(dt)
+    {
+        inf_event.assign(N, 0);
+        out = fopen("bbh_inf.dat", "w");
+    }
+
+    ~Binary_smbh_influence_sphere_output()
+    {
+        fclose(out);
+    }
+    void operator()(const std::vector<int>& ind_act, int n_act, double timesteps, double time_cur);
+private:
+    double factor;
+    const std::vector<double> &pot, &m, &dt;
+    const std::vector<double3> &x, &v;
+    std::vector<int> inf_event;
+    FILE *out;
+};
--- a/cmd.c
+++ b/cmd.c
@ -1,76 +0,0 @@
-
-int cmd(int kstar, double rstar, double lstar, double Rgal, double *abvmag, double *vmag, double *BV, double *Teff, double *dvmag, double *dBV)
-  {
-
-  double lTeff, BC, kb;
-  double bvc[8], bcc[8];
-  double dbmag;
-  double BCsun, abvmagsun;
-  double rand1, rand2, prand;
-
-  kb = 5.6704E-08*0.5*1.3914E9*0.5*1.3914E9/3.846E26;  //Stefan-Boltzmann constant in Lsun Rsun^-2 K^-4
-
-  bvc[0] = -654597.405559323;
-  bvc[1] = 1099118.61158915;
-  bvc[2] = -789665.995692672;
-  bvc[3] = 314714.220932623;
-  bvc[4] = -75148.4728506455;
-  bvc[5] = 10751.803394526;
-  bvc[6] = -853.487897283685;
-  bvc[7] = 28.9988730655392;
-
-  bcc[0] = -4222907.80590972;
-  bcc[1] = 7209333.13326442;
-  bcc[2] = -5267167.04593882;
-  bcc[3] = 2134724.55938336;
-  bcc[4] = -518317.954642773;
-  bcc[5] = 75392.2372207101;
-  bcc[6] = -6082.7301194776;
-  bcc[7] = 209.990478646363;
-
-  BCsun = 0.11;     //sun's bolometric correction
-  abvmagsun = 4.83; //sun's absolute V magnitude
-
-  if( rstar && (kstar<14) )
-    {
-    *Teff = pow(lstar/(4.0*Pi*rstar*rstar*kb),0.25);
-
-    if( (*Teff>3000.0) && (*Teff<55000.0) )
-      {
-      lTeff = log10(*Teff);
-      *BV = bvc[0] + bvc[1]*lTeff + bvc[2]*pow(lTeff,2) + bvc[3]*pow(lTeff,3) + bvc[4]*pow(lTeff,4) + bvc[5]*pow(lTeff,5) + bvc[6]*pow(lTeff,6) + bvc[7]*pow(lTeff,7);
-      BC = bcc[0] + bcc[1]*lTeff + bcc[2]*pow(lTeff,2) + bcc[3]*pow(lTeff,3) + bcc[4]*pow(lTeff,4) + bcc[5]*pow(lTeff,5) + bcc[6]*pow(lTeff,6) + bcc[7]*pow(lTeff,7);
-      if(lstar) *abvmag = -2.5*log10(lstar)-BC+BCsun+abvmagsun;
-      *vmag = *abvmag + 5.0*log10(Rgal) - 5.0;
-
-      do{
-      rand1 = 2.0*drand48()-1.0;
-      rand2 = 2.0*drand48()-1.0;
-      } while (rand1*rand1+rand2*rand2 > 1.0);
-
-      prand = sqrt(-2.0*log(rand1*rand1+rand2*rand2)/(rand1*rand1+rand2*rand2));
-      *dvmag = rand1*prand*sqrt(pow(0.02,2) + pow(0.07*pow(10.0, 0.4*(*vmag-25.0)),2));
-      dbmag = rand2*prand*sqrt(pow(0.02,2) + pow(0.07*pow(10.0, 0.4*(*vmag-25.0)),2));
-      *dBV = *dvmag + dbmag;
-      }
-    else
-      {
-      *vmag = 9999.9;
-      *abvmag = 9999.9;
-      *BV = 9999.9;
-      *dvmag = 0.0;
-      *dBV = 0.0;
-      }
-    }
-  else
-    {
-    *Teff = 0.0;
-    *vmag = 9999.9;
-    *abvmag = 9999.9;
-    *BV = 9999.9;
-    *dvmag = 0.0;
-    *dBV = 0.0;
-    }
-
-  return(0);
-}
--- a/config.cpp
+++ b/config.cpp
@ -0,0 +1,245 @@
+#include "config.h"
+#include <stdexcept>
+#include <cstdio>
+#include <fstream>
+#include <algorithm>
+#include <limits>
+
+// Would be a bit more elegant to do the whole thing with std::variant.
+
+using Dictionary = std::unordered_map<std::string,std::string>;
+
+static constexpr double nix = -std::numeric_limits<double>::max(); // avoid nans
+
+std::string Config::strip(const std::string str)
+{
+    std::string str_new = str;
+    auto pos = str_new.find_first_not_of(" \t");
+    if (pos != std::string::npos) str_new = str_new.substr(pos, str_new.size());
+    pos = str_new.find_last_not_of(" \t");
+    if (pos != std::string::npos) str_new = str_new.substr(0, pos+1);
+    return str_new;
+}
+
+Dictionary Config::read_config_file(const std::string file_name)
+{
+    std::unordered_map<std::string,std::string> dictionary;
+    std::ifstream file(file_name);
+    if (!file.good()) throw std::runtime_error("File not found.");
+    std::string str;
+    int line_number = 0;
+    while (std::getline(file, str)) {
+        line_number++;
+        auto pos = str.find('#');
+        if (pos != std::string::npos) str = str.substr(0, pos);
+        str = strip(str);
+        if (str.size() == 0) continue;
+        pos = str.find_first_of("=");
+        if (pos == std::string::npos) throw std::runtime_error("Error: expected a key-value pair in line " + std::to_string(line_number) + " of file " + file_name);
+        std::string key = strip(str.substr(0, pos));
+        pos = str.find_first_not_of(" \t", pos+1);
+        std::string val = strip(str.substr(pos, str.size()));
+        dictionary[key] = val;
+    }
+    return dictionary;
+}
+
+template<> std::string Config::string_cast(const std::string str)
+{
+    return str;
+}
+
+template<> double Config::string_cast(const std::string str)
+{
+    size_t idx;
+    auto value = std::stod(str, &idx);
+    if (idx == str.size()) return value;
+    else throw std::runtime_error("Cannot convert \"" + str + "\" into a double");
+}
+
+template<> int Config::string_cast(const std::string str)
+{
+    size_t 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");
+}
+
+template<> bool Config::string_cast(const std::string str)
+{
+    if      ((str=="true")  || (str=="True")  || (str=="yes") || (str=="Yes") || (str=="1")) return true;
+    else if ((str=="false") || (str=="False") || (str=="no")  || (str=="No")  || (str=="0")) return false;
+    throw std::runtime_error("Cannot convert \"" + str + "\" into a bool");
+}
+
+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;
+}
+
+template<> std::vector<double> 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<double> result;
+    while (new_str.length() > 0) {
+        size_t idx;
+        auto value = std::stod(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>
+T Config::get_parameter(Dictionary dictionary, std::string name)
+{
+    auto item = dictionary.find(name);
+    if (item==dictionary.end()) throw std::runtime_error("Mandatory parameter " + name + " must be defined");
+    else return string_cast<T>((*item).second);
+}
+
+// For optional parameters
+template<typename T>
+T Config::get_parameter(Dictionary dictionary, std::string name, T default_value)
+{
+    auto item = dictionary.find(name);
+    if (item==dictionary.end()) return default_value;
+    else return string_cast<T>((*item).second);
+}
+
+#include <iostream>
+
+void Config::error_checking()
+{
+#ifndef HAS_HDF5
+    if (output_hdf5)
+        throw std::runtime_error("HDF5 output format (output_hdf5=true) requires the code to be compiled with HAS_HDF5");
+#endif
+    if (live_smbh_count < 0)
+        throw std::runtime_error("The number of live black holes (live_smbh_count) has to be greater than or equals to zero");
+    if (binary_smbh_pn && (live_smbh_count!=2))
+        throw std::runtime_error("Post-Newtonian gravity (binary_smbh_pn=true) requires live_smbh_count=2");
+    if (binary_smbh_pn && (pn_c <= 0))
+        throw std::runtime_error("Post-Newtonian gravity (binary_smbh_pn=true) requires pn_c > 0");
+    if (live_smbh_custom_eps == eps) live_smbh_custom_eps = -1;
+    if (live_smbh_output && (live_smbh_count == 0))
+        throw std::runtime_error("Black hole output (live_smbh_output=true) requires at least one live black hole (live_smbh_count)");
+    if (live_smbh_neighbor_output && (live_smbh_count == 0))
+        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");
+    if (binary_smbh_pn)
+        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 be a 7-component vector filled with ones and zeros only");
+    if ((smbh1_spin.size()!=3) || (smbh2_spin.size()!=3))
+        throw std::runtime_error("Spins must be three-component vectors");
+    if ((pn_usage[6]==1) && ((smbh1_spin[0]==nix) || (smbh2_spin[0]==nix)))
+        throw std::runtime_error("Please define smbh1_spin and smbh2_spin or disable the spin by setting the last component of pn_usage to zero");
+    std::cout << smbh1_spin[0] << std::endl;
+        std::cout << smbh1_spin[1] << std::endl;
+            std::cout << smbh1_spin[2] << std::endl;
+    if ((pn_usage[6]==0) && ((smbh1_spin[0]!=nix) || (smbh2_spin[0]!=nix)))
+        throw std::runtime_error("Spins (smbh1_spin and smbh2_spin) may not be defined if the last element of pn_usage is set to zero");
+
+
+    if (ext_units_physical && ((unit_mass == 0) || (unit_length == 0)))
+        throw std::runtime_error("Physical units for external gravity (ext_units_physical) requires ext_unit_mass and ext_unit_length to be positive numbers");
+    if ((ext_m_bulge > 0) && (ext_b_bulge < 0))
+        throw std::runtime_error("To use external bulge gravity, please specify positive ext_m_bulge and ext_b_bulge");
+    if ((ext_m_halo_plummer > 0) && (ext_b_halo_plummer < 0))
+        throw std::runtime_error("To use external Plummer halo gravity, please specify positive ext_m_halo_plummer and ext_b_halo_plummer");
+    if ((ext_m_disk > 0) && ((ext_a_disk < 0) || (ext_b_disk < 0)))
+        throw std::runtime_error("To use external disk gravity, please specify positive ext_m_disk, ext_a_disk and ext_b_disk");
+    if (((ext_log_halo_r > 0) && (ext_log_halo_v <= 0)) || ((ext_log_halo_r <= 0) && (ext_log_halo_v > 0)))
+        throw std::runtime_error("To use external logarithmic halo gravity, please specify positive ext_log_halo_r and ext_log_halo_v");
+    if ((ext_dehnen_m > 0) && ((ext_dehnen_r <= 0) || (ext_dehnen_gamma <= 0)))
+        throw std::runtime_error("To use external Dehnen model, please specify positive ext_dehnen_r and ext_dehnen_gamma");
+}
+
+Config::Config(std::string file_name)
+{
+    auto dictionary = read_config_file(file_name);
+
+    // TODO check if dt_disk and dt_contr are powers of two
+    eps = get_parameter<double>(dictionary, "eps");
+    t_end = get_parameter<double>(dictionary, "t_end");
+    dt_disk = get_parameter<double>(dictionary, "dt_disk");
+    dt_contr = get_parameter<double>(dictionary, "dt_contr");
+    dt_bh = get_parameter<double>(dictionary, "dt_bh", dt_contr);
+    eta = get_parameter<double>(dictionary, "eta");
+    input_file_name = get_parameter<std::string>(dictionary, "input_file_name", "data.con");
+    devices_per_node = get_parameter<int>(dictionary, "devices_per_node", 0);
+    dt_max_power = get_parameter<int>(dictionary, "dt_max_power", -3);
+    dt_min_power = get_parameter<int>(dictionary, "dt_min_power", -36);
+    eta_s_corr = get_parameter<double>(dictionary, "eta_s_corr", 4);
+    eta_bh_corr = get_parameter<double>(dictionary, "eta_bh_corr", 4);
+
+    output_hdf5 = get_parameter<bool>(dictionary, "output_hdf5", false);
+    output_hdf5_double_precision = get_parameter<bool>(dictionary, "output_hdf5_double_precision", true);
+    output_ascii_precision = get_parameter<int>(dictionary, "output_ascii_precision", 10);
+    output_extra_mode = get_parameter<int>(dictionary, "output_extra_mode", 10);
+    dt_min_warning = get_parameter<bool>(dictionary, "dt_min_warning", false);
+
+    live_smbh_count = get_parameter<int>(dictionary, "live_smbh_count", 0);
+    live_smbh_custom_eps = get_parameter<double>(dictionary, "live_smbh_custom_eps", -1);
+    live_smbh_output = get_parameter<bool>(dictionary, "live_smbh_output", false);
+    live_smbh_neighbor_output = get_parameter<bool>(dictionary, "live_smbh_neighbor_output", false);
+    live_smbh_neighbor_number = get_parameter<int>(dictionary, "live_smbh_neighbor_number", 10);
+
+    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.);
+
+    binary_smbh_pn = get_parameter<bool>(dictionary, "binary_smbh_pn", false);
+    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", 0);
+    smbh1_spin = get_parameter<std::vector<double>>(dictionary, "smbh1_spin", std::vector<double>({nix,nix,nix}));
+    smbh2_spin = get_parameter<std::vector<double>>(dictionary, "smbh2_spin", std::vector<double>({nix,nix,nix}));
+
+    ext_units_physical = get_parameter<bool>(dictionary, "ext_units_physical", false);
+    unit_mass = get_parameter<double>(dictionary, "unit_mass", !ext_units_physical);
+    unit_length = get_parameter<double>(dictionary, "unit_length", !ext_units_physical);
+    ext_m_bulge = get_parameter<double>(dictionary, "ext_m_bulge", 0);
+    ext_b_bulge = get_parameter<double>(dictionary, "ext_b_bulge", -1);
+    ext_m_disk = get_parameter<double>(dictionary, "ext_m_disk", 0);
+    ext_a_disk = get_parameter<double>(dictionary, "ext_a_disk", -1);
+    ext_b_disk = get_parameter<double>(dictionary, "ext_b_disk", -1);
+    ext_m_halo_plummer = get_parameter<double>(dictionary, "ext_m_halo_plummer", 0);
+    ext_b_halo_plummer = get_parameter<double>(dictionary, "ext_b_halo_plummer", -1);
+    ext_log_halo_v = get_parameter<double>(dictionary, "ext_log_halo_v", 0);
+    ext_log_halo_r = get_parameter<double>(dictionary, "ext_log_halo_r", 0);
+    ext_dehnen_m = get_parameter<double>(dictionary, "ext_dehnen_m", 0);
+    ext_dehnen_r = get_parameter<double>(dictionary, "ext_dehnen_r", -1);
+    ext_dehnen_gamma = get_parameter<double>(dictionary, "ext_dehnen_gamma", -1);
+
+#ifdef ETICS
+    dt_scf = get_parameter<double>(dictionary, "dt_scf");
+    grapite_mask_file_name = get_parameter<std::string>(dictionary, "grapite_mask_file_name", "grapite.mask");
+    etics_dump_coeffs = get_parameter<bool>(dictionary, "etics_dump_coeffs", false);
+    grapite_active_search = get_parameter<bool>(dictionary, "grapite_active_search", false);
+    grapite_smbh_star_eps = get_parameter<double>(dictionary, "grapite_smbh_star_eps", -1);
+    grapite_dev_exec_threshold = get_parameter<int>(dictionary, "grapite_dev_exec_threshold", 32);
+#endif
+
+    error_checking();
+}
--- a/config.h
+++ b/config.h
@ -0,0 +1,76 @@
+#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;
+    int                 devices_per_node;
+    int                 dt_max_power;
+    int                 dt_min_power;
+    double              eta_s_corr;
+    double              eta_bh_corr;
+   
+    bool                output_hdf5;
+    bool                output_hdf5_double_precision;
+    int                 output_ascii_precision;
+    int                 output_extra_mode;
+    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;
+    std::vector<double> smbh1_spin;
+    std::vector<double> smbh2_spin;
+
+    bool                ext_units_physical;
+    double              unit_mass;
+    double              unit_length;
+    double              ext_m_bulge;
+    double              ext_b_bulge;
+    double              ext_m_disk;
+    double              ext_a_disk;
+    double              ext_b_disk;
+    double              ext_m_halo_plummer;
+    double              ext_b_halo_plummer;
+    double              ext_log_halo_r;
+    double              ext_log_halo_v;
+    double              ext_dehnen_m;
+    double              ext_dehnen_r;
+    double              ext_dehnen_gamma;
+
+#ifdef ETICS
+    double              dt_scf;
+    std::string         grapite_mask_file_name;
+    bool                etics_dump_coeffs;
+    bool                grapite_active_search;
+    double              grapite_smbh_star_eps;
+    int                 grapite_dev_exec_threshold;
+#endif
+
+private:
+    using Dictionary = std::unordered_map<std::string,std::string>;
+    Dictionary read_config_file(const std::string file_name);
+    std::string strip(const std::string str);
+    template<typename T> T string_cast(const std::string str);
+    template<typename T> T get_parameter(Dictionary dictionary, std::string name);
+    template<typename T> T get_parameter(Dictionary dictionary, std::string name, T default_value);
+    void error_checking();
+};
--- a/debug.h
+++ b/debug.h
@ -1,143 +0,0 @@
-/***************************************************************/
-void d_swap(double *a, double *b)
-{
-register double tmp;
-tmp = *a; *a = *b; *b = tmp;
-}
-
-void i_swap(int *a, int *b)
-{
-register int tmp;
-tmp = *a; *a = *b; *b = tmp;
-}
-/***************************************************************/
-
-/***************************************************************/
-void my_sort(int l, int r, double *arr, int *ind)
-{
-
-int    i, j, cikl;
-double tmp;
-
-i = l; j = r;
-tmp = arr[(l+r)/2];
-
-cikl = 1;
-
-while (cikl)
-  {
-    while (arr[i]<tmp) i++;
-    while (tmp<arr[j]) j--;
-
-    if (i<=j)
-     {
-	d_swap(&arr[i],&arr[j]);
-	i_swap(&ind[i],&ind[j]);
-	i++; j--;
-     }
-    else
-     {
-	cikl = 0;
-     }
-
-  }
-
-  if (l<j) my_sort(l, j, arr, ind);
-  if (i<r) my_sort(i, r, arr, ind);
-}
-/***************************************************************/
-
-/***************************************************************/
-double my_select(int n_1, int n_2, int k, double *arr, int *ind)
-{
-
-int   i, ir, j, l, mid, a_ind;
-double a;
-
-l  = n_1;
-ir = n_2;
-
-for(;;)
-  {
-
-  if (ir <= l+1) 
-    {
-
-    if (ir == l+1 && arr[ir] < arr[l]) 
-      {
-      d_swap(&arr[l],&arr[ir]);
-      i_swap(&ind[l],&ind[ir]);
-      }
-
-    return(arr[k]);
-
-    } 
-    else 
-    {
-
-    mid=(l+ir) >> 1;
-
-    d_swap(&arr[mid],&arr[l+1]);
-    i_swap(&ind[mid],&ind[l+1]);
-
-    if (arr[l+1] > arr[ir]) 
-      {
-      d_swap(&arr[l+1],&arr[ir]);
-      i_swap(&ind[l+1],&ind[ir]);
-      }
-
-    if (arr[l] > arr[ir]) 
-      {
-      d_swap(&arr[l],&arr[ir]);
-      i_swap(&ind[l],&ind[ir]);
-      }
-
-    if (arr[l+1] > arr[l])
-      {
-      d_swap(&arr[l+1],&arr[l]);
-      i_swap(&ind[l+1],&ind[l]);
-      }
-
-    i = l+1;
-    j = ir;
-    a = arr[l];
-    a_ind = ind[l];
-
-    for (;;) 
-     {
-     do i++; while (arr[i] < a);
-     do j--; while (arr[j] > a);
-     if (j < i) break;
-     d_swap(&arr[i],&arr[j]);
-     i_swap(&ind[i],&ind[j]);
-     }
-
-    arr[l] = arr[j];
-    ind[l] = ind[j];
-    arr[j] = a;
-    ind[j] = a_ind;
-    if (j >= k) ir = j-1;
-    if (j <= k) l = i;
-
-    }
-  }
-}
-/***************************************************************/
-
-/***************************************************************/
-/*
-double lagrange_radius(double percent)
-{ 
-int    Nb;
-double tmp;
-
-Nb = (int)(percent * N);
-
-my_sort(0, N-1, d, ind);
-
-tmp = my_select(0, N-1, Nb-1, d, ind); d[Nb-1] = tmp;
-
-return(tmp);
-}
-*/
-/***************************************************************/
--- a/def_DEN.c
+++ b/def_DEN.c
@ -1,38 +0,0 @@
-
-/*************************************************************************/
-/*************************************************************************/
-/*************************************************************************/
-
-void DEF_DEN(double x_gas[][3], double m_gas[], double h_gas[], int sosed_gas[][NB], double DEN_gas[])
-{
-
-double Xij, Yij, Zij, Rij, tmp;
-
-/* SPH style */
-
-for(i=0;i<N_GAS;i++)
-  {
-  tmp = 0.0;
-
-  for(k=0;k<NB;k++) 
-    {
-    j = sosed_gas[i][k]; 
-
-    Xij = x_gas[i][0]-x_gas[j][0]; 
-    Yij = x_gas[i][1]-x_gas[j][1]; 
-    Zij = x_gas[i][2]-x_gas[j][2];
-
-    Rij = sqrt( Xij*Xij + Yij*Yij + Zij*Zij );
-
-    tmp = tmp + 0.5 * m_gas[j] * ( W(Rij,h_gas[i]) + W(Rij,h_gas[j]) );
-    } /* j */
-
-  DEN_gas[i] = tmp;
-
-  } /* i */
-
-}
-
-/*************************************************************************/
-/*************************************************************************/
-/*************************************************************************/
--- a/def_DEN_STAR.c
+++ b/def_DEN_STAR.c
@ -1,74 +0,0 @@
-/*************************************************************************/
-/*************************************************************************/
-/*************************************************************************/
-
-void DEF_DEN_STAR(double x_star, double y_star, double z_star, double x_gas[][3], double m_gas[], double h_gas[], double *DEN_star)
-{
-
-//int    NB_STAR=0;
-//int    sosed_tmp[N_GAS_MAX];
-
-double Xij, Yij, Zij, Rij, tmp=0.0, tmp2;
-
-
-/* DEF DEN dlja odnoj STAR chastici */
-
-//NB_STAR = 0;
-tmp = 0.0;
-tmp2 = 0.0;
-
-for(j=0;j<N_GAS;j++)	// GAS
-  {
-
-  tmp2 = SQR(x_star) + SQR(y_star);
-
-  if( (ABS(z_star) < 1e-2) && (tmp2 < 0.2) )		// zamenity na h & (2*R0)^2 = 0.1936
-    {
-
-    Xij = x_star - x_gas[j][0];
-    Yij = y_star - x_gas[j][1];
-    Zij = z_star - x_gas[j][2];
-
-    Rij = sqrt(Xij*Xij + Yij*Yij + Zij*Zij);
-
-  if( Rij < (2.0*h_gas[j]) )
-    {
-//    sosed_tmp[NB_STAR] = j;
-//    NB_STAR++;
-
-    tmp = tmp + m_gas[j] * W(Rij,h_gas[j]);
-    }
-    }
-
-//  sosed_tmp[j] = j;
-  } /* j */
-
-//printf("%06d \t %.8E \n",NB_STAR, tmp);
-
-
-//  my_sort(0, N-1, d, sosed_tmp); 
-
-//  tmp = my_select(0, N_GAS-1, NB-1, d, sosed_tmp); d[NB-1] = tmp;
-
-/* SPH style */
-
-/*
-  tmp = 0.0;
-
-  for(k=0;k<NB_STAR;k++) 
-    {
-    j = sosed_tmp[k];
-
-    Rij = sqrt(d[k]);
-    
-    tmp = tmp + m[j] * W(Rij,h_gas[j]);
-    } 
-*/
-
-  *DEN_star = tmp;
-
-}
-
-/*************************************************************************/
-/*************************************************************************/
-/*************************************************************************/
--- a/def_H.c
+++ b/def_H.c
@ -1,85 +0,0 @@
-/*************************************************************************/
-/*************************************************************************/
-/*************************************************************************/
-
-double W(double Rij, double h)
-{
-
-double v, v2, v3, tmp=0.0;
-
-  v = Rij/h;  v2 = v*v;  v3 = v*v*v;
-
-  if( (v >= 0.0) && (v <= 1.0) ) 
-    {
-    tmp = 1.0 - 1.5*v2 + 0.75*v3; 
-    }
-  else 
-    {  
-    if( (v > 1.0) && (v < 2.0) ) 
-      {
-      tmp = 0.25*(2.0-v)*(2.0-v)*(2.0-v); 
-      }
-    else 
-      { 
-      tmp = 0.0; /* if ( v >= 2.0 ) */ 
-      } 
-    }
-
-  tmp = 1.0/(Pi*h*h*h)*tmp;
-
-  return(tmp);
-
-}
-
-/*************************************************************************/
-
-void DEF_H(double x_gas[][3], double h_gas[], int sosed_gas[][NB])
-{
-
-int    sosed_tmp[N_GAS_MAX];
-double Xij, Yij, Zij;
-
-for(i=0;i<N_GAS;i++)
-  {
-
-for(j=0;j<N_GAS;j++)
-  {
-  Xij = x_gas[i][0]-x_gas[j][0]; 
-  Yij = x_gas[i][1]-x_gas[j][1]; 
-  Zij = x_gas[i][2]-x_gas[j][2];
-
-  d[j] = Xij*Xij + Yij*Yij + Zij*Zij;
-
-  sosed_tmp[j] = ind_gas[j];
-  } /* j */
-
-  my_sort(0, N_GAS-1, d, sosed_tmp); 
-
-//  tmp = my_select(0, N_GAS-1, NB-1, d, sosed_tmp); d[NB-1] = tmp;
-
-  h_gas[i] = sqrt( d[NB-1] )*0.5;
-
-  for(k=0;k<NB;k++) sosed_gas[i][k] = sosed_tmp[k];
-  
-/*
-  h_gas[i] = MAX(h_gas[i], h_min); 
-  h_gas[i] = MIN(h_gas[i], h_max);
-*/
-
-/*
-  tmp_l = ldiv(i, N/64);
-
-  if(tmp_l.rem == 0) 
-    {
-    printf("."); 
-    fflush(stdout);
-    }
-*/
-
-  } /* i */
-
-}
-
-/*************************************************************************/
-/*************************************************************************/
-/*************************************************************************/
--- a/double3.h
+++ b/double3.h
@ -0,0 +1,81 @@
+#pragma once
+#include <cmath>
+
+struct double3 {
+    union {
+        struct {double x, y, z;};
+        double data[3]; // For legacy access
+    };
+    double3() {}
+    double3(const double x, const double y, const double z)
+        : x(x), y(y), z(z) {}
+    double& operator[](int i) {return data[i];}
+    const double operator[](int i) const {return data[i];}
+    double3& operator=(const double3& a)
+    {
+        x = a.x;
+        y = a.y;
+        z = a.z;
+        return *this;
+    }
+    double3& operator+=(const double3& a)
+    {
+        x += a.x;
+        y += a.y;
+        z += a.z;
+        return *this;
+    }
+    double3& operator-=(const double3& a)
+    {
+        x -= a.x;
+        y -= a.y;
+        z -= a.z;
+        return *this;
+    }
+    double3& operator/=(const double& c)
+    {
+        x /= c;
+        y /= c;
+        z /= c;
+        return *this;
+    }
+    double norm2() const
+    {
+        return x*x + y*y + z*z;
+    }
+    double norm() const
+    {
+        return sqrt(x*x + y*y + z*z);
+    }
+    operator double*() {return data;}
+};
+
+inline double3 operator*(const double& c, const double3& a)
+{
+    return double3(a.x*c, a.y*c, a.z*c);
+}
+
+inline double3 operator*(const double3& a, const double& c)
+{
+    return double3(a.x*c, a.y*c, a.z*c);
+}
+
+inline double operator*(const double3& a, const double3& b)
+{
+    return a.x*b.x+a.y*b.y+a.z*b.z;
+}
+
+inline double3 operator/(const double3& a, const double& c)
+{
+    return double3(a.x/c, a.y/c, a.z/c);
+}
+
+inline double3 operator+(const double3& a, const double3& b)
+{
+    return double3(a.x+b.x, a.y+b.y, a.z+b.z);
+}
+
+inline double3 operator-(const double3& a, const double3& b)
+{
+    return double3(a.x-b.x, a.y-b.y, a.z-b.z);
+}
--- a/drag_force.c
+++ b/drag_force.c
@ -1,185 +0,0 @@
-/*****************************************************************************
-  File Name      : "Drag Force.c"
-  Contents       : Drag force calculus.
-                 : Converted & optimized from Chingis Omarov Fortran code.
-  Coded by       : Taras Panamarev, 
-                 : Denis Yurin, Maxim Makukov & Peter Berczik :-)
-  Last redaction : 2012-11-15 15:33
-*****************************************************************************/
-
-// main parameters
-
-const double Md = 0.01;
-const double s=4, beta_s=0.7, alpha=0.75, h=0.001;
-const double r_lim2=(0.5*0.5), h_lim2=(5e-3*5e-3); 
-const double gradP_rho = 0.90;	// account for -grad(P)/rho in gas disk... 
-
-const double Qtot = 0.01;	//  16e3
-
-/*
-const double Qtot = 0.00547000;	//   8e3
-const double Qtot = 0.00296975;	//  16e3
-const double Qtot = 0.00193372;	//  32e3
-const double Qtot = 0.00054215;	// 128KB
-*/
-
-int    K; 
-
-double RDISC, //projection of radius vector to the disc plane
-       h_z, // disc profile
-       RDISC2, Z2, RDISC_3over2, R, RR0, RR02, RR0_s, RDISC_5over2, R_zeta, R0_zeta, RRcrit, RRcrit_zeta,Rcrit_zeta, R_zetaR, 
-       inner_hole, temp, sqrt_m_bh, pot_drag, 
-       sigma0, R02, hR02, Rcrit, R0, zeta, 
-       R_DOT, //the first derivative of absolute value of projection
-       COEFDENS, 
-       RHO, RHODOT, //disc density and it's firs derivative with respect to time 
-       VXDISC, VYDISC,//disc particle velocity components 
-       VXRE, VYRE, VZRE, VRE,//the absolute value of relative (star-disc) velocity and it's components  
-       coef, 
-       FDRAGX, FDRAGY, FDRAGZ,//drag force components 
-       DVXRE, DVYRE, DVZRE, DVRE,//first derivative of relative velocity and it's absolute value 
-       KFDOT;
-
-void init_drag_force() 
-  {
-  R0 = R_0;
-  Rcrit = R_CRIT;
-
-  sigma0 = (2-alpha)*Md/(TWOPi*sqrt_TWOPi*h*R0);
-  R02 = SQR(R0);
-  hR02 = SQR(h*R0);
-  }
-
-void calc_drag_force() 
-  {
-
-  if (min_t < t_diss_on) return;
-//  if (min_t == t_diss_on) init_drag_force();
-  init_drag_force();
-
-  for (i=0; i<n_act; i++) 
-    {
-  
-    K = ind_act[i];
-
-    RDISC2 = SQR(x_act_new[i][0]) + SQR(x_act_new[i][1]);
-    Z2 = SQR(x_act_new[i][2]);
-
-    if( (RDISC2 < r_lim2) && (Z2 < h_lim2) )
-      {
-
-      RDISC2 = SQR(x_act_new[i][0]) + SQR(x_act_new[i][1]);
-      R = sqrt(RDISC2 + SQR(x_act_new[i][2]));
-        
-      RDISC = sqrt(RDISC2);
-      R_DOT = ( x_act_new[i][0]*v_act_new[i][0] + x_act_new[i][1]*v_act_new[i][1] )/RDISC;
-
-      RR0 = RDISC/R0;
-      RR02 = SQR(RR0);
-//    RR0_s = pow(RR0, s);
-      RR0_s = SQR(RR02);
-      RRcrit = RDISC/Rcrit;
-
-/*
-      if (RDISC < Rcrit) zeta = 1; else zeta = 0;
-
-      RRcrit_zeta = pow(RRcrit, zeta);
-      R_zeta = pow(RDISC, -2*zeta);
-      R0_zeta = pow(R0, -2*zeta);
-      Rcrit_zeta = pow(Rcrit, -2*zeta);
-*/
-
-      if (RDISC < Rcrit)
-        {
-        zeta = 1;
-        RRcrit_zeta = RRcrit;
-        R_zeta = 1.0/RDISC2;
-        R0_zeta = 1.0/R02;
-        Rcrit_zeta = 1.0/SQR(Rcrit);
-        }
-      else 
-        {
-        zeta = 0; 
-        RRcrit_zeta = 1.0;
-        R_zeta = 1.0;
-        R0_zeta = 1.0;
-        Rcrit_zeta = 1.0;
-        }
-
-      R_zetaR = R_zeta/R;
-
-      COEFDENS = pow(RR0, -alpha) * exp(-beta_s*RR0_s)/RRcrit_zeta;
-      RHO = sigma0*COEFDENS*exp(-Z2/(2.0*hR02*SQR(RRcrit_zeta)));
-
-      RHODOT = - R_DOT*(zeta+alpha+beta_s*s*RR0_s )/RDISC 
-               - (v_act_new[i][2]*x_act_new[i][2]*R_zeta-zeta*R_zetaR*Z2*R_DOT)/(hR02*Rcrit_zeta);
-
-      RDISC_3over2 = sqrt(RDISC2*RDISC);
-      RDISC_5over2 = RDISC_3over2*RDISC;
-
-      sqrt_m_bh = sqrt(m_bh);
-
-      sqrt_m_bh *= gradP_rho;			// account for -grad(P)/rho in gas disk... 
-
-      VXDISC = -sqrt_m_bh * x_act_new[i][1]/RDISC_3over2;
-      VYDISC = +sqrt_m_bh * x_act_new[i][0]/RDISC_3over2;
-
-      VXRE = VXDISC - v_act_new[i][0];
-      VYRE = VYDISC - v_act_new[i][1];
-      VZRE =        - v_act_new[i][2];
-
-      VRE = sqrt( SQR(VXRE) + SQR(VYRE) + SQR(VZRE) );
-
-      coef = Qtot*RHO*VRE;
-
-      FDRAGX = coef*VXRE;
-      FDRAGY = coef*VYRE;
-      FDRAGZ = coef*VZRE;
-
-      DVXRE = -sqrt_m_bh * v_act_new[i][1]/RDISC_3over2
-      	      +sqrt_m_bh * 1.5*x_act_new[i][1]*R_DOT/RDISC_5over2
-      	      -a_act_new[i][0];
-
-      DVYRE = +sqrt_m_bh * v_act_new[i][0]/RDISC_3over2
-      	      -sqrt_m_bh * 1.5*x_act_new[i][0]*R_DOT/RDISC_5over2
-      	      -a_act_new[i][1];
-
-      DVZRE = -a_act_new[i][2];
-
-      DVRE  = (VXRE*DVXRE + VYRE*DVYRE + VZRE*DVZRE)/VRE;
-
-      KFDOT = Qtot*RHO;
-
-      adot_drag[K][0] = KFDOT*(RHODOT*VXRE*VRE + DVRE*VXRE + DVXRE*VRE);
-      adot_drag[K][1] = KFDOT*(RHODOT*VYRE*VRE + DVRE*VYRE + DVYRE*VRE);
-      adot_drag[K][2] = KFDOT*(RHODOT*VZRE*VRE + DVRE*VZRE + DVZRE*VRE);
-
-      pot_drag = -0.5*( (a_drag[K][0]+FDRAGX)*(x_act_new[i][0]-x_act[i][0]) +
-                (a_drag[K][1]+FDRAGY)*(x_act_new[i][1]-x_act[i][1]) +
-                (a_drag[K][2]+FDRAGZ)*(x_act_new[i][2]-x_act[i][2]) );
-
-      E_sd += m[K]*pot_drag;
-
-      a_drag[K][0] = FDRAGX;
-      a_drag[K][1] = FDRAGY;
-      a_drag[K][2] = FDRAGZ;
-
-      }
-    else
-      {
-      adot_drag[K][0] = 0.0;
-      adot_drag[K][1] = 0.0;
-      adot_drag[K][2] = 0.0;
-
-      a_drag[K][0] = 0.0;
-      a_drag[K][1] = 0.0;
-      a_drag[K][2] = 0.0;
-      } /* if( (RDISC2 < r_lim2) && (Z2 < h_lim2) ) */
-
-  } /* i */
-
-
-}
-/***************************************************************/
-/***************************************************************/
-/***************************************************************/
--- a/external.cpp
+++ b/external.cpp
@ -0,0 +1,85 @@
+#include "external.h"
+
+inline double square(double x) {return x*x;}
+
+void Plummer::calc_gravity()
+{
+    double r2 = square(b);
+    r2 += x.norm2();
+    double r = sqrt(r2);
+    double rv_ij = v*x;
+    double tmp = m / r;
+    potential = -tmp;
+    tmp /= r2;
+    acceleration = - tmp * x;
+    jerk = - tmp * (v - 3*rv_ij * x/r2);
+}
+
+void Miyamoto_Nagai::calc_gravity()
+{
+    double  x_ij=x[0],  y_ij=x[1],  z_ij=x[2];
+    double vx_ij=v[0], vy_ij=v[1], vz_ij=v[2];
+    auto z2_tmp = z_ij*z_ij + square(b);
+    auto z_tmp  = sqrt(z2_tmp);
+    auto r2_tmp = x_ij*x_ij + y_ij*y_ij + square(z_tmp + a);
+    auto r_tmp  = sqrt(r2_tmp);
+    potential = - m / r_tmp;
+    auto tmp = m / (r2_tmp*r_tmp);
+    acceleration[0] = - tmp * x_ij;
+    acceleration[1] = - tmp * y_ij;
+    acceleration[2] = - tmp * z_ij * (z_tmp + a)/z_tmp;
+    tmp = m / (z_tmp*r2_tmp*r2_tmp*r_tmp);
+    jerk[0] = tmp * (- vx_ij*z_tmp*r2_tmp
+            + 3.0*x_ij*vx_ij*x_ij*z_tmp
+            + 3.0*x_ij*vy_ij*y_ij*z_tmp
+            + 3.0*x_ij*vz_ij*z_ij*square(z_tmp + a));
+    jerk[1] = tmp * (- vy_ij*z_tmp*r2_tmp
+            + 3.0*y_ij*vx_ij*x_ij*z_tmp
+            + 3.0*y_ij*vy_ij*y_ij*z_tmp
+            + 3.0*y_ij*vz_ij*z_ij*square(z_tmp + a));
+    jerk[2] = tmp * (- vz_ij*(z_tmp + a)*(x_ij*x_ij*(z2_tmp*z_tmp + a*square(b)) +
+            y_ij*y_ij*(z2_tmp*z_tmp + a*square(b)) -
+            (2.0*a*(square(z_ij) - square(b))*z_tmp +
+            2.0*square(z_ij*z_ij) +
+            square(b)*square(z_ij) -
+            square(b)*(square(a) + square(b))))
+            + 3.0*vx_ij*x_ij*z_ij*z2_tmp*(z_tmp + a)
+            + 3.0*vy_ij*y_ij*z_ij*z2_tmp*(z_tmp + a)) / z2_tmp;
+}
+
+void Logarithmic_halo::calc_gravity()
+{
+    auto r2 = x.norm2();
+    auto rv_ij = 2.0*(v*x);
+    auto r2_r2_halo = (r2 + r2_halo);
+    potential = - 0.5*v2_halo * log(1.0 + r2/r2_halo);
+    auto tmp = v2_halo/r2_r2_halo;
+    acceleration = - tmp * x;
+    tmp /= (r2 + r2_halo);
+    jerk = tmp * (rv_ij * x - r2_r2_halo * v);
+}
+
+void Dehnen::calc_gravity()
+{
+
+    auto tmp_r = x.norm();
+    auto tmp_r2 = tmp_r*tmp_r;
+    auto tmp_r3 = tmp_r2*tmp_r;
+
+    auto dum   = tmp_r/(tmp_r + r);
+    auto dum2  = dum*dum;
+    auto dum3  = dum2*dum;
+    auto dum_g = pow(dum, -gamma);
+
+    potential = - ( (m/r) / (2.0-gamma) ) * ( 1.0 - dum2 * dum_g );
+
+    auto  tmp = (m/tmp_r3) * dum3 * dum_g;
+
+    acceleration = - tmp * x;
+
+    auto rv_ij = v*x;
+    tmp = ( m/((tmp_r+r)*(tmp_r+r)*(tmp_r+r)) ) * dum_g;
+    auto tmp_g = ( (r*gamma + 3.0*tmp_r)/(tmp_r2*(tmp_r+r)) ) * rv_ij;
+
+    jerk = tmp * (tmp_g * x - v);
+}
--- a/external.h
+++ b/external.h
@ -0,0 +1,87 @@
+#pragma once
+#include <string>
+#include <vector>
+#include "double3.h"
+
+class External_gravity {
+public:
+    void apply(const int n_act, const std::vector<double3> &x, const std::vector<double3> &v, std::vector<double> &pot, std::vector<double3> &a, std::vector<double3> &adot)
+    {
+        for (int i=0; i<n_act; i++) {
+            this->set_coordinates(x[i], v[i]);
+            this->calc_gravity();
+            pot[i]  += potential;
+            a[i]    += acceleration;
+            adot[i] += jerk;
+        }
+    }
+    virtual void calc_gravity() = 0;
+    virtual void print_info() {}
+    void set_name(const std::string &name)
+    {
+        this->name = name;
+    }
+    bool is_active = false;
+protected:
+    double potential;
+    double3 acceleration, jerk;
+    double3 x, v;
+    std::string name = "ext";
+    void set_coordinates(double3 x, double3 v)
+    {
+        this->x = x;
+        this->v = v;
+    }
+};
+
+class Plummer : public External_gravity {
+public:
+    Plummer(double m, double b) : m(m), b(b) {is_active=(m>0);}
+    void calc_gravity() override;
+    void print_info() override
+    {
+        if (!is_active) return;
+        printf("m_%-5s = %.4E                        b_%-5s = %.4E\n", name.c_str(), m, name.c_str(), b);
+    }
+private:
+    double m, b;
+};
+
+class Miyamoto_Nagai : public External_gravity {
+public:
+    Miyamoto_Nagai(double m, double a, double b) : m(m), a(a), b(b) {is_active=(m>0); this->set_name("disk");}
+    void calc_gravity();
+    void print_info() override
+    {
+        if (!is_active) return;
+        printf("m_%-5s = %.4E  a_%-5s = %.4E  b_%-5s = %.4E\n", name.c_str(), m, name.c_str(), a, name.c_str(), b);
+    }
+private:
+    double m, a, b;
+};
+
+class Logarithmic_halo : public External_gravity {
+public:
+    Logarithmic_halo(double v_halo, double r_halo) : v2_halo(v_halo*v_halo), r2_halo(r_halo*r_halo) {is_active=(r_halo>0); this->set_name("halo");}
+    void calc_gravity() override;
+    void print_info() override
+    {
+        if (!is_active) return;
+        printf("v_%-4s = %.6E    r_%-4s = %.4E\n", name.c_str(), sqrt(v2_halo), name.c_str(), sqrt(r2_halo));
+    }
+private:
+    double v2_halo, r2_halo;
+};
+
+class Dehnen : public External_gravity {
+public:
+    Dehnen(double m, double r, double gamma) : m(m), r(r), gamma(gamma) {is_active=(m>0);}
+    void calc_gravity() override;
+    void print_info() override
+    {
+        if (!is_active) return;
+        printf("m_%-5s = %.4E  r_%-5s = %.4E  g_%-5s = %.4E\n", name.c_str(), m, name.c_str(), r, name.c_str(), gamma);
+    }
+private:
+    double m, r, gamma;
+};
--- a/init.py
+++ b/init.py
@ -13,7 +13,7 @@ def gen_plum(N, seed=None, RMAX=10):
        X = np.sqrt(R**2 - Z**2) * np.cos(2*np.pi*X3)
        Y = np.sqrt(R**2 - Z**2) * np.sin(2*np.pi*X3)

-        Ve = np.sqrt(2)*(1.0 + R**2)**(-0.25);
+        Ve = np.sqrt(2)*(1.0 + R**2)**(-0.25)

        X4, X5 = 0, 0
        while 0.1*X5 >= X4**2*(1-X4**2)**3.5:
@ -22,9 +22,9 @@ def gen_plum(N, seed=None, RMAX=10):
        V = Ve*X4
        X6, X7 = np.random.random(2)

-        Vz = (1 - 2*X6)*V;
-        Vx = np.sqrt(V**2 - Vz**2) * np.cos(2*np.pi*X7);
-        Vy = np.sqrt(V**2 - Vz**2) * np.sin(2*np.pi*X7);
+        Vz = (1 - 2*X6)*V
+        Vx = np.sqrt(V**2 - Vz**2) * np.cos(2*np.pi*X7)
+        Vy = np.sqrt(V**2 - Vz**2) * np.sin(2*np.pi*X7)

        X, Y, Z    = np.array([X,  Y,  Z])*3*np.pi/16
        Vx, Vy, Vz = np.array([Vx, Vy, Vz])/np.sqrt(3*np.pi/16)
@ -33,6 +33,38 @@ def gen_plum(N, seed=None, RMAX=10):
        i += 1
    return particle_list

+def kepler_to_cartesian(a, e, i, Omega, w, nu, G=1.0, M=1.0):
+    def to_arrays(*args):
+        result = []
+        for arg in args: result.append(np.atleast_1d(arg))
+        return result
+    a, e, i, Omega, w, nu = to_arrays(a, e, i, Omega, w, nu) # pylint: disable=unbalanced-tuple-unpacking
+    P = [np.cos(w)*np.cos(Omega) - np.sin(w)*np.cos(i)*np.sin(Omega),
+         np.cos(w)*np.sin(Omega) + np.sin(w)*np.cos(i)*np.cos(Omega),
+         np.sin(w)*np.sin(i)]
+    Q = [-np.sin(w)*np.cos(Omega) - np.cos(w)*np.cos(i)*np.sin(Omega),
+         -np.sin(w)*np.sin(Omega) + np.cos(w)*np.cos(i)*np.cos(Omega),
+          np.sin(i)*np.cos(w)]
+    cosnu = np.cos(nu)
+    cosE = (e+cosnu)/(1+e*cosnu)
+    E = np.arccos(cosE)
+    E[nu > np.pi] = 2*np.pi - E[nu > np.pi]
+    X = a*((np.cos(E)-e)*P + np.sqrt(1-e**2)*np.sin(E)*Q)
+    V = (np.sqrt(G*M/a)/(1-e*np.cos(E)))*(-np.sin(E)*P + np.sqrt(1-e**2)*np.cos(E)*Q)
+    if X.shape[1]==1:
+        X=X[:,0]
+        V=V[:,0]
+    return X.T, V.T
+
+def generate_binary(a, e, i, Omega, w, nu, m1, m2):
+    X, V = kepler_to_cartesian(a, e, i, Omega, w, nu, M=m1+m2)
+    q = np.double(m2)/np.double(m1)
+    X1 = -q/(q+1)*X
+    V1 = -q/(q+1)*V
+    X2 =  1/(q+1)*X
+    V2 =  1/(q+1)*V
+    return X1, V1, X2, V2
+
 def write_phi_grape_config(**kargs):
    if 'file_name' in kargs: file_name = kargs['file_name']
    else: file_name = 'phi-GRAPE.cfg'
@ -59,6 +91,8 @@ def gen_mask(particle_list, frac):
        mask = np.ones(N, dtype=int)
    elif frac==1:
        mask = np.zeros(N, dtype=int)
+    elif (frac < 0) or (1 < frac):
+        raise RuntimeError('Fraction has to be between 0 and 1')
    else:
        X = particle_list[:,:3]
        V = particle_list[:,3:]
@ -86,6 +120,7 @@ if __name__=='__main__':
    parser.add_argument('--dt_bh',    type=np.double, default=.125, help='interval for BH output (bh.dat & bh_nb.dat)')
    parser.add_argument('--eta',      type=np.double, default=.01,  help='parameter for timestep determination (0.02 or 0.01)')
    parser.add_argument('--frac',     type=np.double, default=0,    help='fraction of collisional particles (by angular momentum)')
+    parser.add_argument('--bsmbh',    type=bool,      default=0,    help='generate a binary supermassive black hole (parameters hardcoded in the script)')
    args = parser.parse_args()

    try:
@ -97,9 +132,23 @@ if __name__=='__main__':
    write_phi_grape_config(**vars(args))

    particle_list = gen_plum(N, seed=args.seed)
+    m = np.ones(N)/N

-    write_particles(particle_list)
+    if args.bsmbh:
+        m1, m2 = 0.075, 0.025
+        a, e, i, Omega, w, nu = 0.001, 0.5, 0, 0, 0, 0
+        X1, V1, X2, V2 = generate_binary(a, e, i, Omega, w, nu, m1, m2)
+        m[:2]  = m1, m2
+        m[2:] = 1/(N-2)
+        particle_list[0,:3] = X1
+        particle_list[0,3:] = V1
+        particle_list[1,:3] = X2
+        particle_list[1,3:] = V2
+
+    write_particles(particle_list, m=m)

    mask = gen_mask(particle_list, args.frac)
+    if args.bsmbh:
+        mask[:2] = 3

    write_mask(mask)
--- a/io.cpp
+++ b/io.cpp
@ -0,0 +1,204 @@
+#ifdef HAS_HDF5
+#include "hdf5.h"
+#endif
+
+#include "double3.h"
+#include "io.h"
+#include <fstream>
+#include <cstring>
+#include <cmath>
+#include <vector>
+#include <algorithm>
+#include <string>
+#include <stdexcept>
+
+
+bool is_hdf5(std::string file_name)
+{
+    std::ifstream file(file_name, std::ifstream::binary);
+    const char hdf5_magic[] = "\x89HDF\x0d\x0a\x1a\x0a";
+    char buffer[8];
+    file.read(buffer, 8);
+    if (!file) throw std::runtime_error("Unable to read file " + file_name);
+    bool result = (memcmp(buffer, hdf5_magic, 8)==0);
+    file.close();
+    return result;
+}
+
+Input_data ascii_read(const std::string &file_name)
+{
+    Input_data input_data;
+    char rest[512];
+    int result;
+
+    std::ifstream file(file_name);
+    if (!file.good()) throw std::runtime_error("File " + file_name + " not found.");
+    std::string str;
+
+    std::getline(file, str);
+    result = sscanf(str.c_str(), "%d%s", &input_data.step_num, rest);
+    if (result!=1) throw std::runtime_error("Error parsing line 1: expected one integer");
+
+    std::getline(file, str);
+    result = sscanf(str.c_str(), "%d%s", &input_data.N, rest);
+    if (result!=1) throw std::runtime_error("Error parsing line 2: expected one integer");
+
+    std::getline(file, str);
+    result = sscanf(str.c_str(), "%lf%s", &input_data.t, rest);
+    if (result!=1) throw std::runtime_error("Error parsing line 3: expected one real number");
+
+    input_data.m.resize(input_data.N);
+    input_data.x.resize(input_data.N);
+    input_data.v.resize(input_data.N);
+
+    int i = -1;
+    while (std::getline(file, str)) {
+        if (++i > input_data.N) throw std::runtime_error("Error parsing line " + std::to_string(i+4) + ": particle out of range");
+        result = sscanf(str.c_str(), "%*s %lf %lf %lf %lf %lf %lf %lf%s", &input_data.m[i], &(input_data.x[i].x), &(input_data.x[i].y), &(input_data.x[i].z), &(input_data.v[i].x), &(input_data.v[i].y), &(input_data.v[i].z), rest);
+    }
+    file.close();
+    return input_data;
+}
+
+void ascii_write(const std::string file_name, const int step_num, const int N, const double t, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, int precision)
+{
+    auto file = std::ofstream(file_name);
+    if (!file.is_open()) throw std::runtime_error("Cannot open file for output");
+    int id_width = (int)log10(N-1) + 1;
+    char string_template[256], output_string[256];
+    file << step_num << '\n';
+    file << N << '\n';
+    file.precision(16);
+    file << std::scientific << t << '\n';
+    sprintf(string_template, "%%0%dd%%%d.%dE%%%d.%dE%%%d.%dE%%%d.%dE%%%d.%dE%%%d.%dE%%%d.%dE\n", id_width, precision+7, precision, precision+8, precision, precision+8, precision, precision+8, precision, precision+8, precision, precision+8, precision, precision+8, precision);
+    for (int i=0; i<N; i++) {
+        sprintf(output_string, string_template, i, m[i], x[i][0], x[i][1], x[i][2], v[i][0], v[i][1], v[i][2]);
+        file << output_string;
+    }
+    file.close();
+}
+
+Input_data h5_read(const std::string &file_name)
+{
+#ifdef HAS_HDF5
+    Input_data input_data;
+
+    // Open file and root group; count number of top level objects
+    hid_t file_id = H5Fopen(file_name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+    hid_t group_id = H5Gopen2(file_id, "/", H5P_DEFAULT);
+    H5G_info_t object_info;
+    H5Gget_info(group_id, &object_info);
+    // Iterate over objects and add the number of each "step" group into a vector
+    std::vector<int> step_num_arr;
+    for (unsigned int i=0; i<object_info.nlinks; i++) {
+        char name_cstr[256];
+        H5Gget_objname_by_idx(group_id, i, name_cstr, 256);
+        std::string name(name_cstr);
+        if (name.substr(0, 5) == "Step#") step_num_arr.push_back(std::stoi(name.substr(5, std::string::npos)));
+    }
+    // Find the highest step in the file
+    input_data.step_num = *max_element(step_num_arr.begin(), step_num_arr.end());
+
+    // Prepare to read the data sets dimensionality data.
+    char path[256];
+    hid_t attr_id, dataset_id, dataspace_id;
+    int ndims;
+    hsize_t dims[2];
+
+    sprintf(path, "/Step#%d", input_data.step_num);
+    attr_id = H5Aopen_by_name(file_id, path, "Time", H5P_DEFAULT, H5P_DEFAULT);
+    H5Aread(attr_id, H5T_NATIVE_DOUBLE, &input_data.t);
+    H5Aclose(attr_id);
+
+    sprintf(path, "/Step#%d/MASS", input_data.step_num);
+    dataset_id = H5Dopen2(file_id, path, H5P_DEFAULT);
+    dataspace_id = H5Dget_space(dataset_id);
+    ndims = H5Sget_simple_extent_ndims(dataspace_id);
+    if (ndims != 1)
+        throw std::runtime_error("Dataset MASS in Step#" + std::to_string(input_data.step_num) + " of file " + file_name + " is " + std::to_string(ndims) + "-dimensional (expected 1-dimensional)");
+    H5Sget_simple_extent_dims(dataspace_id, dims, NULL);
+    input_data.N = dims[0];
+
+    input_data.m.resize(input_data.N);
+    H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, input_data.m.data());
+    H5Sclose(dataspace_id);
+    H5Dclose(dataset_id);
+
+    sprintf(path, "/Step#%d/X", input_data.step_num);
+    dataset_id = H5Dopen2(file_id, path, H5P_DEFAULT);
+    dataspace_id = H5Dget_space(dataset_id);
+    ndims = H5Sget_simple_extent_ndims(dataspace_id);
+    if (ndims != 2) throw std::runtime_error("Bad dimensionality in " + std::string(path));
+    H5Sget_simple_extent_dims(dataspace_id, dims, NULL);
+    if ((dims[0] != input_data.N) || (dims[1] != 3))  throw std::runtime_error("Bad dimensionality");
+    input_data.x.resize(input_data.N);
+    H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, (double*)input_data.x.data());
+    H5Sclose(dataspace_id);
+    H5Dclose(dataset_id);
+
+    sprintf(path, "/Step#%d/V", input_data.step_num);
+    dataset_id = H5Dopen2(file_id, path, H5P_DEFAULT);
+    dataspace_id = H5Dget_space(dataset_id);
+    ndims = H5Sget_simple_extent_ndims(dataspace_id);
+    if (ndims != 2) throw std::runtime_error("Bad dimensionality in" + std::string(path));
+    H5Sget_simple_extent_dims(dataspace_id, dims, NULL);
+    if ((dims[0] != input_data.N) || (dims[1] != 3)) throw std::runtime_error("Bad dimensionality");
+    input_data.v.resize(input_data.N);
+    H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, (double*)input_data.v.data());
+    H5Sclose(dataspace_id);
+    H5Dclose(dataset_id);
+
+    H5Gclose(group_id);
+    H5Fclose(file_id);
+    return input_data;
+#else
+    throw std::runtime_error("h5_read was called but compiled without HDF5 support");
+#endif
+}
+
+void h5_write(const std::string file_name, const int step_num, const int N, const double t, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, const std::vector<double3> &acc, const std::vector<double3> &jrk, const int extra_mode, const bool use_double_precision)
+{
+#ifdef HAS_HDF5
+    hid_t file_id, group_id, attribute_id, dataspace_id;
+    hsize_t dims[2] = {(hsize_t)N, 3};
+
+    hid_t h5_float_type;
+    if (use_double_precision) h5_float_type = H5T_IEEE_F64LE;
+    else h5_float_type = H5T_IEEE_F32LE;
+    file_id = H5Fcreate(file_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+    char group_name[32];
+    sprintf(group_name, "/Step#%d", step_num);
+    group_id = H5Gcreate2(file_id, group_name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+    dataspace_id = H5Screate(H5S_SCALAR);
+    attribute_id = H5Acreate2 (group_id, "Time", H5T_IEEE_F64LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT);
+    H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &t);
+    H5Sclose(dataspace_id);
+
+    auto write_dataset = [&](const char dataset_name[], int ndims, double *data) {
+        hid_t dataspace_id = H5Screate_simple(ndims, dims, NULL);
+        char dataset_path[32];
+        sprintf(dataset_path, "%s/%s", group_name, dataset_name);
+        hid_t dataset_id = H5Dcreate2(file_id, dataset_path, h5_float_type, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+        H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, data);
+        H5Dclose(dataset_id);
+        H5Sclose(dataspace_id);
+    };
+
+    write_dataset("MASS", 1, (double*)m.data()); // casting away const...
+    write_dataset("X",    2, (double*)x.data());
+    write_dataset("V",    2, (double*)v.data());
+
+    bool write_pot = (extra_mode     ) & 1;
+    bool write_acc = (extra_mode >> 1) & 1;
+    bool write_jrk = (extra_mode >> 2) & 1;
+    if (write_pot) write_dataset("POT", 1, (double*)pot.data());
+    if (write_acc) write_dataset("ACC", 2, (double*)acc.data());
+    if (write_jrk) write_dataset("JRK", 2, (double*)jrk.data());
+
+    H5Gclose(group_id);
+    H5Fclose(file_id);
+    H5close(); // If we don't do that (HDF5 1.10.5) then the file isn't really closed...
+#else
+    throw std::runtime_error("h5_write was called but compiled without HDF5 support");
+#endif
+}
--- a/io.h
+++ b/io.h
@ -0,0 +1,24 @@
+#pragma once
+#include <string>
+#include <vector>
+#include "double3.h"
+
+struct Input_data {
+    int N, step_num;
+    double t;
+    std::vector<double> m;
+    std::vector<double3> x, v;
+};
+
+bool is_hdf5(std::string file_name);
+// This function is implemented independently of the HDF5 library
+
+Input_data ascii_read(const std::string &file_name);
+
+void ascii_write(const std::string file_name, const int step_num, const int N, const double t, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, int precision=10);
+
+Input_data h5_read(const std::string &file_name);
+// In case the code is compiled without HDF5 support, the implementation of this function just throws an error
+
+void h5_write(const std::string file_name, const int step_num, const int N, const double t, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, const std::vector<double3> &acc, const std::vector<double3> &jrk, const int extra_mode=0, const bool use_double_precision=true);
+// In case the code is compiled without HDF5 support, the implementation of this function just throws an error
--- a/n_bh.c
+++ b/n_bh.c
@ -1,76 +0,0 @@
-/***************************************************************************/
-/*
-  Coded by       : Peter Berczik
-  Version number : 1.0
-  Last redaction : 2011.II.20. 13:00
-*/
-
-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
-*/
-
-
-int k;
-
-double r, r2, r3, r4, RP, x[3], v[3];
-
-
-for(k=0;k<3;k++)
-  {
-  x[k] = xx1[k] - xx2[k];
-  v[k] = vv1[k] - vv2[k];
-  }
-
-r2 = SQR(x[0]) + SQR(x[1]) + SQR(x[2]) + SQR(eps_BH);
-
-r  = sqrt(r2);
-r3 = r2*r;
-r4 = r2*r2;
-
-RP = 3.0*(x[0]*v[0]+x[1]*v[1]+x[2]*v[2])/r;
-
-// Newton pot + acc & jerks
-
-*pot_n1 = -m2/r;
-*pot_n2 = -m1/r;
-
-for(k=0;k<3;k++)
-  {
-  a_n1[k] = -m2*x[k]/r3;
-  a_n2[k] =  m1*x[k]/r3;
-
-  adot_n1[k] = -m2*(v[k]/r3 - RP*x[k]/r4);
-  adot_n2[k] =  m1*(v[k]/r3 - RP*x[k]/r4);
-  }
-
-return(0);
-
-}
-/***************************************************************************/
--- a/phi-GRAPE.c
+++ b/phi-GRAPE.c
--- a/phigrape.conf
+++ b/phigrape.conf
@ -0,0 +1,260 @@
+######################
+# GENERAL PARAMETERS #
+######################
+
+# Plummer softening parameter (can be even 0)
+eps =                                   1E-4
+
+# End time of the calculation
+t_end =                                 1
+
+# Interval of snapshot files output
+dt_disk =                               0.125
+
+# Interval for the energy control output (contr.dat)
+dt_contr =                              0.125
+
+# Interval for SMBH output (bh.dat, bh_neighbors.dat, and bh_inf.dat)
+dt_bh =                                 0.125
+
+# Parameter for timestep determination
+eta =                                   0.01
+
+# Name of the input file; use "data.con" in most cases [default: data.con]
+#input_file_name =                       data.con
+
+
+# Number of devices (GRAPEs or GPUs per node) [default: 0]
+# By default, each MPI process will attempt to bind to one device. To change
+# this behaviour, set the value to the number of available devices in the
+# system. For example, if from some reason you want to run multiple MPI
+# processes on a machine with a single device, set the value to 1 and use the
+# mpirun utility (or whatever is used in your job scheduler) to launch as many
+# processes as you like.
+devices_per_node =                      1
+
+
+# The power of 2 of the maximum timestep. For example, if dt_max_power=-3 and
+# the adaptive timestep algotirhm suggests anything bigger than 0.125, the
+# timestep will be 0.125. If the 2^dt_max_power > dt_contr, the maximum timestep
+# is dt_contr. [default: -3]
+#dt_max_power =                          -3
+
+# The power of 2 of the minimum timestep. [default: -36]
+#dt_min_power =                          -36
+
+# Eta correction factor for the first step (the timestep will be calculated with
+# eta divided by eta_s_corr) [default: 4]
+#eta_s_corr =                            4.0
+
+# Eta correction factor for the black holes (the timestep for the black holes
+# will be calculated with eta divided by eta_bh_corr). [default: 4]
+#eta_bh_corr =                           4.0
+
+
+##########
+# OUTPUT #
+##########
+
+# Whether to use HDF5 format for snapshot and restart; regular ASCII snapshorts
+# are saved if false [default: false]
+#output_hdf5 =                           true
+
+# If using HDF5 output, use double precision or not [default: true] Consider
+# setting to false to save disk space. Restart file is always saved in double
+# precision.
+#output_hdf5_double_precision =          true
+
+# If using ASCII output, the number of digits after the decimal point [default: 10]
+# Restart file is saved with 16 digits after the decimal point.
+#output_ascii_precision =                6
+
+# Extra output: optionally save potential, acceleration and jerk in snapshot
+# files [default: 0]
+# This is a number between 0 and 7 that encodes the output options in the
+# following way:
+# [value] = [save jerk]*4 + [save acceleration]*2 + [save potential]
+# Example: choose 5 if it is needed to save the jerk and the potential, but not
+# the acceleration for some reason.
+# Currently implemented in HDF5 output only.
+#output_extra_mode =                     7
+
+# Whether to output a warning on the screen when the minimum time step is
+# encountered. [default: false]
+#dt_min_warning =                        false
+
+
+####################
+# EXTERNAL GRAVITY #
+####################
+
+# Remember that external gravity models are applied at the same coordinate
+# system as the particles. If the idea is to simulate a globular cluster
+# orbiting in an external field, be sure to set the initial conditions
+# appropriately (applying a shift to the coordinates and velocities).
+
+# Whether the parameters for the external gravitational field given below are in
+# physical units or Hénon units. If true, the system used is {kiloparsec, solar
+# mass, kilometre per second} [default: false]
+#ext_units_physical =                    true
+
+# If Physical units were selected, specify the simulation's unit mass (is solar
+# masses) and unit lenght (in parsec; not kiloparsec)
+# TODO: add the option to normalize using other units.
+#unit_mass   =                           4E5 # MSun
+#unit_length =                           15  # pc
+
+
+# The bulge is a Plummer potential with the following total mass and radius.
+#ext_m_bulge =                           5E9 # MSun
+#ext_b_bulge =                           1.9 # kpc
+
+
+# The disk is a Miyamoto-Nagai potential with the following total mass, scale
+# length, and scale height
+#ext_m_disk =                            6.8E10 # MSun
+#ext_a_disk =                            3.00   # kpc
+#ext_b_disk =                            0.28   # kpc
+
+
+# This halo option is yet another Plummer potential with the following total
+# mass and radius.
+#ext_m_halo_plummer =                    8E11 # MSun
+#ext_b_halo_plummer =                    245  # kpc
+
+
+# This halo option is a logarithmic potential with the following velocity and
+# radius parameters.
+#ext_log_halo_v =                        240 # km/s
+#ext_log_halo_r =                        1   # kpc
+
+# This is a spherical Dehnen model.
+#ext_dehnen_m =                          1E11 # MSun
+#ext_dehnen_r =                          2    # kpc
+#ext_dehnen_gamma =                      0.5
+
+
+###################################
+# LIVE SUPERMASSIVE BLACK HOLE(S) #
+###################################
+
+# There is special treatment for particles representing supermassive black holes
+# (SMBHs): they are integrated at every time step, they can have custom
+# softening in SMBH-SMBH interactions, and post Newtonian terms can be added to
+# the gravity.
+
+# The number of SMBH particles. Can be 0 (no SMBH), 1, or 2. [default: 0]
+#live_smbh_count =                       2
+
+# Custom softening length for SMBH-SMBH interactions (can also be zero). If
+# non-negative, the custom softening is applied. [default: -1]
+#live_smbh_custom_eps =                  0
+#TODO this is actually related only to BINARY smbh!
+
+# Output additional diagnostics about live SMBHs. [default: false]
+#live_smbh_output =                      true
+
+# Output additional diagnostics about the SMBH's (or SMBHs') nearest neighbours
+# (number could be set as shown below). [default: false]
+#live_smbh_neighbor_output =             true
+
+# Number of nearest neighbours to the SMBH (or SMBHs) to include in output. [default: 10]
+#live_smbh_neighbor_number =             10
+
+
+##################################
+# BINARY SUPERMASSIVE BLACK HOLE #
+##################################
+
+# The following parameters can be set when live_smbh_count is 2.
+
+# Output additional diagnostics about the sphere of influence (size could be set
+# as shown below). [default: false]
+#binary_smbh_influence_sphere_output =   true
+
+# The influence sphere is centred at the binary SMBH's centre of mass, and its
+# radius is the semi-major axis of the binary times the factor below. [default: 10]
+#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 {Newtonian, 1, 2, 2.5, 3, 3.5, spin}
+# Note: the first element in the array has no effect, the Newtonian force is
+# always included.
+#pn_usage =                              {1, 1, 1, 1, 0, 0, 0}
+
+# The speed of light in N-body units
+#pn_c =                                  477.12
+
+# The spin vectors of the two SMBHs. Only define these if the last component of
+# pn_usage is set to one.
+#smbh1_spin =                            {0, 0, 1}
+#smbh2_spin =                            {0, 0, 1}
+
+
+###############
+# HYBRID CODE #
+###############
+
+# The hybridization with the SCF code is enabled if the ETICS preprocessor flag
+# is defined in the when compiling.
+
+# Time intervals to calculate the SCF series expansion.
+dt_scf =                                0.015625
+
+# Name of the mask file for GRAPite [default: grapite.mask]
+#grapite_mask_file_name =                grapite.mask
+
+# Whether to write to disk a list of SCF coefficients at every dt_disk. [default: false]
+#etics_dump_coeffs =                     true
+
+# Whether to use an alternative procedure for active particle search that is
+# available in the GRAPite library. This requires the number of particles in
+# each MPI process to be exactly divisible by 32. This can substantially
+# accelerate the calculation in some circumstances [default: false]
+#grapite_active_search =                 true
+
+# Custom softening length for SMBH-star interactions in the hybrid scheme only.
+# This value (can also be zero) is used in the direct gravity calculation
+# between SMBHs (tag=3) and both core (tag=0) and halo (tag=1) stars. If
+# negative, the Plummer softening parameter (`eps`) is used in these
+# interactions. Do not confuse with `live_smbh_custom_eps`, which is the
+# softening length for SMBH-SMBH interactions, and works both in the normal and
+# hybrid schemes. [default: -1]
+#grapite_smbh_star_eps =                 1E-6
+
+# If the number of active particles in a particular bunch is bigger than this
+# threshold, then the execution is on the GPU, otherwise on the CPU. When the
+# active bunch is small, the overhead of calculating the SCF gravity on the GPU
+# makes the operation more expensive than if it is done on the CPU. [default: 32]
+#grapite_dev_exec_threshold =            512
+
+# TODO
+########
+# etics dump mode
+# scaling parameter override
+
+
+####################################
+# Negative powers of two           #
+####################################
+#  -1  1/2      0.5                #
+#  -2  1/4      0.25               #
+#  -3  1/8      0.125              #
+#  -4  1/16     0.0625             #
+#  -5  1/32     0.03125            #
+#  -6  1/64     0.015625           #
+#  -7  1/128    0.0078125          #
+#  -8  1/256    0.00390625         #
+#  -9  1/512    0.001953125        #
+# -10  1/1024   0.0009765625       #
+# -11  1/2048   0.00048828125      #
+# -12  1/4096   0.000244140625     #
+# -13  1/8192   0.0001220703125    #
+# -14  1/16384  0.00006103515625   #
+# -15  1/32768  0.000030517578125  #
+# -16  1/65536  0.0000152587890625 #
+####################################
--- a/phigrape.cpp
+++ b/phigrape.cpp
@ -0,0 +1,731 @@
+#include <algorithm>
+#include <chrono>
+#include <mpi.h>
+#include <numeric>
+#include <stdexcept>
+
+#include "black_holes.h"
+#include "config.h"
+#include "double3.h"
+#include "external.h"
+#include "grape6.h"
+#include "io.h"
+
+#ifdef ETICS
+#include "grapite.h"
+#endif
+
+namespace std::chrono {
+    struct Timer {
+        void start()
+        {
+            t_start = steady_clock::now();
+        }
+        void stop()
+        {
+            t_stop = steady_clock::now();
+            time = duration_cast<nanoseconds>(t_stop - t_start).count()*1E-9;
+        }
+        double time; // seconds
+        steady_clock::time_point t_start, t_stop;
+    };
+}
+std::chrono::Timer timer;
+
+class Calc_self_grav {
+public:
+    Calc_self_grav(const int N, const int n_loc, const int clusterid, const int npipe, const double eps)
+        : g6_calls(0), n_loc(n_loc), clusterid(clusterid), npipe(npipe), eps2(eps*eps)
+    {
+        h2.assign(N, eps2); 
+        pot_loc.resize(N);
+        acc_loc.resize(N);
+        jrk_loc.resize(N);
+    }
+    void operator()(const double t, const int n_act, std::vector<int> &ind_act, std::vector<double3> &x_act, std::vector<double3> &v_act,
+              std::vector<double>& pot, std::vector<double3> &acc, std::vector<double3> &jrk)
+    {
+        g6_set_ti(clusterid, t);
+        for (int i=0; i<n_act; i+=npipe) {
+            int nn = npipe;
+            if (n_act-i < npipe) nn = n_act - i;
+            g6calc_firsthalf(clusterid, n_loc, nn, ind_act.data()+i, (double(*)[3])&x_act[i], (double(*)[3])&v_act[i], (double(*)[3])&acc_loc[i], (double(*)[3])&jrk_loc[i], &pot_loc[i], eps2, h2.data());
+            g6calc_lasthalf( clusterid, n_loc, nn, ind_act.data()+i, (double(*)[3])&x_act[i], (double(*)[3])&v_act[i], eps2, h2.data(), (double(*)[3])&acc_loc[i], (double(*)[3])&jrk_loc[i], &pot_loc[i]);
+            g6_calls++;
+        } /* i */
+        /* Reduce the "global" vectors from "local" on all the nodes */
+        MPI_Allreduce(pot_loc.data(), pot.data(), n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+        MPI_Allreduce(acc_loc.data(), acc.data(), 3*n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+        MPI_Allreduce(jrk_loc.data(), jrk.data(), 3*n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    }
+    double g6_calls;
+private:
+    int n_loc, clusterid, npipe;
+    double eps2;
+    std::vector<double> h2;
+    std::vector<double> pot_loc; // the _loc variables are for this node only.
+    std::vector<double3> acc_loc, jrk_loc;
+};
+
+class Calc_ext_grav {
+public:
+    void add_component(External_gravity &component)
+    {
+        components.push_back(&component);
+        if (component.is_active) any_active = true;
+    }
+    void operator()(int n, const std::vector<double3> &x, const std::vector<double3> &v, std::vector<double> &pot, std::vector<double3> &acc, std::vector<double3> &jrk)
+    {
+        for (auto component : components) {
+            if (component->is_active)
+                component->apply(n, x, v, pot, acc, jrk);
+        }
+    }
+    void print_info()
+    {
+        for (auto component : components) {
+                component->print_info();
+        }
+        fflush(stdout);
+    }
+    bool any_active = false;
+private:
+    std::vector<External_gravity*> components;
+};
+
+void energy_contr(const double time_cur, const double timesteps, const double n_act_sum, const double g6_calls, int N, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, const std::vector<double> &pot_ext)
+{
+    double E_pot = 0;
+    for (int i=0; i<N; i++) E_pot += m[i]*pot[i];
+    E_pot *= 0.5;
+
+    double E_kin = 0;
+    for (int i=0; i<N; i++) E_kin += m[i]*v[i].norm2();
+    E_kin *= 0.5;
+
+    double E_pot_ext = 0;
+    for (int i=0; i<N; i++) E_pot_ext += m[i]*pot_ext[i];
+
+    double m_tot = 0;
+    double3 xcm = {0, 0, 0};
+    double3 vcm = {0, 0, 0};
+    for (int i=0; i<N; i++) {
+        m_tot += m[i];
+        xcm += m[i] * x[i];
+        vcm += m[i] * v[i];
+    }
+    xcm /= m_tot;
+    vcm /= m_tot;
+    double rcm_mod = xcm.norm();
+    double vcm_mod = vcm.norm();
+
+    double3 mom = {0, 0, 0};
+    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]);
+    }
+
+    timer.stop();
+
+    double E_tot = E_pot + E_kin + E_pot_ext;
+    
+    static double E_tot_prev;
+    if (timesteps == 0.0) E_tot_prev = E_tot;
+
+    double DE_tot = E_tot - E_tot_prev;
+
+    /* This is the only output to screen */
+    printf("%.3E %.3E  % .4E %.4E % .4E  % .4E % .4E  %.2E\n",
+            time_cur, timesteps,
+            E_pot, E_kin, E_pot_ext, E_tot, DE_tot,
+            timer.time);
+
+    fflush(stdout);
+
+    auto out = fopen("contr.dat", "a");
+    fprintf(out,"%.8E \t %.8E %.8E %.8E \t % .8E % .8E % .8E % .8E % .8E \t % .8E % .8E \t % .8E % .8E % .8E \t %.8E\n",
+            time_cur, timesteps, n_act_sum, g6_calls,
+            E_pot, E_kin, E_pot_ext,
+            E_tot, DE_tot,
+            rcm_mod, vcm_mod,
+            mom[0], mom[1], mom[2],
+            timer.time);
+    fclose(out);
+
+    E_tot_prev = E_tot;
+}
+
+class Active_search {
+    // TODO you can add pointers to t and dt at the constructor, no point giving them at get_minimum_time but without the size.
+public:
+    Active_search(const int myRank, const int n_proc, const int n_loc, const int N, bool grapite_active_search_flag)
+        : myRank(myRank), n_proc(n_proc), n_loc(n_loc), N(N), grapite_active_search_flag(grapite_active_search_flag)
+    {
+        ind_act_loc.resize(n_loc);
+    }
+    double get_minimum_time(const std::vector<double> &t, const std::vector<double> &dt)
+    {
+        double min_t_loc, min_t;
+#ifdef ETICS
+        if (grapite_active_search_flag) {
+            min_t_loc = grapite_get_minimum_time();
+        } else
+#endif
+        {
+            min_t_loc = t[myRank*n_loc]+dt[myRank*n_loc];
+            for (int j=myRank*n_loc+1; j<(myRank+1)*n_loc; j++) {
+                double tmp = t[j] + dt[j];
+                if (tmp < min_t_loc) min_t_loc = tmp;
+            }
+        }
+        /* Reduce the "global" min_t from min_t_loc "local" on all processors) */
+        MPI_Allreduce(&min_t_loc, &min_t, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+        return min_t;
+    }
+    void get_active_indices(const double min_t, const std::vector<double> &t, const std::vector<double> &dt, std::vector<int> &ind_act, int &n_act)
+    {
+#ifdef ETICS
+        if (grapite_active_search_flag) {
+            int n_act_loc;
+            grapite_active_search(min_t, ind_act_loc.data(), &n_act_loc);
+            if (myRank > 0)
+                for (int i=0; i<n_act_loc; i++)
+                    ind_act_loc[i] += myRank*n_loc;
+            int n_act_arr[256], displs[256]; // Assuming maximum of 256 processes... seems safe.
+            MPI_Allgather(&n_act_loc, 1, MPI_INT, n_act_arr, 1, MPI_INT, MPI_COMM_WORLD);
+            n_act = n_act_arr[0];
+            for (int i=1; i<n_proc; i++)
+                n_act += n_act_arr[i];
+            displs[0] = 0;
+            for (int i=1; i<n_proc; i++)
+                displs[i]=displs[i-1]+n_act_arr[i-1];
+            MPI_Allgatherv(ind_act_loc.data(), n_act_loc, MPI_INT, ind_act.data(), n_act_arr, displs, MPI_INT, MPI_COMM_WORLD);
+        } else
+#endif
+        {
+            n_act = 0;
+            for (int i=0; i<N; i++) {
+                if (t[i]+dt[i] == min_t) ind_act[n_act++] = i;
+            } /* i */
+        }
+    }
+private:
+    int myRank, n_proc, n_loc, N;
+    std::vector<int> ind_act_loc;
+    bool grapite_active_search_flag;
+};
+
+inline void calc_high_derivatives(const double dt_tmp, const double3 &a_old, const double3 &a_new, const double3 &a1_old, const double3 &a1_new, double3 &a2, double3 &a3)
+{
+    double dtinv = 1/dt_tmp;
+    double dt2inv = dtinv*dtinv;
+    double dt3inv = dt2inv*dtinv;
+
+    double3 a0mia1 = a_old-a_new;
+    double3 ad04plad12 = 4*a1_old + 2*a1_new;
+    double3 ad0plad1 = a1_old + a1_new;
+
+    a2 = -6*a0mia1*dt2inv - ad04plad12*dtinv;
+    a3 = 12*a0mia1*dt3inv + 6*ad0plad1*dt2inv;
+}
+
+inline void corrector(const double dt_tmp, const double3 &a2, const double3 &a3, double3 &x, double3 &v)
+{
+    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 += dt4over24*a2 + dt5over120*a3;
+    v += dt3over6*a2 + dt4over24*a3;
+}
+
+inline double aarseth_step(const double eta, const double dt, const double3 &a, const double3 &a1, const double3 &a2, const double3 &a3)
+{
+    double a1abs = a.norm();
+    double adot1abs = a1.norm();
+    double3 a2dot1 = a2 + dt*a3;
+    double a2dot1abs = a2dot1.norm();
+    double a3dot1abs = a3.norm();
+    return sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
+}
+
+inline double blockize_step(double dt, double dt_prev, double min_t, double dt_min, double dt_max)
+{
+    double dt_new = dt_prev;
+    if (dt < dt_min) dt_prev = dt_min;
+    if ((dt < dt_prev) && (dt > dt_min)) {
+        int power = log(dt)/M_LN2 - 1;
+        dt_new = pow(2.0, power);
+    }
+    if ((dt > 2*dt_new) && (fmod(min_t, 2*dt_new) == 0) && (2*dt_new <= dt_max)) dt_new *= 2;
+    return dt_new;
+}
+
+inline void predictor(double min_t, const int n_act, const std::vector<int> &ind_act, const std::vector<double> &t, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double3> &a, const std::vector<double3> &adot, std::vector<double3> &x_act_new, std::vector<double3> &v_act_new)
+{
+    for (int i=0; i<n_act; i++) {
+        int j_act = ind_act[i];
+        double dt = min_t - t[j_act];
+        double dt2half = 0.5*dt*dt;
+        double dt3over6 = (1./3.)*dt*dt2half;
+        x_act_new[i] = x[j_act] + v[j_act]*dt + a[j_act]*dt2half + adot[j_act]*dt3over6;
+        v_act_new[i] = v[j_act] + a[j_act]*dt + adot[j_act]*dt2half;
+    } /* i */
+}
+
+int main(int argc, char *argv[])
+{
+    timer.start();
+
+    /* INIT the rand() !!! */
+    srand(19640916);                 /* it is just my birthday :-) */
+
+    /* Init MPI */
+    MPI_Init(&argc, &argv);
+
+    /* Define the total number of processors and the Rank of each processors */
+    int n_proc, myRank;
+    const int rootRank = 0;
+    MPI_Comm_size(MPI_COMM_WORLD, &n_proc);
+    MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
+
+    /* Define the processors names */
+    int name_proc;
+    char processor_name[MPI_MAX_PROCESSOR_NAME];
+    MPI_Get_processor_name(processor_name, &name_proc);
+
+    /* Print the Rank and the names of processors */
+    printf("Rank of the processor %03d on %s \n", myRank, processor_name);
+
+    const Config config("phigrape.conf");
+    Input_data input_data;
+    if (is_hdf5(config.input_file_name)) {
+#ifndef HAS_HDF5
+        fprintf(stderr, "ERROR: input file is in HDF5 format, but the code was compiled without HDF5 support\n");
+        return -1;
+#endif
+        input_data = h5_read(config.input_file_name);
+    }
+    else
+        input_data = ascii_read(config.input_file_name);
+
+    int N = input_data.N;
+    int diskstep = input_data.step_num;
+    double time_cur = input_data.t;
+    auto &m = input_data.m;
+    auto &x = input_data.x;
+    auto &v = input_data.v;
+
+    double t_disk  = config.dt_disk*(1+floor(time_cur/config.dt_disk));
+    double t_contr = config.dt_contr*(1+floor(time_cur/config.dt_contr));
+    double t_bh    = config.dt_bh*(1+floor(time_cur/config.dt_bh));
+
+    if (myRank == rootRank) {
+        printf("\n");
+        printf("Begin the calculation of phi-GRAPE program on %03d processors\n", n_proc);
+        printf("\n");
+        printf("N       = %07d \t eps      = %.6E\n", N, config.eps);
+        printf("t_beg   = %.6E \t t_end    = %.6E\n", time_cur, config.t_end);
+        printf("dt_disk = %.6E \t dt_contr = %.6E\n", config.dt_disk, config.dt_contr);
+        printf("dt_bh   = %.6E \n", config.dt_bh);
+        printf("eta     = %.6E\n\n", config.eta);
+        printf("t_disk = %.6E   t_contr = %.6E   t_bh = %.6E\n\n", t_disk, t_contr, t_bh);
+
+        if ((diskstep == 0) && (time_cur == 0)) {
+            FILE *out = fopen("contr.dat", "w");
+            fclose(out);
+            if (config.live_smbh_output && (config.live_smbh_count > 0)) {
+                out = fopen("bh.dat", "w");
+                fclose(out);
+            }
+            if ((config.live_smbh_neighbor_output) && (config.live_smbh_count > 0)) {
+                out = fopen("bh_neighbors.dat", "w");
+                fclose(out);
+            }
+        }
+    } /* if (myRank == rootRank) */
+
+    double normalization_mass=1, normalization_length=1, normalization_velocity=1;
+    if (config.ext_units_physical) {
+        normalization_mass = 1/config.unit_mass;
+        normalization_length = 1000/config.unit_length;
+        normalization_velocity = 1.52484071426404437233e+01*sqrt(config.unit_length/config.unit_mass);
+    }
+    Calc_ext_grav calc_ext_grav;
+    Plummer ext_bulge(config.ext_m_bulge*normalization_mass, config.ext_b_bulge*normalization_length);
+    ext_bulge.set_name("bulge");
+    calc_ext_grav.add_component(ext_bulge);
+    Miyamoto_Nagai ext_disk(config.ext_m_disk*normalization_mass, config.ext_a_disk*normalization_length, config.ext_b_disk*normalization_length);
+    calc_ext_grav.add_component(ext_disk);
+    Plummer ext_halo_plummer(config.ext_m_halo_plummer*normalization_mass, config.ext_b_halo_plummer*normalization_length);
+    ext_halo_plummer.set_name("halo");
+    calc_ext_grav.add_component(ext_halo_plummer);
+    Logarithmic_halo ext_log_halo(config.ext_log_halo_v*normalization_velocity, config.ext_log_halo_r*normalization_length);
+    calc_ext_grav.add_component(ext_log_halo);
+    Dehnen ext_dehnen(config.ext_dehnen_m*normalization_mass, config.ext_dehnen_r*normalization_length, config.ext_dehnen_gamma);
+    calc_ext_grav.add_component(ext_dehnen);
+    if (myRank == rootRank) calc_ext_grav.print_info();
+
+    /* some local settings for G6a boards */
+    int clusterid, numGPU;
+    if (config.devices_per_node==0) {
+        MPI_Comm shmcomm;
+        MPI_Comm_split_type(MPI_COMM_WORLD, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL, &shmcomm);
+        MPI_Comm_size(shmcomm, &numGPU);
+        MPI_Comm_rank(shmcomm, &clusterid);
+    } else {
+        numGPU = config.devices_per_node;
+        clusterid = myRank % numGPU;
+    }
+    printf("Rank of the processor %03d : Number of GPUs %01d : Cluster ID %01d \n", myRank, numGPU, clusterid);
+    fflush(stdout);
+
+    /* init the local GRAPEs */
+    g6_open(clusterid);
+    int npipe = g6_npipes();
+    g6_set_tunit(51);
+    g6_set_xunit(51);
+
+    bool grapite_active_search_flag = false;
+#ifdef ETICS
+    grapite_set_dev_exec_threshold(config.grapite_dev_exec_threshold);
+    grapite_active_search_flag = config.grapite_active_search;
+#endif
+
+    int n_loc = N/n_proc;
+#ifdef ETICS
+    grapite_read_particle_tags(N, config.grapite_mask_file_name.c_str(), myRank, n_loc);
+    grapite_set_dt_exp(config.dt_scf);
+    grapite_set_t_exp(time_cur);
+#endif
+
+    const double dt_min = pow(2, config.dt_min_power);
+    std::vector<int> ind(N);
+    std::iota(begin(ind), end(ind), 0);
+    /* load the nj particles to the G6 */
+    double3 zeros = {0, 0, 0}; // Dummy; can't really be const because of the GRAPE interface.
+    for (int k=0; k<n_loc; k++) {
+        int j = k + myRank*n_loc;
+        g6_set_j_particle(clusterid, k, ind[j], time_cur, dt_min, m[j], zeros, zeros, zeros, v[j], x[j]);
+    } /* k */
+
+#ifdef ETICS
+    double etics_length_scale;
+    if (myRank == rootRank) etics_length_scale = grapite_get_length_scale(N, m.data(), (double(*)[3])x.data(), (double(*)[3])v.data()); // We don't want all ranks to do it, because they need to write a file and might confuse each other
+    MPI_Bcast(&etics_length_scale, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
+    grapite_set_length_scale(etics_length_scale);
+
+    int grapite_cep_index = grapite_get_cep_index();
+    if (grapite_cep_index >= 0) {
+        double3 xcm, vcm, xdc, vdc;
+        grapite_calc_center(N, m.data(), (double(*)[3])x.data(), (double(*)[3])v.data(), xcm, vcm, xdc, vdc);
+        x[grapite_cep_index] = xdc;
+        v[grapite_cep_index] = vdc;
+        grapite_update_cep(time_cur, xdc, vdc, zeros, zeros);
+    }
+
+    if (config.grapite_smbh_star_eps >= 0) grapite_set_eps_bh(config.grapite_smbh_star_eps);
+#endif
+
+    std::vector<double3> a(N), adot(N);
+    std::vector<double> pot(N);
+
+    /* define the all particles as a active on all the processors for the first time grav calc. */
+    Calc_self_grav calc_self_grav(N, n_loc, clusterid, npipe, config.eps);
+    calc_self_grav(time_cur, N, ind, x, v, pot, a, adot);
+
+    Black_hole_physics black_hole_physics;
+    std::vector<Particle_ref> smbh_list;
+    if (config.live_smbh_count >= 1)
+        black_hole_physics = Black_hole_physics(config.live_smbh_count, myRank, rootRank);
+    else if (config.live_smbh_count >= 2) {
+        if (config.live_smbh_custom_eps >= 0) {
+#ifdef ETICS
+            double eps = (config.grapite_smbh_star_eps >= 0)?config.grapite_smbh_star_eps:config.eps;
+#else
+            double eps = config.eps;
+#endif
+            black_hole_physics.set_softening(eps, config.live_smbh_custom_eps);
+            for (int i = 0; i < config.live_smbh_count; i++)
+                smbh_list.emplace_back(m[i], x[i], v[i], pot[i], a[i], adot[i]);
+            black_hole_physics.adjust_softening(smbh_list);
+        }
+    }
+    if (config.binary_smbh_pn) {
+        throw std::runtime_error("This is the triple+ SMBH version, it cannot do PN yet!");
+        #if 0
+        black_hole_physics.set_post_newtonian(config.pn_c, config.pn_usage.data());
+        if (config.pn_usage[6]) black_hole_physics.set_spins(config.smbh1_spin.data(), config.smbh2_spin.data());
+        black_hole_physics.adjust_post_newtonian(dt_min, a[0], a[1], adot[0], adot[1]);
+        #endif
+    }
+
+    std::vector<double> pot_ext(N, 0.);
+    calc_ext_grav(N, x, v, pot_ext, a, adot);
+
+    double timesteps=0, n_act_sum=0;
+    /* Energy control... */
+    if (myRank == rootRank) energy_contr(time_cur, timesteps, n_act_sum, calc_self_grav.g6_calls, N, m, x, v, pot, pot_ext);
+
+#ifdef ETICS
+    if (config.etics_dump_coeffs && (diskstep==0)) {
+        char out_fname[256];
+        sprintf(out_fname, "coeffs.%06d.%02d.dat", 0, myRank);
+        grapite_dump(out_fname, 2);
+    }
+
+    if (grapite_cep_index >= 0) {
+        double3 xcm, vcm, xdc, vdc;
+        grapite_calc_center(N, m.data(), (double(*)[3])x.data(), (double(*)[3])v.data(), xcm, vcm, xdc, vdc);
+        x[grapite_cep_index] = xdc;
+        v[grapite_cep_index] = vdc;
+        grapite_update_cep(time_cur, xdc, vdc, a[grapite_cep_index], adot[grapite_cep_index]);
+    }
+#endif
+
+    const double dt_max = std::min({config.dt_disk, config.dt_contr, pow(2, config.dt_max_power)});
+    std::vector<double> dt(N);
+    /* Define initial timestep for all particles on all nodes */
+    for (int j=0; j<N; j++) {
+        double a2_mod = a[j].norm2();
+        double adot2_mod = adot[j].norm2();
+
+        double dt_tmp, eta_s = config.eta/config.eta_s_corr;
+        if (adot2_mod==0) dt_tmp = eta_s; // That's weird, when will we have such a case?
+        else              dt_tmp = eta_s*sqrt(a2_mod/adot2_mod);
+
+        int power = log(dt_tmp)/log(2.0) - 1;
+
+        dt_tmp = pow(2.0, (double)power);
+
+        if (dt_tmp < dt_min) dt_tmp = dt_min;
+        if (dt_tmp > dt_max) dt_tmp = dt_max;
+
+        dt[j] = dt_tmp;
+
+        if (config.dt_min_warning && (myRank == 0)) {
+            if (dt[j] == dt_min) {
+                printf("!!! Warning0: dt = dt_min = %.6E \t ind = %07d \n", dt[j], ind[j]);
+                fflush(stdout);
+            }
+        }
+    } /* j */
+
+    if (config.live_smbh_count > 0) {
+        double min_dt = *std::min_element(begin(dt), end(dt));
+        for (int i=0; i<config.live_smbh_count; i++) dt[i] = min_dt;
+    }
+
+    /* load the new values for particles to the local GRAPEs */
+    for (int k=0; k<n_loc; k++) {
+        int j = k + myRank*n_loc;
+        g6_set_j_particle(clusterid, k, ind[j], time_cur, dt[j], m[j], zeros, adot[j]*(1./6.), a[j]*0.5, v[j], x[j]);
+    } /* k */
+
+    timesteps = 0.0; // Why won't those two be long long instead of double + should include the zeroth step
+    n_act_sum = 0.0;
+
+
+    std::vector<int> ind_act(N);
+    std::vector<double3> x_act_new(N), v_act_new(N), a_act_new(N), adot_act_new(N);
+    std::vector<double> t(N, time_cur), pot_act_new(N);
+    std::vector<double> pot_act_ext(N, 0.);
+
+    // Functors for the main integration loop
+    Active_search active_search(myRank, n_proc, n_loc, N, grapite_active_search_flag);
+    Binary_smbh_influence_sphere_output binary_smbh_influence_sphere_output(config.binary_smbh_influence_radius_factor, N, m, x, v, pot, dt);
+    Write_bh_nb_data write_bh_nb_data(config.live_smbh_neighbor_number, config.live_smbh_count, N, m, x, v);
+    if (myRank == rootRank) {
+        if (config.live_smbh_output) black_hole_physics.write_bh_data(time_cur, config.live_smbh_count, m, x, v, pot, a, adot, dt);
+        if (config.live_smbh_neighbor_output) write_bh_nb_data(time_cur);
+    } /* if (myRank == rootRank) */
+
+    /* The main integration loop */
+    while (time_cur <= config.t_end) {
+
+        /* Define the minimal time and the active particles on all the nodes */
+        double min_t = active_search.get_minimum_time(t, dt);
+
+        /* Get indices of all particles that will be active in this bunch */
+        int n_act;
+        active_search.get_active_indices(min_t, t, dt, ind_act, n_act);
+
+        /* Find the BH(s) indices in the active list */
+        smbh_list.clear();
+#ifdef ETICS
+        /* Unlike with the simple active search, with GPU accelerated GRAPite
+        active search, the list of active indices is not sorted. */
+        int n_bh = config.live_smbh_count;
+        if (config.grapite_active_search && (n_bh>0)) {
+            int act_def_grapite_bh_count = 0;
+            for (int i=0; i<n_act; i++) {
+                if (ind_act[i]<n_bh) {
+                    smbh_list.emplace_back(m[ind_act[i]], x_act_new[i], v_act_new[i], pot_act_new[i], a_act_new[i], adot_act_new[i]);
+                    if (act_def_grapite_bh_count++ == n_bh) break;
+                }
+            }
+        }
+#else
+        for (int i = 0; i < config.live_smbh_count; i++)
+            smbh_list.emplace_back(m[ind_act[i]], x_act_new[i], v_act_new[i], pot_act_new[i], a_act_new[i], adot_act_new[i]);
+#endif
+
+        /* predict the active particles positions etc... on all the nodes */
+        predictor(min_t, n_act, ind_act, t, x, v, a, adot, x_act_new, v_act_new);
+
+        /* Calculate gravity on active particles */
+        calc_self_grav(min_t, n_act, ind_act, x_act_new, v_act_new, pot_act_new, a_act_new, adot_act_new);
+
+        if (config.live_smbh_count >= 2) {
+            if (config.live_smbh_custom_eps >= 0) black_hole_physics.adjust_softening(smbh_list);
+            #if 0
+            if (config.binary_smbh_pn) black_hole_physics.adjust_post_newtonian(dt[i_bh1], a_act_new[i_bh1], a_act_new[i_bh2], adot_act_new[i_bh1], adot_act_new[i_bh2]);
+            #endif
+        }
+
+        /* Calculate gravity on active particles due to external forces */
+        if (calc_ext_grav.any_active) {
+            std::fill_n(begin(pot_act_ext), n_act, 0);
+            calc_ext_grav(n_act, x_act_new, v_act_new, pot_act_ext, a_act_new, adot_act_new);
+        }
+
+        /* correct the active particles positions etc... on all the nodes */
+        double min_dt = dt_max; // notice that min_dt is not the same as dt_min; this one is to store the minimum timestep among currently active particles
+        for (int i=0; i<n_act; i++) {
+            int j_act = ind_act[i];
+            double dt_cur = dt[j_act];
+
+            double3 a2, a3;
+            calc_high_derivatives(dt_cur, a[j_act], a_act_new[i], adot[j_act], adot_act_new[i], a2, a3);
+
+            corrector(dt_cur, a2, a3, x_act_new[i], v_act_new[i]);
+
+            //TODO make beautiful
+            double eta_curr;
+            if ((config.live_smbh_count > 0) && (ind_act[i] < config.live_smbh_count)) eta_curr = config.eta/config.eta_bh_corr;
+            else eta_curr = config.eta;
+
+            double dt_new = aarseth_step(eta_curr, dt_cur, a_act_new[i], adot_act_new[i], a2, a3);
+
+            dt_new = blockize_step(dt_new, dt_cur, min_t, dt_min, dt_max);
+
+            if (config.dt_min_warning && (myRank == 0)) {
+                if (dt_new == dt_min) {
+                    printf("!!! Warning1: dt_act = dt_min = %.6E \t ind_act = %07d time_cur=%.16E\n", dt_cur, ind_act[i], time_cur);
+                    fflush(stdout);
+                }
+            }
+            if (dt_new < min_dt) min_dt = dt_new;
+
+            x[j_act] = x_act_new[i];
+            v[j_act] = v_act_new[i];
+            t[j_act] = min_t;
+            dt[j_act] = dt_new;
+            pot[j_act] = pot_act_new[i];
+            pot_ext[j_act] = pot_act_ext[i];
+            a[j_act] = a_act_new[i];
+            adot[j_act] = adot_act_new[i];
+        } /* i */
+
+        /* define the min. dt over all the act. part. and set it also for the BH... */
+        for (int i=0; i < config.live_smbh_count; i++) dt[i] = min_dt;
+
+        if (config.binary_smbh_influence_sphere_output && (myRank == rootRank))
+            binary_smbh_influence_sphere_output(ind_act, n_act, timesteps, time_cur);
+
+        /* load the new values for active particles to the local GRAPE's */
+        for (int i=0; i<n_act; i++) {
+#ifdef ETICS
+            if (ind_act[i] == grapite_cep_index) grapite_update_cep(t[grapite_cep_index], x[grapite_cep_index], v[grapite_cep_index], a[grapite_cep_index], adot[grapite_cep_index]); // All ranks should do it.
+#endif
+            int cur_rank = ind_act[i]/n_loc;
+            if (myRank == cur_rank) {
+                int j_act = ind_act[i];
+                int address = ind_act[i] - myRank*n_loc;
+                g6_set_j_particle(clusterid, address, ind_act[i], t[j_act], dt[j_act], m[j_act], zeros, adot[j_act]*(1./6.), a[j_act]*0.5, v[j_act], x[j_act]);
+            } /* if (myRank == cur_rank) */
+        } /* i */
+
+        /* Current time set to min_t */
+        time_cur = min_t;
+        timesteps += 1.0;
+        n_act_sum += n_act;
+
+        if (time_cur >= t_bh) {
+            if (myRank == rootRank) {
+                /* Write BH data... */
+                if (config.live_smbh_output) black_hole_physics.write_bh_data(time_cur, config.live_smbh_count, m, x, v, pot, a, adot, dt);
+
+                /* Write BH NB data... */
+                if (config.live_smbh_neighbor_output) write_bh_nb_data(time_cur);
+
+            } /* if (myRank == rootRank) */
+
+            t_bh += config.dt_bh;
+        } /* if (time_cur >= t_bh) */
+
+        if (time_cur >= t_contr) {
+            if (myRank == rootRank) {
+                energy_contr(time_cur, timesteps, n_act_sum, calc_self_grav.g6_calls, N, m, x, v, pot, pot_ext);
+                /* write cont data */
+                if (config.output_hdf5) h5_write("data.con", diskstep, N, time_cur, m, x, v, pot, a, adot, 0, true);
+                else ascii_write("data.con", diskstep, N, time_cur, m, x, v, 16);
+            } /* if (myRank == rootRank) */
+
+#ifdef ETICS
+            // We are /inside/ a control step, so all particles must be
+            // synchronized; we can safely calculate their density centre. The
+            // acceleration and jerk currently in the memory are for the
+            // predicted position of the CEP, by calling grapite_calc_center we
+            // "correct" the position and velocity, but not the gravity at that
+            // point.
+            if (grapite_cep_index >= 0) {
+                double3 xcm, vcm, xdc, vdc;
+                grapite_calc_center(N, m.data(), (double(*)[3])x.data(), (double(*)[3])v.data(), xcm, vcm, xdc, vdc);
+                x[grapite_cep_index] = xdc;
+                v[grapite_cep_index] = vdc;
+                grapite_update_cep(time_cur, xdc, vdc, a[grapite_cep_index], adot[grapite_cep_index]);
+            }
+#endif
+
+            t_contr += config.dt_contr;
+        } /* if (time_cur >= t_contr) */
+
+        if (time_cur >= t_disk) {
+            char out_fname[256];
+            diskstep++;
+            if (myRank == rootRank) {
+                sprintf(out_fname, "%06d", diskstep);
+                if (config.output_hdf5) h5_write(std::string(out_fname) + ".h5", diskstep, N, time_cur, m, x, v, pot, a, adot, config.output_extra_mode, config.output_hdf5_double_precision);
+                else ascii_write(std::string(out_fname) + ".dat", diskstep, N, time_cur, m, x, v, config.output_ascii_precision);
+            } /* if (myRank == rootRank) */
+
+#ifdef ETICS
+            if (config.etics_dump_coeffs) {
+                sprintf(out_fname, "coeffs.%06d.%02d.dat", diskstep, myRank);
+                grapite_dump(out_fname, 2);
+            }
+#endif
+            t_disk += config.dt_disk;
+        } /* if (time_cur >= t_disk) */
+    } /* while (time_cur < t_end) */
+
+    /* close the local GRAPEs */
+    timer.stop();
+    g6_close(clusterid);
+
+    double g6_calls_sum;
+    MPI_Reduce(&calc_self_grav.g6_calls, &g6_calls_sum, 1, MPI_DOUBLE, MPI_SUM, rootRank, MPI_COMM_WORLD);
+    if (myRank == rootRank) {
+        /* Write some output for the timestep annalize... */
+        printf("\n");
+        printf("timesteps = %.0f   Total sum of integrated part. = %.0f   g6_calls on all nodes = %.0f \n", timesteps, n_act_sum, g6_calls_sum);
+        printf("\n");
+        printf("Real Speed = %.3f GFlops \n", 57.0*N*n_act_sum/(timer.time)/1.0E+09);
+        fflush(stdout);
+    } /* if (myRank == rootRank) */
+
+    /* Finalize the MPI work */
+    MPI_Finalize();
+}
--- a/pn_bh.c
+++ b/pn_bh.c
@ -1,765 +0,0 @@
-/***************************************************************************/
-/*
-  Coded by       : Peter Berczik (on the base of Gabor Kupi original PN code)
-  Version number : 2.0 SPIN
-  Last redaction : 2012.V.07. 11:16
-*/
-
-int calc_force_pn_BH(double m1, double xx1[], double vv1[], double spin1[],
-                     double m2, double xx2[], double vv2[], double spin2[],
-                     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)
-*/
-
-
-
-int j, k;
-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 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;
-
-double A3, B3, A3_5, B3_5, AK6, BK6, AK7, BK7;
-double A3D, A3_5D, B3D, B3_5D, ADK6, BDK6, ADK7, BDK7;
-
-int Van_Spin=0;
-int Van_QM=0;
-
-double DM, S1[3], SPIN[3][2], S2[3], KSS[3], KSSIG[3], XS[3], XA[3], NCV[3], NCS[3], NCSIG[3], 
-       VCS[3], VCSIG[3], SDNCV, SIGDNCV, NDV, XS2, XA2, NXA, NXS, VDS, VDSIG, NDS, NDSIG, 
-       C1_5[3], C2[3], C2_5[3];
-double LABS, LU[3], S1DLU, S2DLU, SU1[3],  SV1[3],  SS1[3], SS2[3],  SU2[3],  SV2[3]; 
-double AT[3], NDOT[3], NVDOT, NDOTCV[3], NCA[3];
-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], 
-       S2DirNew[3], rS1p, rS2p, S1p[3], S2p[3], Np[3]; 
-
-
-Van_Spin = usedOrNot[6];		// Van vagy nincs SPIN szamolas...
-
-
-for(k=0;k<3;k++)
-  {
-  SPIN[k][0] = spin1[k];
-  SPIN[k][1] = spin2[k];
-  }
-
-
-for(j=0;j<7;j++)
-  {
-for(k=0;k<3;k++)
-  {
-  a_pn1[j][k] = 0.0;   adot_pn1[j][k] = 0.0;
-  a_pn2[j][k] = 0.0;   adot_pn2[j][k] = 0.0;
-  }
-  }
-
-
-// Speed of light "c" and its powers
-
-c_1 = CCC_NB;
-
-c_2 = SQR(c_1);
-c_4 = SQR(c_2);
-c_5 = c_4*c_1;
-c_6 = c_5*c_1;
-c_7 = c_6*c_1;
-
-
-// Mass parameters
-
-M   = m1+m2;
-eta = m1*m2/(M*M);
-nu  = m1/m2;
-
-for(k=0;k<3;k++)
-  {
-  x[k] = xx1[k] - xx2[k];
-  v[k] = vv1[k] - vv2[k];
-  }
-
-
-r2 = SQR(x[0]) + SQR(x[1]) + SQR(x[2]);
-r  = sqrt(r2);
-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;
-
-
-// Newton accelerations
-
-for(k=0;k<3;k++) N[k] = x[k]/r;
-
-// PN accelerations
-
-AK2 = 0.0; BK2 = 0.0;
-AK4 = 0.0; BK4 = 0.0;
-AK5 = 0.0; BK5 = 0.0;
-AK6 = 0.0; BK6 = 0.0;
-AK7 = 0.0; BK7 = 0.0;
-
-for(k=0;k<3;k++)
-  {
-  C1_5[k] = 0.0; 
-  C2[k]   = 0.0; 
-  C2_5[k] = 0.0;
-  QM[k]   = 0.0;
-  }
-
-
-if(usedOrNot[1] == 1)		// PN1		~1/c^2
-  {
-  A1 = 2.0*(2.0+eta)*MOR-(1.0+3.0*eta)*V1_V22 +1.5*eta*RP*RP;
-  B1 = 2.0*(2.0-eta)*RP;
-
-  AK2 = A1/c_2;
-  BK2 = B1/c_2;
-  }
-
-if(usedOrNot[2] == 1)		// PN2		~1/c^4
-  {
-  A2 = -0.75*(12.0+29.0*eta)*MOR*MOR-eta*(3.0-4.0*eta)*V1_V22*V1_V22-1.875*eta*(1.0-3.0*eta)*RP*RP*RP*RP+0.5*eta*(13.0-4.0*eta)*MOR*V1_V22+(2.0+25.0*eta+2.0*eta*eta)*MOR*RP*RP+1.5*eta*(3.0-4.0*eta)*V1_V22*RP*RP;
-  B2 = -0.5*RP*((4.0+41.0*eta+8.0*eta*eta)*MOR-eta*(15.0+4.0*eta)*V1_V22+3.0*eta*(3.0+2.0*eta)*RP*RP);
-
-  AK4 = A2/c_4;
-  BK4 = B2/c_4;
-  }
-
-if(usedOrNot[3] == 1)		// PN2.5	~1/c^5
-  {
-  A2_5 = 1.6*eta*MOR*RP*(17.0*MOR/3.0+3.0*V1_V22);
-  B2_5 = -1.6*eta*MOR*(3.0*MOR+V1_V22);
-
-  AK5 = A2_5/c_5;
-  BK5 = B2_5/c_5;
-  }
-
-if(usedOrNot[4] == 1)		// PN3  	~1/c^6
-  {
-  A3 = MOR*MOR*MOR*(16.0+(1399.0/12.0-41.0*PI2/16.0)*eta+
-       71.0*eta*eta/2.0)+eta*(20827.0/840.0+123.0*PI2/64.0-eta*eta)
-       *MOR*MOR*V1_V22-(1.0+(22717.0/168.0+615.0*PI2/64.0)*eta+
-       11.0*eta*eta/8.0-7.0*eta*eta*eta)*MOR*MOR*RP*RP-
-       0.25*eta*(11.0-49.0*eta+52.0*eta*eta)*V1_V22*V1_V22*V1_V22+
-       35.0*eta*(1.0-5.0*eta+5.0*eta*eta)*RP*RP*RP*RP*RP*RP/16.0-
-       0.25*eta*(75.0+32.0*eta-40.0*eta*eta)*MOR*V1_V22*V1_V22-
-       0.5*eta*(158.0-69.0*eta-60.0*eta*eta)*MOR*RP*RP*RP*RP+
-       eta*(121.0-16.0*eta-20.0*eta*eta)*MOR*V1_V22*RP*RP+
-       3.0*eta*(20.0-79.0*eta+60.0*eta*eta)*V1_V22*V1_V22*RP*RP/8.0-
-       15.0*eta*(4.0-18.0*eta+17.0*eta*eta)*V1_V22*RP*RP*RP*RP/8.0;
-
-  B3 = RP*((4.0+(5849.0/840.0+123.0*PI2/32.0)*eta-25.0*eta*eta-
-       8.0*eta*eta*eta)*MOR*MOR+eta*(65.0-152.0*eta-48.0*eta*eta)*
-       V1_V22*V1_V22/8.0+15.0*eta*(3.0-8.0*eta-2.0*eta*eta)*RP*RP*RP*RP/8.0+
-       eta*(15.0+27.0*eta+10.0*eta*eta)*MOR*V1_V22-eta*(329.0+177.0*eta+
-       108.0*eta*eta)*MOR*RP*RP/6.0-
-       3.0*eta*(16.0-37.0*eta-16.0*eta*eta)*V1_V22*RP*RP/4.0);
-
-  AK6 = A3/c_6;
-  BK6 = B3/c_6;
-  }
-
-if(usedOrNot[5] == 1)		// PN3.5  	~1/c^7
-  {
-  A3_5 = MOR*eta*(V1_V22*V1_V22*(-366.0/35.0-12.0*eta)+V1_V22*RP*RP*(114.0+12.0*eta)-112.0*RP*RP*RP*RP+MOR*(V1_V22*(-692.0/35.0+724.0*eta/15.0)+RP*RP*(-294.0/5.0-376.0*eta/5.0)+MOR*(-3956.0/35.0-184.0*eta/5.0)));
-  B3_5 = 8.0*eta*MOR*((1325.0+546.0*eta)*MOR*MOR/42.0+(313.0+42.0*eta)*V1_V22*V1_V22/28.0+75.0*RP*RP*RP*RP-(205.0+777.0*eta)*MOR*V1_V22/42.0+(205.0+424.0*eta)*MOR*RP*RP/12.0-3.0*(113.0+2.0*eta)*V1_V22*RP*RP/4.0)/5.0;
-
-  AK7 = A3_5/c_7;
-  BK7 = B3_5/c_7;
-  }
-
-
-// Spin accelerations
-
-  if(Van_Spin==1)
-    {
-    DM = m1 - m2;
-
-    Spin1Abs2 = 0.0;
-    Spin2Abs2 = 0.0;
-    rS1 = 0.0;
-    rS2 = 0.0;
-
-    for(k=0;k<3;k++)
-      {
-      Spin1Abs2 += SPIN[k][0]*SPIN[k][0];		// normalizalt spin
-      Spin2Abs2 += SPIN[k][1]*SPIN[k][1];
-      rS1 += N[k]*S1Dir[k];
-      rS2 += N[k]*S2Dir[k];
-
-      S1[k] = SPIN[k][0]*m1*m1/c_1;			// fizikai spin 
-      S2[k] = SPIN[k][1]*m2*m2/c_1;
-      KSS[k] = S1[k]+S2[k];
-      KSSIG[k] = M*(S2[k]/m2-S1[k]/m1);
-      XS[k]  = 0.5*(SPIN[k][0]+SPIN[k][1]);
-      XA[k]  = 0.5*(SPIN[k][0]-SPIN[k][1]);
-      }
-
-   Spin1Abs = sqrt(Spin1Abs2);
-   Spin2Abs = sqrt(Spin2Abs2);
-
-   for(k=0;k<3;k++)
-     {
-     S1Dir[k] = SPIN[k][0]/Spin1Abs;
-     S2Dir[k] = SPIN[k][1]/Spin2Abs;
-     }
-
-              
-  //NCV crossproduct of N[k] and relative v = N[k]Xv[j]
-  NCV[0] = N[1]*v[2] - N[2]*v[1];  
-  NCV[1] = N[2]*v[0] - N[0]*v[2];  
-  NCV[2] = N[0]*v[1] - N[1]*v[0];  
-       
-  //NCS crossproduct of N[k] and KSS = N[k]XKSS
-  NCS[0] = N[1]*KSS[2] - N[2]*KSS[1];  
-  NCS[1] = N[2]*KSS[0] - N[0]*KSS[2];  
-  NCS[2] = N[0]*KSS[1] - N[1]*KSS[0];  
-       
-  //NCSIG crossproduct of N[k] and KSSIG = N[k]XKSSIG
-  NCSIG[0] = N[1]*KSSIG[2] - N[2]*KSSIG[1];  
-  NCSIG[1] = N[2]*KSSIG[0] - N[0]*KSSIG[2];  
-  NCSIG[2] = N[0]*KSSIG[1] - N[1]*KSSIG[0];  
-       
-  //VCS crossproduct of v[k] and KSS = v[k]XKSS
-  VCS[0] = v[1]*KSS[2] - v[2]*KSS[1];  
-  VCS[1] = v[2]*KSS[0] - v[0]*KSS[2];  
-  VCS[2] = v[0]*KSS[1] - v[1]*KSS[0];  
-       
-  //VCSIG crossproduct of v[k] and KSSIG = v[k]XKSSIG
-  VCSIG[0] = v[1]*KSSIG[2] - v[2]*KSSIG[1];  
-  VCSIG[1] = v[2]*KSSIG[0] - v[0]*KSSIG[2];  
-  VCSIG[2] = v[0]*KSSIG[1] - v[1]*KSSIG[0];  
-       
-  SDNCV = KSS[0]*NCV[0]+KSS[1]*NCV[1]+KSS[2]*NCV[2];
-       
-  SIGDNCV = KSSIG[0]*NCV[0]+KSSIG[1]*NCV[1]+KSSIG[2]*NCV[2];
-       
-  NDV = N[0]*v[0] + N[1]*v[1] + N[2]*v[2];
-       
-  XS2 = XS[0]*XS[0]+XS[1]*XS[1]+XS[2]*XS[2];
-  XA2 = XA[0]*XA[0]+XA[1]*XA[1]+XA[2]*XA[2];
-       
-  NXA = N[0]*XA[0]+N[1]*XA[1]+N[2]*XA[2];
-  NXS = N[0]*XS[0]+N[1]*XS[1]+N[2]*XS[2];
-       
-  VDS = v[0]*KSS[0]+v[1]*KSS[1]+v[2]*KSS[2];
-  VDSIG = v[0]*KSSIG[0]+v[1]*KSSIG[1]+v[2]*KSSIG[2];
-       
-  NDS = N[0]*KSS[0]+N[1]*KSS[1]+N[2]*KSS[2];
-  NDSIG = N[0]*KSSIG[0]+N[1]*KSSIG[1]+N[2]*KSSIG[2];
-       
-  for(k=0;k<3;k++)
-    {
-    C1_5[k] = (N[k]*(12.0*SDNCV+6.0*DM*SIGDNCV/M)+9.0*NDV*NCS[k]+3.0*DM*NDV*NCSIG[k]/M -7.0*VCS[k]-3.0*DM*VCSIG[k]/M)/r3;
-    C2[k] = -MOR*MOR*MOR/r*3.0*eta*(N[k]*(XS2-XA2-5.0*NXS*NXS+5.0*NXA*NXA)+2.0*(XS[k]*NXS-XA[k]*NXA));
-    C2_5[k] = (N[k]*(SDNCV*(-30.0*eta*NDV*NDV+24.0*eta*V1_V22-MOR*(38.0+25.0*eta))+DM/M*SIGDNCV*(-15.0*eta*NDV*NDV+12.0*eta*V1_V22
-              -MOR*(18.0+14.5*eta)))+NDV*v[k]*(SDNCV*(-9.0+9.0*eta)+DM/M* SIGDNCV*(-3.0+6.0*eta))+NCV[k]*(NDV*VDS*(-3.0+3.0*eta)
-              -8.0*MOR*eta*NDS-DM/M*(4.0*MOR*eta*NDSIG+3.0*NDV*VDSIG))+NDV*NCS[k]*(-22.5*eta*NDV*NDV+21.0*eta*V1_V22-MOR*(25.0+15.0*eta))
-              +DM/M*NDV*NCSIG[k]*(-15.0*eta*NDV*NDV+12.0*eta*V1_V22-MOR*(9.0+8.5*eta))+VCS[k]*(16.5*eta*NDV*NDV+MOR*(21.0+9.0*eta)
-              -14.0*eta*V1_V22)+DM/M*VCSIG[k]*(9.0*eta*NDV*NDV-7.0*eta*V1_V22+MOR*(9.0+4.5*eta)))/r3;
-
-  if(Van_QM==1)
-    {
-
-    if(m1>m2)
-      {
-      QMAux2_1[k] = (1.0-5.0*rS1*rS1)*N[k]+2.0*rS1*S1Dir[k];
-      QMAux2_2[k] = (1.0-5.0*rS2*rS2)*N[k]+2.0*rS2*S2Dir[k];
-      QMAux1[k] = Spin1Abs2*QMAux2_1[k]/nu+Spin2Abs2*QMAux2_2[k]*nu;
-      QM[k] = -1.5*MOR*MOR*MOR*eta*QMAux1[k]/r;
-      }
-    else
-      {
-      QMAux2_1[k] = (1.0-5.0*rS2*rS2)*N[k]+2.0*rS2*S2Dir[k];
-      QMAux2_2[k] = (1.0-5.0*rS1*rS1)*N[k]+2.0*rS1*S1Dir[k];
-      QMAux1[k] = Spin2Abs2*QMAux2_1[k]/nu+Spin1Abs2*QMAux2_2[k]*nu;
-      QM[k] = -1.5*MOR*MOR*MOR*eta*QMAux1[k]/r;
-      }
-
-    } /* if(Van_QM==1) */
-
-    } /* k */
-
-  } /* if(Van_Spin==1) */   
-
-
-
-
-
-for(k=0;k<3;k++)
-  {
-
-if(usedOrNot[0] == 1)		// PN0 (Newton)	~1/c^0
-  {
-  a_pn1[0][k] = -m2*x[k]/r3;
-  a_pn2[0][k] =  m1*x[k]/r3;
-  }
-
-if(usedOrNot[1] == 1)		// PN1		~1/c^2
-  {
-  a_pn1[1][k] =  ((AK2*N[k] + BK2*v[k])/r2)*m2;
-  a_pn2[1][k] = -((AK2*N[k] + BK2*v[k])/r2)*m1;
-  }
-
-if(usedOrNot[2] == 1)		// PN2		~1/c^4
-  {
-  a_pn1[2][k] =  ((AK4*N[k] + BK4*v[k])/r2)*m2;
-  a_pn2[2][k] = -((AK4*N[k] + BK4*v[k])/r2)*m1;
-  }
-
-if(usedOrNot[3] == 1)		// PN2.5	~1/c^5
-  {
-  a_pn1[3][k] =  ((AK5*N[k] + BK5*v[k])/r2)*m2;
-  a_pn2[3][k] = -((AK5*N[k] + BK5*v[k])/r2)*m1;
-  }
-
-if(usedOrNot[4] == 1)		// PN3  	~1/c^6
-  {
-  a_pn1[4][k] =  ((AK6*N[k] + BK6*v[k])/r2)*m2;
-  a_pn2[4][k] = -((AK6*N[k] + BK6*v[k])/r2)*m1;
-  }
-
-if(usedOrNot[5] == 1)		// PN3.5	~1/c^7
-  {
-  a_pn1[5][k] =  ((AK7*N[k] + BK7*v[k])/r2)*m2;
-  a_pn2[5][k] = -((AK7*N[k] + BK7*v[k])/r2)*m1;
-  }
-
-if(Van_Spin == 1)		// All the SPIN terms
-  {
-  a_pn1[6][k] +=  (C1_5[k]/c_2 + C2[k]/c_4 + C2_5[k]/c_4 + QM[k]/c_4)*m2/M;
-  a_pn2[6][k] += -(C1_5[k]/c_2 + C2[k]/c_4 + C2_5[k]/c_4 + QM[k]/c_4)*m1/M;
-  }
-
-  A[k] = MOR*((AK2+AK4+AK5+AK6+AK7)*N[k] + (BK2+BK4+BK5+BK6+BK7)*v[k])/r + C1_5[k]/c_2 + C2[k]/c_4 + C2_5[k]/c_4 + QM[k]/c_4;
-  }
-
-// PN accelerations
-
-
-
-
-// PN jerks
-
-for(k=0;k<3;k++)
-  {
-  AT[k] = A[k] - MOR*N[k]/r;			// miert van AT - ? 
-  }
-
-/*
-AT[0] = A[0];
-AT[1] = A[1];
-AT[2] = A[2];
-*/
-        
-RPP = V1_V22/r + AT[0]*N[0]+AT[1]*N[1] + AT[2]*N[2] - RP*RP/r;
-VA = AT[0]*v[0] + AT[1]*v[1] + AT[2]*v[2];
-
-for(k=0;k<3;k++) NDOT[k] = (v[k]-N[k]*RP)/r;
-
-NVDOT = NDOT[0]*v[0]+NDOT[1]*v[1]+NDOT[2]*v[2]+N[0]*AT[0]+N[1]*AT[1]+N[2]*AT[2];
-
-//NDOTCV crossproduct of NDOT[k] and relative v = NDOT[k]Xv[j]
-NDOTCV[0] = NDOT[1]*v[2] - NDOT[2]*v[1];  
-NDOTCV[1] = NDOT[2]*v[0] - NDOT[0]*v[2];  
-NDOTCV[2] = NDOT[0]*v[1] - NDOT[1]*v[0];  
-
-//NCA crossproduct of N and AT = N[k]XAT[j]
-NCA[0] = N[1]*AT[2] - N[2]*AT[1];  
-NCA[1] = N[2]*AT[0] - N[0]*AT[2];  
-NCA[2] = N[0]*AT[1] - N[1]*AT[0];  
-
-ADK2 = 0.0; BDK2 = 0.0;
-ADK4 = 0.0; BDK4 = 0.0;
-ADK5 = 0.0; BDK5 = 0.0;
-ADK6 = 0.0; BDK6 = 0.0;
-ADK7 = 0.0; BDK7 = 0.0;
-
-
-for(k=0;k<3;k++)
-  {
-  C1_5D[k] = 0.0;
-  C2D[k]   = 0.0;
-  C2_5D[k] = 0.0;
-  QMD[k]   = 0.0;
-  }
-
-if(usedOrNot[1] == 1)		// PN1		~1/c^2
-  {
-  A1D = -2.0*(2.0+eta)*MOR*RP/r - 2.0*(1.0+3.0*eta)*VA + 3.0*eta*RP*RPP;
-  B1D =  2.0*(2.0-eta)*RPP;
-
-  ADK2 = A1D/c_2;
-  BDK2 = B1D/c_2;
-  }
-
-if(usedOrNot[2] == 1)		// PN2		~1/c^4
-  {
-  A2D =  1.5*(12.0+29.0*eta)*MOR*MOR*RP/r -eta*(3.0-4.0*eta)*4.0*V1_V22*VA - 7.5*eta*(1.0-3.0*eta)*RPP -0.5*eta*(13.0-4.0*eta)*MOR*RP*V1_V22/r+eta*(13.0-4.0*eta)*MOR*VA -(2.0+25.0*eta+2.0*eta*eta)*MOR*RP*RP*RP/r+2.0*(2.0+25.0*eta+2.0*eta*eta)*MOR*RP*RPP + 3.0*eta*(3.0-4.0*eta)*VA*RP*RP + 3.0*eta*(3.0-4.0*eta)*V1_V22*RP*RPP;
-  B2D = -0.5*RPP*((4.0+41.0*eta+8.0*eta*eta)*MOR - eta*(15.0+4.0*eta)*V1_V22+3.0*eta*(3.0+2.0*eta)*RP*RP) - 0.5*RP*(-(4.0+41.0*eta+8.0*eta*eta)*MOR*RP/r - 2.0*eta*(15.0+4.0*eta)*VA + 6.0*eta*(3.0+2.0*eta)*RP*RPP);
-
-  ADK4 = A2D/c_4;
-  BDK4 = B2D/c_4;
-  }
-
-if(usedOrNot[3] == 1)		// PN2.5	~1/c^5
-  {
-  A2_5D = -1.6*eta*MOR*RP*RP*(17.0/3.0*MOR+3.0*V1_V22)/r +1.6*eta*MOR*RPP*(17.0/3.0*MOR+3.0*V1_V22)+1.6*eta*MOR*RP*(-17.0*MOR*RP/3.0/r+6.0*VA);
-  B2_5D = 1.6*eta*MOR*RP*(3.0*MOR+V1_V22)/r - 1.6*eta*MOR*(-3.0*MOR*RP/r+2.0*VA);
-
-  ADK5 = A2_5D/c_5;
-  BDK5 = B2_5D/c_5;
-  }
-
-if(usedOrNot[4] == 1)		// PN3  	~1/c^6
-  {
-  A3D =  6.0*eta*RP*RP*RP*RP*RP*RPP*(35.0-175.0*eta+175.0*eta*eta)/16.0 + eta*(4.0*RP*RP*RP*RPP*V1_V22 + 2.0*RP*RP*RP*RP*VA)*(-15.0+135.0*eta/2.0-255.0*eta*eta/4.0)/2.0 + eta*(2.0*RP*RPP*V1_V22*V1_V22+4.0*RP*RP*V1_V22*VA)/2.0*(15.0-237.0*eta/2.0+45.0*eta*eta) + 6.0*V1_V22*V1_V22*VA*eta*(-11.0/4.0-49.0*eta/4.0-13.0*eta*eta) + MOR*(4.0*RP*RP*RP*RPP*eta*(-79.0+69.0/2.0*eta+30.0*eta*eta) + eta*(2.0*RP*RPP*V1_V22+2.0*RP*RP*VA)*(121.0-16.0*eta-20.0*eta*eta)+4.0*V1_V22*VA*eta*(-75.0/4.0-8.0*eta+10.0*eta*eta)) - MOR*RP*((-79.0+69.0*eta/2.0+30.0*eta*eta)*RP*RP*RP*RP*eta+eta*RP*RP*V1_V22*(121.0-16.0*eta-20.0*eta*eta)+eta*V1_V22*V1_V22*(-75.0/4.0-8.0*eta+10.0*eta*eta))/r - 2.0*MOR*MOR*RP*(RP*RP*((-1.0-615.0*PI2*eta/64.0)-22717.0*eta/168.0-11.0*eta*eta/8.0+7.0*eta*eta*eta)+eta*V1_V22*((20827.0/840.0+123.0*PI2/64.0)-eta*eta))/r + MOR*MOR*(2.0*RP*RPP*((-1.0-615*PI2*eta/64.0)-22717.0*eta/168.0-11.0*eta*eta/8.0+7*eta*eta*eta)+2.0*eta*VA*((20827.0/840.0 +123.0*PI2/64.0)-eta*eta)) - 3.0*MOR*MOR*MOR*RP*(16.0+(1399.0/12.0-41.0*PI2/16.0)*eta+71.0*eta*eta/2.0)/r;
-  B3D = 75.0*RP*RP*RP*RP*RPP*eta*(3.0/8.0-eta-.25*eta*eta)+eta*(3.0*RP*RP*RPP*V1_V22+2.0*RP*RP*RP*VA)*(-12.0+111.0*eta/4.0+12.0*eta*eta)+eta*(RPP*V1_V22*V1_V22+4.0*RP*V1_V22*VA)*(65.0/8.0-19.0*eta-6.0*eta*eta)-MOR*RP*(RP*RP*RP*eta*(-329.0/6.0-59.0*eta/2.0-18.0*eta*eta)+RP*V1_V22*eta*(15.0+27.0*eta+10.0*eta*eta))/r+MOR*(3.0*RP*RP*RPP*eta*(-329.0/6.0-59.0*eta/2.0-18.0*eta*eta)+eta*(RPP*V1_V22+2.0*RP*VA)*(15.0+27.0*eta+10.0*eta*eta))-2.0*MOR*MOR*RP*(RP*((4.0+123.0*PI2*eta/32.0)+5849.0*eta/840.0-25.0*eta*eta-8.0*eta*eta*eta))/r+MOR*MOR*(RPP*((4.0+123.0*PI2*eta/32.0)+5849.0/840.0*eta-25.0*eta*eta-8.0*eta*eta*eta));
-
-  ADK6 = A3D/c_6;
-  BDK6 = B3D/c_6;
-  }
-
-if(usedOrNot[5] == 1)		// PN3.5	~1/c^7
-  {
-  A3_5D = MOR*eta*(-RP*(V1_V22*V1_V22*(-366.0/35.0-12.0*eta)+V1_V22*RP*RP*(114.0+12.0*eta)+RP*RP*RP*RP*(-112.0))/r+4.0*V1_V22*VA*(-366.0/35.0-12.0*eta)+2.0*(VA*RP*RP+RP*RPP*V1_V22)*(114.0+12.0*eta)+4.0*RP*RP*RP*RPP*(-112.0)+MOR*(2.0*VA*(-692.0/35.0+724.0*eta/15.0)+2.0*RP*RPP*(-294.0/5.0-376.0*eta/5.0)-2.0*RP*(V1_V22*(-692.0/35.0+724.0*eta/15.0)+RP*RP*(-294.0/5.0-376.0*eta/5.0))/r-3.0*MOR*RP*(-3956.0/35.0-184.0*eta/5.0)/r));
-  B3_5D = MOR*eta*(4.0*V1_V22*VA*(626.0/35.0+12.0*eta/5.0)+2.0*(VA*RP*RP+V1_V22*RP*RPP)*(-678.0/5.0-12.0*eta/5.0)+4.0*RP*RP*RP*RPP*120.0-RP*(V1_V22*V1_V22*(626.0/35.0+12.0*eta/5.0)+V1_V22*RP*RP*(-678.0/5.0-12.0*eta/5.0)+120.0*RP*RP*RP*RP)/r+MOR*(2.0*VA*(-164.0/21.0-148.0*eta/5.0)+2*RP*RPP*(82.0/3.0+848.0*eta/15.0)-2.0*RP*(V1_V22*(-164.0/21-148.0*eta/5.0)+RP*RP*(82.0/3.0+848.0*eta/15.0))/r-3.0*MOR*RP*(1060.0/21.0+104.0*eta/5.0)/r));
-
-  ADK7 = A3_5D/c_7;
-  BDK7 = B3_5D/c_7;
-  }
-
-
-
-
-
-  if(Van_Spin==1)
-    {
-
-    //L crossproduct of x[k] and relative v = x[k]Xv[j]
-    L[0] = x[1]*v[2] - x[2]*v[1];  
-    L[1] = x[2]*v[0] - x[0]*v[2];  
-    L[2] = x[0]*v[1] - x[1]*v[0];  
-
-    LABS = sqrt(L[0]*L[0]+L[1]*L[1]+L[2]*L[2]);
-
-    LU[0] = L[0]/LABS;
-    LU[1] = L[1]/LABS;
-    LU[2] = L[2]/LABS;
-
-    S1DLU = S1[0]*LU[0]+S1[1]*LU[1]+S1[2]*LU[2];
-    S2DLU = S2[0]*LU[0]+S2[1]*LU[1]+S2[2]*LU[2];
-
-    for(k=0;k<3;k++)
-      {
-      SU1[k] = MOR*eta*(N[k]*(-4.0*VDS-2.0*DM/M*VDSIG)+ v[k]*(3.0*NDS+DM/M*NDSIG)+NDV*(2.0*KSS[k]+DM/M*KSSIG[k])) /r;
-      SV1[k] = MOR*(N[k]*(VDSIG*(-2.0+4.0*eta)-2.0*DM/M*VDS)+ v[k]*(NDSIG*(1.0-eta)+DM/M*NDS)+NDV*(KSSIG[k]*(1.0- 2.0*eta)+ DM/M*KSS[k]))/r;
-
-      SS1[k] = 0.5*(L[k]*(4.0+3.0*(m2/m1))+ (S2[k]-3.0*S2DLU*LU[k]))/r3;
-      SS2[k] = 0.5*(L[k]*(4.0+3.0*(m1/m2))+ (S1[k]-3.0*S1DLU*LU[k]))/r3;
-         
-      SU2[k] = MOR*eta/r*(N[k]*(VDS*(-2.0*V1_V22+3.0*NDV*NDV- 6.0*eta*NDV*NDV+7.0*MOR-8.0*eta*MOR)-14.0*MOR*NDS*NDV+ DM/M*VDSIG*eta*(-3.0*NDV*NDV-4.0*MOR)+DM/M*MOR*NDSIG*NDV* (2.0-eta/2.))+v[k]*(NDS*(2.0*V1_V22-4.0*eta*V1_V22-3.0*NDV* NDV+7.5*eta*NDV*NDV+4.0*MOR-6.0*eta*MOR)+VDS*NDV*(2.0- 6.0*eta)+ DM/M*NDSIG*(-1.5*eta*V1_V22+3.0*eta*NDV*NDV-MOR-3.5*eta* MOR)-3.0*DM/M*VDSIG*NDV*eta)+KSS[k]*NDV*(V1_V22-2.0*eta* V1_V22-1.5*NDV*NDV+3.0*eta*NDV*NDV-MOR+2.0*eta*MOR)+ DM/M*KSSIG[k]*NDV*(-eta*V1_V22+1.5*eta*NDV*NDV+ (eta-1.)*MOR));
-      SV2[k] = MOR/r*(N[k]*(VDSIG*eta*(-2.0*V1_V22+6.0*eta*NDV* NDV+(3.0+8.0*eta)*MOR)+MOR*NDSIG*NDV*(2.0-22.5*eta+2.0* eta*eta)+ DM/M*VDS*eta*(-3.0*NDV*NDV-4.0*MOR)+DM/M*MOR*NDS*NDV*(2.0- 0.5*eta))+v[k]*(NDSIG*(0.5*eta*V1_V22+2.0*eta*eta*V1_V22- 4.5*eta*eta*NDV*NDV+(4.5*eta-1.0+8.0*eta*eta)*MOR)+VDSIG*NDV* eta*(6.0*eta-1.)-3.0*DM/M*VDS*NDV*eta+DM/M*NDS*(-1.5* eta*V1_V22+ 3.0*eta*NDV*NDV-(1.0+3.5*eta)*MOR))+KSSIG[k]*NDV*(2.0*eta*eta* V1_V22-3.0*eta*eta*NDV*NDV+(-1.0+4.0*eta-2.0*eta*eta)*MOR)+ DM/M*KSS[k]*NDV*(-eta*V1_V22+1.5*eta*NDV*NDV+(-1.0+eta)* MOR));
-      }
-
-    //SS1 crossproduct of SS1 and S1 = SS1[k]XS1[j]
-    SS1aux[0] =   SS1[1]*S1[2] - SS1[2]*S1[1];  
-    SS1aux[1] =   SS1[2]*S1[0] - SS1[0]*S1[2];  
-    SS1aux[2] =   SS1[0]*S1[1] - SS1[1]*S1[0];  
-
-    SS1[0] = SS1aux[0];  
-    SS1[1] = SS1aux[1];  
-    SS1[2] = SS1aux[2];  
-
-    //SS2 crossproduct of SS2 and S2 = SS2[k]XS2[j]
-    SS2aux[0] =   SS2[1]*S2[2] - SS2[2]*S2[1];  
-    SS2aux[1] =   SS2[2]*S2[0] - SS2[0]*S2[2];  
-    SS2aux[2] =   SS2[0]*S2[1] - SS2[1]*S2[0];  
-
-    SS2[0] = SS2aux[0];  
-    SS2[1] = SS2aux[1];  
-    SS2[2] = SS2aux[2];  
-
-    SPINPrev[0][0] = SPIN[0][0];
-    SPINPrev[1][0] = SPIN[1][0];
-    SPINPrev[2][0] = SPIN[2][0];
-
-    SPINPrev[0][1] = SPIN[0][1];
-    SPINPrev[1][1] = SPIN[1][1];
-    SPINPrev[2][1] = SPIN[2][1];
-
-    SpinPrev2_1 = SPINPrev[0][0]*SPINPrev[0][0] + SPINPrev[1][0]*SPINPrev[1][0] + SPINPrev[2][0]*SPINPrev[2][0];
-    SpinPrev2_2 = SPINPrev[0][1]*SPINPrev[0][1] + SPINPrev[1][1]*SPINPrev[1][1] + SPINPrev[2][1]*SPINPrev[2][1];
-
-    SPSPP1 = 0.0;
-    SPSPP2 = 0.0;
-
-    Spin1AbsNew2 = 0.0;
-    Spin2AbsNew2 = 0.0;
-
-    for(k=0;k<3;k++)
-      {
-      SU[k] = SU1[k]/c_2 + SU2[k]/c_4 + (SS1[k] + SS2[k])/c_2;
-      SV[k] = SV1[k]/c_2 + SV2[k]/c_4+M*(SS2[k]/m2-SS1[k]/ m1)/c_2;
-      
-      KSS[k] = KSS[k] + SU[k]*dt_bh;						// integrate for dt_bh timestep 
-      KSSIG[k] = KSSIG[k] + SV[k]*dt_bh;
-      
-      SPIN[k][0] = m1*(M*KSS[k]-m2*KSSIG[k])/M/M/m1/m1*c_1;       
-      SPIN[k][1] = m2*(M*KSS[k]+m1*KSSIG[k])/M/M/m2/m2*c_1;
-      Spin1AbsNew2 += SPIN[k][0]*SPIN[k][0];
-      Spin2AbsNew2 += SPIN[k][1]*SPIN[k][1];
-      XAD[k] = 0.5/(M*M*m1*m2)*(-SU[k]*M*DM-SV[k]*(m1*m1+m2*m2));
-      XSD[k] = 0.5/(M*M*m1*m2)*(SU[k]*M*M+SV[k]*(m1*m1-m2*m2));
- 
-      if(m1>m2)
-        {
-        SPSPP1 += SPINPrev[k][0]*(SPIN[k][0]-SPINPrev[k][0])/dt_bh;
-        SPSPP2 += SPINPrev[k][1]*(SPIN[k][1]-SPINPrev[k][1])/dt_bh;
-        } 
-      else
-        {
-        SPSPP1 += SPINPrev[k][1]*(SPIN[k][1]-SPINPrev[k][1])/dt_bh;
-        SPSPP2 += SPINPrev[k][0]*(SPIN[k][0]-SPINPrev[k][0])/dt_bh;
-        }
-      }
-
-    Spin1AbsNew = sqrt(Spin1AbsNew2);
-    Spin2AbsNew = sqrt(Spin2AbsNew2);
-
-    for(k=0;k<3;k++)
-      {
-      S1DirNew[k] = SPIN[k][0]/Spin1AbsNew;
-      S2DirNew[k] = SPIN[k][1]/Spin2AbsNew;
-      }
-
-
-    //NDOTCS crossproduct of NDOT and KSS = NDOT[k]XKSS[j]
-    NDOTCS[0] =   NDOT[1]*KSS[2] - NDOT[2]*KSS[1];  
-    NDOTCS[1] =   NDOT[2]*KSS[0] - NDOT[0]*KSS[2];  
-    NDOTCS[2] =   NDOT[0]*KSS[1] - NDOT[1]*KSS[0];  
-    //NCSU crossproduct of N and SU = N[k]XSU[j]
-    NCSU[0] =   N[1]*SU[2] - N[2]*SU[1];  
-    NCSU[1] =   N[2]*SU[0] - N[0]*SU[2];  
-    NCSU[2] =   N[0]*SU[1] - N[1]*SU[0];  
-    //NDOTCSIG crossproduct of NDOT and KSSIG = NDOT[k]XKSSIG[j]
-    NDOTCSIG[0] =   NDOT[1]*KSSIG[2] - NDOT[2]*KSSIG[1];  
-    NDOTCSIG[1] =   NDOT[2]*KSSIG[0] - NDOT[0]*KSSIG[2];  
-    NDOTCSIG[2] =   NDOT[0]*KSSIG[1] - NDOT[1]*KSSIG[0];  
-    //NCSV crossproduct of N and SV = N[k]XSV[j]
-    NCSV[0] =   N[1]*SV[2] - N[2]*SV[1];  
-    NCSV[1] =   N[2]*SV[0] - N[0]*SV[2];  
-    NCSV[2] =   N[0]*SV[1] - N[1]*SV[0];  
-    //ACS crossproduct of AT and KSS = AT[k]XKSS[j]
-    ACS[0] =   AT[1]*KSS[2] - AT[2]*KSS[1];  
-    ACS[1] =   AT[2]*KSS[0] - AT[0]*KSS[2];  
-    ACS[2] =   AT[0]*KSS[1] - AT[1]*KSS[0];  
-    //VCSU crossproduct of relative v and SU = v[k]XSU[j]
-    VCSU[0] =   v[1]*SU[2] - v[2]*SU[1];  
-    VCSU[1] =   v[2]*SU[0] - v[0]*SU[2];  
-    VCSU[2] =   v[0]*SU[1] - v[1]*SU[0];  
-    //ACSIG crossproduct of AT and KSSIG = AT[k]XKSSIG[j]
-    ACSIG[0] =   AT[1]*KSSIG[2] - AT[2]*KSSIG[1];  
-    ACSIG[1] =   AT[2]*KSSIG[0] - AT[0]*KSSIG[2];  
-    ACSIG[2] =   AT[0]*KSSIG[1] - AT[1]*KSSIG[0];  
-    //VCSV crossproduct of relative v and SV = v[k]XSV[j]
-    VCSV[0] =   v[1]*SV[2] - v[2]*SV[1];  
-    VCSV[1] =   v[2]*SV[0] - v[0]*SV[2];  
-    VCSV[2] =   v[0]*SV[1] - v[1]*SV[0];  
-
-    SNVDOT = SU[0]*NCV[0]+SU[1]*NCV[1]+SU[2]*NCV[2]+ KSS[0]*NDOTCV[0]+KSS[1]*NDOTCV[1]+KSS[2]*NDOTCV[2]+ KSS[0]*NCA[0]+KSS[1]*NCA[1]+KSS[2]*NCA[2];
-
-    SIGNVDOT = SV[0]*NCV[0]+SV[1]*NCV[1]+SV[2]*NCV[2]+ KSSIG[0]*NDOTCV[0]+KSSIG[1]*NDOTCV[1]+KSSIG[2]*NDOTCV[2]+ KSSIG[0]*NCA[0]+KSSIG[1]*NCA[1]+KSSIG[2]*NCA[2];
-
-    NSDOT = NDOT[0]*KSS[0]+NDOT[1]*KSS[1]+NDOT[2]*KSS[2]+ N[0]*SU[0]+N[1]*SU[1]+N[2]*SU[2];
-    NSIGDOT = NDOT[0]*KSSIG[0]+NDOT[1]*KSSIG[1]+NDOT[2]*KSSIG[2]+ N[0]*SV[0]+N[1]*SV[1]+N[2]*SV[2];
-    VSDOT = AT[0]*KSS[0]+AT[1]*KSS[1]+AT[2]*KSS[2]+ v[0]*SU[0]+v[1]*SU[1]+v[2]*SU[2];
-    VSIGDOT = AT[0]*KSSIG[0]+AT[1]*KSSIG[1]+AT[2]*KSSIG[2]+ v[0]*SV[0]+v[1]*SV[1]+v[2]*SV[2];
-
-    NXSDOT = NDOT[0]*XS[0]+NDOT[1]*XS[1]+NDOT[2]*XS[2]+ N[0]*XSD[0]+N[1]*XSD[1]+N[2]*XSD[2];
-    NXADOT = NDOT[0]*XA[0]+NDOT[1]*XA[1]+NDOT[2]*XA[2]+ N[0]*XAD[0]+N[1]*XAD[1]+N[2]*XAD[2];
-
-    rS1p = -rS1*NDV/r;
-    rS2p = -rS2*NDV/r;
-
-    for(k=0;k<3;k++)
-      {
-      S1p[k] = (S1DirNew[k] - S1Dir[k])/dt_bh;
-      S2p[k] = (S2DirNew[k] - S2Dir[k])/dt_bh;
-
-      rS1p += v[k]*S1Dir[k]/r + N[k]*S1p[k];
-      rS2p += v[k]*S2Dir[k]/r + N[k]*S2p[k];
-
-      Np[k] = (v[k] - N[k]*NDV)/r;
-      }
-
-  for(k=0;k<3;k++)
-    { 
-    C1_5D[k] = -3.0*RP/r*C1_5[k]+(NDOT[k]*(12.0*SDNCV+6.0*DM/M* SIGDNCV)+N[k]*(12.0*SNVDOT+6.0*DM/M*SIGNVDOT)+9.0*NVDOT* NCS[k]+9.0*NDV*(NDOTCS[k]+NCSU[k])+3.0*DM/M*(NVDOT*NCSIG[k]+ NDV*(NDOTCSIG[k]+NCSV[k]))-7.0*(ACS[k]+VCSU[k])-3.0*DM/M* (ACSIG[k]+VCSV[k]))/(r3);
-    C2D[k] = -4.0*RP/r*C2[k]-MOR*MOR*MOR*3.0*eta/r*(NDOT[k]* (XS2-XA2-5.0*NXS*NXS+5.0*NXA*NXA)+N[k]*(2.0*(XS[0]*XSD[0]+ XS[1]*XSD[1]+XS[2]*XSD[2]-XA[0]*XAD[0]-XA[1]*XAD[1]- XA[2]*XAD[2])-10.0*NXS*NXSDOT+10.0*NXA*NXADOT)+2.0*(XSD[k]* NXS+XS[k]*NXSDOT-XAD[k]*NXA-XA[k]*NXADOT));
-    C2_5D[k] = -3.0*RP/r*C2_5[k]+(NDOT[k]*(SDNCV*(-30.0*eta* NDV*NDV+24.0*eta*V1_V22-MOR*(38.0+25.0*eta))+DM/M*SIGDNCV* (-15.0*eta*NDV*NDV+12.0*eta*V1_V22-MOR*(18.0+14.5*eta)))+ N[k]*(SNVDOT*(-30.0*eta*NDV*NDV+24.0*eta*V1_V22-MOR* (38.0+25.0*eta))+SDNCV*(-60.0*eta*NDV*NVDOT+48.0*eta*VA+ MOR*RP/r*(38.0+25.0*eta))+DM/M*SIGNVDOT*(-15.0*eta*NDV* NDV+12.0*eta*V1_V22-MOR*(18.0+14.5*eta))+DM/M*SIGDNCV* (-30.0*eta*NDV*NVDOT+24.0*eta*VA+MOR*RP/r*(18.0+14.5*eta)))+ (NVDOT*v[k]+NDV*AT[k])*(SDNCV*(-9.0+9.0*eta)+DM/M*SIGDNCV* (-3.0+6.0*eta))+NDV*v[k]*(SNVDOT*(-9.0+9.0*eta)+DM/M* SIGNVDOT*(-3.0+6.0*eta))+(NDOTCV[k]+NCA[k])*(NDV*VDS*(-3.0+ 3.0*eta)-8.0*MOR*eta*NDS-DM/M*(4.0*MOR*eta*NDSIG+3.0*NDV*VDSIG) )+NCV[k]*((NVDOT*VDS+NDV*VSDOT)*(-3.0+3.0*eta)-8.0*eta*MOR* (NSDOT-RP/r*NDS)-DM/M*(4.0*eta*MOR*(NSIGDOT-RP/r*NDSIG)+ 3.0*(NVDOT*VDSIG+NDV*VSIGDOT)))+(NVDOT*NCS[k]+NDV* (NDOTCS[k]+NCSU[k]))*(-22.5*eta*NDV*NDV+21.0*eta*V1_V22- MOR*(25.0+15.0*eta))+NDV*NCS[k]*(-45.0*eta*NDV*NVDOT+42.0*eta* VA+MOR*RP/r*(25.0+15.0*eta))+DM/M*(NVDOT*NCSIG[k]+NDV* (NDOTCSIG[k]+NCSV[k]))*(-15.0*eta*NDV*NDV+12.0*eta*V1_V22- MOR*(9.0+8.5*eta))+DM/M*NDV*NCSIG[k]*(-30.0*eta*NDV*NVDOT+ 24.0*eta*VA+MOR*RP/r*(9.0+8.5*eta))+(ACS[k]+VCSU[k])* (16.5*eta*NDV*NDV+MOR*(21.0+9.0*eta)-14.0*eta*V1_V22)+ VCS[k]*(33.0*eta*NDV*NVDOT-MOR*RP/r*(21.0+9.0*eta)- 28.0*eta*VA)+DM/M*(ACSIG[k]+VCSV[k])*(9.0*eta*NDV*NDV- 7.0*eta*V1_V22+MOR*(9.0+4.5*eta))+DM/M*VCSIG[k]*(18.0* eta*NDV*NVDOT-14.0*eta*VA-MOR*RP/r*(9.0+4.5*eta)))/ (r3);
-
-
-  if(Van_QM==1)
-    {
-
-    if(m1>m2)
-      {
-      QMD[k] = -1.5*MOR*MOR*MOR*eta*(-4.0*RP*QMAux1[k]/r2+(  2.0*(SPSPP1*QMAux2_1[k]/nu+SPSPP2*QMAux2_2[k]*nu) + SpinPrev2_1*(-10.0*rS1*rS1p*N[k]+(1.0-5.0*rS1*rS1)*Np[k]+2.0*rS1p*S1Dir[k]+2.0*rS1*S1p[k])/nu + SpinPrev2_2*(-10.0*rS2*rS2p*N[k]+(1.0-5.0*rS2*rS2)*Np[k]+2.0*rS2p*S2Dir[k]+2.0*rS2*S2p[k])*nu )/r);
-      } 
-    else
-      {
-      QMD[k] = -1.5*MOR*MOR*MOR*eta*(-4.0*RP*QMAux1[k]/r2+(  2.0*(SPSPP2*QMAux2_1[k]/nu+SPSPP1*QMAux2_2[k]*nu) + SpinPrev2_2*(-10.0*rS2*rS2p*N[k]+(1.0-5.0*rS2*rS2)*Np[k]+2.0*rS2p*S2Dir[k]+2.0*rS2*S2p[k])/nu + SpinPrev2_1*(-10.0*rS1*rS1p*N[k]+(1.0-5.0*rS1*rS1)*Np[k]+2.0*rS1p*S1Dir[k]+2.0*rS1*S1p[k])*nu )/r);
-      }
-
-    } /* if(Van_QM==1) */
-
-    } /* k */
-
-
-    } /* 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];
-  spin2[k] = SPIN[k][1];
-  }
-
-
-
-
-
-
-
-for(k=0;k<3;k++)
-  {
-
-if(usedOrNot[0] == 1)		// PN0 (Newton)	~1/c^0
-  {
-  adot_pn1[0][k] = -m2*(v[k]/r3 - 3.0*RP*x[k]/r2/r2);
-  adot_pn2[0][k] =  m1*(v[k]/r3 - 3.0*RP*x[k]/r2/r2);
-  }
- 
-if(usedOrNot[1] == 1)		// PN1		~1/c^2
-  {                  
-  adot_pn1[1][k] =  (-2.0*MOR*RP*(AK2*N[k]+BK2*v[k])/r2 + MOR*(ADK2*N[k]+BDK2*v[k])/r + MOR*(AK2*(v[k]-N[k]*RP)/r+BK2*A[k])/r)*m2/M;
-  adot_pn2[1][k] = -(-2.0*MOR*RP*(AK2*N[k]+BK2*v[k])/r2 + MOR*(ADK2*N[k]+BDK2*v[k])/r + MOR*(AK2*(v[k]-N[k]*RP)/r+BK2*A[k])/r)*m1/M;
-  }
- 
-if(usedOrNot[2] == 1)		// PN2		~1/c^4
-  {
-  adot_pn1[2][k] =  (-2.0*MOR*RP*(AK4*N[k]+BK4*v[k])/r2 + MOR*(ADK4*N[k]+BDK4*v[k])/r + MOR*(AK4*(v[k]-N[k]*RP)/r+BK4*A[k])/r)*m2/M;
-  adot_pn2[2][k] = -(-2.0*MOR*RP*(AK4*N[k]+BK4*v[k])/r2 + MOR*(ADK4*N[k]+BDK4*v[k])/r + MOR*(AK4*(v[k]-N[k]*RP)/r+BK4*A[k])/r)*m1/M;
-  }
- 
-if(usedOrNot[3] == 1)		// PN2.5	~1/c^5
-  {
-  adot_pn1[3][k] =  (-2.0*MOR*RP*(AK5*N[k]+BK5*v[k])/r2 + MOR*(ADK5*N[k]+BDK5*v[k])/r + MOR*(AK5*(v[k]-N[k]*RP)/r+BK5*A[k])/r)*m2/M;
-  adot_pn2[3][k] = -(-2.0*MOR*RP*(AK5*N[k]+BK5*v[k])/r2 + MOR*(ADK5*N[k]+BDK5*v[k])/r + MOR*(AK5*(v[k]-N[k]*RP)/r+BK5*A[k])/r)*m1/M;
-  }
- 
-if(usedOrNot[4] == 1)		// PN3  	~1/c^6
-  {
-  adot_pn1[4][k] =  (-2.0*MOR*RP*(AK6*N[k]+BK6*v[k])/r2 + MOR*(ADK6*N[k]+BDK6*v[k])/r + MOR*(AK6*(v[k]-N[k]*RP)/r+BK6*A[k])/r)*m2/M;
-  adot_pn2[4][k] = -(-2.0*MOR*RP*(AK6*N[k]+BK6*v[k])/r2 + MOR*(ADK6*N[k]+BDK6*v[k])/r + MOR*(AK6*(v[k]-N[k]*RP)/r+BK6*A[k])/r)*m1/M;
-  }
-
-if(usedOrNot[5] == 1)		// PN3.5  	~1/c^7
-  {
-  adot_pn1[5][k] =  (-2.0*MOR*RP*(AK7*N[k]+BK7*v[k])/r2 + MOR*(ADK7*N[k]+BDK7*v[k])/r + MOR*(AK7*(v[k]-N[k]*RP)/r+BK7*A[k])/r)*m2/M;
-  adot_pn2[5][k] = -(-2.0*MOR*RP*(AK7*N[k]+BK7*v[k])/r2 + MOR*(ADK7*N[k]+BDK7*v[k])/r + MOR*(AK7*(v[k]-N[k]*RP)/r+BK7*A[k])/r)*m1/M;
-  }
-
-
-if(Van_Spin == 1)		// All the SPIN terms
-  {
-  adot_pn1[6][k] +=  (C1_5D[k]/c_2 + C2D[k]/c_4 +C2_5D[k]/c_4 + QMD[k]/c_4)*m2/M;
-  adot_pn2[6][k] += -(C1_5D[k]/c_2 + C2D[k]/c_4 +C2_5D[k]/c_4 + QMD[k]/c_4)*m1/M;
-  }
-
-  
-  }
-
-// PN jerks
-
-
-
-
-// Check RS_DIST conditions !!!
-
-RS_DIST = 4.0*(2.0*m1/c_2 + 2.0*m2/c_2);
-
-if(r < RS_DIST) 
-  {
-  if(myRank == rootRank)
-    {
-    fprintf(stdout,"PN RSDIST: r  = %.8E \t RS = %.8E \n", r, RS_DIST);
-    fflush(stdout);
-    }
-  return(505);
-  }
-else 
-  {
-  return(0);
-  }
-
-
-}
-/***************************************************************************/
--- 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];
--- a/sse_sse.f
+++ b/sse_sse.f
--- a/star_destr.c
+++ b/star_destr.c
@ -1,275 +0,0 @@
-/*****************************************************************************
-  File Name      : "Star Destr.c"
-  Contents       : star "destruction" by tidal field of "live" BH (1 or 2)
-  Coded by       : Peter Berczik
-  Last redaction : 2010.IX.14 1:43PM
-*****************************************************************************/
-
-void star_destr(double time,
-                int n,
-                int ind[],
-                double m[],
-    	        double x[][3],
-		double v[][3],
-		double pot[],
-		double a[][3], 
-		double adot[][3], 
-		double t[],
-		double dt[],
-		int N,
-		double m_N[], 
-		double x_N[][3], 
-		double v_N[][3], 
-		double a_N[][3], 
-		double adot_N[][3], 
-		double t_N[], 
-		double *m_bh,
-		int *num_bh,
-                int i_bh)
-{
-
-int    n_end, k_act, N_end;
-
-double R_t, e_kin=0.0, e_pot_BH=0.0, e_pot=0.0, e_corr=0.0;
-
-double eps_bh, eps2, eps_bh2, rsb, rkb2, rks2, xp[3], v_bh[3]; 
-
-//double vir = 0.6, gamma = 1000.0;
-
-double x_max = 1.0E+03, v_max = 1.0E-08, 
-       a_max = 1.0E-08, adot_max = 1.0E-08;
-
-
-if( time < t_diss_on ) return;
-
-
-eps_bh = eps;
-
-eps2    = SQR(eps);
-eps_bh2 = SQR(eps_bh);
-
-/*
-m_s = 1.0/N;
-R_s = 2.52E-08/2.507328103;
-R_t = gamma * R_s * pow( 2.0*(*m_bh/m_s), over3 );
-*/
-
-R_t = R_TIDAL;
-
-/*
-if(myRank == rootRank)
-  {
-  printf("%.6E \t %06d %06d %06d \t %.6E \t %.6E %06d \n", time, n, i_bh, ind[i_bh], R_t, *m_bh, *num_bh);
-  fflush(stdout);
-  }
-*/
-
-#ifdef ADD_BH2
-  n_end = n-2;
-  N_end = N-2;
-#else
-#ifdef ADD_BH1
-  n_end = n-1;
-  N_end = N-1;
-#endif		// ADD_BH1
-#endif		// ADD_BH2
-
-/*
-if(myRank == rootRank)
-  {
-  printf("%.6E \t %06d %06d %06d %06d \t %.6E \t %.6E %06d \n", time, n, n_end, i_bh, ind[i_bh], R_t, *m_bh, *num_bh);
-  fflush(stdout);
-  }
-*/
-
-for(i=0; i<n_end; i++)
-  {
-
-  rsb = sqrt( SQR(x[i][0]-x[i_bh][0]) + SQR(x[i][1]-x[i_bh][1]) + SQR(x[i][2]-x[i_bh][2]) );
-  
-//  R_t = (R_s/5.848035) * pow( (2.0*(*m_bh)/m_s), over3 );
-//  if( (r < R_t) && (e_kin < ABS(vir*e_pot)) )
-
-  if( (rsb < R_t) && (m[i] != 0.0) )
-    {
-
-
-    if(myRank == rootRank)
-      {
-
-      out = fopen("accr-data.dat","a");
-
-      fprintf(out,"%04d %04d %04d %04d \t %.6E %.6E %.6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \n",
-      i, ind[i], i_bh, ind[i_bh], 
-      time,        
-      m[i],        *m_bh, 
-      x[i][0],     x[i][1],     x[i][2], 
-      x[i_bh][0],  x[i_bh][1],  x[i_bh][2], 
-      v[i][0],     v[i][1],     v[i][2],
-      v[i_bh][0],  v[i_bh][1],  v[i_bh][2],
-      pot[i], 
-      a[i][0],     a[i][1],     a[i][2],
-      adot[i][0],  adot[i][1],  adot[i][2]);
-
-      fclose(out);
-
-      } /* if(myRank == rootRank) */
-
-
-    v_bh[0] = (m[i] * v[i][0] + *m_bh * v[i_bh][0]) / (m[i] + *m_bh);
-    v_bh[1] = (m[i] * v[i][1] + *m_bh * v[i_bh][1]) / (m[i] + *m_bh);
-    v_bh[2] = (m[i] * v[i][2] + *m_bh * v[i_bh][2]) / (m[i] + *m_bh);
-
-    e_kin  =  m[i] * ( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
-    e_kin += *m_bh * ( SQR(v[i_bh][0]) + SQR(v[i_bh][1]) + SQR(v[i_bh][2]) );
-    e_kin -= (m[i] + *m_bh) * ( SQR(v_bh[0]) + SQR(v_bh[1]) + SQR(v_bh[2]) );
-    e_kin *= 0.5;
-
-//    e_kin = 0.5*m[i] * ( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) - SQR(v[i_bh][0]) - SQR(v[i_bh][1]) - SQR(v[i_bh][2]) );
-
-    e_pot_BH = - *m_bh * m[i] / sqrt( SQR(rsb) + eps_bh2 );
-
-
-    tmp = 0.0;
-
-    for(k=0; k<N_end; k++)
-      {
-    if( k != ind[i] )
-      {
-
-      /* define if k prinadlezhit ACT ili net */
-
-      k_act = 0;
-      for(ii=0; ii<n_end; ii++) 
-        {
-        if( k == ind[ii] ) 
-          {
-          k_act = 1; break;
-          } 
-        } /* ii */
-
-/*
-    if(myRank == rootRank)
-      {
-      printf("0: \t %04d %04d \t %04d %04d \n", i, ind[i], k, k_act);
-      fflush(stdout);
-      }
-*/
-
-      if(k_act == 1)	/* k prinadlezhit ACT */
-        {
-
-        rkb2 = SQR(x_N[k][0]-x[i_bh][0]) + SQR(x_N[k][1]-x[i_bh][1]) + SQR(x_N[k][2]-x[i_bh][2]);
-        tmp -= m_N[k]/sqrt( rkb2 + eps_bh2 );
-
-        rks2 = SQR(x_N[k][0]-x[i][0]) + SQR(x_N[k][1]-x[i][1]) + SQR(x_N[k][2]-x[i][2]);
-        tmp += m_N[k]/sqrt( rks2 + eps2 );
-
-        }
-      else 		/* k ne prinadlezhit ACT */
-        {
-
-        dt_tmp = time - t_N[k];
-        dt2half = over2*dt_tmp*dt_tmp;
-        dt3over6 = over3*dt_tmp*dt2half;
-
-        xp[0] = x_N[k][0] + v_N[k][0]*dt_tmp + a_N[k][0]*dt2half + adot_N[k][0]*dt3over6;
-        xp[1] = x_N[k][1] + v_N[k][1]*dt_tmp + a_N[k][1]*dt2half + adot_N[k][1]*dt3over6;
-        xp[2] = x_N[k][2] + v_N[k][2]*dt_tmp + a_N[k][2]*dt2half + adot_N[k][2]*dt3over6;
-
-        rkb2 = SQR(xp[0]-x[i_bh][0]) + SQR(xp[1]-x[i_bh][1]) + SQR(xp[2]-x[i_bh][2]);
-        tmp -= m_N[k]/sqrt( rkb2 + eps_bh2 );
-
-        rks2 = SQR(xp[0]-x[i][0]) + SQR(xp[1]-x[i][1]) + SQR(xp[2]-x[i][2]);
-        tmp += m_N[k]/sqrt( rks2 + eps2 );
-
-        } /* if(k_act == 1) */
-
-/*
-    if(myRank == rootRank)
-      {
-      printf("\t %04d %04d \t %04d \t % .8E \n", k, ind[i], k_act, tmp);
-      fflush(stdout);
-      }
-*/
-
-      } /* if( k != ind[i] ) */
-
-      } /* k */
-
-    e_pot = - m[i] * tmp; 
-
-    e_corr = e_kin + e_pot_BH + e_pot;
-
-    E_corr += e_corr;
-
-    *m_bh += m[i];		// add the star mass to the BH mass
-
-    *num_bh = *num_bh + 1;
-
-
-    if(myRank == rootRank)
-      {
-      printf("ACCR: %.6E %04d %06d   %.4E %.4E %04d %06d   %.4E % .4E % .4E % .4E % .4E   %.4E %04d \n", 
-            time, i_bh, ind[i_bh], rsb, R_t, i, ind[i], 
-	    e_kin, e_pot_BH, e_pot, e_corr, E_corr, *m_bh, *num_bh);
-      fflush(stdout);
-      }
-
-
-    v[i_bh][0] = v_bh[0];
-    v[i_bh][1] = v_bh[1];
-    v[i_bh][2] = v_bh[2];
-
-
-//    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]) );
-
-
-    /* for part "i" */
-
-    x[i][0] = x_max + 1.0*my_rand2();
-    x[i][1] = x_max + 1.0*my_rand2();
-    x[i][2] = x_max + 1.0*my_rand2();
-
-    v[i][0] = v_max*my_rand2();
-    v[i][1] = v_max*my_rand2();
-    v[i][2] = v_max*my_rand2();
-
-    pot[i] = 0.0;
-
-    a[i][0] = a_max*my_rand2();
-    a[i][1] = a_max*my_rand2();
-    a[i][2] = a_max*my_rand2();
-
-    adot[i][0] = adot_max*my_rand2();
-    adot[i][1] = adot_max*my_rand2();
-    adot[i][2] = adot_max*my_rand2();
-
-    m[i] = 0.0;
-
-    t[i] = t[i] + 0.125;
-
-    dt[i] = 0.125;
-
-    } /* if( (r < R_t) ) */
-
-  } /* i */
-
-
-  m[i_bh] = *m_bh;
-
-/*
-  if(myRank == rootRank)
-    {
-    printf("%06d \t %.6E \n", i_bh, m[i_bh]);
-    fflush(stdout);
-    }
-*/
-
-}
-
-/***************************************************************/
-/***************************************************************/
-/***************************************************************/
--- a/star_destr_ext.c
+++ b/star_destr_ext.c
@ -1,286 +0,0 @@
-/*****************************************************************************
-  File Name      : "Star Destr EXT.c"
-  Contents       : star "destruction" by tidal field of EXT fixed BH
-  Coded by       : Peter Berczik
-  Last redaction : 2010.IX.14 1:43PM
-*****************************************************************************/
-
-void star_destr_ext(double time,
-                int n,
-                int ind[],
-                double m[],
-    	        double x[][3],
-		double v[][3],
-		double pot[],
-		double a[][3], 
-		double adot[][3], 
-		double t[],
-		double dt[],
-		int N,
-		double m_N[], 
-		double x_N[][3], 
-		double v_N[][3], 
-		double a_N[][3], 
-		double adot_N[][3], 
-		double t_N[], 
-		double *m_bh,
-		int *num_bh)
-{
-
-int    n_end, k_act, N_end;
-
-double R_t, e_kin=0.0, e_pot_BH=0.0, e_pot=0.0, e_corr=0.0;
-
-double eps_bh, eps_bh2, eps2, rsb, rkb2, rks2, xp[3]; 
-
-//double vir = 0.6, gamma = 1000.0;
-
-double x_max = 1.0E+03, v_max = 1.0E-08, 
-       a_max = 1.0E-08, adot_max = 1.0E-08;
-
-double a_ext_i[3], adot_ext_i[3];
-
-
-if(time < t_diss_on) return;
-
-
-eps2 = SQR(eps);
-
-eps_bh = b_bh;
-eps_bh2 = SQR(eps_bh);
-
-
-/*
-m_s = 1.0/N;
-
-R_s = 2.52E-08/2.507328103;
-R_t = gamma * R_s * pow( 2.0*(*m_bh/m_s), over3 );
-*/
-
-R_t = R_TIDAL;
-
-
-/*
-if(myRank == rootRank)
-  {
-  printf("%.6E \t %06d %06d %06d \t %.6E \t %.6E %06d \n", time, n, i_bh, ind[i_bh], R_t, *m_bh, *num_bh);
-  fflush(stdout);
-  }
-*/
-
-
-  n_end = n;
-  N_end = N;
-
-/*
-if(myRank == rootRank)
-  {
-  printf("%.6E \t %06d %06d %06d %06d \t %.6E \t %.6E %06d \n", time, n, n_end, i_bh, ind[i_bh], R_t, *m_bh, *num_bh);
-  fflush(stdout);
-  }
-*/
-
-
-for(i=0; i<n_end; i++)
-  {
-
-  rsb = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) );
-
-//  R_t = (R_s/5.848035) * pow( (2.0*(*m_bh)/m_s), over3 );
-//  if( (r < R_t) && (e_kin < ABS(vir*e_pot)) )
-
-  if( (rsb < R_t) )
-    {
-
-
-    if(myRank == rootRank)
-      {
-
-      x_ij = x[i][0];
-      y_ij = x[i][1];
-      z_ij = x[i][2];
-
-      vx_ij = v[i][0];
-      vy_ij = v[i][1];
-      vz_ij = v[i][2];
-
-      r2 = SQR(x_ij) + SQR(y_ij) + SQR(z_ij) + eps_bh2; 
-      r = sqrt(r2);
-
-      rv_ij = vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij;
-
-      tmp = *m_bh / (r * r2);
-
-      a_ext_i[0] = -tmp * x_ij;
-      a_ext_i[1] = -tmp * y_ij;
-      a_ext_i[2] = -tmp * z_ij;
-
-      adot_ext_i[0] = -tmp * (vx_ij - 3.0*rv_ij * x_ij/r2);
-      adot_ext_i[1] = -tmp * (vy_ij - 3.0*rv_ij * y_ij/r2);
-      adot_ext_i[2] = -tmp * (vz_ij - 3.0*rv_ij * z_ij/r2);
-
-      out = fopen("accr-data.dat","a");
-
-      ii = ind[i];
-
-#ifndef STARDISK
-      a_drag[ii][0] = 0.0;
-      a_drag[ii][1] = 0.0;
-      a_drag[ii][2] = 0.0;
-
-      adot_drag[ii][0] = 0.0;
-      adot_drag[ii][1] = 0.0;
-      adot_drag[ii][2] = 0.0;
-#endif // no STARDISK
-
-      fprintf(out, "%04d %04d \t %.6E %.6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \n",
-      i, ii, 
-      time_cur, *m_bh, 
-      x[i][0],     x[i][1],     x[i][2], 
-      v[i][0],     v[i][1],     v[i][2],
-      a[i][0],     a[i][1],     a[i][2],
-      adot[i][0],  adot[i][1],  adot[i][2],
-      a_ext_i[0],       a_ext_i[1],       a_ext_i[2],
-      adot_ext_i[0],    adot_ext_i[1],    adot_ext_i[2],
-      a_drag[ii][0],    a_drag[ii][1],    a_drag[ii][2],
-      adot_drag[ii][0], adot_drag[ii][1], adot_drag[ii][2]);
-
-      fclose(out);
-
-      } /* if(myRank == rootRank) */
-
-
-
-    e_kin = 0.5*m[i] * ( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
-
-    e_pot_BH = - *m_bh * m[i] / sqrt( SQR(rsb) + eps_bh2 );
-
-
-    tmp = 0.0;
-
-    for(k=0; k<N_end; k++)
-      {
-    if( k != ind[i] )
-      {
-
-      /* define if k prinadlezhit ACT ili net */
-
-      k_act = 0;
-      for(ii=0; ii<n_end; ii++) 
-        {
-        if(k==ind[ii]) 
-          {
-          k_act=1; break;
-          } 
-        } /* ii */
-
-/*
-    if(myRank == rootRank)
-      {
-      printf("0: \t %04d %04d \t %04d %04d \n", i, ind[i], k, k_act);
-      fflush(stdout);
-      }
-*/
-
-      if(k_act == 1)	/* k prinadlezhit ACT */
-        {
-
-        rkb2 = SQR(x_N[k][0]) + SQR(x_N[k][1]) + SQR(x_N[k][2]);
-        tmp -= m_N[k]/sqrt( rkb2 + eps_bh2 );
-
-        rks2 = SQR(x_N[k][0]-x[i][0]) + SQR(x_N[k][1]-x[i][1]) + SQR(x_N[k][2]-x[i][2]);
-        tmp += m_N[k]/sqrt( rks2 + eps2 );
-
-        }
-      else 		/* k ne prinadlezhit ACT */
-        {
-
-        dt_tmp = time - t_N[k];
-        dt2half = over2*dt_tmp*dt_tmp;
-        dt3over6 = over3*dt_tmp*dt2half;
-
-        xp[0] = x_N[k][0] + v_N[k][0]*dt_tmp + a_N[k][0]*dt2half + adot_N[k][0]*dt3over6;
-        xp[1] = x_N[k][1] + v_N[k][1]*dt_tmp + a_N[k][1]*dt2half + adot_N[k][1]*dt3over6;
-        xp[2] = x_N[k][2] + v_N[k][2]*dt_tmp + a_N[k][2]*dt2half + adot_N[k][2]*dt3over6;
-
-        rkb2 = SQR(xp[0]) + SQR(xp[1]) + SQR(xp[2]);
-        tmp -= m_N[k]/sqrt( rkb2 + eps_bh2 );
-
-        rks2 = SQR(xp[0]-x[i][0]) + SQR(xp[1]-x[i][1]) + SQR(xp[2]-x[i][2]);
-        tmp += m_N[k]/sqrt( rks2 + eps2 );
-
-        } /* if(k_act == 1) */
-
-/*
-    if(myRank == rootRank)
-      {
-      printf("\t %04d %04d \t %04d \t % .8E \n", k, ind[i], k_act, tmp);
-      fflush(stdout);
-      }
-*/
-
-      } /* if( k != ind[i] ) */
-
-      } /* k */
-
-
-    e_pot = - m[i] * tmp; 
-
-    e_corr = e_kin + e_pot_BH + e_pot;
-
-    E_corr += e_corr;
-
-    *m_bh += m[i];		// add the star mass to the BH mass
-
-    *num_bh = *num_bh + 1;
-
-    
-    if(myRank == rootRank)
-      {
-      printf("ACCR: %.4E %.4E %.4E %04d %06d %.4E % .4E % .4E % .4E % .4E %.4E %04d \n", 
-            time, rsb, R_t, i, ind[i], e_kin, e_pot_BH, e_pot, e_corr, E_corr, *m_bh, *num_bh);
-      fflush(stdout);
-      }
-
-
-//    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]) );
-
-
-    /* for part "i" */
-
-    x[i][0] = x_max + 1.0*my_rand2();
-    x[i][1] = x_max + 1.0*my_rand2();
-    x[i][2] = x_max + 1.0*my_rand2();
-
-    v[i][0] = v_max*my_rand2();
-    v[i][1] = v_max*my_rand2();
-    v[i][2] = v_max*my_rand2();
-
-    pot[i] = 0.0;
-
-    a[i][0] = a_max*my_rand2();
-    a[i][1] = a_max*my_rand2();
-    a[i][2] = a_max*my_rand2();
-
-    adot[i][0] = adot_max*my_rand2();
-    adot[i][1] = adot_max*my_rand2();
-    adot[i][2] = adot_max*my_rand2();
-
-    m[i] = 0.0;
-
-    t[i] = t[i] + 0.125;
-
-    dt[i] = 0.125;
-
-    } /* if( (r < R_t) ) */
-
-  } /* i */
-
-
-}
-
-/***************************************************************/
-/***************************************************************/
-/***************************************************************/
Author	SHA1	Message	Date
Yohai Meiron	6351f8075b	Added #include <stdexcept>	2022-05-30 10:28:12 -04:00
Yohai Meiron	633c82b917	Enabled triple (plus) SMBHs	2022-05-19 20:35:24 -04:00
Yohai Meiron	fcdafd94a6	Fixed PN force adjustment before PN parameters set	2022-01-14 17:54:38 -05:00
Yohai Meiron	1fe119fdb2	Fixed BH softening correction when grapite_smbh_star_eps is defined	2021-12-26 21:40:50 -05:00
Yohai Meiron	8b71f4fe92	Cleaned up corrector loop	2021-12-26 17:20:20 -05:00
Yohai Meiron	15227a9eee	exposed grapite_set_eps_bh in conf	2021-08-30 09:31:53 -04:00
Yohai Meiron	6048c7748e	special softening set before first acc adjustment	2021-08-27 19:08:57 -04:00
Yohai Meiron	064eb4e3c5	Removed init. of pot_act_ext from main loop; zero filling of pot_act_ext now only if external potential active; improved warning on minimum timestep; changed CUDAHOME to CUDA_HOME	2021-08-23 19:40:55 -04:00
Yohai Meiron	0a2344563d	run directory now separate from source directory	2021-03-17 20:37:12 -04:00
Yohai Meiron	11c0db23a3	Compatibility with Sapporo	2020-11-03 21:41:21 -05:00
Yohai Meiron	dd31b63215	Improved printing of external potential	2020-10-31 15:08:48 -04:00
Yohai Meiron	8adb4ac813	Fixed memory issue with reading of HDF5 input file and improved style	2020-10-28 22:46:27 -04:00
Yohai Meiron	b6c55a30da	Made ETICS blocks work with the vector structures	2020-10-28 21:40:46 -04:00
Yohai Meiron	54231b8537	Removed remaining macros and unused includes	2020-10-28 21:15:51 -04:00
Yohai Meiron	d9ed8e5131	Got rid of remaining manual allocations, moved some variable declarations to where they are needed	2020-10-28 19:18:37 -04:00
Yohai Meiron	aef0f26878	Improved external gravity calls style	2020-10-28 16:31:45 -04:00
Yohai Meiron	747f9f9d89	Moved config to inside main and changed several arrays into vectors	2020-10-27 23:55:55 -04:00
Yohai Meiron	329dd2ca4d	New timing, moved MPI reduction into function call, started converting some of the pointers into std::vectors	2020-05-09 18:21:01 -04:00
Yohai Meiron	1a438449a8	Now fast active search works, and can be controlled from config file	2020-05-01 20:41:50 -04:00
Yohai Meiron	2f8f8c582c	Attempting to get fast active search working	2020-05-01 00:27:41 -04:00
Yohai Meiron	d046c189b3	Cleaned up energy control as well as many other things	2020-04-29 21:51:05 -04:00
Yohai Meiron	df5d89a0c8	Finished moving the BH stuff into their compilation unit	2020-04-19 22:10:19 -04:00
Yohai Meiron	19ce85da88	Moved out almost all BH-related stuff to a different compilation unit	2020-04-19 14:16:03 -04:00
Yohai Meiron	2a50f0fc9a	Lots of cleanup and started moving the PN out of the main	2020-04-18 21:48:46 -04:00
Yohai Meiron	30ae8631a9	Moved calc_self_grav to a class and cleaned up a little	2020-04-14 20:20:39 -04:00
Yohai Meiron	b51613695f	Now working as the hybrid code Removed the need for the annoying ETICS_CEP flag; now dt_scf is read from the config file.	2020-04-13 23:08:42 -04:00
Yohai Meiron	75ad0f0e89	Added spins to config file and got rid of more _act arrays	2020-04-08 20:38:36 -04:00
Yohai Meiron	ebec3280f8	Cleaned unused variables and fixed a bug with PN (also tested PN against old version)	2020-04-07 22:40:54 -04:00
Yohai Meiron	d9e3aea243	Moved SMBH de-softening to a function (still not pretty) and improved init script	2020-04-05 19:54:10 -04:00
Yohai Meiron	56abe820c3	Number of devices can be overriden in config file	2020-04-05 11:17:21 -04:00
Yohai Meiron	7fae35a2d9	Moved some tasks from the correction+timestep adjustment loop outside	2020-04-04 23:03:06 -04:00
Yohai Meiron	99ec3064ed	Started moving toward dynamical allocation; removed the _loc variables which were only used in the zeroth step and were pretty unnecessary anyway	2020-04-03 21:40:11 -04:00
Yohai Meiron	ea94dbb626	Moved the active search outside of the main loop	2020-04-02 21:49:52 -04:00
Yohai Meiron	c79cef895a	Added Dehnen potential and moved external potentials to a new compilation unit	2020-04-01 15:47:13 -04:00
Yohai Meiron	62b0d7e491	Improved output and config	2020-03-31 22:06:04 -04:00
Yohai Meiron	5cb4282be4	Completed HDF5 support	2020-03-31 19:04:18 -04:00
Yohai Meiron	b2943be7c1	Improved HDF5 reader and added ASCII writer with custom precision	2020-03-28 01:34:37 -04:00
Yohai Meiron	ba5546534f	Continued working on HDF5 plugin	2020-03-27 16:53:01 -04:00
Yohai Meiron	0b9fa6e46d	Working on HDF5 reader	2020-03-26 19:56:12 -04:00
Yohai Meiron	f065566c6b	Added an HDF5 writer	2020-03-26 15:09:15 -04:00
Yohai Meiron	131642aa71	Added logarithmic external potential	2020-03-25 21:42:29 -04:00
Yohai Meiron	8cba5b5e1d	Got rid of more global variables	2020-03-25 16:58:46 -04:00
Yohai Meiron	23d32c6dcf	New approach to external potential	2020-03-25 16:37:55 -04:00
Yohai Meiron	562ea9618a	External potential (EXTPOT_GAL) cleaned up	2020-03-24 20:37:23 -04:00
Yohai Meiron	d7332be188	Added post-Newtonian parameters to the config file	2020-03-21 22:52:13 -04:00
Yohai Meiron	9b900eff6e	Prettied up the black hole output and removed some redundant variables	2020-03-21 12:46:26 -04:00
Yohai Meiron	8f984b5158	Moved more variables from global to main and removed some unused variables	2020-03-19 14:01:17 -04:00
Yohai Meiron	c4830423a3	Fixed the control output problem	2020-03-18 22:50:00 -04:00
Yohai Meiron	22f983cf17	Moved some variables from global to main and removed some unused variables	2020-03-18 11:59:28 -04:00
Yohai Meiron	9296db0609	Now config file handles most parameters (but not external potential)	2020-03-17 14:09:30 -04:00
Yohai Meiron	bb9a343763	Improved parameter parsing from config file	2020-03-17 10:43:58 -04:00
Yohai Meiron	cd282cc406	Improved the new-style config file and its reader	2020-03-15 21:35:53 -04:00
Yohai Meiron	90f070d1fb	Moved much of the black hole stuff from ifdef blocks to regular if clauses	2020-03-15 18:43:13 -04:00
Yohai Meiron	7721bc1830	Minor edits	2020-03-14 16:50:40 -04:00
Yohai Meiron	fd33819a7a	Got rid of most k-loops	2020-03-14 16:29:00 -04:00
Yohai Meiron	7177148b6c	Fixed typo in Makefile	2020-03-13 20:59:53 -04:00
Yohai Meiron	34030c06d0	Introduced double3 structure in preparation of dynamic allocation and getting rid of k-loops	2020-03-13 20:42:09 -04:00
Yohai Meiron	4b8ef0ca61	Removed redundant copying from calc_self_grav and got rid of calc_self_grav_zero	2020-03-13 19:36:18 -04:00
Yohai Meiron	2c05355eaa	Removed more ifdef blocks	2020-03-13 16:59:09 -04:00
Yohai Meiron	510b0ff02d	started working on simple config file parser	2020-03-08 00:59:44 -05:00
Yohai Meiron	63a39cc9c8	Removed ifdef blocks specific to non-Yebisu libraries	2020-03-06 15:18:32 -05:00
Yohai Meiron	b5071b792e	Added TODO	2020-03-06 15:04:55 -05:00
Yohai Meiron	86240127fd	mostly fixed indentation in main	2020-03-05 23:48:08 -05:00
Yohai Meiron	dd77ecde75	removed many ifdefs blocks and associated included c files	2020-03-05 17:00:13 -05:00
Yohai Meiron	e9455c0a0e	Started cleanup process	2020-03-04 13:46:04 -05:00