diff --git a/docs/gadfly/data/huber2011_sample.ecsv b/docs/gadfly/huber2011_sample.ecsv similarity index 100% rename from docs/gadfly/data/huber2011_sample.ecsv rename to docs/gadfly/huber2011_sample.ecsv diff --git a/docs/gadfly/start.rst b/docs/gadfly/start.rst index fca1cef..8b33abb 100644 --- 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 index d768074..e3424bd 100644 --- 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 index f268a1f..040008e 100644 --- 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,12 +247,19 @@ 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) @@ -250,52 +267,48 @@ def for_star( 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,7 +582,7 @@ 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)) @@ -577,7 +590,10 @@ def __init__(self, identifier_or_filter, download=False): 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 index a33a094..901e7e4 100644 --- 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, diff --git a/gadfly/scale.py b/gadfly/scale.py index 416b133..fb59da6 100644 --- a/gadfly/scale.py +++ b/gadfly/scale.py @@ -81,7 +81,7 @@ def _amplitudes_huber(mass, temperature, luminosity): @u.quantity_input(mass=u.g, temperature=u.K, luminosity=u.L_sun) -def _p_mode_amplitudes_huber(mass, temperature, luminosity): +def p_mode_amplitudes(mass, temperature, luminosity): """ p-mode oscillation power amplitudes. @@ -107,29 +107,39 @@ def _p_mode_amplitudes_huber(mass, temperature, luminosity): ) -def tau_osc(temperature): +def tau_osc(temperature, nu, nu_max, nu_solar): # mode lifetime - return 1 / (np.pi * _lifetimes_lund(temperature.value)) + return 1 / (np.pi * _lifetimes_lund(temperature, nu, nu_max, nu_solar)) -def _p_mode_amplitudes_kb2011( - mass, radius, temperature, luminosity, nu, wavelength, r=2.0 +def _p_mode_amplitudes_kallinger2014( + mass, radius, temperature, filter=None ): - # Kjeldsen & Bedding (2011), Eqn 20 - return ( - luminosity * tau_osc(temperature) ** 0.5 / - (wavelength * mass ** 1.5 * temperature ** (2.25 - r)) + # Kallinger et al. 2014, page 13 + g = mass / radius ** 2 + + bolometric_amplitude_scaling = ( + g ** -0.66 * mass ** -0.35 * temperature ** 0.04 ) + # scale with approach from Morris 2020 + if filter is None: + alpha = 1 + else: + alpha = amplitude_with_wavelength(filter, temperature) + + scaled_amplitude = alpha * bolometric_amplitude_scaling + + return scaled_amplitude + @u.quantity_input( mass=u.g, radius=u.m, temperature=u.K, luminosity=u.L_sun, wavelength=u.nm, nu=u.Hz, nu_solar=u.Hz ) -def p_mode_amplitudes( - mass, radius, temperature, luminosity, - nu, nu_solar, wavelength=550*u.nm +def _p_mode_amplitudes_kallinger( + mass, radius, temperature, filter ): """ p-mode oscillation power amplitudes. @@ -140,24 +150,27 @@ def p_mode_amplitudes( ---------- mass : ~astropy.units.Quantity Stellar mass + radius : ~astropy.units.Quantity + Stellar radius temperature : ~astropy.units.Quantity Effective temperature - luminosity : ~astropy.units.Quantity - Stellar luminosity + filter : str or ~tynt.Filter + Either the SVO FPS name for a filter bandpass available + via the ``tynt`` package, or an instance of the + :py:class:`~tynt.Filter` object itself. References ---------- .. [1] `Huber et al. (2011) `_ """ - num = _p_mode_amplitudes_kb2011( - mass, radius, temperature, luminosity, nu, wavelength + num = _p_mode_amplitudes_kallinger2014( + mass, radius, temperature, filter ) - denom = _p_mode_amplitudes_kb2011( - _solar_mass, _solar_radius, _solar_temperature, - _solar_luminosity, nu_solar, 550 * u.nm + denom = _p_mode_amplitudes_kallinger2014( + _solar_mass, _solar_radius, _solar_temperature ) - return (num / denom).to(u.dimensionless_unscaled).value + return (num / denom).to_value(u.dimensionless_unscaled) @u.quantity_input(mass=u.g, radius=u.m) @@ -288,29 +301,56 @@ def _fwhm_deprecated(temperature, quiet=False): ) -def _lifetimes_lund(temperature): +def _lifetimes_at_numax_lund(temperature): # Lund 2017 Eqn 32 Gamma_0 = 0.07 alpha = 0.91 beta = 15.3 return ( - Gamma_0 + alpha * (temperature / 5777) ** beta + Gamma_0 + alpha * (temperature.value / 5777) ** beta ) * u.uHz -def quality(nu_max, temperature): +def _lifetimes_with_nu_lund(nu, nu_max): + # parameter pairs correspond to columns "b" and "a" in Table 4 + # of Lund 2017 + def f(param): + # apply linear fit for parameter value as function + # of temperature + a, b = param + return b + a * float(nu_max / (3090 * u.uHz)) + + alpha = f([2.95, 0.39]) + Gamma_alpha = f([3.08, 3.32]) + + return np.exp( + # power law trend: + alpha * np.log(nu / nu_max) + np.log(Gamma_alpha) + 0 + ) + + +def _scale_lifetimes_with_nu_lund(nu, nu_max, nu_solar): + return ( + _lifetimes_with_nu_lund(nu, nu_max) / + _lifetimes_with_nu_lund(nu_solar, _solar_nu_max) + ) + + +def _lifetimes_lund(temperature, nu, nu_max, nu_solar): + return _lifetimes_with_nu_lund(nu, nu_max) + + +def quality(nu_solar, nu_stellar, nu_max, temperature): """ Scale the mode lifetime of p-mode oscillation peaks in power. Gets the mode lifetime scaling as a function of the p-mode central frequency before and after scaling, as shown in Figure 20 of Lund et al. (2017) [1]_, and the - scaling with frequency described by Eqn 30 of the same paper. + scaling with frequency described by Eqn 32 of the same paper. Parameters ---------- - nu_max : ~astropy.units.Quantity - Frequency of the maximum p-mode oscillations. temperature : ~astropy.units.Quantity Stellar temperature @@ -319,15 +359,23 @@ def quality(nu_max, temperature): .. [1] `Lund et al. (2017) `_ """ - Gamma_sun = _lifetimes_lund(_solar_temperature.value) - solar_Q = (_solar_nu_max / Gamma_sun).to(u.dimensionless_unscaled).value - Gamma_star = _lifetimes_lund(temperature.value) - star_Q = (nu_max / Gamma_star).to(u.dimensionless_unscaled).value - scale_factor_Q = star_Q / solar_Q + Gamma_sun = _lifetimes_lund(_solar_temperature, nu_solar, _solar_nu_max, nu_solar) + solar_Q = nu_solar / Gamma_sun + Gamma_star = _lifetimes_lund(temperature, nu_stellar, nu_max, nu_solar) + star_Q = nu_stellar / Gamma_star + scale_factor_Q = (star_Q / solar_Q).to_value(u.dimensionless_unscaled) return scale_factor_Q def _tau_gran(mass, temperature, luminosity): + """ + Granulation timescale scaling. + + Kjeldsen & Bedding (2011) Eqn 9 [1]_. + + .. [1] `Kjeldsen & Bedding (2011) + `_ + """ return luminosity / (mass * temperature ** 3.5) @@ -486,9 +534,9 @@ def _v_osc_kiefer_scaled(freq, nu_max, delta_nu): S = -0.1 # asymmetry parameter a = 3299 * 1e4 * u.m**2 / u.s**2 / u.Hz # height factor b = -581 * u.m**2 / u.s**2 / u.Hz # offset factor - A = 1 / np.pi * (np.arctan(S * (freq - nu_max) / Sigma).to(u.rad).value + 0.5) + A = 1 / np.pi * (np.arctan(S * (freq - nu_max) / Sigma).to_value(u.rad) + 0.5) voigt_profile = Voigt1D(x_0=nu_max, amplitude_L=a, fwhm_L=2*gamma, fwhm_G=2.355 * sigma) - return (A * (b + voigt_profile(freq)) * u.uHz).to(u.m**2/u.s**2).value + return (A * (b + voigt_profile(freq)) * u.uHz).to_value(u.m**2/u.s**2) @u.quantity_input(freq=u.uHz, wavelength=u.nm, temperature=u.K) @@ -579,7 +627,7 @@ def p_mode_intensity( # normalized so that at nu_max the relative amplitude is unity: relative_intensity = ( intensity_stellar_freq / intensity_stellar_numax - ).to(u.dimensionless_unscaled).value + ).to_value(u.dimensionless_unscaled) return relative_intensity diff --git a/setup.cfg b/setup.cfg index 3620e19..9f5e3d3 100644 --- a/setup.cfg +++ b/setup.cfg @@ -35,7 +35,6 @@ docs = sphinx-automodapi scipy lightkurve - scipy sphinx-book-theme>=0.3.3 numpydoc