diff --git a/src/scope/input.txt b/src/scope/input.txt index 487a52c..3dde4dd 100644 --- a/src/scope/input.txt +++ b/src/scope/input.txt @@ -19,10 +19,11 @@ seed 150 # seed for the random number generator log_level DEBUG # logging level above which statements are printed. Supported levels are ``DEBUG``, ``INFO``, ``WARNING``, ``ERROR``, and ``CRITICAL``. # Filepaths -planet_spectrum_path /home/asavel/scratch/retrieval_stuff/fm_wasp_77ab_noisy/fm_wasp_77ab_noisy/spectrum_tran_WASP-166b_0.3_c_to_o.pic -star_spectrum_path /home/asavel/scratch/scope/src/scope/data/PHOENIX_5605_4.33.txt -data_cube_path /home/asavel/scratch/scope/src/scope/data/data_RAW_20201214_wl_algn_03.pic -snr_path /home/asavel/scratch/scope/src/scope/data/snr_IGRINS_LTT-9779b.json # only required for CRIRES+ simulations +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 bfe9125..3b140ae 100644 --- a/src/scope/run_simulation.py +++ b/src/scope/run_simulation.py @@ -4,7 +4,6 @@ calculating the log likelihood and cross-correlation function of the data given the model parameters. """ - import logging import os @@ -26,14 +25,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, @@ -58,6 +57,7 @@ def make_data( out_of_transit_dur=0.1, v_sys_measured=0.0, vary_throughput=True, + combined_spectrum=np.nan, instrument="IGRINS", pca_removal="subtract", ): @@ -106,27 +106,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") @@ -443,6 +447,7 @@ def simulate_observation( seed=42, LD=True, vary_throughput=True, + combined_planet_star_path=np.nan, instrument="IGRINS", snr_path=None, planet_name="yourfirstplanet", @@ -501,8 +506,9 @@ def simulate_observation( # need to read in the filepath n_order, n_pixel, wl_cube_model, SNR = read_crires_data(snr_path) - wl_model, Fp, Fstar = np.load(planet_spectrum_path, allow_pickle=True) + # 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) @@ -512,21 +518,36 @@ def simulate_observation( # instrument profile convolution instrument_kernel = get_instrument_kernel(instrument) Fp_conv = np.convolve(Fp_conv_rot, instrument_kernel, mode="same") - print("here spectrum") - print(Fp_conv) - 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 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, + ) + + 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) - # todo: swap out - Fstar_conv = get_star_spline( - star_wave, star_flux, wl_model, instrument_kernel, smooth=False - ) # fill with NaNs because we're checkpointing. 0 would be a valid value, NaNs indicate that we haven't calculated it yet. lls, ccfs = np.zeros((n_kp, n_vsys)) * np.nan, np.zeros((n_kp, n_vsys)) * np.nan @@ -536,14 +557,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, @@ -562,6 +583,7 @@ def simulate_observation( divide_out_of_transit=divide_out_of_transit, out_of_transit_dur=0.1, v_sys_measured=v_sys, + combined_spectrum=True, LD=LD, instrument=instrument, pca_removal=pca_removal,