diff --git a/etl/prs/prs_analytical_representations.py b/etl/prs/prs_analytical_representations.py index bec3f562c30dc044314455e3f0aba767653b2b98..e780377e52a3dc4fe79a7ab78e4e0724e5e9b173 100644 --- a/etl/prs/prs_analytical_representations.py +++ b/etl/prs/prs_analytical_representations.py @@ -12,7 +12,7 @@ from typing import Union, List def approx(x, x_shift, b, exponent, scale): - arg = 2 * (x - x_shift) / scale + arg = (x - x_shift) / scale clean_arg = np.where(np.abs(arg) > 0.99999, np.sign(arg) * 0.99999, arg) return (np.arctanh(clean_arg) + b) ** exponent @@ -28,11 +28,11 @@ def main(ratios_to_fit: Union[List[str], None] = None): logger.info(f'Data contains {data.shape[0]} rows') x0 = { - '87-86': ((0.55, 6.75, 0.9, 1),), - '88-87': ((0.59, 6.9, 1.05, 0.96),), - '88-86': ((0.55, 6.3, 0.83, 1.08),), - '257-256': ((5, 4.45, 0.75, 8), (4.5, 8.5, 0.79, -7)), - '381-380': ((2.04, 5.1, 0.6, 2.1), (2.18, 5.7, 0.79, -2.4)), + '87-86': ((0.55, 6.75, 0.9, 0.5),), + '88-87': ((0.59, 6.9, 1.05, 0.48),), + '88-86': ((0.55, 6.3, 0.83, 0.54),), + '257-256': ((5, 4.45, 0.75, 4), (4.5, 8.5, 0.79, -3.5)), + '381-380': ((2.04, 5.1, 0.6, 1.05), (2.18, 5.7, 0.79, -1.2)), } min_ratio = { '87-86': 0.08, @@ -78,7 +78,7 @@ def main(ratios_to_fit: Union[List[str], None] = None): kde_max[1][clean_mask], np.log10(kde_max[0][clean_mask]), p0=guesses, - bounds=((0.1, 2, 0.3, -9), (7, 10, 1.5, 9))) + bounds=((0.1, 2, 0.3, -4.5), (7, 10, 1.5, 4.5))) approx_density = 10 ** approx(x_reg, *best_fit) print(best_fit) best_fit_params[ratio_string].append({param: float(element) for element, param in zip(best_fit, ['a', 'b', 'nu', 's'])}) diff --git a/etl/prs/prs_make_profiles_figure.py b/etl/prs/prs_make_profiles_figure.py index a6b8b57f865359c6a89c02cda0463a0f8b01d077..0757ec40d7991b77a2e6890f2d0d63e5a967cd38 100644 --- a/etl/prs/prs_make_profiles_figure.py +++ b/etl/prs/prs_make_profiles_figure.py @@ -9,11 +9,11 @@ q = -0.5 n_0 = 1e5 T_0 = 25 -radius = np.linspace(0, 0.9, 100) +radius = np.linspace(0, 0.9, 200) density = n_0 * (radius / 0.5) ** p -density[0] = 1e8 +density = np.where(density > 1e8, 1e8, density) temperature = T_0 * (radius / 0.5) ** q -temperature[0] = 2000 +temperature[0] = temperature[1] plt.clf() fig, ax = plt.subplots(layout='constrained') @@ -26,6 +26,7 @@ ln2 = ax2.plot(radius, temperature, color='red', label='T') lns = ln1 + ln2 labs = [l.get_label() for l in lns] ax.legend(lns, labs, loc=0) +ax2.set_ylim(10, 500) ax2.semilogy() ax2.set_ylabel('Temperature [K]') plt.savefig(os.path.join('prs', 'output', 'density_temperature_profiles.png'))