Fixed bug with the bbh influence output when grapite active search is enabled

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'
@@ -86,6 +118,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 +130,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)

phi-GRAPE.c

@@ -356,7 +356,7 @@ double m_act[N_MAX],
        pot_act_tmp[N_MAX], a_act_tmp[N_MAX][3], adot_act_tmp[N_MAX][3],
        pot_act_tmp_loc[N_MAX], a_act_tmp_loc[N_MAX][3], adot_act_tmp_loc[N_MAX][3];
 
-FILE   *inp, *out, *tmp_file, *dbg;
+FILE   *inp, *out, *tmp_file, *dbg, *out_bbhinf;
 
 double CPU_time_real0, CPU_time_user0, CPU_time_syst0;
 double CPU_time_real,  CPU_time_user,  CPU_time_syst;
@@ -4528,8 +4528,7 @@ int main(int argc, char *argv[])
 
 
 #ifdef BBH_INF
-  out = fopen("bbh.inf","w");
-  fclose(out);
+  out_bbhinf = fopen("bbh.inf","w");
 #endif
 
 
@@ -6295,8 +6294,6 @@ for(i=0; i<N; i++)
   if(myRank == rootRank)
     {
 
-    out = fopen("bbh.inf","a");
-
     m_bh1 = m_act[i_bh1];
     m_bh2 = m_act[i_bh2];
 
@@ -6335,12 +6332,13 @@ for(i=0; i<N; i++)
 
 
 
-    for(i=2; i<n_act; i++)
+    for(i=0; i<n_act; i++)
       {
 
-      tmp_r2 = SQR(x_act[i][0] - x_bbhc[0]) + SQR(x_act[i][1] - x_bbhc[1]) + SQR(x_act[i][2] - x_bbhc[2]);
-
       iii = ind_act[i];
+      if (iii < 2) continue;
+
+      tmp_r2 = SQR(x_act[i][0] - x_bbhc[0]) + SQR(x_act[i][1] - x_bbhc[1]) + SQR(x_act[i][2] - x_bbhc[2]);
 
 //      if( (tmp_r2 < DR2*R_INF2) )
       if( (tmp_r2 < SEMI_a2*R_INF2) )
@@ -6361,8 +6359,8 @@ for(i=0; i<N; i++)
         fprintf(out,"INF1 %.6E %.16E %07d %07d   %.6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E \n", 
                 Timesteps, time_cur, i, ind_act[i], 
                 sqrt(DR2), x_bbhc[0], x_bbhc[1], x_bbhc[2], v_bbhc[0], v_bbhc[1], v_bbhc[2], 
-                m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_act[0], 
-                m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_act[1], 
+                m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_act[i_bh1], 
+                m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_act[i_bh2], 
                 sqrt(tmp_r2), 
                 m_act[i], x_act[i][0], x_act[i][1], x_act[i][2], v_act[i][0], v_act[i][1], v_act[i][2], pot_act[i], 
                 dt_act[i]);
@@ -6389,8 +6387,8 @@ for(i=0; i<N; i++)
         fprintf(out,"INF2 %.6E %.16E %07d %07d   %.6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E   %.6E   %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E  %.6E \n", 
                 Timesteps, time_cur, i, ind_act[i], 
                 sqrt(DR2), x_bbhc[0], x_bbhc[1], x_bbhc[2], v_bbhc[0], v_bbhc[1], v_bbhc[2], 
-                m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_act[0], 
-                m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_act[1], 
+                m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_act[i_bh1], 
+                m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_act[i_bh2], 
                 sqrt(tmp_r2), 
                 m_act[i], x_act[i][0], x_act[i][1], x_act[i][2], v_act[i][0], v_act[i][1], v_act[i][2], pot_act[i], 
                 dt_act[i]);
@@ -6402,8 +6400,6 @@ for(i=0; i<N; i++)
 
       } /* i */
 
-    fclose(out);
-
     } /* if(myRank == rootRank) */
 #endif
 
@@ -7273,7 +7269,9 @@ for(i=0; i<N; i++)
     fclose(tmp_file);
 #endif
 
-
+#ifdef BBH_INF
+  fclose(out_bbhinf);
+#endif
   } /* if(myRank == rootRank) */