From 08d1c155f551bbd5a904a2a4a3923e45f3eba8be Mon Sep 17 00:00:00 2001
From: Yohai Meiron <sahmes@gmail.com>
Date: Wed, 4 Mar 2020 12:47:27 -0500
Subject: [PATCH] ACT_DEF_GRAPITE should now work with SMBHs

---
 phi-GRAPE.c | 67 +++++++++++++++++++++++++++++++++--------------------
 1 file changed, 42 insertions(+), 25 deletions(-)

diff --git a/phi-GRAPE.c b/phi-GRAPE.c
index df15d03..236ee7e 100644
--- a/phi-GRAPE.c
+++ b/phi-GRAPE.c
@@ -67,8 +67,8 @@
 
 //#define ADD_BH1		// add the Single BH
 
-// #define ADD_BH2			// add the Binary BH's
-// #define ADD_N_BH		// eps_BH = 0.0, but added only the Newtonian forces
+#define ADD_BH2			// add the Binary BH's
+#define ADD_N_BH		// eps_BH = 0.0, but added only the Newtonian forces
 // #define ADD_PN_BH		// extra - added also the Post-Newton forces
 // #define ADD_SPIN_BH		// extra - added the SPIN for the BH's - DEFAULT !!!
 
@@ -5701,6 +5701,35 @@ for(i=0; i<N; i++)
       } /* i */
 #endif      // ACT_DEF_GRAPITE
 
+#if defined(ACT_DEF_GRAPITE) && (defined(ADD_BH1) || defined(ADD_BH2))
+#ifdef ADD_BH1
+#define ACT_DEF_GRAPITE_NUMBH 1
+#else
+#define ACT_DEF_GRAPITE_NUMBH 2
+#endif
+    int acf_def_grapite_bh_count = 0;
+    for (i=0; i<n_act; i++) {
+      if (ind_act[i]==0) {
+        i_bh1 = i;
+        acf_def_grapite_bh_count++;
+      }
+#ifdef ADD_BH2
+      else if (ind_act[i]==1) {
+        i_bh2 = i;
+        acf_def_grapite_bh_count++;
+      }
+#endif
+      if (acf_def_grapite_bh_count==ACT_DEF_GRAPITE_NUMBH) break;
+    }
+    if (i==n_act) {
+      fprintf(stderr, "ERROR: black holes were not found in the active particle list");
+      exit(1);
+    }
+#elif defined(ADD_BH1) || defined(ADD_BH2)
+    i_bh1 = 0;
+    i_bh2 = 1;
+#endif
+
 #ifdef TIMING
   get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
   DT_ACT_DEF2 += (CPU_tmp_user - CPU_tmp_user0);
@@ -5854,9 +5883,6 @@ for(i=0; i<N; i++)
 
 #ifdef ADD_BH2
 
-  i_bh1 = 0;
-  i_bh2 = 1;
-
 #ifdef ADD_N_BH
 
   m_bh1 = m_act[i_bh1];
@@ -6194,8 +6220,8 @@ for(i=0; i<N; i++)
       if( dt_act[i] < min_dt ) min_dt = dt_act[i];
       } /* i */
 
-    dt_act[0] = min_dt;
-    dt_act[1] = min_dt;
+    dt_act[i_bh1] = min_dt;
+    dt_act[i_bh2] = min_dt;
 
 #else
 
@@ -6208,7 +6234,7 @@ for(i=0; i<N; i++)
       if( dt_act[i] < min_dt ) min_dt = dt_act[i];
       } /* i */
 
-    dt_act[0] = min_dt;
+    dt_act[i_bh1] = min_dt;
 
 #endif		// ADD_BH1
 
@@ -6271,9 +6297,6 @@ for(i=0; i<N; i++)
 
     out = fopen("bbh.inf","a");
 
-    i_bh1 = 0;
-    i_bh2 = 1;
-
     m_bh1 = m_act[i_bh1];
     m_bh2 = m_act[i_bh2];
 
@@ -6452,30 +6475,24 @@ for(i=0; i<N; i++)
 
 #ifdef ADD_BH2
 
-  i_bh = 0;
+  star_destr(min_t, n_act, ind_act, m_act, x_act, v_act, pot_act, a_act, adot_act, t_act, dt_act,
+                    N, ind, m, x, v, pot, a, adot, t, &m_bh1, &num_bh1, i_bh1);
+
+  m_act[i_bh1] = m_bh1;
 
   star_destr(min_t, n_act, ind_act, m_act, x_act, v_act, pot_act, a_act, adot_act, t_act, dt_act,
-                    N, ind, m, x, v, pot, a, adot, t, &m_bh1, &num_bh1, i_bh);
+                    N, ind, m, x, v, pot, a, adot, t, &m_bh2, &num_bh2, i_bh2);
 
-  m_act[i_bh] = m_bh1;
-
-  i_bh = 1;
-
-  star_destr(min_t, n_act, ind_act, m_act, x_act, v_act, pot_act, a_act, adot_act, t_act, dt_act,
-                    N, ind, m, x, v, pot, a, adot, t, &m_bh2, &num_bh2, i_bh);
-
-  m_act[i_bh] = m_bh2;
+  m_act[i_bh2] = m_bh2;
 
 #else
 
 #ifdef ADD_BH1
 
-  i_bh = 0;
-
   star_destr(min_t, n_act, ind_act, m_act, x_act, v_act, pot_act, a_act, adot_act, t_act, dt_act,
-                    N, ind, m, x, v, pot, a, adot, t, &m_bh1, &num_bh1, i_bh);
+                    N, ind, m, x, v, pot, a, adot, t, &m_bh1, &num_bh1, i_bh1);
 
-  m_act[i_bh] = m_bh1;
+  m_act[i_bh1] = m_bh1;
 
 
 #endif		// ADD_BH1