Amplitude fixes

.github/workflows/ci_tests.yml

@@ -26,17 +26,15 @@ jobs:
             python: 3.x
             toxenv: codestyle
 
-          - name: Python 3.9 with minimal dependencies
+          - name: Python 3.11 with minimal dependencies
             os: ubuntu-latest
-            python: 3.9
-            toxenv: py39-test
-            toxposargs: --remote-data=none
+            python: 3.11
+            toxenv: py311-test
 
-          - name: Python 3.9 with latest dev versions of key dependencies
+          - name: Python 3.11 with latest dev versions of key dependencies
             os: ubuntu-latest
-            python: 3.9
-            toxenv: py39-test-devdeps
-            toxposargs: --remote-data=none
+            python: 3.11
+            toxenv: py311-test-devdeps
 
           - name: Test building of Sphinx docs
             os: ubuntu-latest

.readthedocs.yml

@@ -3,11 +3,8 @@ version: 2
 build:
   image: latest
 
-conda:
-  environment: environment.yml
-
 python:
-  version: 3.8
+  version: 3.11
   install:
     - method: pip
       path: .

MANIFEST.in

@@ -1,6 +1,5 @@
 include README.rst
 include CHANGES.rst
-include setup.cfg
 include LICENSE.rst
 include pyproject.toml
 

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

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)

environment.yml

@@ -1,13 +1,2 @@
 name: gadfly
-channels:
-  - astropy
-dependencies:
-  - astropy
-  - pip:
-    - pytest-astropy
-    - celerite2
-    - matplotlib
-    - PyYAML
-    - scipy
-    - lightkurve
-    - tynt
+

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)