Fix discrepancies between the old and new spectra

examples/plot_old_new.py

@@ -9,15 +9,15 @@
 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
 from pttools.bubble import fluid_bag
 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()
 

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)