adds files for making gyre transfer functions

evanhanders · Mar 18, 2023 · 9d5d65d · 9d5d65d
1 parent 7dda3ef
commit 9d5d65d
Show file tree

Hide file tree

Showing 4 changed files with 234 additions and 0 deletions.
diff --git a/GYRE/03msol_Zsolar/transfer_gyre.py b/GYRE/03msol_Zsolar/transfer_gyre.py
@@ -0,0 +1,53 @@
+"""
+This file reads in gyre eigenfunctions, calculates the velocity and velocity dual basis, and outputs in a clean format so that it's ready to be fed into the transfer function calculation.
+"""
+import os
+import numpy as np
+import pygyre as pg
+
+from compstar.tools.mesa import find_core_cz_radius
+from compstar.waves.clean_gyre_eig import GyreMSGPostProcessor, solar_z
+
+plot = True
+use_delta_L = False
+Lmin = 1
+Lmax = 16
+ell_list = np.arange(Lmin, Lmax+1)
+folder = 'gyre_output'
+for ell in ell_list:
+    om_list = np.logspace(-8, -2, 1000) #Hz * 2pi
+
+    mesa_LOG   = '../../MESA/03msol_Zsolar/LOGS/profile43.data'
+    pulse_file = '{}.GYRE'.format(mesa_LOG)
+    mode_base = './gyre_output/mode_id{:05d}_ell{:03d}_m+00_n{:06d}.h5'
+    files = []
+
+    max_cond = 1e12
+    summary_file='gyre_output/summary_ell{:02d}.txt'.format(ell)
+    summary = pg.read_output(summary_file)
+
+    #sort eigenvalues by 1/freq
+    sorting = np.argsort(summary['freq'].real**(-1))
+    summary = summary[sorting]
+
+    good_freqs = []
+    for row in summary:
+        this_ell = row['l']
+        this_id = row['id']
+        n_pg = row['n_pg']
+        if complex(row['freq']).real < 0:
+            continue
+        if n_pg >= 0: continue
+        if this_ell != ell: continue
+        files.append(mode_base.format(this_id, ell, n_pg))
+        good_freqs.append(complex(row['freq']))
+
+    post = GyreMSGPostProcessor(ell, summary_file, files, pulse_file, mesa_LOG,
+                  filters=['Red',], initial_z=solar_z,
+                  MSG_DIR = os.environ['MSG_DIR'],
+                  GRID_DIR=os.path.join('..','..','data','MSG','specgrid'),
+                  PASS_DIR=os.path.join('..','..','data','MSG','passbands'))
+    post.sort_eigenfunctions()
+    data_dicts = post.evaluate_magnitudes()
+    data_dict = post.calculate_duals(max_cond=max_cond)
+    post.calculate_transfer(plot=plot, use_delta_L=use_delta_L, N_om=3000)
diff --git a/GYRE/15msol_ZLMC/transfer_gyre.py b/GYRE/15msol_ZLMC/transfer_gyre.py
@@ -0,0 +1,53 @@
+"""
+This file reads in gyre eigenfunctions, calculates the velocity and velocity dual basis, and outputs in a clean format so that it's ready to be fed into the transfer function calculation.
+"""
+import os
+import numpy as np
+import pygyre as pg
+
+from compstar.tools.mesa import find_core_cz_radius
+from compstar.waves.clean_gyre_eig import GyreMSGPostProcessor, solar_z
+
+plot = True
+use_delta_L = False
+Lmin = 1
+Lmax = 16
+ell_list = np.arange(Lmin, Lmax+1)
+folder = 'gyre_output'
+for ell in ell_list:
+    om_list = np.logspace(-8, -2, 1000) #Hz * 2pi
+
+    mesa_LOG   = '../../MESA/15msol_ZLMC/LOGS/profile47.data'
+    pulse_file = '{}.GYRE'.format(mesa_LOG)
+    mode_base = './gyre_output/mode_id{:05d}_ell{:03d}_m+00_n{:06d}.h5'
+    files = []
+
+    max_cond = 1e12
+    summary_file='gyre_output/summary_ell{:02d}.txt'.format(ell)
+    summary = pg.read_output(summary_file)
+
+    #sort eigenvalues by 1/freq
+    sorting = np.argsort(summary['freq'].real**(-1))
+    summary = summary[sorting]
+
+    good_freqs = []
+    for row in summary:
+        this_ell = row['l']
+        this_id = row['id']
+        n_pg = row['n_pg']
+        if complex(row['freq']).real < 0:
+            continue
+        if n_pg >= 0: continue
+        if this_ell != ell: continue
+        files.append(mode_base.format(this_id, ell, n_pg))
+        good_freqs.append(complex(row['freq']))
+
+    post = GyreMSGPostProcessor(ell, summary_file, files, pulse_file, mesa_LOG,
+                  filters=['Red',], initial_z=0.006, specgrid='OSTAR2002',
+                  MSG_DIR = os.environ['MSG_DIR'],
+                  GRID_DIR=os.path.join('..','..','data','MSG','specgrid'),
+                  PASS_DIR=os.path.join('..','..','data','MSG','passbands'))
+    post.sort_eigenfunctions()
+    data_dicts = post.evaluate_magnitudes()
+    data_dict = post.calculate_duals(max_cond=max_cond)
+    post.calculate_transfer(plot=plot, use_delta_L=use_delta_L, N_om=3000)
diff --git a/GYRE/40msol_Zsolar/transfer_gyre.py b/GYRE/40msol_Zsolar/transfer_gyre.py
@@ -0,0 +1,53 @@
+"""
+This file reads in gyre eigenfunctions, calculates the velocity and velocity dual basis, and outputs in a clean format so that it's ready to be fed into the transfer function calculation.
+"""
+import os
+import numpy as np
+import pygyre as pg
+
+from compstar.tools.mesa import find_core_cz_radius
+from compstar.waves.clean_gyre_eig import GyreMSGPostProcessor, solar_z
+
+plot = True
+use_delta_L = False
+Lmin = 1
+Lmax = 16
+ell_list = np.arange(Lmin, Lmax+1)
+folder = 'gyre_output'
+for ell in ell_list:
+    om_list = np.logspace(-8, -2, 1000) #Hz * 2pi
+
+    mesa_LOG   = '../../MESA/40msol_Zsolar/LOGS/profile53.data'
+    pulse_file = '{}.GYRE'.format(mesa_LOG)
+    mode_base = './gyre_output/mode_id{:05d}_ell{:03d}_m+00_n{:06d}.h5'
+    files = []
+
+    max_cond = 1e12
+    summary_file='gyre_output/summary_ell{:02d}.txt'.format(ell)
+    summary = pg.read_output(summary_file)
+
+    #sort eigenvalues by 1/freq
+    sorting = np.argsort(summary['freq'].real**(-1))
+    summary = summary[sorting]
+
+    good_freqs = []
+    for row in summary:
+        this_ell = row['l']
+        this_id = row['id']
+        n_pg = row['n_pg']
+        if complex(row['freq']).real < 0:
+            continue
+        if n_pg >= 0: continue
+        if this_ell != ell: continue
+        files.append(mode_base.format(this_id, ell, n_pg))
+        good_freqs.append(complex(row['freq']))
+
+    post = GyreMSGPostProcessor(ell, summary_file, files, pulse_file, mesa_LOG,
+                  filters=['Red',], initial_z=solar_z, specgrid='OSTAR2002',
+                  MSG_DIR = os.environ['MSG_DIR'],
+                  GRID_DIR=os.path.join('..','..','data','MSG','specgrid'),
+                  PASS_DIR=os.path.join('..','..','data','MSG','passbands'))
+    post.sort_eigenfunctions()
+    data_dicts = post.evaluate_magnitudes()
+    data_dict = post.calculate_duals(max_cond=max_cond)
+    post.calculate_transfer(plot=plot, use_delta_L=use_delta_L, N_om=3000)
diff --git a/GYRE/generate_magnitude_spectra.py b/GYRE/generate_magnitude_spectra.py
@@ -0,0 +1,75 @@
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+import h5py
+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
+eig_dir = 'gyre_output'
+output_file = 'magnitude_spectra.h5'
+
+star_dirs = ['03msol_Zsolar', '40msol_Zsolar', '15msol_ZLMC']
+luminosity_amplitudes = [7.34e-15, 5.3e-10, 2.33e-11]
+Lmax = [15, 15, 15]
+obs_length_days = 365
+obs_length_sec  = obs_length_days*24*60*60
+obs_cadence = 30*60 #30 min
+df = 1/obs_length_sec
+N_data = int(obs_length_sec/obs_cadence)
+freqs = np.arange(N_data)*df
+
+out_f = h5py.File(output_file, 'w')
+out_f['frequencies'] = freqs
+
+signals = []
+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
+    transfer_oms = []
+    transfer_signal = []
+    pure_transfers  = []
+
+    ell_list = np.arange(1, Lmax[i]+1)
+    for ell in ell_list:
+        print(sdir, " ell = %i" % ell)
+
+
+        with h5py.File('{:s}/{:s}/ell{:03d}_eigenvalues.h5'.format(sdir, eig_dir, ell), 'r') as f:
+            om = f['transfer_om'][()]
+            transfer_root_lum = f['transfer_root_lum'][()].real
+        micromag = transfer_root_lum*np.sqrt(np.abs(wave_luminosity(om/(2*np.pi), ell)))
+        transfer_oms.append(om[np.isfinite(micromag)])
+        transfer_signal.append(micromag[np.isfinite(micromag)])
+        pure_transfers.append(transfer_root_lum[np.isfinite(micromag)])
+
+
+
+    transfers  = np.zeros((Lmax[i], N_data))
+    magnitudes = np.zeros((Lmax[i], N_data))
+    for j in range(freqs.size-1):
+        for k, oms, signal, transfer in zip(range(Lmax[i]), transfer_oms, transfer_signal, pure_transfers):
+            good = (2*np.pi*freqs[j+1] >= oms)*(2*np.pi*freqs[j] < oms)
+            if np.sum(good) > 0:
+                magnitudes[k,j] = np.max(signal[good])
+                transfers[k,j] = np.max(transfer[good])
+            elif 2*np.pi*freqs[j] > oms.min() and 2*np.pi*freqs[j] <= oms.max():
+                magnitudes[k,j] = signal[np.argmin(np.abs(2*np.pi*freqs[j] - oms))]
+                transfers[k,j]  = transfer[np.argmin(np.abs(2*np.pi*freqs[j] - oms))]
+
+    total_signal = np.sqrt(np.sum(magnitudes**2, axis=0))
+    out_f['{}_magnitude_cube'.format(sdir)] = magnitudes
+    out_f['{}_transfer_cube'.format(sdir)] = transfers
+    out_f['{}_magnitude_sum'.format(sdir)] = total_signal
+    plt.loglog(freqs, total_signal, c='k')
+    plt.legend()
+    plt.xlim(1e-7, 1e-3)
+    plt.ylim(1e-6, 1)
+    plt.xlabel('freq (Hz)')
+    plt.ylabel(r'$\Delta m\,(\mu\rm{mag})$')
+    plt.savefig('obs_ell_contributions_{}.png'.format(sdir), bbox_inches='tight')
+    plt.clf()
+out_f.close()