diff --git a/GYRE/03msol_Zsolar/transfer_gyre.py b/GYRE/03msol_Zsolar/transfer_gyre.py new file mode 100644 index 0000000..d2315e5 --- /dev/null +++ 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 new file mode 100644 index 0000000..eccbd5b --- /dev/null +++ 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 new file mode 100644 index 0000000..dfea312 --- /dev/null +++ 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 new file mode 100644 index 0000000..7f9c140 --- /dev/null +++ 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()