From d52ac08abb38ebe4e6e9d624a7699abecc875c5f Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Wed, 26 Feb 2025 13:21:34 -0500
Subject: [PATCH 01/22] attempt crires+ sim

---
 src/scope/input.txt         |  4 ++-
 src/scope/input_output.py   |  7 ++++-
 src/scope/noise.py          | 21 +++++++++++++++
 src/scope/run_simulation.py | 19 +++++++++++---
 src/scope/utils.py          | 51 ++++++++++++++++++++++++++++++++-----
 5 files changed, 89 insertions(+), 13 deletions(-)

diff --git a/src/scope/input.txt b/src/scope/input.txt
index a758926..3d462d6 100644
--- a/src/scope/input.txt
+++ b/src/scope/input.txt
@@ -21,6 +21,7 @@ seed                   42               # seed for the random number generator w
 planet_spectrum_path    /home/asavel/scratch/retrieval_stuff/fm_wasp_77ab_noisy/fm_wasp_77ab_noisy/spectrum_emission_LTT-9779b.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
 
 # Astrophysical Parameters. Can specify DATABASE to pull based on name for parameters with [DB]
 Rp                     DATABASE          # planetary radius in Jupiter radii. [DB]
@@ -40,6 +41,7 @@ lambda_misalign        0                 # misalignment angle of the planet's or
 inc                    90.0              # inclination of the planet's orbit with respect to the line of sight, in degrees. only matters if include_rm is set to True. and observation is set to transmission.
 
 # Instrument Parameters
+spectrograph           IGRINS        # name of the spectrograph to simulate. currently spectrographs are ``IGRINS`` and ``CRIRES+``.
 blaze                  True          # whether to include a blaze function or not.
 wav_error              False         # whether to include wavelength solution errors or not.
 order_dep_throughput   True          # whether to include order-dependent throughput variations.
@@ -53,7 +55,7 @@ phase_end              0.45           # phase of the end of the observations. 0
 n_exposures            10             # number of exposures to simulate. sets the phases of the exposures.
 star                   True          # whether to include the star in the simulation. In general, you'd like to!
 telluric               True          # whether to include tellurics in the simulation. In general, you'd like to!
-SNR                    300           # signal-to-noise ratio of the observations, per pixel. I.e., what sets the photon counts at the detector.
+SNR                    300           # signal-to-noise ratio of the observations, per pixel. I.e., what sets the photon counts at the detector. Ignored for CRIRES+ simulations.
 tell_type              data-driven   # type of telluric simulation. supported modes are ``ATRAN`` and ``data-driven``.
 time_dep_tell          False         # whether the tellurics are time-dependent or not.
 
diff --git a/src/scope/input_output.py b/src/scope/input_output.py
index ea527fc..1a3da61 100644
--- a/src/scope/input_output.py
+++ b/src/scope/input_output.py
@@ -228,7 +228,12 @@ def write_input_file(data, output_file_path="input.txt"):
     # Define the order and categories of parameters
 
     categories = {
-        "Filepaths": ["planet_spectrum_path", "star_spectrum_path", "data_cube_path"],
+        "Filepaths": [
+            "planet_spectrum_path",
+            "star_spectrum_path",
+            "data_cube_path",
+            "snr_path",
+        ],
         "Astrophysical Parameters": ["Rp", "Rp_solar", "Rstar", "kp", "v_rot", "v_sys"],
         "Instrument Parameters": ["SNR"],
         "Observation Parameters": [
diff --git a/src/scope/noise.py b/src/scope/noise.py
index 1bb0101..760b62e 100644
--- a/src/scope/noise.py
+++ b/src/scope/noise.py
@@ -112,6 +112,26 @@ def add_igrins_noise(flux_cube_model, wl_grid, SNR):
     return add_quadratic_noise(flux_cube_model, wl_grid, SNR, IGRINS=True)
 
 
+def add_crires_plus_noise(flux_cube_model, wl_grid, SNR):
+    """
+    Adds CRIRES+ noise to the flux cube model.
+    """
+    noisy_flux = np.zeros_like(flux_cube_model)
+    for exposure in range(flux_cube_model.shape[1]):
+        for order in range(flux_cube_model.shape[0]):
+            flux_level = (
+                SNR**2
+                * flux_cube_model[order][exposure]
+                / np.nanmax(flux_cube_model[order][exposure])
+            )
+            noisy_flux[order][exposure] = np.random.poisson(flux_level)
+
+            # a few quick checks to make sure that nothing has gone wrong with adding noise
+            noisy_flux[order][exposure][noisy_flux[order][exposure] < 0.0] = 0.0
+            noisy_flux[order][exposure][~np.isfinite(noisy_flux[order][exposure])] = 0.0
+    return noisy_flux
+
+
 def add_custom_noise(SNR):
     """
     Adds custom noise to the flux cube model.
@@ -142,6 +162,7 @@ def add_noise_cube(flux_cube_model, wl_grid, SNR, noise_model="constant", **kwar
     noise_models = {
         "constant": add_constant_noise,
         "IGRINS": add_igrins_noise,
+        "CRIRES+": add_crires_plus_noise,
         "custom_quadratic": add_quadratic_noise,
         "custom": add_custom_noise,
     }
diff --git a/src/scope/run_simulation.py b/src/scope/run_simulation.py
index ffe2943..915963a 100644
--- a/src/scope/run_simulation.py
+++ b/src/scope/run_simulation.py
@@ -55,6 +55,7 @@ def make_data(
     out_of_transit_dur=0.1,
     v_sys_measured=0.0,
     vary_throughput=True,
+    instrument="IGRINS",
 ):
     """
     Creates a simulated HRCCS dataset. Main function.
@@ -168,7 +169,7 @@ def make_data(
     flux_cube[np.isnan(flux_cube)] = 0.0
     if SNR > 0:  # 0 means don't add noise!
         if order_dep_throughput:
-            noise_model = "IGRINS"
+            noise_model = instrument
         else:
             noise_model = "constant"
         flux_cube = add_noise_cube(flux_cube, wlgrid, SNR, noise_model=noise_model)
@@ -420,6 +421,8 @@ def simulate_observation(
     inc=90.0,
     seed=42,
     vary_throughput=True,
+    instrument="IGRINS",
+    snr_path=None,
     **kwargs,
 ):
     """
@@ -459,10 +462,17 @@ def simulate_observation(
     Rp_solar = Rp * rjup_rsun  # convert from jupiter radii to solar radii
     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")
 
-    wl_cube_model = mike_wave.copy().astype(np.float64)
+    if instrument == "IGRINS":
+        n_order, n_pixel = (44, 1848)
+        mike_wave, mike_cube = pickle.load(
+            open(data_cube_path, "rb"), encoding="latin1"
+        )
+
+        wl_cube_model = mike_wave.copy().astype(np.float64)
+    elif instrument == "CRIRES+":
+        # 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)
 
@@ -523,6 +533,7 @@ def simulate_observation(
         divide_out_of_transit=False,
         out_of_transit_dur=0.1,
         v_sys_measured=v_sys,
+        instrument=instrument,
     )
 
     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}"
diff --git a/src/scope/utils.py b/src/scope/utils.py
index 0ae2b30..41e0e29 100644
--- a/src/scope/utils.py
+++ b/src/scope/utils.py
@@ -2,6 +2,7 @@
 Utility functions for simulating HRCCS data.
 """
 
+import json
 import os
 import pickle
 
@@ -24,6 +25,44 @@
 end_clip = 100
 
 
+def read_crires_data(data_path):
+    """
+    Reads in CRIRES data.
+
+    Inputs
+    ------
+        :data_path: (str) path to the data
+
+    Outputs
+    -------
+        n_orders: (int) number of orders
+        n_pixel: (int) number of pixels
+        wl_cube_model: (array) wavelength cube model
+        snrs: (array) signal-to-noise ratios
+    """
+    with open(data_path, "r") as file:
+        data = json.load(file)
+
+    n_orders = 0  # an integer :)
+    for i in range(len(data["data"]["orders"])):
+        order_len = len(data["data"]["orders"][i]["detectors"])
+        n_orders += order_len
+
+    n_wavs = len(data["data"]["orders"][i]["detectors"][0]["wavelength"])
+
+    wl_grid = np.zeros((n_orders, n_wavs))
+    snr_grid = np.zeros((n_orders, n_wavs))
+
+    for i in range(len(data["data"]["orders"])):
+        order_len = len(data["data"]["orders"][i]["detectors"])
+        for j in range(order_len):
+            wl_grid[i * order_len + j] = data["data"]["orders"][i]["detectors"][j][
+                "wavelength"
+            ]
+
+    return n_orders, n_wavs, wl_grid, snr_grid
+
+
 @njit
 def doppler_shift_planet_star(
     model_flux_cube,
@@ -63,7 +102,6 @@ def doppler_shift_planet_star(
             if not reprocessing:
                 model_flux_cube[exposure,] *= flux_star * Rstar**2
         elif observation == "emission":
-
             model_flux_cube[exposure,] = flux_planet * (Rp_solar * rsun) ** 2
         elif (
             observation == "transmission"
@@ -410,14 +448,13 @@ def calc_rvs(v_sys, v_sys_measured, Kp, Kstar, phases):
 
     """
     v_sys_tot = v_sys + v_sys_measured  # total offset
-    rv_planet = (
-        v_sys_tot + Kp * np.sin(2.0 * np.pi * phases)
+    rv_planet = v_sys_tot + Kp * np.sin(
+        2.0 * np.pi * phases
     )  # input in km/s, convert to m/s
 
-
-    rv_star = (
-        v_sys_measured - Kstar * np.sin(2.0 * np.pi * phases)
-    )   # measured in m/s. note opposite sign!
+    rv_star = v_sys_measured - Kstar * np.sin(
+        2.0 * np.pi * phases
+    )  # measured in m/s. note opposite sign!
     return rv_planet * 1e3, rv_star * 1e3
 
 

From a6415972222c4513bde30119a08ee84244484bbf Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Wed, 26 Feb 2025 13:24:25 -0500
Subject: [PATCH 02/22] ah fix the version thing

---
 pyproject.toml        |  1 +
 src/scope/_version.py | 17 +++++++++++++++++
 2 files changed, 18 insertions(+)
 create mode 100644 src/scope/_version.py

diff --git a/pyproject.toml b/pyproject.toml
index c4ca0ad..5dc41a9 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -6,6 +6,7 @@ build-backend = "setuptools.build_meta"
 name = "scope-astr"
 authors = [
     {name = "Arjun Savel", email = "asavel@umd.edu"},
+    {name = "Arjun Savel"},
     {name= "Megan Bedell"},
     {name= "Eliza M.-R. Kempton"},
     {name= "Peter Smith"},
diff --git a/src/scope/_version.py b/src/scope/_version.py
new file mode 100644
index 0000000..4556cfd
--- /dev/null
+++ b/src/scope/_version.py
@@ -0,0 +1,17 @@
+# file generated by setuptools_scm
+# don't change, don't track in version control
+TYPE_CHECKING = False
+if TYPE_CHECKING:
+    from typing import Tuple, Union
+
+    VERSION_TUPLE = Tuple[Union[int, str], ...]
+else:
+    VERSION_TUPLE = object
+
+version: str
+__version__: str
+__version_tuple__: VERSION_TUPLE
+version_tuple: VERSION_TUPLE
+
+__version__ = version = "0.4.2b0"
+__version_tuple__ = version_tuple = (0, 4, 2)

From ad0b408202f397de269c57b886b7bd972b7bad33 Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Wed, 26 Feb 2025 13:34:36 -0500
Subject: [PATCH 03/22] correctly do doppler shift now

---
 src/scope/tests/test_utils.py | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/src/scope/tests/test_utils.py b/src/scope/tests/test_utils.py
index 5522db4..dd2fbbf 100644
--- a/src/scope/tests/test_utils.py
+++ b/src/scope/tests/test_utils.py
@@ -131,7 +131,9 @@ def test_same_after_shift(self):
         interped_flux = np.interp(eval_wave, template_wave, template_flux)
         v = 1e-6  # m/s. should not change much
 
-        shifted_flux = calc_doppler_shift(eval_wave, template_wave, template_flux, v)
+        shifted_wave, shifted_flux = calc_doppler_shift(
+            eval_wave, template_wave, template_flux, v
+        )
         # test that shifted flux and interpred flux are very similar
         assert (
             np.testing.assert_allclose(shifted_flux, interped_flux, rtol=1e-2) == None
@@ -148,7 +150,9 @@ def test_diff_after_shift(self):
         interped_flux = np.interp(eval_wave, template_wave, template_flux)
         v = 2e3  # m/s. should not change much
 
-        shifted_flux = calc_doppler_shift(eval_wave, template_wave, template_flux, v)
+        shifted_wave, shifted_flux = calc_doppler_shift(
+            eval_wave, template_wave, template_flux, v
+        )
         # test that shifted flux and interpred flux are very similar
         np.testing.assert_raises(
             AssertionError, np.testing.assert_array_equal, shifted_flux, interped_flux
@@ -165,7 +169,9 @@ def test_gaussian_shifted_blue(self):
         interped_flux = np.interp(eval_wave, template_wave, template_flux)
 
         v = -5e4  # m/s. should not change much
-        shifted_flux = calc_doppler_shift(eval_wave, template_wave, template_flux, v)
+        shifted_wave, shifted_flux = calc_doppler_shift(
+            eval_wave, template_wave, template_flux, v
+        )
         wav_max_shifted = eval_wave[np.argmax(shifted_flux)]
         wav_max = eval_wave[np.argmax(interped_flux)]
         assert wav_max_shifted < wav_max

From c09d540c5808c43a259f6c6f0a9720f9c68288ce Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Wed, 26 Feb 2025 15:45:14 -0500
Subject: [PATCH 04/22] missed a test!

---
 src/scope/tests/test_utils.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/scope/tests/test_utils.py b/src/scope/tests/test_utils.py
index dd2fbbf..ac6d9a7 100644
--- a/src/scope/tests/test_utils.py
+++ b/src/scope/tests/test_utils.py
@@ -186,7 +186,9 @@ def test_gaussian_shifted_red(self):
         template_flux = np.exp(-0.5 * ((template_wave - 1.5e-6) / 0.01) ** 2)
 
         v = 5e4  # m/s. should not change much
-        shifted_flux = calc_doppler_shift(eval_wave, template_wave, template_flux, v)
+        shifted_wave, shifted_flux = calc_doppler_shift(
+            eval_wave, template_wave, template_flux, v
+        )
         interped_flux = np.interp(eval_wave, template_wave, template_flux)
         wav_max_shifted = eval_wave[np.argmax(shifted_flux)]
         wav_max = eval_wave[np.argmax(interped_flux)]

From e1a571ded2193a8f86e47158ada990a387fe5e21 Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Wed, 26 Feb 2025 16:43:27 -0500
Subject: [PATCH 05/22] try a crires test

---
 src/scope/tests/test_run_simulation.py | 49 ++++++++++++++++++++++++++
 1 file changed, 49 insertions(+)

diff --git a/src/scope/tests/test_run_simulation.py b/src/scope/tests/test_run_simulation.py
index 12c4eef..65dc328 100644
--- a/src/scope/tests/test_run_simulation.py
+++ b/src/scope/tests/test_run_simulation.py
@@ -516,3 +516,52 @@ def test_nopca_fluxlevels_deterministic(
 
     # the flux levels must be higher when there's no tellurics.
     np.testing.assert_array_equal(flux_cube_nopca, flux_cube_nopca_b)
+
+
+def test_crires_simulation(test_inputs):
+    data_cube_path = os.path.join(test_data_path, "data_RAW_20201214_wl_algn_03.pic")
+    planet_spectrum_path = os.path.join(
+        test_data_path, "best_fit_spectrum_tran_gj1214_steam_nir.pic"
+    )
+    snr_path = os.path.join(test_data_path, "output1.json")
+    star_spectrum_path = os.path.join(test_data_path, "PHOENIX_5605_4.33.txt")
+    simulate_observation(
+        planet_spectrum_path=planet_spectrum_path,
+        star_spectrum_path=star_spectrum_path,
+        data_cube_path=data_cube_path,
+        phase_start=0,
+        phase_end=1,
+        n_exposures=2,
+        observation="emission",
+        blaze=True,
+        n_princ_comp=4,
+        star=True,
+        SNR=250,
+        telluric=False,
+        tell_type="data-driven",
+        time_dep_tell=False,
+        wav_error=False,
+        rv_semiamp_orbit=0.3229,
+        order_dep_throughput=True,
+        Rp=1.21,  # Jupiter radii,
+        Rstar=0.955,  # solar radii
+        kp=192.02,  # planetary orbital velocity, km/s
+        v_rot=4.5,
+        scale=1.0,
+        v_sys=0.0,
+        modelname="yourfirstsimulation",
+        divide_out_of_transit=False,
+        out_of_transit_dur=0.1,
+        include_rm=False,
+        v_rot_star=3.0,
+        a=0.033,  #
+        lambda_misalign=0.0,
+        inc=90.0,
+        seed=42,
+        vary_throughput=True,
+        snr_path=snr_path,
+    )
+
+    # check that there are nonzeros, basically that it runs
+    ccfs = np.loadtxt("yourfirstsimulation_ccfs.txt")
+    assert np.sum(ccfs) > 0

From 6a0a297f38ecbed206fa271a9d938b44d83ff154 Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Wed, 26 Feb 2025 16:53:54 -0500
Subject: [PATCH 06/22] actually make intermediate directories too!

---
 src/scope/utils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/scope/utils.py b/src/scope/utils.py
index 41e0e29..fbb4422 100644
--- a/src/scope/utils.py
+++ b/src/scope/utils.py
@@ -126,7 +126,7 @@ def save_results(outdir, run_name, lls, ccfs):
 
 def make_outdir(outdir):
     try:
-        os.mkdir(outdir)
+        os.makedirs(outdir)
     except FileExistsError:
         print("Directory already exists. Continuing!")
 

From f6b093fd69aaa609ebcefff5464aeb833906cab4 Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Wed, 26 Feb 2025 17:53:02 -0500
Subject: [PATCH 07/22] write planet name defaulty

---
 src/scope/run_simulation.py            | 2 +-
 src/scope/tests/test_run_simulation.py | 1 +
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/scope/run_simulation.py b/src/scope/run_simulation.py
index 915963a..bf79e1d 100644
--- a/src/scope/run_simulation.py
+++ b/src/scope/run_simulation.py
@@ -423,7 +423,7 @@ def simulate_observation(
     vary_throughput=True,
     instrument="IGRINS",
     snr_path=None,
-    **kwargs,
+    planet_name="yourfirstplanet" ** kwargs,
 ):
     """
     Run a simulation of the data, given a grid index and some paths. Side effects:
diff --git a/src/scope/tests/test_run_simulation.py b/src/scope/tests/test_run_simulation.py
index 65dc328..dc0a503 100644
--- a/src/scope/tests/test_run_simulation.py
+++ b/src/scope/tests/test_run_simulation.py
@@ -550,6 +550,7 @@ def test_crires_simulation(test_inputs):
         scale=1.0,
         v_sys=0.0,
         modelname="yourfirstsimulation",
+        planet_name="GJ1214b",
         divide_out_of_transit=False,
         out_of_transit_dur=0.1,
         include_rm=False,

From 7ac644f2e71f5ddc7d88beed39f5d5dce23aa3ef Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Wed, 26 Feb 2025 17:56:16 -0500
Subject: [PATCH 08/22] forgot a comma!

---
 src/scope/run_simulation.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/scope/run_simulation.py b/src/scope/run_simulation.py
index bf79e1d..451de84 100644
--- a/src/scope/run_simulation.py
+++ b/src/scope/run_simulation.py
@@ -423,7 +423,8 @@ def simulate_observation(
     vary_throughput=True,
     instrument="IGRINS",
     snr_path=None,
-    planet_name="yourfirstplanet" ** kwargs,
+    planet_name="yourfirstplanet",
+    **kwargs,
 ):
     """
     Run a simulation of the data, given a grid index and some paths. Side effects:

From 238f33103084ebb7d6841ce60ae17486e8d43dc9 Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Wed, 26 Feb 2025 17:59:59 -0500
Subject: [PATCH 09/22] sim obs does not need vary throughput

---
 src/scope/run_simulation.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/scope/run_simulation.py b/src/scope/run_simulation.py
index 451de84..aba0ad8 100644
--- a/src/scope/run_simulation.py
+++ b/src/scope/run_simulation.py
@@ -567,7 +567,6 @@ def simulate_observation(
                 star=star,
                 observation=observation,
                 v_sys_measured=v_sys,
-                vary_throughput=vary_throughput,
             )
             lls[l, k], ccfs[l, k] = res
 

From fc183423688f653b9250bf6d25b819dd2e4557d1 Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Fri, 28 Feb 2025 14:29:38 -0500
Subject: [PATCH 10/22] grab correct file

---
 src/scope/tests/test_run_simulation.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/scope/tests/test_run_simulation.py b/src/scope/tests/test_run_simulation.py
index dc0a503..8ade506 100644
--- a/src/scope/tests/test_run_simulation.py
+++ b/src/scope/tests/test_run_simulation.py
@@ -564,5 +564,6 @@ def test_crires_simulation(test_inputs):
     )
 
     # check that there are nonzeros, basically that it runs
-    ccfs = np.loadtxt("yourfirstsimulation_ccfs.txt")
+    filetype = glob("output/yourfirstsimulation/*ccf*")
+    ccfs = np.loadtxt(filetype[0])
     assert np.sum(ccfs) > 0

From d431f62c51f9580da701349b15112840dc21615e Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Sun, 2 Mar 2025 16:08:12 -0500
Subject: [PATCH 11/22] actually test crires+

---
 src/scope/tests/test_run_simulation.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/scope/tests/test_run_simulation.py b/src/scope/tests/test_run_simulation.py
index 8ade506..ffac2be 100644
--- a/src/scope/tests/test_run_simulation.py
+++ b/src/scope/tests/test_run_simulation.py
@@ -554,6 +554,7 @@ def test_crires_simulation(test_inputs):
         divide_out_of_transit=False,
         out_of_transit_dur=0.1,
         include_rm=False,
+        instrument="CRIRES+",
         v_rot_star=3.0,
         a=0.033,  #
         lambda_misalign=0.0,
@@ -564,6 +565,6 @@ def test_crires_simulation(test_inputs):
     )
 
     # check that there are nonzeros, basically that it runs
-    filetype = glob("output/yourfirstsimulation/*ccf*")
+    filetype = glob("src/scope/output/yourfirstsimulation/*ccf*")
     ccfs = np.loadtxt(filetype[0])
     assert np.sum(ccfs) > 0

From 57d5ddfda9612776f291f017d2286b4f9fd7437d Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Sun, 2 Mar 2025 16:10:02 -0500
Subject: [PATCH 12/22] only test a few kp and vsys values!

---
 src/scope/run_simulation.py            | 8 +++++---
 src/scope/tests/test_run_simulation.py | 2 ++
 2 files changed, 7 insertions(+), 3 deletions(-)

diff --git a/src/scope/run_simulation.py b/src/scope/run_simulation.py
index aba0ad8..ebc97f8 100644
--- a/src/scope/run_simulation.py
+++ b/src/scope/run_simulation.py
@@ -424,6 +424,8 @@ def simulate_observation(
     instrument="IGRINS",
     snr_path=None,
     planet_name="yourfirstplanet",
+    n_kp=200,
+    n_vsys=200,
     **kwargs,
 ):
     """
@@ -461,8 +463,8 @@ def simulate_observation(
 
     phases = np.linspace(phase_start, phase_end, n_exposures)
     Rp_solar = Rp * rjup_rsun  # convert from jupiter radii to solar radii
-    Kp_array = np.linspace(kp - 100, kp + 100, 200)
-    v_sys_array = np.arange(v_sys - 100, v_sys + 100)
+    Kp_array = np.linspace(kp - 100, kp + 100, n_kp)
+    v_sys_array = np.linspace(v_sys - 100, v_sys + 100, n_vsys)
 
     if instrument == "IGRINS":
         n_order, n_pixel = (44, 1848)
@@ -501,7 +503,7 @@ def simulate_observation(
         star_wave, star_flux, wl_model, instrument_kernel, smooth=False
     )
 
-    lls, ccfs = np.zeros((200, 200)), np.zeros((200, 200))
+    lls, ccfs = np.zeros((n_kp, n_vsys)), np.zeros((n_kp, n_vsys))
 
     # redoing the grid. how close does PCA get to a tellurics-free signal detection?
     A_noplanet, flux_cube, flux_cube_nopca, just_tellurics = make_data(
diff --git a/src/scope/tests/test_run_simulation.py b/src/scope/tests/test_run_simulation.py
index ffac2be..2ef862b 100644
--- a/src/scope/tests/test_run_simulation.py
+++ b/src/scope/tests/test_run_simulation.py
@@ -559,6 +559,8 @@ def test_crires_simulation(test_inputs):
         a=0.033,  #
         lambda_misalign=0.0,
         inc=90.0,
+        n_kp=20,
+        n_vsys=20,
         seed=42,
         vary_throughput=True,
         snr_path=snr_path,

From 819be3dae8361e21700df3bf4bdf19983be63882 Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Sun, 2 Mar 2025 16:18:42 -0500
Subject: [PATCH 13/22] print what the noise model is

---
 src/scope/run_simulation.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/scope/run_simulation.py b/src/scope/run_simulation.py
index ebc97f8..a2df7b2 100644
--- a/src/scope/run_simulation.py
+++ b/src/scope/run_simulation.py
@@ -172,6 +172,7 @@ def make_data(
             noise_model = instrument
         else:
             noise_model = "constant"
+        print(f"Adding noise with model {noise_model}")
         flux_cube = add_noise_cube(flux_cube, wlgrid, SNR, noise_model=noise_model)
 
     flux_cube = detrend_cube(flux_cube, n_order, n_exposure)

From d1ce8391b8cfa0bb621f81a58f7a367078d5ca7c Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Sun, 2 Mar 2025 16:47:45 -0500
Subject: [PATCH 14/22] do not include blaze function for crires sims

---
 src/scope/run_simulation.py            | 8 ++++++--
 src/scope/tests/test_run_simulation.py | 2 +-
 src/scope/utils.py                     | 8 +++++++-
 3 files changed, 14 insertions(+), 4 deletions(-)

diff --git a/src/scope/run_simulation.py b/src/scope/run_simulation.py
index a2df7b2..2c242fd 100644
--- a/src/scope/run_simulation.py
+++ b/src/scope/run_simulation.py
@@ -172,7 +172,7 @@ def make_data(
             noise_model = instrument
         else:
             noise_model = "constant"
-        print(f"Adding noise with model {noise_model}")
+        # print(f"Adding noise with model {noise_model}")
         flux_cube = add_noise_cube(flux_cube, wlgrid, SNR, noise_model=noise_model)
 
     flux_cube = detrend_cube(flux_cube, n_order, n_exposure)
@@ -220,7 +220,11 @@ def make_data(
 
         # add blaze
         out_of_transit_flux = add_blaze_function(
-            wlgrid, out_of_transit_flux, n_order, n_exposures_baseline
+            wlgrid,
+            out_of_transit_flux,
+            n_order,
+            n_exposures_baseline,
+            instrument=instrument,
         )
         out_of_transit_flux[flux_cube < 0.0] = 0.0
         out_of_transit_flux[np.isnan(flux_cube)] = 0.0
diff --git a/src/scope/tests/test_run_simulation.py b/src/scope/tests/test_run_simulation.py
index 2ef862b..e99a86b 100644
--- a/src/scope/tests/test_run_simulation.py
+++ b/src/scope/tests/test_run_simulation.py
@@ -533,7 +533,7 @@ def test_crires_simulation(test_inputs):
         phase_end=1,
         n_exposures=2,
         observation="emission",
-        blaze=True,
+        blaze=False,
         n_princ_comp=4,
         star=True,
         SNR=250,
diff --git a/src/scope/utils.py b/src/scope/utils.py
index fbb4422..8142081 100644
--- a/src/scope/utils.py
+++ b/src/scope/utils.py
@@ -541,7 +541,9 @@ def change_wavelength_solution(wl_cube_model, flux_cube_model, doppler_shifts):
     return flux_cube_model
 
 
-def add_blaze_function(wl_cube_model, flux_cube_model, n_order, n_exp):
+def add_blaze_function(
+    wl_cube_model, flux_cube_model, n_order, n_exp, instrument="IGRINS"
+):
     """
     Adds the blaze function to the model.
 
@@ -557,6 +559,10 @@ def add_blaze_function(wl_cube_model, flux_cube_model, n_order, n_exp):
         :flux_cube_model: (array) flux cube model with blaze function included.
     """
     # read in...have to somehow match the telluric spectra
+    if instrument != "IGRINS":
+        raise NotImplementedError(
+            "Only the IGRINS blaze function is currently supported."
+        )
 
     with open(abs_path + "/data/K_blaze_spectra.pic", "rb") as f:
         K_blaze_cube = pickle.load(f)

From 31dd8182272671a80addef40ce3c60f945ae63fa Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Sun, 2 Mar 2025 16:52:04 -0500
Subject: [PATCH 15/22] SNR greater than 0?

---
 src/scope/run_simulation.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/scope/run_simulation.py b/src/scope/run_simulation.py
index 2c242fd..9024782 100644
--- a/src/scope/run_simulation.py
+++ b/src/scope/run_simulation.py
@@ -167,7 +167,7 @@ def make_data(
 
     flux_cube = detrend_cube(flux_cube, n_order, n_exposure)
     flux_cube[np.isnan(flux_cube)] = 0.0
-    if SNR > 0:  # 0 means don't add noise!
+    if np.any(SNR > 0):  # 0 means don't add noise!
         if order_dep_throughput:
             noise_model = instrument
         else:

From 647d3a7b1a61725f83770de2758414e7b95997d8 Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Sun, 2 Mar 2025 16:57:28 -0500
Subject: [PATCH 16/22] a little debug

---
 src/scope/run_simulation.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/scope/run_simulation.py b/src/scope/run_simulation.py
index 9024782..f85776d 100644
--- a/src/scope/run_simulation.py
+++ b/src/scope/run_simulation.py
@@ -240,7 +240,8 @@ def make_data(
 
         # divide out the flux cube
         flux_cube /= median_out_of_transit  # todo: check axes work out
-
+    print("mean flux cube", np.mean(flux_cube))
+    print("std flux cube", np.std(flux_cube))
     if do_pca:
         for j in range(n_order):
             flux_cube[j] -= np.mean(flux_cube[j])

From 511be6540494a4cf1742d35f198fc402b63ff7cc Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Mon, 3 Mar 2025 09:47:03 -0500
Subject: [PATCH 17/22] test igrins

---
 src/scope/tests/test_run_simulation.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/scope/tests/test_run_simulation.py b/src/scope/tests/test_run_simulation.py
index e99a86b..2b58f3e 100644
--- a/src/scope/tests/test_run_simulation.py
+++ b/src/scope/tests/test_run_simulation.py
@@ -526,6 +526,7 @@ def test_crires_simulation(test_inputs):
     snr_path = os.path.join(test_data_path, "output1.json")
     star_spectrum_path = os.path.join(test_data_path, "PHOENIX_5605_4.33.txt")
     simulate_observation(
+        # problem: this is all 1s!
         planet_spectrum_path=planet_spectrum_path,
         star_spectrum_path=star_spectrum_path,
         data_cube_path=data_cube_path,
@@ -554,7 +555,7 @@ def test_crires_simulation(test_inputs):
         divide_out_of_transit=False,
         out_of_transit_dur=0.1,
         include_rm=False,
-        instrument="CRIRES+",
+        instrument="IGRINS",
         v_rot_star=3.0,
         a=0.033,  #
         lambda_misalign=0.0,

From a6523b29ea653d2108112e7ccc595b8fe211f4f6 Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Mon, 3 Mar 2025 10:36:54 -0500
Subject: [PATCH 18/22] test igrins with blaze?

---
 src/scope/tests/test_run_simulation.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/scope/tests/test_run_simulation.py b/src/scope/tests/test_run_simulation.py
index 2b58f3e..a324021 100644
--- a/src/scope/tests/test_run_simulation.py
+++ b/src/scope/tests/test_run_simulation.py
@@ -534,7 +534,7 @@ def test_crires_simulation(test_inputs):
         phase_end=1,
         n_exposures=2,
         observation="emission",
-        blaze=False,
+        blaze=True,
         n_princ_comp=4,
         star=True,
         SNR=250,

From d99f9250eb74f5a654c99f3c8638461fa248b1f3 Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Mon, 3 Mar 2025 13:26:08 -0500
Subject: [PATCH 19/22] more debug stuff

---
 src/scope/run_simulation.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/scope/run_simulation.py b/src/scope/run_simulation.py
index f85776d..4a06781 100644
--- a/src/scope/run_simulation.py
+++ b/src/scope/run_simulation.py
@@ -261,6 +261,8 @@ def make_data(
         pca_noise_matrix = np.ones_like(pca_noise_matrix)
     if tellurics:
         return pca_noise_matrix, flux_cube, flux_cube_nopca, just_tellurics
+    print("mean flux cube", np.mean(flux_cube))
+    print("std flux cube", np.std(flux_cube))
     return (
         pca_noise_matrix,
         flux_cube,

From d13b5e0c2f0c2da89acb537af086abb97b168c12 Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Mon, 3 Mar 2025 16:16:21 -0500
Subject: [PATCH 20/22] maybe need some tellurics

---
 src/scope/tests/test_run_simulation.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/scope/tests/test_run_simulation.py b/src/scope/tests/test_run_simulation.py
index a324021..25c53bc 100644
--- a/src/scope/tests/test_run_simulation.py
+++ b/src/scope/tests/test_run_simulation.py
@@ -538,7 +538,7 @@ def test_crires_simulation(test_inputs):
         n_princ_comp=4,
         star=True,
         SNR=250,
-        telluric=False,
+        telluric=True,
         tell_type="data-driven",
         time_dep_tell=False,
         wav_error=False,

From 1cf638473e79e5e5304978a70a29f315f840f642 Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Mon, 3 Mar 2025 16:55:14 -0500
Subject: [PATCH 21/22] fix run simulation tests!

---
 src/scope/run_simulation.py            | 17 ++++---
 src/scope/tests/test_run_simulation.py | 67 +++++++++++++++++++++++---
 src/scope/utils.py                     |  2 +-
 3 files changed, 70 insertions(+), 16 deletions(-)

diff --git a/src/scope/run_simulation.py b/src/scope/run_simulation.py
index 4a06781..baf1621 100644
--- a/src/scope/run_simulation.py
+++ b/src/scope/run_simulation.py
@@ -137,6 +137,11 @@ def make_data(
             flux_cube[:, i, :] *= throughput_factor
 
     if tellurics:
+        if instrument != "IGRINS":
+            raise NotImplementedError(
+                "Only IGRINS data is currently supported for telluric simulations."
+                "Please set tellurics=False, and at minimum ensure that your SNR input includes telluric losses."
+            )
         flux_cube = add_tellurics(
             wlgrid,
             flux_cube,
@@ -172,9 +177,8 @@ def make_data(
             noise_model = instrument
         else:
             noise_model = "constant"
-        # print(f"Adding noise with model {noise_model}")
+        print(f"Adding noise with model {noise_model}")
         flux_cube = add_noise_cube(flux_cube, wlgrid, SNR, noise_model=noise_model)
-
     flux_cube = detrend_cube(flux_cube, n_order, n_exposure)
 
     if wav_error:
@@ -186,7 +190,6 @@ def make_data(
     flux_cube = detrend_cube(flux_cube, n_order, n_exposure)
     flux_cube[np.isnan(flux_cube)] = 0.0
     flux_cube_nopca = flux_cube.copy()
-
     if observation == "transmission" and divide_out_of_transit:
         # generate the out of transit baseline
         n_exposures_baseline = (
@@ -240,8 +243,7 @@ def make_data(
 
         # divide out the flux cube
         flux_cube /= median_out_of_transit  # todo: check axes work out
-    print("mean flux cube", np.mean(flux_cube))
-    print("std flux cube", np.std(flux_cube))
+
     if do_pca:
         for j in range(n_order):
             flux_cube[j] -= np.mean(flux_cube[j])
@@ -257,12 +259,11 @@ def make_data(
                 flux_cube[j][i] -= np.mean(flux_cube[j][i])
 
     if np.all(pca_noise_matrix == 0):
-        print("was all zero")
         pca_noise_matrix = np.ones_like(pca_noise_matrix)
+
     if tellurics:
         return pca_noise_matrix, flux_cube, flux_cube_nopca, just_tellurics
-    print("mean flux cube", np.mean(flux_cube))
-    print("std flux cube", np.std(flux_cube))
+
     return (
         pca_noise_matrix,
         flux_cube,
diff --git a/src/scope/tests/test_run_simulation.py b/src/scope/tests/test_run_simulation.py
index 25c53bc..31f7d74 100644
--- a/src/scope/tests/test_run_simulation.py
+++ b/src/scope/tests/test_run_simulation.py
@@ -518,11 +518,9 @@ def test_nopca_fluxlevels_deterministic(
     np.testing.assert_array_equal(flux_cube_nopca, flux_cube_nopca_b)
 
 
-def test_crires_simulation(test_inputs):
+def test_igrins_emission(test_inputs):
     data_cube_path = os.path.join(test_data_path, "data_RAW_20201214_wl_algn_03.pic")
-    planet_spectrum_path = os.path.join(
-        test_data_path, "best_fit_spectrum_tran_gj1214_steam_nir.pic"
-    )
+    planet_spectrum_path = os.path.join(test_data_path, "best_fit_spectrum.npy")
     snr_path = os.path.join(test_data_path, "output1.json")
     star_spectrum_path = os.path.join(test_data_path, "PHOENIX_5605_4.33.txt")
     simulate_observation(
@@ -530,9 +528,9 @@ def test_crires_simulation(test_inputs):
         planet_spectrum_path=planet_spectrum_path,
         star_spectrum_path=star_spectrum_path,
         data_cube_path=data_cube_path,
-        phase_start=0,
-        phase_end=1,
-        n_exposures=2,
+        phase_start=0.4,
+        phase_end=0.46,
+        n_exposures=5,
         observation="emission",
         blaze=True,
         n_princ_comp=4,
@@ -571,3 +569,58 @@ def test_crires_simulation(test_inputs):
     filetype = glob("src/scope/output/yourfirstsimulation/*ccf*")
     ccfs = np.loadtxt(filetype[0])
     assert np.sum(ccfs) > 0
+
+
+def test_crires_emission(test_inputs):
+    data_cube_path = os.path.join(test_data_path, "data_RAW_20201214_wl_algn_03.pic")
+    planet_spectrum_path = os.path.join(
+        test_data_path, "best_fit_spectrum_tran_gj1214_steam_nir.pic"
+    )
+    snr_path = os.path.join(test_data_path, "output1.json")
+    star_spectrum_path = os.path.join(test_data_path, "PHOENIX_5605_4.33.txt")
+    simulate_observation(
+        # problem: this is all 1s!
+        planet_spectrum_path=planet_spectrum_path,
+        star_spectrum_path=star_spectrum_path,
+        data_cube_path=data_cube_path,
+        phase_start=0.4,
+        phase_end=0.46,
+        n_exposures=5,
+        observation="transmission",
+        blaze=False,
+        n_princ_comp=4,
+        star=True,
+        SNR=250,
+        telluric=False,
+        tell_type="data-driven",
+        time_dep_tell=False,
+        wav_error=False,
+        rv_semiamp_orbit=0.3229,
+        order_dep_throughput=True,
+        Rp=1.21,  # Jupiter radii,
+        Rstar=0.955,  # solar radii
+        kp=192.02,  # planetary orbital velocity, km/s
+        v_rot=4.5,
+        scale=1.0,
+        v_sys=0.0,
+        modelname="yoursecondsimulation",
+        planet_name="GJ1214b",
+        divide_out_of_transit=False,
+        out_of_transit_dur=0.1,
+        include_rm=False,
+        instrument="CRIRES+",
+        v_rot_star=3.0,
+        a=0.033,  #
+        lambda_misalign=0.0,
+        inc=90.0,
+        n_kp=20,
+        n_vsys=20,
+        seed=42,
+        vary_throughput=True,
+        snr_path=snr_path,
+    )
+
+    # check that there are nonzeros, basically that it runs
+    filetype = glob("src/scope/output/yoursecondsimulation/*ccf*")
+    ccfs = np.loadtxt(filetype[0])
+    assert np.sum(ccfs) != 0
diff --git a/src/scope/utils.py b/src/scope/utils.py
index 8142081..7d06f51 100644
--- a/src/scope/utils.py
+++ b/src/scope/utils.py
@@ -60,7 +60,7 @@ def read_crires_data(data_path):
                 "wavelength"
             ]
 
-    return n_orders, n_wavs, wl_grid, snr_grid
+    return n_orders, n_wavs, wl_grid * 1e6, snr_grid
 
 
 @njit

From 451cc8b14bd726c054f989dbe9a48ad16368ca44 Mon Sep 17 00:00:00 2001
From: arjunsavel <asavel@gmail.com>
Date: Mon, 3 Mar 2025 17:07:24 -0500
Subject: [PATCH 22/22] map the instrumental resolutions

---
 src/scope/run_simulation.py            | 2 +-
 src/scope/tests/conftest.py            | 2 +-
 src/scope/tests/test_run_simulation.py | 2 +-
 src/scope/utils.py                     | 8 +++++++-
 4 files changed, 10 insertions(+), 4 deletions(-)

diff --git a/src/scope/run_simulation.py b/src/scope/run_simulation.py
index baf1621..3c6cc65 100644
--- a/src/scope/run_simulation.py
+++ b/src/scope/run_simulation.py
@@ -495,7 +495,7 @@ def simulate_observation(
     Fp_conv_rot = np.convolve(Fp, rot_ker, mode="same")
 
     # instrument profile convolution
-    instrument_kernel = get_instrument_kernel()
+    instrument_kernel = get_instrument_kernel(instrument)
     Fp_conv = np.convolve(Fp_conv_rot, instrument_kernel, mode="same")
 
     star_wave, star_flux = np.loadtxt(
diff --git a/src/scope/tests/conftest.py b/src/scope/tests/conftest.py
index 2d5a501..196b7cf 100644
--- a/src/scope/tests/conftest.py
+++ b/src/scope/tests/conftest.py
@@ -28,7 +28,7 @@ def test_inputs():
     rot_ker = get_rot_ker(v_rot, wl_model)
     Fp_conv_rot = np.convolve(Fp, rot_ker, mode="same")
     # instrument profile convolution
-    instrument_kernel = get_instrument_kernel()
+    instrument_kernel = get_instrument_kernel("IGRINS")
     Fp_conv = np.convolve(Fp_conv_rot, instrument_kernel, mode="same")
     star_spectrum_path = os.path.join(test_data_path, "PHOENIX_5605_4.33.txt")
     star_wave, star_flux = np.loadtxt(
diff --git a/src/scope/tests/test_run_simulation.py b/src/scope/tests/test_run_simulation.py
index 31f7d74..24926a2 100644
--- a/src/scope/tests/test_run_simulation.py
+++ b/src/scope/tests/test_run_simulation.py
@@ -568,7 +568,7 @@ def test_igrins_emission(test_inputs):
     # check that there are nonzeros, basically that it runs
     filetype = glob("src/scope/output/yourfirstsimulation/*ccf*")
     ccfs = np.loadtxt(filetype[0])
-    assert np.sum(ccfs) > 0
+    assert np.sum(ccfs) != 0
 
 
 def test_crires_emission(test_inputs):
diff --git a/src/scope/utils.py b/src/scope/utils.py
index 7d06f51..2c814ac 100644
--- a/src/scope/utils.py
+++ b/src/scope/utils.py
@@ -131,7 +131,7 @@ def make_outdir(outdir):
         print("Directory already exists. Continuing!")
 
 
-def get_instrument_kernel(resolution_ratio=250000 / 45000, kernel_size=41):
+def get_instrument_kernel(instrument, model_resolution=250000, kernel_size=41):
     """
     Creates a Gaussian kernel for instrument response using an alternative implementation.
 
@@ -147,6 +147,12 @@ def get_instrument_kernel(resolution_ratio=250000 / 45000, kernel_size=41):
     numpy.ndarray
         Normalized Gaussian kernel
     """
+    instrument_resolution_dict = {
+        "IGRINS": 45000,
+        "CRIRES+": 145000,
+    }
+    instrument_resolution = instrument_resolution_dict[instrument]
+    resolution_ratio = instrument_resolution / model_resolution
     # Ensure kernel size is odd
     if kernel_size % 2 == 0:
         raise ValueError("Kernel size must be odd")