From 98828f4104aaacdb0bb4de1c4ff76d6171941d38 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mika=20M=C3=A4ki?= Date: Tue, 15 Oct 2024 13:02:14 +0300 Subject: [PATCH] Fix discrepancies between the old and new spectra --- examples/plot_old_new.py | 51 +++++++++++++++++++++--------------- pttools/ssmtools/spectrum.py | 23 +++++++++++++--- 2 files changed, 50 insertions(+), 24 deletions(-) diff --git a/examples/plot_old_new.py b/examples/plot_old_new.py index 3e725a9..22adb80 100644 --- a/examples/plot_old_new.py +++ b/examples/plot_old_new.py @@ -9,7 +9,7 @@ import numpy as np from examples import utils -from pttools.analysis.utils import A4_PAPER_SIZE +from pttools.analysis.utils import A3_PAPER_SIZE from pttools.bubble import boundary from pttools.bubble.boundary import Phase, SolutionType from pttools.bubble.bubble import Bubble @@ -17,7 +17,7 @@ from pttools.bubble import relativity from pttools.models.model import Model from pttools.models.bag import BagModel -from pttools.ssmtools.spectrum import Spectrum, power_gw_scaled_bag, spec_den_v_bag +from pttools.ssmtools.spectrum import Spectrum, power_gw_scaled_bag, spec_den_v_bag, power_v_bag from tests.paper.plane import xiv_plane from tests.paper.plot_plane import plot_plane @@ -88,43 +88,52 @@ def main(): z = spectra[0].z data = xiv_plane(separate_phases=False) - fig: plt.Figure = plt.figure(figsize=A4_PAPER_SIZE) - ax, ax2, ax3 = fig.subplots(1, 3) - plot_plane(ax=ax, data_s=data, selected_solutions=False) + fig: plt.Figure = plt.figure(figsize=A3_PAPER_SIZE) + axs = fig.subplots(2, 2) + ax1 = axs[0, 0] + ax2 = axs[0, 1] + ax3 = axs[1, 0] + ax4 = axs[1, 1] + plot_plane(ax=ax1, data_s=data, selected_solutions=False) print("Solving & plotting old bubbles") for v_wall, alpha_n, sol_type in zip(v_walls, alpha_ns, sol_types): v, w, xi = fluid_bag.sound_shell_bag(v_wall=v_wall, alpha_n=alpha_n) - ax.plot(xi, v, color="blue", label=rf"$v_w={v_wall}, \alpha_n={alpha_n}$") + ax1.plot(xi, v, color="blue", label=rf"$v_w={v_wall}, \alpha_n={alpha_n}$") validate(bag, v, w, xi, sol_type) + label = rf"old, $v_w={v_wall}, \alpha_n={alpha_n}$" sdv = spec_den_v_bag(z, (v_wall, alpha_n)) - ax2.plot(z, sdv, label=rf"old, $v_w={v_wall}, \alpha_n={alpha_n}$") + ax2.plot(z, sdv, label=label) + + pow_v = power_v_bag(z, (v_wall, alpha_n)) + ax3.plot(z, pow_v, label=label) gw = power_gw_scaled_bag(z, (v_wall, alpha_n)) - ax3.plot(z, gw, label=rf"old, $v_w={v_wall}, \alpha_n={alpha_n}$") + ax4.plot(z, gw, label=label) print("Plotting new bubbles") for spectrum in spectra: bubble = spectrum.bubble - ax.plot(bubble.xi, bubble.v, ls=":", color="red") + ax1.plot(bubble.xi, bubble.v, ls=":", color="red") validate(bag, bubble.v, bubble.w, bubble.xi, bubble.sol_type) - ax2.plot(spectrum.z, spectrum.spec_den_v, label=rf"new, $v_w={bubble.v_wall}, \alpha_n={bubble.alpha_n}$") - ax3.plot(spectrum.z, spectrum.pow_gw, label=rf"new, $v_w={bubble.v_wall}, \alpha_n={bubble.alpha_n}$") + label = rf"new, $v_w={bubble.v_wall}, \alpha_n={bubble.alpha_n}$" + ax2.plot(spectrum.z, spectrum.spec_den_v, label=label) + ax3.plot(spectrum.z, spectrum.pow_v, label=label) + ax4.plot(spectrum.z, spectrum.pow_gw, label=label) - ax.legend() + ax2.set_ylabel("spec_den_v") + ax3.set_ylabel("pow_v") + ax4.set_ylabel(r"$\mathcal{P}_{\text{gw}}(z)$") - ax2.set_xscale("log") - ax2.set_yscale("log") - ax2.set_xlabel("$z = kR*$") - ax2.legend() + for ax in (ax2, ax3, ax4): + ax.set_xscale("log") + ax.set_yscale("log") + ax.set_xlabel("$z = kR*$") - ax3.set_xscale("log") - ax3.set_yscale("log") - ax3.set_xlabel("$z = kR*$") - ax3.set_ylabel(r"$\mathcal{P}_{\text{gw}}(z)$") - ax3.legend() + for ax in axs.flat: + ax.legend() fig.tight_layout() diff --git a/pttools/ssmtools/spectrum.py b/pttools/ssmtools/spectrum.py index b945143..cd36d81 100644 --- a/pttools/ssmtools/spectrum.py +++ b/pttools/ssmtools/spectrum.py @@ -87,11 +87,26 @@ def compute(self): nuc_type=self.nuc_type, nt=self.nt, z_st_thresh=self.z_st_thresh, cs=self.cs, return_a2=True ) self.pow_v = pow_spec(self.z, spec_den=self.spec_den_v) - # V2_pow_v = np.trapz(pow_v/self.z, self.z) - self.spec_den_gw, y = spec_den_gw_scaled(self.z, self.spec_den_v, cs=self.cs) + # This gives different spectra than the old code + self.spec_den_gw, y = spec_den_gw_scaled(z_lookup=self.z, P_v_lookup=self.spec_den_v, cs=self.cs) + + # This gives the same results as the old code + eps = 1e-8 # Seems to be needed for max(z) <= 100. Why? + nx = const.NPTDEFAULT[2] + 100 + # Defined on p. 12 between eq. 3.44 and 3.45 + # Todo: in the article the z here is y there ?! And y = k L_f + z_plus = np.max(self.z) * (0.5 * (1. + self.cs) / self.cs) + eps + z_minus = np.min(self.z) * (0.5 * (1. - self.cs) / self.cs) - eps + # The variable to integrate over in eq. 3.44 and 3.47 + x = np.logspace(np.log10(z_minus), np.log10(z_plus), nx) + sdv2, a2_v2 = spec_den_v( + bub=self.bubble, z=x, a=1., + nuc_type=self.nuc_type, nt=self.nt, z_st_thresh=self.z_st_thresh, cs=self.cs, return_a2=True + ) + self.spec_den_gw, y = spec_den_gw_scaled(z_lookup=x, P_v_lookup=sdv2, y=self.z, cs=self.cs) + self.pow_gw = pow_spec(self.z, spec_den=self.spec_den_gw) - # gw_power = np.trapz(self.pow_gw/y, y) # Plotting @@ -328,6 +343,8 @@ def _spec_den_gw_scaled_no_y( zmax = z_lookup.max() / (0.5 * (1. + cs) / cs) zmin = z_lookup.min() / (0.5 * (1. - cs) / cs) new_z = speedup.logspace(np.log10(zmin), np.log10(zmax), z_lookup.size) + # Todo: think if this is right + # return _spec_den_gw_scaled_core(new_z, P_v_lookup, z_lookup, cs, Gamma) return _spec_den_gw_scaled_core(z_lookup, P_v_lookup, new_z, cs, Gamma)