From a3ef8e89e4715b554a2d34e50c8abd87daa6d14a Mon Sep 17 00:00:00 2001 From: Noe Brucy Date: Thu, 2 Jun 2022 10:37:15 +0200 Subject: [PATCH] add ismfeed helper --- ismfeed.py | 130 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 130 insertions(+) create mode 100644 ismfeed.py diff --git a/ismfeed.py b/ismfeed.py new file mode 100644 index 0000000..43a822b --- /dev/null +++ b/ismfeed.py @@ -0,0 +1,130 @@ +from scipy.integrate import solve_ivp +from plotter import U +import select_runs +import numpy as np +import pandas as pd + +ssfr_base = 2.5e-10 + + +def get_sinks(self, overwrite=False, stellar=True, convert_units=False): + + self.sinks_from_log(overwrite=overwrite) + self.coarse_step_from_log(overwrite=overwrite) + self.sinks = self.get_value("/series/sinks_from_log") + + if stellar: + self.stellar_from_log(overwrite=overwrite) + self.stellar = self.get_value("/dataset/stellar_from_log") + self.stellar["sn_time"] = {} + for run in self.runs: + self.stellar["sn_time"][run] = (self.stellar["time"][run] + self.stellar["lifetime"][run] )*self.info["unit_time"].express(U.year) + ind_sort = self.stellar["sn_time"][run].argsort() + for key in self.stellar: + self.stellar[key][run] = self.stellar[key][run][ind_sort] + self.stellar["cum_mass"] = {} + sn = {} + sn["time"] = {} + sn["cum_mass"] = {} + + self.sinks["cum_mass"] = {} + + for run in self.runs: + + if convert_units: + self.sinks["time"][run] *= self.info["unit_time"].express(U.year) + + self.sinks["cum_mass"][run] = self.sinks["mass_sink"][run].copy() + + if stellar: + self.stellar["cum_mass"][run] = np.cumsum(self.stellar["mass"][run]) * self.info["unit_mass"].express(U.Msun) + sn["time"][run], idx, count = np.unique(self.stellar["sn_time"][run], return_index=True, return_counts=True) + idx += count - 1 + sn["cum_mass"][run] = self.stellar["cum_mass"][run][idx] + ind_sn = np.searchsorted(self.sinks["time"][run], sn["time"][run]) + for i, indi in enumerate(ind_sn[ind_sn < ind_sn[-1]]): + self.sinks["cum_mass"][run][indi:ind_sn[i+1]] += sn["cum_mass"][run][i] + + + time_gas, total_mass, mass_gas = self.total_mass() + + sk_ssfr = {} + for run in self.runs: + time_gas[run] *= self.info["unit_time"].express(U.Myr) + sk_ssfr[run] = ssfr_base * 1e6 * (mass_gas[run]/1e6)**1.4 + + def dermass(t, sm, run, fact_sfr=1): + ind = min(np.searchsorted(time_gas[run], t), sk_ssfr[run].size - 1) + return sk_ssfr[run][ind]*fact_sfr + + tspan = np.linspace(0, 200, 100) + + self.time_esm = {} + self.expected_stellar_mass = {} + self.expected_stellar_mass_max = {} + self.expected_stellar_mass_min = {} + + for run in self.runs: + sol = solve_ivp(dermass, (tspan[0], tspan[-1]), [0], t_eval=tspan, args=[run]) + self.time_esm[run], self.expected_stellar_mass[run] = sol["t"], sol["y"][0] + sol_max = solve_ivp(dermass, (tspan[0], tspan[-1]), [0], t_eval=tspan, args=[run, 3]) + self.expected_stellar_mass_max[run] = sol_max["y"][0] + sol_min = solve_ivp(dermass, (tspan[0], tspan[-1]), [0], t_eval=tspan, args=[run, 1./3]) + self.expected_stellar_mass_min[run] = sol_min["y"][0] + +def compute_sfr(self, target_start=0.05, target_end=0.3): + mass = self.sinks["cum_mass"] + time = self.sinks["time"] + sfr = {} + sfr_err = {} + tend = {} + tstart = {} + for run in self.runs: + tend[run] = select_runs.time_mcons(self, run, target=target_end) * 1e6 + tstart[run] = select_runs.time_mcons(self, run, target=target_start) * 1e6 + if tend[run] < tstart[run]: + tend[run] = time[run][-1] + + + for run in self.runs: + idx1 = time[run].searchsorted(tstart[run]) + idx2 = min(time[run].searchsorted(tend[run]), len(time[run]) - 1) + + sfr[run] = (mass[run][idx2] - mass[run][idx1]) / (time[run][idx2] - time[run][idx1]) + sfr_other = [] + for i in range(idx1, idx2-10): + sfr_other.append((mass[run][i+10:idx2+1] - mass[run][i]) / (time[run][i+10:idx2+1] - time[run][i])) + + sfr_err[run] = np.std(np.concatenate(sfr_other)) + self.sfr = np.array(list(sfr.values())) + self.sfr_err = np.array(list(sfr_err.values())) + return self.sfr, self.sfr_err + + +def compute_power(self, target_start=0.05, target_end=0.3): + mass = self.sinks["cum_mass"] + time = self.sinks["time"] + sfr = {} + sfr_err = {} + tend = {} + tstart = {} + for run in self.runs: + tend[run] = select_runs.time_mcons(self, run, target=target_end) * 1e6 + tstart[run] = select_runs.time_mcons(self, run, target=target_start) * 1e6 + if tend[run] < tstart[run]: + tend[run] = time[run][-1] + + + for run in self.runs: + idx1 = time[run].searchsorted(tstart[run]) + idx2 = min(time[run].searchsorted(tend[run]), len(time[run]) - 1) + + sfr[run] = (mass[run][idx2] - mass[run][idx1]) / (time[run][idx2] - time[run][idx1]) + sfr_other = [] + for i in range(idx1, idx2-10): + sfr_other.append((mass[run][i+10:idx2+1] - mass[run][i]) / (time[run][i+10:idx2+1] - time[run][i])) + + sfr_err[run] = np.std(np.concatenate(sfr_other)) + self.sfr = np.array(list(sfr.values())) + self.sfr_err = np.array(list(sfr_err.values())) + return self.sfr, self.sfr_err