mirror of
https://github.com/Findus23/rebound-collisions.git
synced 2024-09-18 14:43:49 +02:00
120 lines
3.8 KiB
Python
120 lines
3.8 KiB
Python
import random
|
|
from sys import argv
|
|
|
|
from matplotlib import pyplot as plt
|
|
from matplotlib.axes import Axes
|
|
from matplotlib.figure import Figure
|
|
from rebound import SimulationArchive, Simulation
|
|
from scipy.constants import mega
|
|
|
|
from extradata import ExtraData, CollisionMeta
|
|
from utils import filename_from_argv, earth_mass, earth_water_mass, plot_settings, is_ci, is_potentially_habitable
|
|
|
|
files = argv[1:]
|
|
multifile = len(files) > 1
|
|
|
|
plot_settings()
|
|
fig: Figure = plt.figure()
|
|
ax_masses: Axes = fig.add_subplot(2, 1, 1)
|
|
ax_wmfs: Axes = fig.add_subplot(2, 1, 2)
|
|
|
|
ax_wmfs.set_ylabel("water mass fraction")
|
|
ax_masses.set_ylabel("masses [kg]")
|
|
for ax in [ax_wmfs, ax_masses]:
|
|
ax.set_xlim(1e4, 200 * mega)
|
|
ax.set_xlabel("time [yr]")
|
|
ax.set_xscale("log")
|
|
ax.set_yscale("log")
|
|
|
|
earth_mass_ax = ax_masses.secondary_yaxis("right",
|
|
functions=(lambda x: x / earth_mass, lambda x: x * earth_mass))
|
|
earth_mass_ax.set_ylabel('masses [$M_\\oplus$]')
|
|
ax_wmfs.axhline(earth_water_mass / earth_mass, linestyle="dotted")
|
|
|
|
num_formed_planets = 0
|
|
num_large_planets = 0
|
|
num_habitable_planets = 0
|
|
num_water_rich_planets = 0
|
|
|
|
random.seed(1)
|
|
random.shuffle(files)
|
|
for file in files:
|
|
fn = filename_from_argv(file)
|
|
|
|
ed = ExtraData.load(fn)
|
|
sa = SimulationArchive(str(fn.with_suffix(".bin")))
|
|
|
|
last_sim: Simulation = sa[-1]
|
|
print([p.hash.value for p in last_sim.particles])
|
|
print(last_sim.t)
|
|
|
|
for particle in last_sim.particles:
|
|
if ed.pd(particle).type in ["sun", "gas giant"]:
|
|
continue
|
|
# if not is_potentially_habitable(particle):
|
|
# continue
|
|
masses = []
|
|
objects = []
|
|
times = []
|
|
hash = particle.hash.value
|
|
objects.append(ed.pdata[hash])
|
|
times.append(ed.meta.current_time)
|
|
|
|
num_formed_planets += 1
|
|
if particle.m > .6 * earth_mass:
|
|
num_large_planets += 1
|
|
if is_potentially_habitable(particle):
|
|
num_habitable_planets += 1
|
|
if ed.pd(particle).water_mass_fraction > 1e-4:
|
|
num_water_rich_planets += 1
|
|
|
|
while True:
|
|
print(f"looking at {hash}")
|
|
try:
|
|
collision = ed.tree.get_tree()[hash]
|
|
except KeyError:
|
|
print("found end of the tree")
|
|
break
|
|
meta: CollisionMeta = collision["meta"]
|
|
parents = collision["parents"]
|
|
print("mass:", ed.pdata[hash].total_mass / earth_mass)
|
|
masses.append(ed.pdata[hash].total_mass)
|
|
objects.append(ed.pdata[hash])
|
|
times.append(meta.time)
|
|
# print(collision)
|
|
if ed.pdata[parents[0]].total_mass > ed.pdata[parents[1]].total_mass:
|
|
hash = parents[0]
|
|
else:
|
|
hash = parents[1]
|
|
objects.append(ed.pdata[hash])
|
|
times.append(0)
|
|
if len(times) < 3:
|
|
continue
|
|
masses = [p.total_mass for p in objects]
|
|
wmfs = [p.water_mass_fraction for p in objects]
|
|
figs = []
|
|
if multifile:
|
|
args = {
|
|
"linewidth": 1,
|
|
"color": "black",
|
|
"alpha": .3
|
|
}
|
|
else:
|
|
args = {}
|
|
ax_masses.step(times, masses, label=particle.hash.value, **args)
|
|
ax_wmfs.step(times, wmfs, label=particle.hash.value, **args)
|
|
|
|
if not multifile:
|
|
for ax in [ax_wmfs, ax_masses]:
|
|
ax.legend()
|
|
|
|
print(num_large_planets / num_formed_planets)
|
|
habitable_large_planet_fraction = num_habitable_planets / num_large_planets
|
|
water_rich_planet_fraction = num_water_rich_planets / num_habitable_planets
|
|
print(habitable_large_planet_fraction)
|
|
print(water_rich_planet_fraction)
|
|
|
|
fig.tight_layout()
|
|
if not is_ci():
|
|
fig.savefig(f"/home/lukas/tmp/collisionhistory_{fn.name}.pdf", transparent=True)
|
|
plt.show()
|