-
Notifications
You must be signed in to change notification settings - Fork 0
/
BH_merger.py
278 lines (193 loc) · 11.1 KB
/
BH_merger.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
import subprocess
import numpy as np
from tqdm import tqdm
import matplotlib.pyplot as plt
import os
import argparse
from BH_BlackBox import *
import time
import shutil
script_description = "This script simulates n number of black hole mergers and stores the physical parameters for each generation of the hierarchical merger "\
"in different star cluster environments: Open/Young cluster, Globular cluster, Nuclear cluster. "
print('\n\n--------------------------------------------------\n'\
'Make sure that the directory in which the script is located contains the "bh_kick_code" directory.\n'\
'--------------------------------------------------\n')
temp1 = input('Press Enter to continue...\nor to exit, type "exit"')
if temp1.lower() == 'exit':
exit()
parser = argparse.ArgumentParser(description=script_description)
# Arguments for the directories
parser.add_argument('-path', type=str, default=os.getcwd(), help='Path to the directory containing the gwkik2 code (Default: Current working directory)')
parser.add_argument('-output', type=str, default=os.getcwd(), help='Path to the output directory (Default: Current working directory)')
parser.add_argument('-n_sim', type=int, help='Number of simulations/iterations to run')
parser.add_argument('-n_max_gen', type=int, default=4, help='Maximum black hole generation in a single simulation (Default: 4)')
# Argument for the system
parser.add_argument('-cluster_env', type=str, help='Star cluster environment: [Open (or Young), Globular, Nuclear]')
args = parser.parse_args()
path = args.path
output = args.output
n_sim = args.n_sim
n_max_gen = args.n_max_gen
cluster_env = args.cluster_env.lower()
# Check if the path exists
if os.path.exists(path):
os.chdir(path)
else:
print('The path does not exist. Please enter a valid path.')
exit()
# Check if the output directory exists
if not os.path.exists(output):
print('The output directory does not exist. Creating the directory...')
os.makedirs(output)
# Check if the number of max generations is valid
if n_max_gen < 1:
print('The number of maximum generations should be greater than 1.')
exit()
elif n_max_gen > 10:
print('The number of maximum generations should be less than 10.')
exit()
# Check if the number of given simulations is valid
if n_sim < 1:
print('The number of maximum generations should be greater than 1.')
exit()
# Check if the inputs are proper
if cluster_env.lower() in ['open', 'young', 'globular', 'nuclear']:
pass
else:
print('Invalid cluster environment. Please enter a valid cluster environment: [Open (or Young), Globular, Nuclear]')
exit()
# Get the operating system
op_sys = get_OS()
print('\n--------------------------------------------------\n'\
'Python script initiated\n'\
'--------------------------------------------------\n')
print("\nThis python script simulates hierarchical merging in different star cluster environements. The workflow of this script is as follows:\n"\
"1. Choose random black hole mass and its corresponding kerr parameter from 'bhlist' file.\n"\
"2. Calculate the post merger parameters using some equations written in the fortran code (gwkik2.f)\n"\
"3. Randomly choose another black hole from the 'bhlist' file and merger it with the remnant black hole from the previous merger.\n"\
"4. Continue this until either of two conditions are met:\n"\
"\t-Kick velocity of the remnant black hole exceeds the escape velocity of the cluster environment.\n"\
"\t-The hierarchical merging reaches the maximum generation given by the user.\n"\
"5. This continues for n number of simulations.\n"\
"6. Plots are produced to:\n"\
"\t-Check the inherent probability of getting upto a specific generation of black hole in different cluster environments.\n"\
"\t-Weight this probability by the observed volume fraction of the merger events from diffferent GW observatories.\n"\
"\t-See how parameters vary with increasing black hole generation.\n"\
"\t-Know the correlations between black hole parameters.\n"\
)
input("Press Enter to continue...\n")
print("\nIt uses grkick code, which is code written in Fortran to simulate the black hole mergers. \n"\
"The references to the equations used in the grkick code is given in the file 'gwkick2.f'.\n\n"\
"--------------------------------------"\
"\nNOTE: Before running the script, make sure that the gwkik2 code is compiled and the executable is in the same directory as this script.\n"\
"--------------------------------------\n\n")
print('******************** Program initiates ********************')
print(f'Operating system: {op_sys}')
print(f'Number of simulations: {n_sim}')
print(f'Number of maximum mergers/generations in one simulation: {n_max_gen}')
print(f'Star cluster environment: {cluster_env}\n')
input('Press Enter to continue...')
# Import the data from the bhlist file
bh_mass, bh_kerr = get_bh_data()
print('\nBlack hole data imported successfully!')
print('\nSelecting random black hole masses and Kerr parameters...')
# Get random black hole masses and Kerr parameters
sample_bh_mass, sample_kerr_param = select_random_bh(bh_mass, bh_kerr, num_samples = 2)
print('\nRandom selection of Black Hole parameters successful!')
##### Define the remaining input parameters #####
theta1, theta2 = 90.0, 0.0
phi1, phi2 = 0.0, 0.0
initial_params = {
'm1': sample_bh_mass[0], # Mass of Black Hole 1
'm2': sample_bh_mass[1], # Mass of Black Hole 2
's1': -20.0, # Spin of Black Hole 1 (with dimension)
'theta1': theta1, # Polar angle of Black Hole 1
'phi1': phi1, # Azimuthal angle of Black Hole 1
's2': -20.0, # Spin of Black Hole 2 (with dimension)
'theta2': theta2, # Polar angle of Black Hole 2
'phi2': phi2, # Azimuthal angle of Black Hole 2
'a1': sample_kerr_param[0], # Kerr parameter of Black Hole 1 (dimensionless spin)
'a2': sample_kerr_param[1] # Kerr parameter of Black Hole 2 (dimensionless spin)
}
################### CHANGE THE ESCAPE CONDITIONS HERE ###################
escape_velocity = {
'Nuclear': 500,
'Young': 30,
'Globular': 100,
}
print('\n********Initial parameters for the first black hole merger********')
print(initial_params)
print('\n\n********Escape velocities of different star cluster environments [km/s]********\n', escape_velocity)
print('\n You can change the magnitudes in the script.')
temp1 = input('Press Enter to continue...\nor to exit, type "exit"')
if temp1.lower() == 'exit':
exit()
print(f'\nStarting the simulation for {cluster_env} cluster environment...')
# Check the cluster environment
nuclear_escape, young_escape, globular_escape = False, False, False
if cluster_env.lower() == 'nuclear':
nuclear_escape = True
elif cluster_env.lower() in ['open', 'young']:
young_escape = True
elif cluster_env.lower() == 'globular':
globular_escape = True
# Create a directory Results to save the files
if os.path.exists(f'{output}/Results_data_{cluster_env}'):
shutil.rmtree(f'{output}/Results_data_{cluster_env}')
os.mkdir(f'{output}/Results_data_{cluster_env}')
else:
os.mkdir(f'{output}/Results_data_{cluster_env}')
# Create a directory in Results to save the plots
if os.path.exists(f'{output}/Plots_{cluster_env}'):
shutil.rmtree(f'{output}/Plots_{cluster_env}')
os.mkdir(f'{output}/Plots_{cluster_env}')
else:
os.mkdir(f'{output}/Plots_{cluster_env}')
if nuclear_escape:
# Run the simulation
result_one_sim, m1_one_sim, q1_one_sim = simulate_hierarchical_merging(op_sys, bh_mass, bh_kerr, n_max_gen, escape_velocity, initial_params, nuclear_escape=nuclear_escape, young_escape=young_escape, globular_escape=globular_escape)
print('One simulation completed successfully!\nYou can see the plots in the Results directory.')
# Save the results
plot_one_sim(result_one_sim, cluster_env, output)
sim_result, sim_result_extended_array, m1_vals, q_vals = run_n_simulations(op_sys, bh_mass, bh_kerr, n_simulations= n_sim, max_generation=n_max_gen,\
escape_velocity=escape_velocity, initial_params=initial_params,\
nuclear_escape=nuclear_escape, young_escape=young_escape, globular_escape=globular_escape)
print('\nSimulation completed successfully!')
os.chdir(f'{output}')
# Save the results
print('Saving the results...')
save_complete_simulation(f'Results_data_{cluster_env}/Simulation_param_evolution_{cluster_env}.pkl', sim_result)
save_concat_simulation_data(f'Results_data_{cluster_env}/Concatenated_param_evol_{cluster_env}.csv', sim_result_extended_array)
save_complete_simulation(f'Results_data_{cluster_env}/M1_vals_{cluster_env}.pkl', m1_vals)
save_complete_simulation(f'Results_data_{cluster_env}/q_vals_{cluster_env}.pkl', q_vals)
print('Results saved successfully!\n')
# Plot the inherent merging probabilities for different generations in the star cluster environment
inherent_merge_probability, gens_count, gens_distrbution = get_inherent_merge_probability(n_sim, sim_result, n_max_gen)
plot_merge_prob(inherent_merge_probability, gens_distrbution, n_max_gen, cluster_env)
plot_weighted_merge_prob(gens_distrbution, n_max_gen, cluster_env, q_vals, m1_vals)
print('\nPlots saved successfully!')
# Get black hole parameters for each generation of the complete simulation
generation_data = each_gen_data(sim_result, n_max_gen)
# Save the parameters for each generation in the simulation
print('\nSaving the parameters for each generation of the complete simulation...')
save_each_gen_params(generation_data, cluster_env)
print('Parameters saved successfully!')
print('\n\n*************************** Understand the results ****************************\n')
print("\nResults are stored in a dictionary format to access them as you wish, and also in .csv format, which contains the concatenated list parameter evolution.\n"\
"For example, the first column represents the mass of all the remnant black holes across the simulation, irrespective of its generation.\n"\
"M1 is the maximum of the two BH masses at each step/merger.\n"\
"Q is the mass ratio of the two merging black hole at every merger.")
print("\n\n*************************** Analyzing the results ****************************\n")
print('Creating plots to analyze the results...')
# Change directory to where the data is stored
os.chdir(f'{output}')
# Extract the data for each generation
simulation_param_data = np.array([np.loadtxt(f'Results_data_{cluster_env}/each_gen_data_Gen_{i}.txt', skiprows=1) for i in range(2, n_max_gen+1)])
concat_data = np.loadtxt(f'Results_data_{cluster_env}/Concatenated_param_evol_{cluster_env}.csv', skiprows=1)
try:
# Plotting correlation plot for the black hole parameters
plot_correlations(simulation_param_data, concat_data, cluster_env, n_max_gen)
except:
print('\nIncrease the number of simulations to plot the correlation plot.')
print('\nPlots saved successfully!')
print('\nExiting the script...')