Adding amplitude scaling for p-modes

docs/index.rst

@@ -51,21 +51,33 @@ Scaling relations for p-modes
 For computation of stellar p-mode oscillation frequencies, we use the scaling
 relations found in
 `Huber et al. (2011) <http://adsabs.harvard.edu/abs/2011ApJ...743..143H>`_ and
-references therein (e.g. Kjeldsen & Bedding 1995), namely:
+references therein (e.g. Kjeldsen & Bedding 1995), namely Equation 1:
 
 .. math::
 
     \nu_\textrm{max} \approx \frac{M / M_\odot (T_\textrm{eff}/
     T_{\textrm{eff},\odot})^{3.5} }{L/L_\odot} \nu_{\textrm{max}, \odot}
 
-(Equation 1), and
+and Equation 2
 
 .. math::
 
     \Delta \nu_\textrm{max} \approx \frac{(M / M_\odot)^{0.5} (T_\textrm{eff}/
-    T_{\textrm{eff},\odot})^{3} }{(L/L_\odot)^{0.75}} \Delta \nu_{\odot}
+    T_{\textrm{eff},\odot})^{3} }{(L/L_\odot)^{0.75}} \Delta \nu_{\odot}.
+
+
+The amplitude scaling of the p-mode oscillations is given by Equation 9:
+
+.. math::
+
+    A \propto \frac{L^s}{M^t T_\textrm{eff}^{r-1} c(T_\textrm{eff})}
+
+where :math:`r = 2`, :math:`s = 0.886`, :math:`t = 1.89` and
+
+.. math::
+
+    c(T_\textrm{eff}) = \left( \frac{T_\textrm{eff}}{5934 \textrm{K}} \right)^{0.8}.
 
-(Equation 2).
 
 Scaling relations for granulation
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

example.py

@@ -12,7 +12,7 @@
 
 fig, ax = plt.subplots(3, 2, figsize=(14, 8))
 ax[0, 0].plot(x, y)
-ax[0, 0].set(xlabel='Time [day]', ylabel='Flux')
+ax[0, 0].set(xlabel='Time [day]', ylabel='Flux', title='Sun')
 
 ftest, Ptest = periodogram(y, fs=1/60)
 
@@ -33,7 +33,7 @@
 x = np.arange(len(y)) / 60 / 24
 
 ax[1, 0].plot(x, y)
-ax[1, 0].set(xlabel='Time [day]', ylabel='Flux')
+ax[1, 0].set(xlabel='Time [day]', ylabel='Flux', title='Kepler-62')
 
 ftest, Ptest = periodogram(y, fs=1/60)
 
@@ -54,7 +54,7 @@
 x = np.arange(len(y)) / 60 / 60 / 24
 
 ax[2, 0].plot(x, y)
-ax[2, 0].set(xlabel='Time [day]', ylabel='Flux')
+ax[2, 0].set(xlabel='Time [day]', ylabel='Flux', title='Kepler-296')
 
 ftest, Ptest = periodogram(y, fs=1)
 

shocksgo/core.py

@@ -116,8 +116,9 @@ def generate_stellar_fluxes(size, M, T_eff, L, cadence=60*u.s):
     # Scale p-mode frequencies
     ##########################
 
+    # Scale frequencies
     tunable_amps = np.exp(parameter_vector[::3][2:])
-    tunable_freqs = np.exp(parameter_vector[2::3][2:]) * 1e6/2/np.pi
+    tunable_freqs = np.exp(parameter_vector[1::3][2:]) * 1e6/2/np.pi
     peak_ind = np.argmax(tunable_amps)
     peak_freq = tunable_freqs[peak_ind]
     delta_freqs = tunable_freqs - peak_freq
@@ -137,7 +138,21 @@ def generate_stellar_fluxes(size, M, T_eff, L, cadence=60*u.s):
 
     new_log_omegas = np.log(2*np.pi*new_freqs*1e-6).value
 
-    parameter_vector[2::3][2:] = new_log_omegas
+    parameter_vector[1::3][2:] = new_log_omegas
+
+    # Scale amplitudes
+    c = (T_eff/(5934 * u.K))**0.8
+    c_sun = ((5777 * u.K)/(5934 * u.K))**0.8
+    s = 0.886
+    r = 2
+    t = 1.89
+    pmode_amp_factor = (L/L_sun)**s / ((M/M_sun)**t * T_eff.value**(r-1) * c)
+    pmode_amp_factor_sun = (L_sun/L_sun)**s / ((M_sun/M_sun)**t * 5777**(r-1)
+                                               * c_sun)
+
+    new_pmode_amps = np.log(np.exp(parameter_vector[0::3][2:]) *
+                            pmode_amp_factor/pmode_amp_factor_sun)
+    parameter_vector[0::3][2:] = new_pmode_amps
 
     #############################
     # Scale granulation frequency
@@ -146,7 +161,8 @@ def generate_stellar_fluxes(size, M, T_eff, L, cadence=60*u.s):
     tau_eff_factor = (new_peak_freq/peak_freq)**-0.89
     granulation_amplitude_factor = (new_peak_freq/peak_freq)**-0.5
     parameter_vector[4] = np.log(np.exp(parameter_vector[4]) * tau_eff_factor)
-    parameter_vector[3] = np.log(np.exp(parameter_vector[3]) * granulation_amplitude_factor)
+    parameter_vector[3] = np.log(np.exp(parameter_vector[3]) *
+                                 granulation_amplitude_factor)
 
     ##########################
     # Assemble celerite kernel