updates paper figures per refereeing

GYRE/generate_magnitude_spectra.py

@@ -12,6 +12,7 @@
 output_file = 'magnitude_spectra.h5'
 
 star_dirs = ['03msol_Zsolar', '40msol_Zsolar', '15msol_ZLMC', '15msol_ZLMC']
+amp_corr = 2/5
 luminosity_amplitudes = [7.34e-15, 5.3e-10, 2.33e-11, 3.004e-10]
 min_trustworthy = np.array([3e-2, 3e-2, 3e-2, 3e-2])/(24*60*60) #1/d -> Hz
 Lmax = 15
@@ -29,7 +30,7 @@
 specLums = []
 logTeffs = []
 for i, sdir in enumerate(star_dirs):
-    wave_luminosity = lambda f, l: luminosity_amplitudes[i]*f**(-6.5)*np.sqrt(l*(l+1))**4
+    wave_luminosity = lambda f, l: (amp_corr**2)*luminosity_amplitudes[i]*f**(-6.5)*np.sqrt(l*(l+1))**4
     transfer_oms = []
     transfer_signal = []
     pure_transfers  = []

dedalus/zams_15Msol_LMC/wave_propagation/nr128/transfer.py

@@ -19,8 +19,6 @@
 from compstar.waves.transfer import calculate_refined_transfer
 
 plot = False
-forcing_width = 0.1
-forcing_center = 1.1
 
 #Grab relevant information about the simulation stratification.
 out_dir, out_file = name_star()
@@ -63,8 +61,8 @@
         else:
             raise ValueError("Cannot find MESA profile file in {} or {}".format(config.star['path'], stock_file_path))
     core_cz_radius = find_core_cz_radius(mesa_file_path)
-    min_forcing_radius = forcing_center - forcing_width / 2
-    max_forcing_radius = forcing_center + forcing_width / 2
+    min_forcing_radius = 1.005
+    max_forcing_radius = 1.20 
     N2_force_max = N2(max_forcing_radius)
     N2_max = N2s.max()
 #    N2_adjust = np.sqrt(N2_max/N2_force_max)

dedalus/zams_15Msol_LMC/wave_propagation/nr256/transfer.py

@@ -19,8 +19,6 @@
 from compstar.waves.transfer import calculate_refined_transfer
 
 plot = False
-forcing_width = 0.1
-forcing_center = 1.1
 
 #Grab relevant information about the simulation stratification.
 out_dir, out_file = name_star()
@@ -63,8 +61,8 @@
         else:
             raise ValueError("Cannot find MESA profile file in {} or {}".format(config.star['path'], stock_file_path))
     core_cz_radius = find_core_cz_radius(mesa_file_path)
-    min_forcing_radius = forcing_center - forcing_width / 2
-    max_forcing_radius = forcing_center + forcing_width / 2
+    min_forcing_radius = 1.005
+    max_forcing_radius = 1.20 
     N2_force_max = N2(max_forcing_radius)
     N2_max = N2s.max()
 #    N2_adjust = np.sqrt(N2_max/N2_force_max)

figures/manuscript/plot_fig02_variability_spectra.py

@@ -176,7 +176,7 @@
 ax3.set_xlabel(r'$\log_{10}\, $T$_{\rm eff}/$K')
 
 plt.axes(ax1)
-ax1.text(0.06, 1.5e-1, 'Predicted wave signal $(15 M_{\odot})$', ha='left', color=cmap.mpl_colors[2])
+ax1.text(0.063, 1e-1, 'Predicted wave signal $(15 M_{\odot})$', ha='left', color=cmap.mpl_colors[2])
 ax1.text(0.25, 3.5e2, 'Observed red noise', ha='left', va='center', color=cmap.mpl_colors[3])
 
 for i in range(len(LogL)):
@@ -200,18 +200,18 @@
 
 #con2 = ConnectionPatch(xyA=(1e1,1e-4), xyB=(4e-2,1e-4), coordsA='data', coordsB='data', axesA=ax1, axesB=ax2, color='grey', lw=0.5)
 #ax1.add_artist(con2)
-con1 = ConnectionPatch(xyA=(1e1,1e0), xyB=(5e-2,1e0), coordsA='data', coordsB='data', axesA=ax1, axesB=ax2, color='grey', lw=1)
+con1 = ConnectionPatch(xyA=(1e1,2e-1), xyB=(6e-2,2e-1), coordsA='data', coordsB='data', axesA=ax1, axesB=ax2, color='grey', lw=1)
 ax1.add_artist(con1)
-ax1.plot([7e0,1e1],[1e0, 1e0], c='grey', lw=1)
+ax1.plot([7e0,1e1],[2e-1, 2e-1], c='grey', lw=1)
 
 #ax1.text(0.12, 0.04, r'15 $M_{\odot}$', color=cmap.mpl_colors[2], ha='center', va='center', size=8)
-ax2.text(8e-2, 3.3e-1, r'40 $M_{\odot}$', color=cmap.mpl_colors[1], ha='center', va='center', size=8)
-ax2.text(8e-2, 1.1e-1, r'15 $M_{\odot}$', color=cmap.mpl_colors[2], ha='center', va='center', size=8)
-ax2.text(1.3e-1, 1.1e-2, r'3 $M_{\odot}$', color=cmap.mpl_colors[0], ha='center', va='center', size=8)
-ax2.set_ylim(5e-4, 1e0)
+ax2.text(6.3e-2, 1.7e-1, r'40 $M_{\odot}$', color=cmap.mpl_colors[1], ha='left', va='center', size=8)
+ax2.text(6.5e-2, 7e-2, r'15 $M_{\odot}$', color=cmap.mpl_colors[2], ha='left', va='center', size=8)
+ax2.text(1.3e-1, 6e-3, r'3 $M_{\odot}$', color=cmap.mpl_colors[0], ha='center', va='center', size=8)
+ax2.set_ylim(5e-4, 2e-1)
 ax2.set_xlabel(r'frequency (d$^{-1}$)')
 for ax in [ax1, ax2]:
-    ax.set_xlim(5e-2, 1e1)
+    ax.set_xlim(6e-2, 1e1)
 
 plt.savefig('fig02_obs_prediction.png', bbox_inches='tight', dpi=300)
 plt.savefig('fig02_obs_prediction.pdf', bbox_inches='tight', dpi=300)

figures/manuscript/plot_fig03_scatter_plot.py

@@ -13,8 +13,8 @@
 #Get simulation details.
 star_dirs = ['03msol_Zsolar', '40msol_Zsolar', '15msol_ZLMC']
 sim_mass = [3, 40, 15]
-sim_alpha = [9e-3, 0.3, 0.12]
-sim_nuchar = [0.26, 0.1, 0.16]
+sim_alpha = [5.5e-3, 0.16, 0.06]
+sim_nuchar = [0.3, 0.13, 0.22]
 sim_specLums = []
 sim_logTeffs = []
 for i, sdir in enumerate(star_dirs):

figures/supplemental/plot_alt_fig06-07_transfer_and_surface.py

@@ -0,0 +1,101 @@
+import os
+import numpy as np
+import matplotlib as mpl 
+import matplotlib.pyplot as plt 
+import h5py
+from scipy import interpolate
+from pathlib import Path
+from scipy.interpolate import interp1d
+
+import mesa_reader as mr
+plt.rcParams['font.family'] = ['Times New Roman']
+plt.rcParams['mathtext.fontset'] = 'dejavusans'
+plt.rcParams['mathtext.fontset'] = 'cm'
+plt.rcParams['mathtext.rm'] = 'serif'
+
+wave_lum = lambda f, ell: (2.33e-11)*f**(-6.5)*np.sqrt(ell*(ell+1))**4 #f in Hz.
+correction = 0.4
+
+root_dir = '../../data/'
+output_file = '{}/dedalus/surface_signals/wave_propagation_power_spectra.h5'.format(root_dir)
+with h5py.File('../../dedalus/zams_15Msol_LMC/wave_propagation/nr128/star/star_128+192+64_bounds0-0.93R_Re4.00e+03_de1.5_cutoff1.0e-10.h5', 'r') as starf:
+    s_nd = starf['s_nd'][()]
+    L_nd = starf['L_nd'][()]
+    m_nd = starf['m_nd'][()]
+    tau_nd = starf['tau_nd'][()]
+    energy_nd = L_nd**2 * m_nd / (tau_nd**2)
+    lum_nd = energy_nd/tau_nd
+
+wave_lums_data = dict()
+with h5py.File('../../data/dedalus/wave_fluxes/zams_15Msol_LMC/re01200/wave_luminosities.h5', 'r') as lum_file:
+    radius_str = '1.25'
+    wave_lums_data['freqs'] = lum_file['cgs_freqs'][()]
+    wave_lums_data['ells'] = lum_file['ells'][()].ravel()
+    wave_lums_data['lum'] = lum_file['cgs_wave_luminosity(r={})'.format(radius_str)][0,:]
+
+
+for use_fit in [False, True]:
+    with h5py.File(output_file, 'r') as out_f:
+        freqs_hz = out_f['freqs'][()] / tau_nd  #Hz
+        freqs = freqs_hz * (24 * 60 * 60) #invday
+        ells = out_f['ells'][()].ravel()
+        s1 = np.sqrt(out_f['shell(s1_S2,r=R)'][0,:,:])*s_nd
+
+        fig = plt.figure(figsize=(7.5, 4)) 
+        ax1 = fig.add_axes([0.04, 0.50, 0.32, 0.45])
+        ax2 = fig.add_axes([0.36, 0.50, 0.32, 0.45])
+        ax3 = fig.add_axes([0.68, 0.50, 0.32, 0.45])
+        ax4 = fig.add_axes([0.04, 0.05, 0.32, 0.45])
+        ax5 = fig.add_axes([0.36, 0.05, 0.32, 0.45])
+        ax6 = fig.add_axes([0.68, 0.05, 0.32, 0.45])
+        axs = [ax1, ax2, ax3, ax4, ax5, ax6]
+        cmap = mpl.cm.plasma
+        Lmax = 6
+
+        for j in range(1,Lmax+1):
+            axs[j-1].loglog(freqs, s1[:,ells==j], color='k', lw=1, label='simulation')
+
+            with h5py.File('../../data/dedalus/transfer/transfer_ell{:03d}_eigenvalues.h5'.format(j), 'r') as ef:
+                #transfer dimensionality is (in simulation units) entropy / sqrt(wave luminosity).
+                transfer_func_root_lum = ef['transfer_root_lum'][()] * s_nd/np.sqrt(lum_nd)
+                transfer_freq_hz = ef['om'][()]/(2*np.pi) / tau_nd #Hz
+                transfer_freq = transfer_freq_hz * (24 * 60 * 60) #invday
+                transfer_interp = lambda f: 10**interp1d(np.log10(transfer_freq_hz), np.log10(transfer_func_root_lum), bounds_error=False, fill_value=-1e99)(np.log10(f))
+
+            if use_fit:
+                surface_s1_amplitude = correction*np.abs(transfer_interp(transfer_freq_hz))*np.sqrt(wave_lum(transfer_freq_hz, j))
+                axs[j-1].loglog(transfer_freq, surface_s1_amplitude, c='orange', label='transfer solution (Eqn. 74 $L_w$)', lw=0.5)
+            else:
+                log_wave_lum = interp1d(np.log10(wave_lums_data['freqs']), np.log10(wave_lums_data['lum'][:,j == wave_lums_data['ells']].ravel()))
+                data_wave_lum = lambda f: 10**(log_wave_lum(np.log10(f)))
+                surface_s1_amplitude = correction*np.abs(transfer_interp(transfer_freq_hz))*np.sqrt(data_wave_lum(transfer_freq_hz))
+                axs[j-1].loglog(transfer_freq, surface_s1_amplitude, c='orange', label='transfer solution (measured $L_w$)', lw=0.5)
+            axs[j-1].text(0.96, 0.92, r'$\ell =$ ' + '{}'.format(j), transform=axs[j-1].transAxes, ha='right', va='center')
+
+        for ax in axs:
+            ax.set_xlim(7e-2,9e0)
+            ax.set_ylim(1e-6, 1e2)
+            ax.set_yticks((1e-5, 1e-3, 1e-1, 1e1))
+        for ax in [ax4, ax5, ax6]:
+            ax.set_xlabel('frequency (d$^{-1}$)')
+        for ax in [ax1, ax4]:
+            ax.set_ylabel('$|s_1|$ (erg g$^{-1}$ K$^{-1}$)')
+        for ax in [ax1, ax2, ax3]:
+            #        ax.set_ylim(2e-5,3e-1)
+            ax.set_xticklabels(())
+        for ax in [ax2, ax3, ax5, ax6]:
+            ax.set_yticklabels(())
+            ax.tick_params(axis="y",direction="in", which='both')
+
+        for ax in [ax1, ax2, ax3]:
+            ax.tick_params(axis="x", direction="in", which='both')
+
+
+        ax1.legend(loc='lower left', frameon=False, borderaxespad=0.1, handletextpad=0.2, fontsize=8)
+        if use_fit:
+            plt.savefig('alt_fig07_wavepropagation_transferVerification_fit.png', bbox_inches='tight', dpi=300)
+            plt.savefig('alt_fig07_wavepropagation_transferVerification_fit.pdf', bbox_inches='tight', dpi=300)
+        else:
+            plt.savefig('alt_fig06_wavepropagation_transferVerification_raw.png', bbox_inches='tight', dpi=300)
+            plt.savefig('alt_fig06_wavepropagation_transferVerification_raw.pdf', bbox_inches='tight', dpi=300)
+        plt.clf()

figures/supplemental/plot_fig05_sbdf2_gamma_eff.py

@@ -12,6 +12,7 @@
 from dedalus.tools.general import natural_sort
 
 
+a_corr = 0.4
 
 files = natural_sort(glob.glob('../../data/dedalus/eigenvalues/dual*ell*.h5'))
 with h5py.File('../../dedalus/zams_15Msol_LMC/wave_propagation/nr256/star/star_256+192+64_bounds0-0.93R_Re1.00e+04_de1.5_cutoff1.0e-10.h5', 'r') as starf:
@@ -68,8 +69,8 @@
         plt.xscale('log')
 
         plt.axes(axbot)
-        plt.loglog(transfer_om/(2*np.pi), transfer, lw=2, c='k')
-        plt.loglog(raw_transfer_om/(2*np.pi), raw_transfer, lw=1, c='orange')
+        plt.loglog(transfer_om/(2*np.pi), a_corr*transfer, lw=2, c='k')
+        plt.loglog(raw_transfer_om/(2*np.pi), a_corr*raw_transfer, lw=1, c='orange')
         plt.xlabel(r'frequency (d$^{-1}$)')
 
         for ax in [axtop, axbot]:

figures/supplemental/plot_fig06-07_transfer_and_surface.py

@@ -14,7 +14,8 @@
 plt.rcParams['mathtext.rm'] = 'serif'
 
 wave_lum = lambda f, ell: (2.33e-11)*f**(-6.5)*np.sqrt(ell*(ell+1))**4 #f in Hz.
-fudge = 1
+correction = 0.4
+min_freqs = [0.17, 0.22, 0.3, 0.35, 0.4, 0.45] #invday
 
 root_dir = '../../data/'
 output_file = '{}/dedalus/surface_signals/wave_propagation_power_spectra.h5'.format(root_dir)
@@ -64,14 +65,17 @@
                 transfer_freq_hz = ef['om'][()]/(2*np.pi) / tau_nd #Hz
                 transfer_freq = transfer_freq_hz * (24 * 60 * 60) #invday
                 transfer_interp = lambda f: 10**interp1d(np.log10(transfer_freq_hz), np.log10(transfer_func_root_lum), bounds_error=False, fill_value=-1e99)(np.log10(f))
+                bad = transfer_freq < min_freqs[j-1]
 
             if use_fit:
-                surface_s1_amplitude = fudge*np.abs(transfer_interp(transfer_freq_hz))*np.sqrt(wave_lum(transfer_freq_hz, j))
+                surface_s1_amplitude = correction*np.abs(transfer_interp(transfer_freq_hz))*np.sqrt(wave_lum(transfer_freq_hz, j))
+                surface_s1_amplitude[bad] = 0 #set low freq to 0
                 axs[j-1].loglog(transfer_freq, surface_s1_amplitude, c='orange', label='transfer solution (Eqn. 74 $L_w$)', lw=0.5)
             else:
                 log_wave_lum = interp1d(np.log10(wave_lums_data['freqs']), np.log10(wave_lums_data['lum'][:,j == wave_lums_data['ells']].ravel()))
                 data_wave_lum = lambda f: 10**(log_wave_lum(np.log10(f)))
-                surface_s1_amplitude = fudge*np.abs(transfer_interp(transfer_freq_hz))*np.sqrt(data_wave_lum(transfer_freq_hz))
+                surface_s1_amplitude = correction*np.abs(transfer_interp(transfer_freq_hz))*np.sqrt(data_wave_lum(transfer_freq_hz))
+                surface_s1_amplitude[bad] = 0 #set low freq to 0
                 axs[j-1].loglog(transfer_freq, surface_s1_amplitude, c='orange', label='transfer solution (measured $L_w$)', lw=0.5)
             axs[j-1].text(0.96, 0.92, r'$\ell =$ ' + '{}'.format(j), transform=axs[j-1].transAxes, ha='right', va='center')
 

figures/supplemental/plot_fig08_transfer_functions.py

@@ -80,8 +80,8 @@
 for ax in ax_rows[-1]:
     ax.set_xlabel('$f$ (d$^{-1}$)')
 
-ax2_3.text(0.25, 0.1, 'star', ha='left',   va='center', transform=ax2_3.transAxes)
-ax2_3.text(0.25, 0.6, 'WP sim', ha='left', va='center',  c='orange', transform=ax2_3.transAxes)
+ax2_3.text(0.25, 0.2, 'star', ha='left',   va='center', transform=ax2_3.transAxes)
+ax2_3.text(0.25, 0.65, 'WP sim', ha='left', va='center',  c='orange', transform=ax2_3.transAxes)
 
 
 plt.savefig('fig08_transfer_functions.png', bbox_inches='tight', dpi=300)

figures/supplemental/plot_fig10_rednoise_fits.py

@@ -21,11 +21,11 @@ def red_noise(nu, alpha0, nu_char, gamma=2):
 
 star_dirs = ['03msol_Zsolar', '15msol_ZLMC', '40msol_Zsolar']
 Lmax = [15, 15, 15]
-alpha0 = [9e-3, 1.2e-1, 3e-1]
-alpha_latex = [ r'$9 \times 10^{-3}$ $\mu$mag',\
-                r'$0.12$ $\mu$mag',\
-                r'$0.3$ $\mu$mag']
-nu_char = [2.6e-1, 0.16, 0.105]
+alpha0 = [5.5e-3, 6e-2, 1.6e-1]
+alpha_latex = [ r'$5.5 \times 10^{-3}$ $\mu$mag',\
+                r'$6 \times 10^{-2}$ $\mu$mag',\
+                r'$0.16$ $\mu$mag']
+nu_char = [3e-1, 0.22, 0.13]
 gamma  = [4.5, 3.9, 4.3]
 out_f = h5py.File(output_file, 'r')
 freqs = out_f['frequencies'][()]*24*60*60
@@ -35,6 +35,9 @@ def red_noise(nu, alpha0, nu_char, gamma=2):
 ax1 = fig.add_axes([0.050, 0.025, 0.275, 0.95])
 ax2 = fig.add_axes([0.375, 0.025, 0.275, 0.95])
 ax3 = fig.add_axes([0.700, 0.025, 0.275, 0.95])
+ax1.text(0.98, 0.98, '3 $M_{\odot}$', ha='right', va='top', transform=ax1.transAxes)
+ax2.text(0.98, 0.98, '15 $M_{\odot}$', ha='right', va='top', transform=ax2.transAxes)
+ax3.text(0.98, 0.98, '40 $M_{\odot}$', ha='right', va='top', transform=ax3.transAxes)
 axs = [ax1, ax2, ax3]
 
 for i, sdir in enumerate(star_dirs):

figures/supplemental/plot_fig11_rotation_figure.py

@@ -356,8 +356,8 @@ def find_tams(center_h1,model):
 
 #Fiducial rotating star
 from palettable.colorbrewer.qualitative import Dark2_5 as cmap
-ax3.scatter(21.7, 0.4, c=sm2.to_rgba(fiducial_T), marker='*', edgecolors='k', linewidth=1, s=200, zorder=100)
-ax4.scatter(fiducial_Rop, 0.4, c=sm2.to_rgba(fiducial_T), marker='*', edgecolors='k', linewidth=1, s=200, zorder=100)
+ax3.scatter(21.7, 0.17, c=sm2.to_rgba(fiducial_T), marker='*', edgecolors='k', linewidth=1, s=200, zorder=100)
+ax4.scatter(fiducial_Rop, 0.17, c=sm2.to_rgba(fiducial_T), marker='*', edgecolors='k', linewidth=1, s=200, zorder=100)
 
 plt.savefig("fig11_rednoise_Ro.pdf",bbox_inches='tight', dpi=300)
 plt.savefig("fig11_rednoise_Ro.png",bbox_inches='tight', dpi=300)

figures/supplemental/plot_fig13_rot_wave_luminosity.py

@@ -135,9 +135,9 @@ def red_noise(nu, alpha0, nu_char, gamma=2):
 out_f = h5py.File(output_file, 'r')
 freqs = out_f['frequencies'][()] * 24 * 60 * 60
 
-alpha = 0.4e0 
-alpha_latex = '0.4 $\mu$mag'
-nu_char = 0.16
+alpha = 0.21e0 
+alpha_latex = '0.21 $\mu$mag'
+nu_char = 0.22
 gamma = 3.9
 
 for i, sdir in enumerate(star_dirs):