adds figs 06-10 of supplemental which had been missed?

.gitignore

@@ -0,0 +1,4 @@
+*.png
+*.pdf
+*.data
+*.data.GYRE

figures/supplemental/plot_fig06-07_transfer_and_surface.py

@@ -0,0 +1,106 @@
+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.
+fudge = 1
+
+root_dir = '../../data/'
+output_file = '{}/dedalus/surface_signals/wave_propagation_power_spectra.h5'.format(root_dir)
+with h5py.File('../../dedalus/zams_15Msol_LMC/wave_generation/re10000/star/star_512+192_bounds0-2L_Re3.00e+04_de1.5_cutoff1.0e-10.h5', 'r') as lumf:
+    s_nd = lumf['s_nd'][()]
+
+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:
+    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/re03200/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
+        norm = mpl.colors.Normalize(vmin=1, vmax=Lmax)
+        sm = mpl.cm.ScalarMappable(norm=norm, cmap=mpl.colors.ListedColormap(Dark2_5.mpl_colors[:Lmax]))
+
+        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 = fudge*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 = fudge*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('fig07_wavepropagation_transferVerification_fit.png', bbox_inches='tight', dpi=300)
+            plt.savefig('fig07_wavepropagation_transferVerification_fit.pdf', bbox_inches='tight', dpi=300)
+        else:
+            plt.savefig('fig06_wavepropagation_transferVerification_raw.png', bbox_inches='tight', dpi=300)
+            plt.savefig('fig06_wavepropagation_transferVerification_raw.pdf', bbox_inches='tight', dpi=300)
+        plt.clf()

figures/supplemental/plot_fig08_transfer_functions.py

@@ -0,0 +1,90 @@
+"""
+Calculate transfer function to get surface response of convective forcing.
+Outputs a function which, when multiplied by sqrt(wave flux), gives you the surface response.
+"""
+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
+
+plt.rcParams['font.family'] = ['Times New Roman']
+plt.rcParams['mathtext.fontset'] = 'dejavusans'
+plt.rcParams['mathtext.fontset'] = 'cm'
+plt.rcParams['mathtext.rm'] = 'serif'
+
+
+#Calculate transfer functions
+output_file = '../../data/dedalus/predictions/magnitude_spectra.h5'
+star_dirs = ['3msol', '15msol', '40msol']
+Lmax = [3,3,3]
+out_f = h5py.File(output_file, 'r')
+freqs = out_f['frequencies'][()] * 24 * 60 * 60 #invday
+
+fig = plt.figure(figsize=(7.5, 5))
+ax1_1 = fig.add_axes([0.05, 0.66, 0.28, 0.33])
+ax2_1 = fig.add_axes([0.385, 0.66, 0.28, 0.33])
+ax3_1 = fig.add_axes([0.72, 0.66, 0.28, 0.33])
+ax1_2 = fig.add_axes([0.05, 0.33, 0.28, 0.33])
+ax2_2 = fig.add_axes([0.385, 0.33, 0.28, 0.33])
+ax3_2 = fig.add_axes([0.72, 0.33, 0.28, 0.33])
+ax1_3 = fig.add_axes([0.05, 0.00, 0.28, 0.33])
+ax2_3 = fig.add_axes([0.385, 0.00, 0.28, 0.33])
+ax3_3 = fig.add_axes([0.72, 0.00, 0.28, 0.33])
+#axs = [ax1, ax2, ax3, ax4, ax5, ax6]
+ax_rows = [[ax1_1, ax2_1, ax3_1], [ax1_2, ax2_2, ax3_2], [ax1_3, ax2_3, ax3_3]]
+
+
+for i, sdir in enumerate(star_dirs):
+    transfer = out_f['{}_transfer_cube'.format(sdir)][()]
+
+    for j in range(Lmax[i]):
+        ax = ax_rows[j][i]
+        ax.loglog(freqs, transfer[j,:], label='star', c='k', lw=0.75)
+        ax.text(0.03, 0.93, r'$M = $ ' + '{}'.format(int(sdir.split('msol')[0])) + r'$M_{\odot}$', transform=ax.transAxes, ha='left', va='center', c='k')
+        ax.text(0.03, 0.85, r'$\ell = $ ' + '{}'.format(j+1), transform=ax.transAxes, ha='left', va='center', c='k')
+        ax.set_xlim(5e-2, 1e1)
+        ax.set_ylim(1e-21, 1e-9)
+        ax.set_yticks((1e-20, 1e-17, 1e-14, 1e-11))
+
+        if i == 0:
+            ax.set_ylabel('Transfer')
+
+#Plot dedalus transfer
+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:
+    s_nd = starf['s_nd'][()]
+    Cp = starf['Cp'][()]*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
+
+for ell in range(3):
+    with h5py.File('../../data/dedalus/transfer/transfer_ell{:03d}_eigenvalues.h5'.format(ell+1), 'r') as ef:
+        #transfer dimensionality is (in simulation units) entropy / sqrt(wave luminosity).
+        transfer_func_root_lum = 1e6*ef['transfer_root_lum'][()] * s_nd/np.sqrt(lum_nd) / Cp #turns sqrt(L) -> 1e6 * s/cp
+        transfer_freq_hz = ef['om'][()]/(2*np.pi) / tau_nd #Hz
+        transfer_freq = transfer_freq_hz * (24 * 60 * 60) #invday
+    ax_rows[ell][1].loglog(transfer_freq, transfer_func_root_lum, label='WP simulation', c='orange', lw=0.75)
+    ax_rows[ell][1].set_yticks((1e-20, 1e-17, 1e-14, 1e-11))
+
+
+for row in ax_rows[:2]:
+    for ax in row:
+        ax.set_xticklabels(())
+        ax.tick_params(axis='x', direction='in', which='both')
+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)
+
+
+plt.savefig('fig08_transfer_functions.png', bbox_inches='tight', dpi=300)
+plt.savefig('fig08_transfer_functions.pdf', bbox_inches='tight', dpi=300)
+
+out_f.close()

figures/supplemental/plot_fig09_ellsum_spectra.py

@@ -0,0 +1,94 @@
+"""
+Calculate transfer function to get surface response of convective forcing.
+Outputs a function which, when multiplied by sqrt(wave flux), gives you the surface response.
+"""
+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
+plt.rcParams['font.family'] = ['Times New Roman']
+plt.rcParams['mathtext.fontset'] = 'dejavusans'
+plt.rcParams['mathtext.fontset'] = 'cm'
+plt.rcParams['mathtext.rm'] = 'serif'
+
+
+#Calculate transfer functions
+output_file = '../../data/dedalus/predictions/magnitude_spectra.h5'
+
+
+star_dirs = ['3msol', '15msol', '40msol']
+Lmax = [15, 15, 15]
+out_f = h5py.File(output_file, 'r')
+freqs = out_f['frequencies'][()] * 24 * 60 * 60
+
+fig = plt.figure(figsize=(7.5, 7))
+ax1 = fig.add_axes([0.00 , 0.60, 0.45, 0.30])#fig.add_subplot(1,2,1)
+ax2 = fig.add_axes([0.575, 0.60, 0.45, 0.30])#fig.add_subplot(1,2,2)
+ax3 = fig.add_axes([0.00 , 0.30, 0.45, 0.30])#fig.add_subplot(1,2,1)
+ax4 = fig.add_axes([0.575, 0.30, 0.45, 0.30])#fig.add_subplot(1,2,2)
+ax5 = fig.add_axes([0.00 , 0.00, 0.45, 0.30])#fig.add_subplot(1,2,1)
+ax6 = fig.add_axes([0.575, 0.00, 0.45, 0.30])#fig.add_subplot(1,2,2)
+cax1 = fig.add_axes([0.650, 0.975, 0.30, 0.025])
+cax2 = fig.add_axes([0.075, 0.975, 0.30, 0.025])
+axs = [ax1, ax2, ax3, ax4, ax5, ax6]
+ax_pairs = [[ax1, ax2], [ax3, ax4], [ax5, ax6]]
+
+even_cmap = mpl.cm.Purples_r
+odds_cmap = mpl.cm.Greens_r
+even_norm = mpl.colors.Normalize(vmin=2, vmax=24)
+odds_norm = mpl.colors.Normalize(vmin=1, vmax=23)
+even_sm = mpl.cm.ScalarMappable(norm=even_norm, cmap=even_cmap)
+odds_sm = mpl.cm.ScalarMappable(norm=odds_norm, cmap=odds_cmap)
+cb1 = plt.colorbar(even_sm, cax=cax1, orientation='horizontal', boundaries=[1, 3, 5, 7, 9, 11, 13, 15], ticks=[2, 4, 6, 8, 10, 12, 14, 16])
+cb2 = plt.colorbar(odds_sm, cax=cax2, orientation='horizontal', ticks=[1, 3, 5, 7, 9, 11, 13, 15], boundaries=[0, 2, 4, 6, 8, 10, 12, 14, 16])
+cb1.set_label(r'$\ell$ (even)')
+cb2.set_label(r'$\ell$ (odd)')
+
+
+
+
+for i, sdir in enumerate(star_dirs):
+    magnitude_cube = out_f['{}_magnitude_cube'.format(sdir)][()]
+    ell_list = np.arange(1, Lmax[i]+1)
+    axl, axr = ax_pairs[i]
+
+    for j in range(Lmax[i]):
+        sum_mag = np.sqrt(np.sum(magnitude_cube[:Lmax[i]-j,:]**2, axis=0))
+        if (j+1) % 2 == 0: 
+            axl.loglog(freqs, magnitude_cube[j,:], color=even_sm.to_rgba(j+1), lw=0.5, zorder=2)
+            axr.loglog(freqs, sum_mag, color=even_sm.to_rgba(Lmax[i]-j), lw=0.5, zorder=2)
+        else:
+            axl.loglog(freqs, magnitude_cube[j,:], color=odds_sm.to_rgba(j+1), lw=0.5, zorder=2)
+            axr.loglog(freqs, sum_mag, color=odds_sm.to_rgba(Lmax[i]-j), lw=0.5, zorder=2)
+    axr.loglog(freqs, np.sqrt(np.sum(magnitude_cube[:Lmax[i],:]**2, axis=0)), c='k', lw=0.25, zorder=2)
+    axl.text(0.04, 0.93, r'$M = $ ' + '{}'.format(int(sdir.split('msol')[0])) + r'$M_{\odot}$', transform=axl.transAxes, ha='left', va='center')
+    axr.text(0.96, 0.93, r'$M = $ ' + '{}'.format(int(sdir.split('msol')[0])) + r'$M_{\odot}$', transform=axr.transAxes, ha='right', va='center')
+
+for ax in axs:
+    ax.set_xlim(3e-2, 3e1)
+
+for ax in [ax1, ax3, ax5]:
+    ax.set_ylabel(r'$\Delta m_{\ell}\,(\mu\rm{mag})$')
+for ax in [ax2, ax4, ax6]:
+    ax.set_ylabel(r'$\sqrt{\sum_{i}^{\ell}\Delta m_{i}^2}\,(\mu\rm{mag})$')
+
+for ax in [ax1, ax2]:
+    ax.set_ylim(1e-7, 1e0)
+    ax.set_xticklabels(())
+    ax.tick_params(axis="x", direction="in", which='both', zorder=5, top=True, bottom=False)
+for ax in [ax3, ax4]:
+    ax.set_ylim(1e-6, 3e-1)
+    ax.set_xticklabels(())
+    ax.tick_params(axis="x", direction="in", which='both', zorder=5, top=True, bottom=False)
+for ax in [ax5, ax6]:
+    ax.set_ylim(3e-6, 3e0)
+    ax.set_xlabel('frequency (d$^{-1}$)')
+plt.savefig('fig09_ellsums.png', bbox_inches='tight', dpi=300)
+plt.savefig('fig09_ellsums.pdf', bbox_inches='tight', dpi=300)
+#    plt.clf()
+
+out_f.close()

figures/supplemental/plot_fig10_rednoise_fits.py

@@ -0,0 +1,68 @@
+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
+plt.rcParams['font.family'] = ['Times New Roman']
+plt.rcParams['mathtext.fontset'] = 'dejavusans'
+plt.rcParams['mathtext.fontset'] = 'cm'
+plt.rcParams['mathtext.rm'] = 'serif'
+
+
+#Calculate transfer functions
+output_file = '../../data/dedalus/predictions/magnitude_spectra.h5'
+
+def red_noise(nu, alpha0, nu_char, gamma=2):
+    return alpha0/(1 + (nu/nu_char)**gamma)
+
+
+star_dirs = ['3msol', '15msol', '40msol']
+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]
+gamma  = [4.5, 3.9, 4.3]
+out_f = h5py.File(output_file, 'r')
+freqs = out_f['frequencies'][()]*24*60*60
+
+
+fig = plt.figure(figsize=(7.5, 2.5))
+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])
+axs = [ax1, ax2, ax3]
+
+for i, sdir in enumerate(star_dirs):
+    plt.axes(axs[i])
+    magnitude_sum = out_f['{}_magnitude_sum'.format(sdir)][()]
+
+    alpha_lbl = r'$\alpha_0 = $' + alpha_latex[i]
+    nu_lbl    = r'$\nu_{\rm char} = $' + '{:.2f}'.format(nu_char[i]) + ' d$^{-1}$'
+    gamma_lbl = r'$\gamma = $' + '{:.1f}'.format(gamma[i])
+    label     = '{}\n{}\n{}'.format(alpha_lbl, nu_lbl, gamma_lbl)
+    plt.loglog(freqs, magnitude_sum, label=sdir, c='k')
+    plt.loglog(freqs, red_noise(freqs, alpha0[i], nu_char[i], gamma=gamma[i]), label=label, c='orange')
+    if i == 0: #3
+        plt.ylim(1e-5, 8e-2)
+        plt.xlim(5e-2, 5e0)
+    if i == 1: #15
+        plt.ylim(1e-4, 1)
+        plt.xlim(3e-2, 2e0)
+    if i == 2: #40
+        plt.ylim(1e-3, 2e0)
+        plt.xlim(2e-2, 2e0)
+    axs[i].text(0.01, 0.99, label, ha='left', va='top', transform=axs[i].transAxes)
+
+#    plt.legend()
+    plt.xlabel('frequency (d$^{-1}$)')
+    if i == 0:
+        plt.ylabel(r'$\Delta m$ ($\mu$mag)')
+fig.savefig('fig10_rednoise_fit.png', bbox_inches='tight', dpi=300)
+fig.savefig('fig10_rednoise_fit.pdf', bbox_inches='tight')
+
+out_f.close()