-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathscript_analysis-Agora_GEARv3.py
128 lines (86 loc) · 4.66 KB
/
script_analysis-Agora_GEARv3.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
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Jan 15 16:22:40 2025
@author: asier
"""
import yt
import numpy as np
import pandas as pd
from astropy.table import Table
from unyt import unyt_quantity, unyt_array
from yt.utilities.logger import ytLogger
import concurrent.futures
import galexquared as gal
from galexquared.class_methods import random_vector_spherical
ytLogger.setLevel(50)
gal.config.code = "GEAR"
files = ["GEAR_satellitesV5_Rvir1.0.csv", "GEAR_satellitesV5_Rvir1.5.csv", "GEAR_satellitesV5_Rvir2.0.csv", "GEAR_satellitesV5_Rvir3.0.csv", "GEAR_satellitesV5_Rvir4.0.csv", "GEAR_satellitesV5_Rvir5.0.csv"]
def analyzer(f):
print(f"Analyzing {f}")
tdir = "/home/asier/StructuralAnalysis/satellite_tables/"
errors = []
file = tdir + f
candidates = pd.read_csv(file)
candidates['delta_rel'] = pd.Series()
candidates['stellar_mass'] = pd.Series()
candidates['dark_mass'] = pd.Series()
candidates['all_gas_mass'] = pd.Series()
candidates['rh_stars_physical'] = pd.Series()
candidates['rh_dm_physical'] = pd.Series()
candidates['e_rh_stars_physical'] = pd.Series()
candidates['e_rh_dm_physical'] = pd.Series()
candidates['rh3D_stars_physical'] = pd.Series()
candidates['rh3D_dm_physical'] = pd.Series()
candidates['Mhl'] = pd.Series()
candidates['sigma*'] = pd.Series()
candidates['e_sigma*'] = pd.Series()
candidates['Mdyn'] = pd.Series()
candidates['e_Mdyn'] = pd.Series()
MT = gal.MergerTree("../merger_trees_csv/thinnedGEAR_CompleteTrees.csv")
MT.set_equivalence("GEAR_equivalence.dat")
pdir = "/media/asier/T7/GEAR/"
for index, halorow in candidates.iterrows():
print(f"##### {halorow['Sub_tree_id']}:snapshot-{halorow['Snapshot']} #####\n")
lines_of_sight = random_vector_spherical(N=20)
try:
halo = MT.load_halo(pdir, sub_tree=halorow["Sub_tree_id"], snapshot=halorow["Snapshot"])
if halo.stars.empty:
errors.append(f"Sub_tree_id {halorow['Sub_tree_id']} in snapshot-{halorow['Snapshot']} has no stars at all.")
continue
halo.compute_energy(refine=False)
halo.compute_stars_in_halo()
halo.switch_to_bound()
if halo.stars["mass"].sum() == 0:
errors.append(f"Sub_tree_id {halorow['Sub_tree_id']} in snapshot-{halorow['Snapshot']} has no stars bound.")
continue
halo.stars.refined_center6d(nmin=40)
halo.darkmatter.refined_center6d(nmin=50)
candidates.loc[index, "all_gas_mass"] = halo.gas["mass"].sum().to("Msun").value
candidates.loc[index, "dark_mass"] = halo.darkmatter["mass"].sum().to("Msun").value
candidates.loc[index, "stellar_mass"] = halo.stars["mass"].sum().to("Msun").value
candidates.loc[index, "delta_rel"] = float(halo.stars.delta_rel)
candidates.loc[index, 'rh3D_stars_physical'] = halo.stars.half_mass_radius()[0].in_units("kpc").value
candidates.loc[index, 'rh3D_dm_physical'] = halo.darkmatter.half_mass_radius()[0].in_units("kpc").value
rh_stars = halo.stars.half_mass_radius(project=True, lines_of_sight=lines_of_sight).to("kpc").value
rh_dm = halo.darkmatter.half_mass_radius(project=True, lines_of_sight=lines_of_sight).to("kpc").value
sigma = halo.stars.los_dispersion(lines_of_sight=lines_of_sight).to("km/s").value
mdyn = 580 * 1E3 * rh_stars * sigma ** 2
candidates.loc[index, "rh_stars_physical"] = rh_stars.mean()
candidates.loc[index, "e_rh_stars_physical"] = rh_stars.std()
candidates.loc[index, "rh_dm_physical"] = rh_dm.mean()
candidates.loc[index, "e_rh_dm_physical"] = rh_dm.std()
candidates.loc[index, "sigma*"] = sigma.mean()
candidates.loc[index, "e_sigma*"] = sigma.std()
candidates.loc[index, 'Mhl'] = halo.enclosed_mass(halo.stars.q["cm"], halo.stars.q["rh3d"], components="all").to("Msun").value
candidates.loc[index, 'Mdyn'] = mdyn.mean()
candidates.loc[index, 'e_Mdyn'] = mdyn.std()
except:
errors.append(f"Strange error ocurred in {halorow['Sub_tree_id']} in snapshot-{halorow['Snapshot']}")
continue
candidates.to_csv(file, index=False)
return f"{file} processed", errors
with concurrent.futures.ProcessPoolExecutor(max_workers=6) as executor:
results = list(executor.map(analyzer, files))
print(results)
print(f"Finished")