From 6fa4cf27f0c90c78538a97fd250ef6d00f32a1ad Mon Sep 17 00:00:00 2001 From: arjunsavel Date: Mon, 4 Nov 2024 14:49:57 -0500 Subject: [PATCH] add combined spectrum technique --- src/scope/input.txt | 1 + src/scope/run_simulation.py | 96 +++++++++++++++++++++++-------------- 2 files changed, 60 insertions(+), 37 deletions(-) diff --git a/src/scope/input.txt b/src/scope/input.txt index 7bc6b88..182debd 100644 --- a/src/scope/input.txt +++ b/src/scope/input.txt @@ -20,6 +20,7 @@ seed 42 # seed for the random number generator w planet_spectrum_path /Users/arjunsavel/Desktop/research/gj1214_hires_sims/best_fit_spectrum_tran_gj1214_steam_nir.pic star_spectrum_path /Users/arjunsavel/Desktop/research/scope/src/scope/data/PHOENIX_5605_4.33.txt data_cube_path /Users/arjunsavel/Desktop/research/scope/src/scope/data/data_RAW_20201214_wl_algn_03.pic +planet_and_star_path NULL # path to *combined* planet and star spectrum, if applicable. If not set to NULL, the planet and stellar spectrum are combined from the individual files specified above. # Astrophysical Parameters. Can specify DATABASE to pull based on name for parameters with [DB] Rp DATABASE # planetary radius in Jupiter radii. [DB] diff --git a/src/scope/run_simulation.py b/src/scope/run_simulation.py index ad57c4b..9a7dff1 100644 --- a/src/scope/run_simulation.py +++ b/src/scope/run_simulation.py @@ -22,14 +22,14 @@ def make_data( scale, wlgrid, wl_model, - Fp_conv, - Fstar_conv, n_order, n_exposure, n_pixel, phases, Rp_solar, Rstar, + Fp_conv=np.ones(3), + Fstar_conv=np.ones(3), seed=42, do_pca=True, blaze=False, @@ -54,6 +54,7 @@ def make_data( out_of_transit_dur=0.1, v_sys_measured=0.0, vary_throughput=True, + combined_spectrum=np.nan, ): """ Creates a simulated HRCCS dataset. Main function. @@ -96,27 +97,31 @@ def make_data( pca_noise_matrix = np.zeros((n_order, n_exposure, n_pixel)) for order in range(n_order): wlgrid_order = np.copy(wlgrid[order,]) # Cropped wavelengths - flux_cube[order] = doppler_shift_planet_star( - flux_cube[order], - n_exposure, - phases, - rv_planet, - rv_star, - wlgrid_order, - wl_model, - Fp_conv, - Rp_solar, - Fstar_conv, - Rstar, - u1, - u2, - a, - b, - LD, - scale, - star, - observation, - ) + + if combined_spectrum is not np.nan: + flux_cube[order] = np.interp(wlgrid_order, wl_model, combined_spectrum) + else: + flux_cube[order] = doppler_shift_planet_star( + flux_cube[order], + n_exposure, + phases, + rv_planet, + rv_star, + wlgrid_order, + wl_model, + Fp_conv, + Rp_solar, + Fstar_conv, + Rstar, + u1, + u2, + a, + b, + LD, + scale, + star, + observation, + ) throughput_baselines = np.loadtxt(abs_path + "/data/throughputs.txt") @@ -424,6 +429,7 @@ def simulate_observation( inc=90.0, seed=42, vary_throughput=True, + combined_planet_star_path=np.nan, **kwargs, ): """ @@ -464,12 +470,13 @@ def simulate_observation( Kp_array = np.linspace(kp - 100, kp + 100, 200) v_sys_array = np.arange(v_sys - 100, v_sys + 100) n_order, n_pixel = (44, 1848) # todo: generalize. - mike_wave, mike_cube = pickle.load(open(data_cube_path, "rb"), encoding="latin1") + # read in the data cube + mike_wave, mike_cube = pickle.load(open(data_cube_path, "rb"), encoding="latin1") wl_cube_model = mike_wave.copy().astype(np.float64) + # read in the planet spectrum and stellar spectrum wl_model, Fp, Fstar = np.load(planet_spectrum_path, allow_pickle=True) - wl_model = wl_model.astype(np.float64) # Fp_conv_rot = broaden_spectrum(wl_model / 1e6, Fp, 0, vl=v_rot) @@ -480,18 +487,32 @@ def simulate_observation( instrument_kernel = get_instrument_kernel() Fp_conv = np.convolve(Fp_conv_rot, instrument_kernel, mode="same") - star_wave, star_flux = np.loadtxt( - star_spectrum_path - ).T # Phoenix stellar model packing + if combined_planet_star_path is np.nan: + star_wave, star_flux = np.loadtxt( + star_spectrum_path + ).T # Phoenix stellar model packing - if include_rm: - star_flux, _ = make_stellar_disk( - star_flux, star_wave, v_rot_star, phases, Rstar, inc, lambda_misalign, a, Rp - ) + if include_rm: + star_flux, _ = make_stellar_disk( + star_flux, + star_wave, + v_rot_star, + phases, + Rstar, + inc, + lambda_misalign, + a, + Rp, + ) - Fstar_conv = get_star_spline( - star_wave, star_flux, wl_model, instrument_kernel, smooth=False - ) + Fstar_conv = get_star_spline( + star_wave, star_flux, wl_model, instrument_kernel, smooth=False + ) + else: + wl_model, combined_spectrum = np.load( + combined_planet_star_path, allow_pickle=True + ) + Fstar_conv = np.ones_like(Fp_conv) lls, ccfs = np.zeros((200, 200)), np.zeros((200, 200)) @@ -500,14 +521,14 @@ def simulate_observation( scale, wl_cube_model, wl_model, - Fp_conv, - Fstar_conv, n_order, n_exposures, n_pixel, phases, Rp_solar, Rstar, + Fp_conv=Fp_conv, + Fstar_conv=Fstar_conv, seed=seed, do_pca=True, blaze=blaze, @@ -526,6 +547,7 @@ def simulate_observation( divide_out_of_transit=False, out_of_transit_dur=0.1, v_sys_measured=v_sys, + combined_spectrum=True, ) run_name = f"{n_princ_comp}_NPC_{blaze}_blaze_{star}_star_{telluric}_telluric_{SNR}_SNR_{tell_type}_{time_dep_tell}_{wav_error}_{order_dep_throughput}"