diff --git a/init.py b/init.py index e5aa76c..649f75f 100644 --- 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' @@ -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) diff --git a/phi-GRAPE.c b/phi-GRAPE.c index 236ee7e..352336b 100644 --- a/phi-GRAPE.c +++ b/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