calc SNR with web queries

environment.yml

@@ -3,6 +3,7 @@ channels:
   - conda-forge
   - defaults
 dependencies:
+  - selenium
   - bzip2=1.0.8=h99b78c6_7
   - ca-certificates=2024.7.4=hf0a4a13_0
   - colorama=0.4.6=pyhd8ed1ab_0

pyproject.toml

@@ -46,7 +46,10 @@ dependencies = [
     "exoplanet-core",
     "pymc>=4",
     "schwimmbad",
-    "astroquery"
+    "selenium",
+    "requests",
+    "astroquery",
+    "webdriver_manager"
 ]
 dynamic = ["version"]
 

src/scope/calc_quantities.py

@@ -1,6 +1,7 @@
 """
 This module contains functions to calculate derived quantities from the input data.
 """
+import os
 from math import floor
 
 import astropy.constants as const
@@ -9,6 +10,11 @@
 import numpy as np
 from exoplanet.orbits.keplerian import KeplerianOrbit
 
+from scope.logger import get_logger
+from scope.scrape_igrins_etc import scrape_crires_plus_etc, scrape_igrins_etc
+
+logger = get_logger()
+
 # refactor the below two functions to be a single function
 
 
@@ -72,6 +78,34 @@ def calc_param_boilerplate(param, args, data, distance_unit):
     return data
 
 
+def calc_snr(param, args, data):
+    detector_reset_times = {
+        "IGRINS": 60,
+        "CRIRES+": 0,
+    }  # these are in seconds. crires+ considers this as part of DIT in the ETC.
+    if data[param] == -1:  # only calculate if set to -1
+        kmag = data["Kmag"]
+        n_exposures = data["n_exposures"]
+        instrument = data["instrument"]
+        detector_reset_time = detector_reset_times[instrument]
+        duration = (
+            data["P_rot"] * (data["phase_end"] - data["phase_start"]) * 24 * 3600
+        )  # in seconds
+        exptime = duration / n_exposures - detector_reset_time
+        # ok this is pretty good lol
+
+        if instrument == "IGRINS":
+            snr = scrape_igrins_etc(kmag, exptime)
+        elif instrument == "CRIRES+":
+            snr = scrape_crires_plus_etc(kmag, exptime)
+        else:
+            logger.error(f"Unknown instrument {instrument}!")
+        # need to construct a JSON object to send to the CRRIES+ SNR calculator.
+        # ok need to submit the job automatically to the IGRINS SNR calculator.
+        data[param] = snr
+        data["snr_path"] = os.getcwd() + "/snr.json"
+
+
 def calc_max_npix(param, args, data):
     # todo: refactor the instrument data to be somwhere!
     instrument_resolutions_dict = {
@@ -80,7 +114,7 @@ def calc_max_npix(param, args, data):
     }
     pixel_per_res = {"IGRINS": 3.3, "CRIRES+": 3.3}
 
-    if data[param] == 0:
+    if data[param] == -1:
         check_args(args, data, param)
         period = data["P_rot"]
         mstar = data["Mstar"]
@@ -262,7 +296,7 @@ def calc_crossing_time(
     period = period * u.day
     dphase_per_exp = (np.min(res_crossing_time_phases) / period).si
 
-    phasedur = phase_end - phase_start  # degrees. todo: calculate.
+    phasedur = phase_end - phase_start
     n_exp = phasedur / dphase_per_exp
     # todo: actually get phase_start and phase_end in there
 

src/scope/input_output.py

@@ -4,6 +4,7 @@
 
 import argparse
 import io
+import json
 import os
 import warnings
 from datetime import datetime
@@ -40,6 +41,7 @@ def __init__(self, message="scope input file error:"):
     "e": "pl_orbeccen",
     "peri": "pl_orblper",
     "v_rot_star": "st_vsin",
+    "Kmag": "sy_kmag",
 }
 
 
@@ -428,6 +430,48 @@ def parse_arguments():
     return parser.parse_args()
 
 
+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"
+            ]
+
+            snr_grid[i * order_len + j] = data["data"]["orders"][i]["detectors"][j][
+                "plots"
+            ]["snr"]["snr"]
+
+    return n_orders, n_wavs, wl_grid * 1e6, snr_grid
+
+
 def refresh_db():
     """
     Refresh the database with the latest exoplanet data.

src/scope/noise.py

@@ -81,7 +81,7 @@ def add_quadratic_noise(flux_cube_model, wl_grid, SNR, IGRINS=False, **kwargs):
                     * 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
@@ -121,7 +121,7 @@ def add_crires_plus_noise(flux_cube_model, wl_grid, SNR):
     for exposure in range(flux_cube_model.shape[1]):
         for order in range(flux_cube_model.shape[0]):
             flux_level = (
-                SNR**2
+                SNR[order] ** 2
                 * flux_cube_model[order][exposure]
                 / np.nanmax(flux_cube_model[order][exposure])
             )

src/scope/run_simulation.py

@@ -567,7 +567,7 @@ def simulate_observation(
         pca_removal=pca_removal,
     )
 
-    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}_{seed}"
+    run_name = f"{n_princ_comp}_NPC_{blaze}_blaze_{star}_star_{telluric}_telluric_{round(np.mean(SNR))}_SNR_{tell_type}_{time_dep_tell}_{wav_error}_{order_dep_throughput}_{seed}"
 
     save_data(outdir, run_name, flux_cube, flux_cube_nopca, A_noplanet, just_tellurics)