From 1a8db42d29536c29cc7e7b5a4067466a0f7b16bb Mon Sep 17 00:00:00 2001 From: Noe Brucy Date: Mon, 30 Jan 2023 12:12:14 +0100 Subject: [PATCH] black (still without mortimer) --- .flake8 | 2 + fragdisk/gui.py | 10 +- galactica/fragdisk_galactica.py | 4 +- galturb/galaxy.py | 2 +- ismfeed/ismfeed.py | 506 +++++++++++++++++++++++++++----- scripts/mngseg.py | 6 +- scripts/pspec.py | 10 +- turbox/__init__.py | 9 +- turbox/turbox.py | 6 +- utils/params.py | 6 +- utils/units.py | 12 +- 11 files changed, 474 insertions(+), 99 deletions(-) diff --git a/.flake8 b/.flake8 index 3551fe6..8f299df 100644 --- a/.flake8 +++ b/.flake8 @@ -7,3 +7,5 @@ ignore = W503 # Invalid escape sequence '\[' W605 + # Imported but unused + F401 diff --git a/fragdisk/gui.py b/fragdisk/gui.py index 52055ec..8862e9d 100644 --- a/fragdisk/gui.py +++ b/fragdisk/gui.py @@ -320,14 +320,14 @@ class InteractiveGUI: 10 ** (a * centers + b), "--", color="navy", - label=r"a = {:.3g}, $R^2$ = {:.3g}".format(a, r ** 2), + label=r"a = {:.3g}, $R^2$ = {:.3g}".format(a, r**2), ) plt.yscale("log") plt.title("PDF") else: self.step.set_ydata(values) self.fit.set_ydata(10 ** (a * centers + b)) - self.fit.set_label(r"a = {:.3g}, $R^2$ = {:.3g}".format(a, r ** 2)) + self.fit.set_label(r"a = {:.3g}, $R^2$ = {:.3g}".format(a, r**2)) plt.legend() # PROFILE @@ -382,7 +382,7 @@ class InteractiveGUI: a * np.log10(x[mask_fit_prof]) + b, "--", color="navy", - label=r"a = {:.3g}, $R^2$ = {:.3g}".format(a, r ** 2), + label=r"a = {:.3g}, $R^2$ = {:.3g}".format(a, r**2), ) self.ax_profile.set_xlabel(r"$r$") self.ax_profile.set_ylabel(r"$\log(\rho)$") @@ -390,7 +390,7 @@ class InteractiveGUI: self.prof.set_data(x, rho_prof) self.prof_z.set_data(z_vert, rho_vert) self.fit_prof.set_data(x[mask_fit_prof], a * np.log10(x[mask_fit_prof]) + b) - self.fit_prof.set_label(r"a = {:.3g}, $R^2$ = {:.3g}".format(a, r ** 2)) + self.fit_prof.set_label(r"a = {:.3g}, $R^2$ = {:.3g}".format(a, r**2)) plt.legend() plt.draw() @@ -486,7 +486,7 @@ class InteractiveGUI: y = np.linspace(self.im_extent[2], self.im_extent[3], self.shape[1]) self.xx, self.yy = np.meshgrid(x, y) self.xy = np.column_stack((self.xx.flatten(), self.yy.flatten())) - self.rr = np.sqrt(self.xx ** 2 + self.yy ** 2) + self.rr = np.sqrt(self.xx**2 + self.yy**2) self.rhocs = np.column_stack((self.frho_map.flatten(), self.fcs_map.flatten())) self.rmin = self.pp.params.disk.rmin_pdf / 2.0 diff --git a/galactica/fragdisk_galactica.py b/galactica/fragdisk_galactica.py index 1848e9f..ff36c80 100644 --- a/galactica/fragdisk_galactica.py +++ b/galactica/fragdisk_galactica.py @@ -186,7 +186,7 @@ rho_u = U.Unit.create_unit( ) orp = U.Unit.create_unit( "ORP", - base_unit=2 * np.pi * np.sqrt(0.25 ** 3) * info["unit_time"], + base_unit=2 * np.pi * np.sqrt(0.25**3) * info["unit_time"], ) vkd = U.Unit.create_unit( "vkd", @@ -200,7 +200,7 @@ Pd = U.Unit.create_unit( P_u = U.Unit.create_unit( "P_u", latex="v_\\mathrm{k,d}^2.\\mathrm{M}_\\star.r_\\mathrm{d}^{-3}", - base_unit=vkd ** 2 * rho_u, + base_unit=vkd**2 * rho_u, ) Sigmad = U.Unit.create_unit( diff --git a/galturb/galaxy.py b/galturb/galaxy.py index 28a24a1..2cd8ac3 100644 --- a/galturb/galaxy.py +++ b/galturb/galaxy.py @@ -118,7 +118,7 @@ def get_coldens(pp): xx, yy = np.meshgrid(x, y) # Physical radius - pp.rr_map = np.sqrt(xx ** 2 + yy ** 2) + pp.rr_map = np.sqrt(xx**2 + yy**2) pp.phi_map = np.angle(xx + yy * 1j) diff --git a/ismfeed/ismfeed.py b/ismfeed/ismfeed.py index b8a06b4..bd34478 100644 --- a/ismfeed/ismfeed.py +++ b/ismfeed/ismfeed.py @@ -1,24 +1,45 @@ from scipy.integrate import solve_ivp from plotter import U -import utils.snapshotselector +from snapshotprocessor import mean_by_bins +from utils import snapshotselector as select_snapshot + import numpy as np -import pandas as pd +import os +import matplotlib as mpl -ssfr_base = 2.5e-10 +import matplotlib.ticker as tick +import matplotlib.pyplot as plt +from matplotlib import FuncFormatter + +from utils.units import convert_exp + +ssfr_sun = 2.5e-3 + +mp = 1.4 * 1.66 * 10 ** (-24) * U.g +z0 = 150 * U.pc -def get_sinks(self, overwrite=False, stellar=True, convert_units=False): - +def convert_coldens_s(n0, z0=z0): + return (np.sqrt(2 * np.pi) * mp * z0 * (n0 * U.cm ** (-3))).express(U.coldens) + + +convert_coldens = np.vectorize(convert_coldens_s) + + +def get_sinks(self, overwrite=False, stellar=True, convert_units=False, sk=True): + 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) + 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] @@ -26,105 +47,452 @@ def get_sinks(self, overwrite=False, stellar=True, convert_units=False): 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() - + self.sinks["time"][run] *= self.info["unit_time"].express(U.year) # year + + self.sinks["cum_mass"][run] = self.sinks["mass_sink"][run].copy() # Msun + 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) + 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] - + 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 + if sk: + time_gas, total_mass, mass_gas = self.total_mass() # year, Msun, Msun - 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 + sk_ssfr = {} + for run in self.runs: + coldens = mass_gas[run] / 1e6 # in Msun/pc² + time_gas[run] *= self.info["unit_time"].express(U.Myr) # Myr + sk_ssfr[run] = ssfr_sun * (coldens / 10) ** 1.4 # Msun/kpc²/yr-1 - tspan = np.linspace(0, 200, 100) + 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 * 1e6 # Msun/Myr-1 + + tspan = np.linspace(0, 200, 100) # Myr + + 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], + ) # Myr, Msun + 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] # Myr, Msun + sol_min = solve_ivp( + dermass, (tspan[0], tspan[-1]), [0], t_eval=tspan, args=[run, 1.0 / 3] + ) + self.expected_stellar_mass_min[run] = sol_min["y"][0] # Myr, Msun - 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"] + mass = self.sinks["cum_mass"] # Msun + time = self.sinks["time"] # year 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 + tend[run] = ( + select_snapshot.time_mcons(self, run, target=target_end) * 1e6 + ) # year + tstart[run] = ( + select_snapshot.time_mcons(self, run, target=target_start) * 1e6 + ) # year 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[run] = (mass[run][idx2] - mass[run][idx1]) / ( + time[run][idx2] - time[run][idx1] + ) # Msun/year 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])) + 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]) + ) # Msun/year - 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())) + sfr_err[run] = np.std(np.concatenate(sfr_other)) # Msun/year + self.sfr = np.array(list(sfr.values())) # Msun/year + self.sfr_err = np.array(list(sfr_err.values())) # Msun/year 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 = {} +def get_nml_array( + pl, + xkey, + cmap=None, + logcmap=False, + put_colorbar=True, + cbar_fmt=tick.FormatStrFormatter("%.2g"), +): + pl.study.nml(xkey, overwrite=True) + x = pl.study.get_value(f"/comp/nml_{xkey}") + xname = os.path.basename(xkey) + + if xname in pl.value_convert: + x = np.array([pl.value_convert[xname](x_v) for x_v in x]) + if xname in pl.label_convert: + xlabel = pl.label_convert[xname] + else: + xlabel = xname + + colors = None + + if cmap is not None: + if logcmap: + sm = plt.cm.ScalarMappable( + cmap=cmap, norm=mpl.colors.LogNorm(vmin=min(x), vmax=max(x)) + ) + else: + sm = plt.cm.ScalarMappable( + cmap=cmap, norm=mpl.colors.Normalize(vmin=min(x), vmax=max(x)) + ) + print(x) + if put_colorbar: + cb = plt.colorbar(sm) + cb.set_label(xlabel) + cb.ax.yaxis.set_major_formatter(cbar_fmt) + + def colors(xi): + return sm.cmap(sm.norm(xi)) + + return x, xname, xlabel, colors + + +def histo_speed(pli, redo=False): + fig, axes = plt.subplots(1, 4, figsize=(16, 5), sharey=True) + fig2, axes2 = plt.subplots(1, 4, figsize=(16, 5), sharey=True) + fig3, axes3 = plt.subplots(1, 4, figsize=(16, 5), sharey=True) + + direction = ["x", "y", "z"] + + for pp, ax, ax2, ax3 in zip( + pli.get_snap_list( + select={"filter_nml": (pli.nml_key[0], "in", [1.5, 3, 6, 12])} + ), + axes, + axes2, + axes3, + ): + + pp.load_cells() + vel = pp.cells["vel"][:] * pp.info["unit_velocity"].express(U.km_s) + pos = pp.cells["pos"][:] * pp.info["unit_length"].express(U.pc) + + if redo: + pp.vmean_z = [] + pp.vmean_x = [] + + mass = ( + pp.cells["rho"] + * pp.info["unit_density"].express(U.Msun / U.pc**3) + * (pp.cells["dx"] * pp.info["unit_length"].express(U.pc)) ** 3 + ) + for i in range(3): + ax.hist( + vel[:, i], + range=[-100, 100], + bins=100, + histtype="step", + density=True, + ls="-", + label=f"$v_{direction[i]}$", + lw=2, + ) + ax.set_yscale("log") + ax.set_ylim(1e-3, 2e-1) + ax.set_xlabel("$v$ [km/s]") + ax.set_title(pli.get_label_run(nml_key=pli.nml_key[:2], run=pp.run)) + ax2.set_title(pli.get_label_run(nml_key=pli.nml_key[:2], run=pp.run)) + ax3.set_title(pli.get_label_run(nml_key=pli.nml_key[:2], run=pp.run)) + + if redo: + zbin, vmean_z = mean_by_bins(pos[:, 2], vel[:, i], weights=mass) + pp.zbin = zbin + pp.vmean_z.append(vmean_z) + + xbin, vmean_x = mean_by_bins(pos[:, 0], vel[:, i], weights=mass) + pp.xbin = zbin + pp.vmean_x.append(vmean_x) + else: + zbin = pp.zbin + xbin = pp.xbin + + vmean_x = pp.vmean_x[i] + vmean_z = pp.vmean_z[i] + + ax2.plot(zbin, vmean_z, label=f"$v_{direction[i]}$", lw=2) + ax2.set_xlabel("z [pc]") + + ax3.plot(xbin, vmean_x, label=f"$v_{direction[i]}$", lw=2) + ax3.set_xlabel("x [pc]") + + fig.suptitle(pli.name) + fig2.suptitle(pli.name) + fig3.suptitle(pli.name) + axes[0].set_ylabel("Mass weighted PDF") + axes[0].legend() + axes2[0].set_ylabel("v [km/s]") + axes2[0].legend() + axes3[0].set_ylabel("v [km/s]") + axes3[0].legend() + + +# def make_clump_hop(self, threshold_density=10): +# """ +# Apply HOP algorithm + +# Args: +# threshold_density (float): select only cells over threshold +# """ + +# # Selection of cells +# mask = ( +# self.cells["rho"] * self.info["unit_density"].express(U.H_cc) +# > threshold_density +# ) +# ncells = np.sum(mask) + +# # fill the matrice with ID, x,y,z and masses of particles +# cells_group = np.zeros((ncells, 14)) +# cells_group[:, 0] = np.arange(ncells) # index +# position = self.cells["pos"][mask] * self.info["unit_length"].express(U.pc) +# cells_group[:, 1:4] = position # position +# density = self.cells["rho"][mask] * self.info["unit_density"].express(U.H_cc) +# size = self.cells["dx"][mask] * self.info["unit_length"].express(U.pc) +# cells_group[:, 4] = ( +# density * size ** 3 * (U.H_cc * U.pc ** 3).express(U.Msun) +# ) # mass +# cells_group[:, 6] = density +# velocity = self.cells["vel"][mask] * self.info["unit_velocity"].express(U.km_s) +# pressure = self.cells["P"][mask] * self.info["unit_pressure"].express(U.bar) +# temperature = (pressure / density) * ( +# (U.bar / U.H_cc) * (1.4 * U.mH) / U.kB +# ).express(U.K) +# cells_group[:, 7] = temperature +# cells_group[:, 8] = pressure +# cells_group[:, 9:12] = velocity +# cells_group[:, 12] = self.cells["phi"][mask] +# cells_group[:, 13] = size ** 3 + +# # save file.txt +# head = str(ncells) +# np.savetxt( +# self.filename[:-3] + "_hop.txt", +# cells_group[:, :5], +# fmt="%10d %.10e %.10e %.10e %.10e", +# header=head, +# delimiter=" ", +# comments=" ", +# ) + +# # save file.den +# f = open(self.filename[:-3] + "_hop.den", "wb") +# f.write(pack("I", ncells)) +# self.cells["rho"][mask].astype("f").tofile(f) +# f.close() + +# # exec HOP algo +# h = HOP(self.filename[:-3] + "_hop.txt", os.path.dirname(self.filename)) +# h.process_hop(force=True) + +# # get the igroup array +# group_ids = h.get_group_ids() + +# # sort it and apply the sorting to the coordinates +# # this means that the particules of group 1 are written first then of group 2 etc... +# ind_sort = np.argsort(group_ids) +# cells_group = cells_group[ind_sort] +# cells_group[:, 5] = group_ids[ind_sort] + +# # Make it a pandas' DataFrame +# cells_group = pd.DataFrame( +# cells_group, +# columns=[ +# "id", +# "x", +# "y", +# "z", +# "mass", +# "group", +# "density", +# "temperature", +# "pressure", +# "vx", +# "vy", +# "vz", +# "phi", +# "volume", +# ], +# ) + +# self.clumps = cells_group + +# return cells_group + + +def plot_mass_sfr( + pl, + xkey="turb_params/comp_frac", + start=0.01, + end=0.4, + cmap=None, + logcmap=False, + redo=False, + fig_height=4, + ax_sfr=None, + logsfr=True, + marker="o", + color_sfr="k", + scale=0.9, + do_sfr=True, + plot_dir=".", +): + plt.figure(figsize=np.array([5, fig_height]) * scale) + get_sinks(pl.study, overwrite=redo) + + pl.value_convert["Bx"] = lambda x: x * 7.6189439 + pl.label_convert["Bx"] = r"B$_0$ [$\mu$G]" + pl.value_convert["bx_bound"] = lambda x: x * 7.6189439 + pl.label_convert["bx_bound"] = r"B$_0$ [$\mu$G]" + pl.label_convert["comp_frac"] = r"$\chi$" + + x_fmt = FuncFormatter(lambda x, p: f"{convert_exp(x,2)}") + x, xname, xlabel, colors = get_nml_array(pl, xkey, cmap, logcmap, cbar_fmt=x_fmt) + + try: + n0 = pl.study.get_nml("galbox_params/dens0", pl.runs[0]) + except KeyError(): + n0 = pl.study.get_nml("cloud_params/dens0", pl.runs[0]) + + mass = pl.study.sinks["cum_mass"] + time = pl.study.sinks["time"] 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 + + for i, run in enumerate(pl.runs): + tend[run] = (select_snapshot.time_mcons(pl.study, run, target=end) * 1e6,) + tstart[run] = select_snapshot.time_mcons(pl.study, run, target=start) * 1e6 if tend[run] < tstart[run]: tend[run] = time[run][-1] + idx1 = time[run].searchsorted(tend[run]) + idx2 = time[run].searchsorted(tstart[run]) + mask = time[run] < tend[run] - for run in self.runs: - idx1 = time[run].searchsorted(tstart[run]) - idx2 = min(time[run].searchsorted(tend[run]), len(time[run]) - 1) + if cmap is None: + (p,) = plt.plot( + time[run][mask] * 1e-6, + mass[run][mask] * 1e-6, + label=f"{xlabel} = {x[i]:.2f}", + ) + else: + color = colors(x[i]) + (p,) = plt.plot(time[run][mask] * 1e-6, mass[run][mask] * 1e-6, color=color) + plt.scatter( + [time[run][idx1] * 1e-6], + [mass[run][idx1] * 1e-6], + marker="<", + color=p.get_color(), + ) + plt.scatter( + [time[run][idx2] * 1e-6], + [mass[run][idx2] * 1e-6], + marker=">", + color=p.get_color(), + ) - 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])) + plt.title(f"$\Sigma_0$ = {convert_coldens(n0):.1f}" + " M$_\odot$.pc$^{-1}$") + plt.xlabel("Time [Myr]") + plt.ylabel("Mass accreted [$10^6$ M$_\odot$]") + plt.savefig(f"{plot_dir}/{xname}_n{n0}_mass.pdf") - 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 + ################################################## + if do_sfr: + label = f"$\Sigma_0$ = {convert_coldens(n0):.1f}" + " M$_\odot$.pc$^{-2}$" + + if ax_sfr is None: + plt.figure(figsize=(5, fig_height)) + else: + plt.sca(ax_sfr) + if redo or not hasattr(pl.study, "sfr"): + sfr, sfr_err = compute_sfr(pl.study, target_start=start, target_end=end) + else: + sfr, sfr_err = pl.study.sfr, pl.study.sfr_err + if logsfr: + plt.errorbar( + x=x, + y=np.log10(sfr), + yerr=sfr_err / sfr, + color=color_sfr, + marker=marker, + ls=":", + lw=0.5, + label=label, + elinewidth=1.5, + ) + plt.hlines( + np.log10( + ssfr_sun * (convert_coldens(n0) * (1 - end / 2.0) / 10) ** 1.4 + ), + xmin=min(x), + xmax=max(x), + color=color_sfr, + ls="--", + ) + plt.ylabel("log$(\Sigma_{\mathrm{SFR}})$ " + rf"[${U.ssfrK.latex}$]") + else: + plt.errorbar( + x=x, + y=sfr, + yerr=sfr_err, + color=color_sfr, + marker=marker, + ls=":", + lw=0.5, + label=label, + ) + plt.hlines( + ssfr_sun * (convert_coldens(n0) * (1 - end / 2.0) / 10) ** 1.4, + xmin=min(x), + xmax=max(x), + color=color_sfr, + ls="--", + ) + plt.ylabel("$\Sigma_{\mathrm{SFR}}$ " + rf"[${U.ssfrK.latex}$]") + + plt.xlabel(f"{xlabel}") + plt.gca().xaxis.set_major_formatter(x_fmt) + plt.legend() + + plt.savefig(f"{plot_dir}/{xname}_n{n0}_logsfr.pdf") diff --git a/scripts/mngseg.py b/scripts/mngseg.py index f3ed002..fb7d37d 100644 --- a/scripts/mngseg.py +++ b/scripts/mngseg.py @@ -99,7 +99,7 @@ def plot_each_scale(S11a, wav_k, q, label, coherent=False, reso=1): for scl in range(0, M): plt.figure(figsize=(6, 6)) # determine bins (large scales should have less bins) - nbins = np.int(nsize ** 2.0 * (wav_k[scl] * reso) ** 2.0) + nbins = np.int(nsize**2.0 * (wav_k[scl] * reso) ** 2.0) nbins = max(9, nbins) nbins = min(500, nbins) # calc histogram gaussian component w.r.t. its mean value (easier to compare) @@ -192,8 +192,8 @@ def plot_components_power_spectrum( # show resolution limits # Gaussian part not accurate below levelmin due to the way AMR works - sim_res_eff = 2 ** lvlmin / 10 - sim_res_lvl_min = 2 ** lvlmin + sim_res_eff = 2**lvlmin / 10 + sim_res_lvl_min = 2**lvlmin plt.plot( [sim_res_lvl_min / nsize, sim_res_lvl_min / nsize], [1e-6, 1e8], diff --git a/scripts/pspec.py b/scripts/pspec.py index e7f358b..be5b22b 100644 --- a/scripts/pspec.py +++ b/scripts/pspec.py @@ -111,7 +111,7 @@ def calc_k(n, nbinsk, nbig, dkbig, dim=3, saxis=2): a[saxis] = np.zeros_like(k) kx, ky, kz = np.meshgrid(a[0], a[1], a[2], indexing="ij") - cube_k = np.sqrt(kx ** 2 + ky ** 2 + kz ** 2) + cube_k = np.sqrt(kx**2 + ky**2 + kz**2) cubes_k = {"kx": kx, "ky": ky, "kz": kz, "k": cube_k} @@ -293,7 +293,7 @@ def proj_B(cubes_k, kbins, vec, var="", dim=3, saxis=2, update=False): vind = tuple(vind) vec_z = vec_z[vind + (slice(None),)] - vnorm = np.sqrt(np.sum(vec_z ** 2, axis=-1)) + vnorm = np.sqrt(np.sum(vec_z**2, axis=-1)) kpar = np.zeros_like(cubes_k["k"]) for i, d in enumerate(["x", "y", "z"]): @@ -572,7 +572,7 @@ def main(arg): distance=arg.size / 2.0, far_cut_depth=arg.size / 2.0, up_vector="y", - map_max_size=2 ** clvl, + map_max_size=2**clvl, ) cubes = {} @@ -584,7 +584,7 @@ def main(arg): operator = ScalarOperator(cube_vars[i], cube_units[i]) extractor = cube3d.CubeExtractor(amr, operator) cubes[v] = extractor.process( - cam, cube_size=arg.size, resolution=2 ** clvl + cam, cube_size=arg.size, resolution=2**clvl ).data else: h5f = T.open_file("cube.hdf5", "r") @@ -833,7 +833,7 @@ def main(arg): print("Compute 2D power spectra") # 2D power spectra ------------------------------------------------------------- - ns = 2 ** clvl + ns = 2**clvl f = "_%%(i)0%dd" % (np.floor(np.log10(ns)) + 1) for v in list(pcubes2.keys()): for i in range(ns): diff --git a/turbox/__init__.py b/turbox/__init__.py index 293dbe8..fdd7cf8 100644 --- a/turbox/__init__.py +++ b/turbox/__init__.py @@ -1,2 +1,7 @@ -from .turbox import get_pspec, get_pspec_slope, build_suite, apply_rule_pdf, span_resolution - +from .turbox import ( + get_pspec, + get_pspec_slope, + build_suite, + apply_rule_pdf, + span_resolution, +) diff --git a/turbox/turbox.py b/turbox/turbox.py index eee3af0..3d31d56 100644 --- a/turbox/turbox.py +++ b/turbox/turbox.py @@ -81,14 +81,14 @@ def get_pspec_slope(pp, field: str, resol: int, plotdebug: bool = False): f"Fit results in get_slope({field}, {resol}): slope:{results.slope:.2f}" + f", b:{results.intercept:.2f}, R2:{results.rvalue**2:.2f}" ) - if results.rvalue ** 2 < 0.8: + if results.rvalue**2 < 0.8: pp.logger.warning( f"Bad fit in get_slope({field}, {resol}) with {logkmin} <= logk < {logkmax}" ) pp.logger.warning(f"log(k) is \n {logk[mask]}") pp.logger.warning(f"log(power) is \n {logpower[mask]}") - return results.slope, results.intercept, results.rvalue ** 2 + return results.slope, results.intercept, results.rvalue**2 def build_suite(pl, redo=False, cs0=0.28834810480560674): @@ -155,7 +155,7 @@ def build_suite(pl, redo=False, cs0=0.28834810480560674): ) ) - df["sigma"] = list(map(lambda l: np.mean(l), df["sigma_all"].values)) + df["sigma"] = list(map(lambda sig_list: np.mean(sig_list), df["sigma_all"].values)) df["Mach"] = df["sigma"] / cs0 df["turnover"] = (df["L"] * U.pc.express(U.km) / (2 * df["sigma"])) * U.s.express( U.Myr diff --git a/utils/params.py b/utils/params.py index 3f7ca68..9b3e96a 100644 --- a/utils/params.py +++ b/utils/params.py @@ -11,9 +11,9 @@ _dir_path = os.path.dirname(os.path.realpath(__file__)) # Add support for '1e3' kind of float _loader = yaml.SafeLoader _loader.add_implicit_resolver( - u"tag:yaml.org,2002:float", + "tag:yaml.org,2002:float", re.compile( - u"""^(?: + """^(?: [-+]?(?:[0-9][0-9_]*)\\.[0-9_]*(?:[eE][-+]?[0-9]+)? |[-+]?(?:[0-9][0-9_]*)(?:[eE][-+]?[0-9]+) |\\.[0-9_]+(?:[eE][-+][0-9]+)? @@ -22,7 +22,7 @@ _loader.add_implicit_resolver( |\\.(?:nan|NaN|NAN))$""", re.X, ), - list(u"-+0123456789."), + list("-+0123456789."), ) diff --git a/utils/units.py b/utils/units.py index 7898020..bfd1a4d 100644 --- a/utils/units.py +++ b/utils/units.py @@ -62,25 +62,25 @@ def unit_str(unit, base=None, prefix="", format=" [{unit}]"): U.coldens = create_unit( - "Msun.pc^-2", base_unit=U.Msun / U.pc ** 2, descr="Column density" + "Msun.pc^-2", base_unit=U.Msun / U.pc**2, descr="Column density" ) U.km_s = create_unit("km.s^-1", base_unit=U.km / U.s, descr="Speed") -U.Msun_pc3 = create_unit("Msun.pc^-3", base_unit=U.Msun / U.pc ** 3, descr="Density") +U.Msun_pc3 = create_unit("Msun.pc^-3", base_unit=U.Msun / U.pc**3, descr="Density") -U.kg_m3 = create_unit("kg.m^-3", base_unit=U.kg / U.m ** 3, descr="Density") +U.kg_m3 = create_unit("kg.m^-3", base_unit=U.kg / U.m**3, descr="Density") U.ssfr = create_unit( "Msun.year^-1.pc^-2", - base_unit=U.Msun / U.year / U.pc ** 2, + base_unit=U.Msun / U.year / U.pc**2, descr="Surfacic SFR", ) # latex='M$_{\odot}$.yr$^{-1}$.pc$^{-2}$') U.ssfrG = create_unit( "Msun.Gyr^-1.pc^-2", - base_unit=1e-9 * U.Msun / U.year / U.pc ** 2, + base_unit=1e-9 * U.Msun / U.year / U.pc**2, descr="Surfacic SFR", latex="\mathrm{M}_{\odot}.\mathrm{Gyr}^{-1}.\mathrm{pc}^{-2}", ) @@ -92,7 +92,7 @@ U.uG = create_unit( U.ssfrK = create_unit( "Msun.year^-1.kpc^-2", - base_unit=U.Msun / U.year / U.kpc ** 2, + base_unit=U.Msun / U.year / U.kpc**2, descr="Surfacic SFR", latex="\mathrm{M}_{\odot}.\mathrm{yr}^{-1}.\mathrm{kpc}^{-2}", )