amplitude fixes

bmorris3 · Jun 12, 2024 · 416255d · 416255d
1 parent c249174
commit 416255d
Show file tree

Hide file tree

Showing 7 changed files with 185 additions and 108 deletions.
diff --git a/docs/gadfly/data/huber2011_sample.ecsv → docs/gadfly/huber2011_sample.ecsv b/docs/gadfly/data/huber2011_sample.ecsv → docs/gadfly/huber2011_sample.ecsv
diff --git a/docs/gadfly/start.rst b/docs/gadfly/start.rst
@@ -50,10 +50,22 @@ observed power spectra.
     )
 
     # Plot the kernel's PSD, and the observed (binned) solar PSD:
-    fig, ax = kernel.plot(
+    fig, ax = plt.subplots(1, 2, figsize=(10, 3.5))
+
+    # Plot one PSD with coarse binning:
+    kernel.plot(
+        ax=ax[0],
+        # plot the observed power spectrum
+        obs=ps.bin(bins=50)
+    )
+
+    # Plot another PSD with finer binning and p-modes:
+    kernel.plot(
+        ax=ax[1],
         p_mode_inset=True,
         # also plot the observed power spectrum
-        obs=ps.bin(bins=100)
+        obs=ps.bin(bins=1_500),
+        obs_kwargs=dict(marker='.', mfc='C1', ms=1)
     )
 
 The high power at low frequencies corresponds to super-granulation, and

diff --git a/docs/gadfly/validation.rst b/docs/gadfly/validation.rst
@@ -117,7 +117,7 @@ Now we'll call a big loop to do most of the work:
         # plots:
         obs_kw = dict(color='k', marker='.', lw=0)
 
-        ps.bin(600).plot(
+        ps.bin(1_500).plot(
             ax=axis,
             kernel=kernel,
             freq=ps.frequency,
@@ -202,7 +202,7 @@ Ok, let's see the output:
         obs_kw = dict(color='k', marker='.', lw=0)
 
         # Plot the binned PSD of the light curve:
-        ps.bin(600).plot(
+        ps.bin(1_500).plot(
             ax=axis,
             kernel=kernel,
             freq=ps.frequency,
@@ -320,7 +320,7 @@ to see the code that generates this plot):
         ps = PowerSpectrum.from_light_curve(
             light_curve, name=name,
             detrend_poly_order=1
-        ).bin(200)
+        ).bin(1_500)
 
         kernel_kw = dict(color=f"C{i}", alpha=0.9)
         obs_kw = dict(color='k', marker='.', lw=0)
@@ -374,7 +374,7 @@ Kepler would observe for each target.
 
     # read a curated table of a few targets from Huber et al. (2011)
     sequence = Table.read(
-        'data/huber2011_sample.ecsv', format='ascii.ecsv'
+        'huber2011_sample.ecsv', format='ascii.ecsv'
     )[['KIC', 'mass', 'rad', 'teff', 'lum', 'cad_2']]
 
     # create a frequency grid for PSD estimates:
@@ -407,7 +407,7 @@ Kepler would observe for each target.
         ps = PowerSpectrum.from_light_curve(
             light_curve, name=name,
             detrend_poly_order=3,
-        ).bin(50)
+        ).bin(500)
 
         # adjust some plot settings:
         kernel_kw = dict(color=f"C{i}", alpha=0.9)

diff --git a/gadfly/core.py b/gadfly/core.py
@@ -177,6 +177,20 @@ def for_star(
         # get access to the astronomical filter bandpass via tynt:
         filt = Filter(bandpass)
 
+        # scale the amplitudes for the observing bandpass:
+        amp_with_wavelength = scale.amplitude_with_wavelength(
+            filt, temperature
+        )
+
+        # scale the amplitudes by a term for granulation
+        granulation_amp = scale.granulation_amplitude(
+            mass, temperature, luminosity
+        )
+
+        granulation_timescale = scale.tau_gran(
+            mass, temperature, luminosity
+        )
+
         # scale the hyperparameters for each of the granulation components
         scaled_hyperparameters = []
         for item in granulation_hyperparams:
@@ -185,19 +199,13 @@ def for_star(
             scale_S0 = (
                 params['S0'] *
                 # scale the amplitudes by a term for granulation:
-                scale.granulation_amplitude(
-                    mass, temperature, luminosity
-                ) *
+                granulation_amp *
                 # also scale the amplitudes for the observing bandpass:
-                scale.amplitude_with_wavelength(
-                    filt, temperature
-                )
+                amp_with_wavelength
             )
 
             # scale the timescales:
-            scaled_w0 = params['w0'] / scale.tau_gran(
-                mass, temperature, luminosity
-            )
+            scaled_w0 = params['w0'] / granulation_timescale
 
             if scaled_w0 > 0:
                 scaled_hyperparameters.append(
@@ -224,11 +232,13 @@ def for_star(
         # amplitude of granulation at the frequency where the
         # oscillations occur:
         solar_nu = solar_w0 / (2 * np.pi) * u.uHz
+
         granulation_background_solar = _sho_psd(
-            2 * np.pi * solar_nu.value[:, None],
-            solar_gran_S0[None, :],
-            solar_gran_w0[None, :], solar_gran_Q[None, :]
-        )
+            2 * np.pi * solar_nu[:, None],
+            solar_gran_S0[None, :] * u.cds.ppm**2 / u.uHz,
+            solar_gran_w0[None, :] * u.uHz,
+            solar_gran_Q[None, :]
+        ) * amp_with_wavelength
 
         scale_delta_nu = scale.delta_nu(mass, radius)
         solar_delta_nu = solar_nu - solar_nu_max
@@ -237,65 +247,68 @@ def for_star(
         scaled_w0 = 2 * np.pi * scaled_nu.to(u.uHz).value
 
         # limit to positive scaled frequencies:
-        # only_positive_omega = scaled_w0 > 0
-        # solar_nu = solar_nu[only_positive_omega]
-        # S0_fit = S0_fit[only_positive_omega]
-        # Q_fit = Q_fit[only_positive_omega]
-        # scaled_nu = scaled_nu[only_positive_omega]
-        # scaled_w0 = scaled_w0[only_positive_omega]
+        only_positive_omega = scaled_w0 > 0
+        solar_nu = solar_nu[only_positive_omega]
+        S0_fit = S0_fit[only_positive_omega]
+        Q_fit = Q_fit[only_positive_omega]
+        scaled_nu = scaled_nu[only_positive_omega]
+        scaled_w0 = scaled_w0[only_positive_omega]
+
+        # if bandpass isn't SOHO, use mean wavelength:
+        wavelength = (
+            filt.mean_wavelength
+            if filt.mean_wavelength is not None
+            else 550 * u.nm
+        )
 
         p_mode_scale_factor = (
             # scale p-mode "heights" like Kiefer et al. (2018)
             # as a function of frequency (scaled_nu)
             scale.p_mode_intensity(
                 temperature, scaled_nu, scaled_nu_max,
                 scale._solar_delta_nu * scale_delta_nu,
-                filt.mean_wavelength
+                wavelength
             ) *
-            # scale the p-mode amplitudes like Kjeldsen & Bedding (2011)
-            # according to stellar spectroscopic parameters:
+            # scale the p-mode amplitudes according to
+            # stellar spectroscopic parameters:
             scale.p_mode_amplitudes(
-                mass, radius, temperature, luminosity,
-                scaled_nu, solar_nu,
-                filt.mean_wavelength
+                mass, temperature, luminosity
             )
         )
+
         # scale the quality factors:
-        scale_factor_Q = scale.quality(
-            scaled_nu_max, temperature
+        scaled_Gamma = 1.02 * np.exp(
+            (temperature - scale._solar_temperature) / (436 * u.K)
         )
-        scaled_Q = Q_fit * scale_factor_Q
+        solar_Gamma = solar_nu.to_value(u.uHz) / Q_fit / 2  # [uHz]
+        scaled_Q = Q_fit * scaled_Gamma / solar_Gamma
 
         solar_psd_at_p_mode_peaks = _sho_psd(
-            2 * np.pi * solar_nu.value, S0_fit,
-            2 * np.pi * solar_nu.value, Q_fit
-        ) * u.cds.ppm ** 2 / u.uHz
+            2 * np.pi * solar_nu,
+            S0_fit * u.cds.ppm**2 / u.uHz,
+            solar_w0[only_positive_omega] * u.uHz,
+            Q_fit
+        )
 
         # Following Chaplin 2008 Eqn 3:
         A = 2 * np.sqrt(4 * np.pi * solar_nu * solar_psd_at_p_mode_peaks)
-        solar_mode_width = scale._lifetimes_lund(5777)
-        scaled_mode_width = scale._lifetimes_lund(temperature.value)
-        unscaled_height = 2 * A ** 2 / (np.pi * solar_mode_width)
+        unscaled_height = 2 * A ** 2 / (np.pi * solar_Gamma)
 
         scaled_height = (
-            unscaled_height *
-            # scale to trace the envelope in power with frequency:
-            p_mode_scale_factor *
-            # scale to the correct bandpass:
-            scale.amplitude_with_wavelength(
-                filt, temperature
-            )
+            unscaled_height * p_mode_scale_factor
         )
-        scaled_A = np.sqrt(np.pi * scaled_mode_width * scaled_height / 2)
+        scaled_A = np.sqrt(np.pi * scaled_Gamma * scaled_height / 2)
         scaled_psd_at_p_mode_peaks = (
             (scaled_A / 2)**2 / (4 * np.pi * scaled_nu)
-        ).to(u.cds.ppm**2/u.uHz).value
-        gran_background_solar = granulation_background_solar.sum(1)
+        ).to_value(u.cds.ppm**2/u.uHz)
 
         scaled_S0 = (
-            np.pi * scaled_psd_at_p_mode_peaks /
+            0.5 *
+            (np.pi / 2)**0.5 *
+            scaled_psd_at_p_mode_peaks /
             scaled_Q ** 2
-        ) * gran_background_solar
+        ) * granulation_background_solar.sum(1)[only_positive_omega].value
+
         scaled_w0 = np.ravel(
             np.repeat(scaled_w0[None, :], len(S0_fit), 0)
         )
@@ -569,15 +582,18 @@ def __init__(self, identifier_or_filter, download=False):
             identifier_or_filter = self.default_filter
 
         # Define a "bandpass" for SOHO, which is bolometric
-        if identifier_or_filter.upper() == 'SOHO VIRGO':
+        if isinstance(identifier_or_filter, str) and identifier_or_filter.upper() == 'SOHO VIRGO':
             wavelength = np.logspace(-1.5, 1.5, 1000) * u.um
             super().__init__(wavelength, np.ones_like(wavelength.value))
 
         else:
             if isinstance(identifier_or_filter, (tynt.Filter, Filter)):
                 # this happens if the user supplies a filter, we just
                 # return the same filter:
-                return identifier_or_filter
+                super().__init__(
+                    identifier_or_filter.wavelength,
+                    identifier_or_filter.transmittance
+                )
             else:
                 # otherwise try to reconstruct the bandpass transmittance from
                 # the tynt FFT parameterization:
@@ -600,4 +616,6 @@ def mean_wavelength(self):
         """
         Transmittance-weighted mean wavelength of the filter bandpass.
         """
+        if np.all(self.transmittance == 1):
+            return None
         return np.average(self.wavelength, weights=self.transmittance)
diff --git a/gadfly/data/hyperparameters.json b/gadfly/data/hyperparameters.json
@@ -1,8 +1,8 @@
 [
     {
         "hyperparameters": {
-            "S0": 18340.83781221723,
-            "w0": 5.353217604472537,
+            "S0": 18340.835205958254,
+            "w0": 5.35321843959353,
             "Q": 0.6
         },
         "metadata": {
@@ -11,8 +11,8 @@
     },
     {
         "hyperparameters": {
-            "S0": 17.898584214330256,
-            "w0": 114.55253839253504,
+            "S0": 17.89841172794779,
+            "w0": 114.55520319850163,
             "Q": 0.6
         },
         "metadata": {
@@ -21,8 +21,8 @@
     },
     {
         "hyperparameters": {
-            "S0": 0.5660643649317958,
-            "w0": 575.7650056817263,
+            "S0": 0.5659662607164637,
+            "w0": 575.893886406892,
             "Q": 0.6
         },
         "metadata": {
@@ -31,8 +31,8 @@
     },
     {
         "hyperparameters": {
-            "S0": 0.2959520205501044,
-            "w0": 5048.227233604119,
+            "S0": 0.2960462107291726,
+            "w0": 5050.277958149664,
             "Q": 0.6
         },
         "metadata": {
@@ -41,8 +41,8 @@
     },
     {
         "hyperparameters": {
-            "S0": 0.07174195850535738,
-            "w0": 20287.111882131547,
+            "S0": 0.07160842899044823,
+            "w0": 20343.609557336833,
             "Q": 0.6
         },
         "metadata": {
@@ -51,8 +51,8 @@
     },
     {
         "hyperparameters": {
-            "S0": 2.67678237950369e-06,
-            "Q": 2650.177949940339
+            "S0": 7.160744534863139e-06,
+            "Q": 2638.6200296239595
         },
         "metadata": {
             "degree": 0,
@@ -61,8 +61,8 @@
     },
     {
         "hyperparameters": {
-            "S0": 6.569552590518916e-06,
-            "Q": 1742.4758327254947
+            "S0": 1.7501924007383615e-05,
+            "Q": 1739.8359964306403
         },
         "metadata": {
             "degree": 1,
@@ -71,8 +71,8 @@
     },
     {
         "hyperparameters": {
-            "S0": 3.5515872723177826e-06,
-            "Q": 1284.9828290484409
+            "S0": 9.435829574482945e-06,
+            "Q": 1286.0000231757188
         },
         "metadata": {
             "degree": 2,
@@ -81,8 +81,8 @@
     },
     {
         "hyperparameters": {
-            "S0": 3.159395314588797e-07,
-            "Q": 1418.0862074222077
+            "S0": 8.392743868118509e-07,
+            "Q": 1424.2738543042017
         },
         "metadata": {
             "degree": 3,