calc n_max exposures self-consistently

src/scope/calc_quantities.py

@@ -1,9 +1,13 @@
 """
 This module contains functions to calculate derived quantities from the input data.
 """
+from math import floor
 
+import astropy.constants as const
 import astropy.units as u
+import matplotlib.pyplot as plt
 import numpy as np
+from exoplanet.orbits.keplerian import KeplerianOrbit
 
 # refactor the below two functions to be a single function
 
@@ -47,9 +51,66 @@ def calc_param_boilerplate(param, args, data, distance_unit):
     # calculate planetary orbital velocity
     calc_param_boilerplate("kp", ["a", "P_rot"], data, u.AU)
 
+    # calculate max npix per resolution element
+    calc_max_npix(
+        "n_exposures",
+        [
+            "phase_start",
+            "phase_end",
+            "P_rot",
+            "Mstar",
+            "e",
+            "Mp",
+            "Rstar",
+            "omega",
+            "b",
+            "instrument",
+        ],
+        data,
+    )
+
     return data
 
 
+def calc_max_npix(param, args, data):
+    # todo: refactor the instrument data to be somwhere!
+    instrument_resolutions_dict = {
+        "IGRINS": 45000,
+        "CRIRES+": 145000,
+    }
+    pixel_per_res = {"IGRINS": 3.3, "CRIRES+": 3.3}
+
+    if data[param] == 0:
+        check_args(args, data, param)
+        period = data["P_rot"]
+        mstar = data["Mstar"]
+        e = data["e"]  # need
+        mplanet = data["Mp"]
+        rstar = data["Rstar"]
+        peri = np.radians(data["omega"])
+        b = data["b"]  # need
+        R = instrument_resolutions_dict[data["instrument"]]
+        pixel_per_res = pixel_per_res[data["instrument"]]
+        plot = False
+        phase_start = data["phase_start"]
+        phase_end = data["phase_end"]
+        max_time, max_time_per_pixel, dphase_per_exp, n_exp = calc_crossing_time(
+            period=period,
+            mstar=mstar,
+            e=e,
+            mplanet=mplanet,
+            rstar=rstar,
+            peri=peri,
+            b=b,
+            R=R,
+            pix_per_res=pixel_per_res,
+            phase_start=phase_start,
+            phase_end=phase_end,
+            plot=plot,
+        )
+        data[param] = floor(n_exp.value)
+
+
 def calculate_velocity(distance, P, distance_unit=u.AU, time_unit=u.day):
     """
     Calculate a velocity over some period over some distance.
@@ -77,3 +138,137 @@ def convert_tdur_to_phase(tdur, period):
         The phase of the transit duration.
     """
     return ((tdur * u.hour) / (period * u.day)).si.value
+
+
+def calc_crossing_time(
+    period=1.80988198,
+    mstar=1.458,
+    e=0.000,
+    mplanet=0.894,
+    rstar=1.756,
+    peri=0,
+    b=0.027,
+    R=45000,
+    pix_per_res=3.3,
+    phase_start=0.9668567402328337,
+    phase_end=1.0331432597671664,
+    plot=False,
+):
+    """
+
+    todo: refactor this into separate functions, maybe?
+
+    Inputs
+    -------
+    autofilled for WASP-76b and IGRINS.
+
+    R: (int) resolution of spectrograph.
+    pix_per_res: (float) pixels per resolution element
+
+
+    Outputs
+    -------
+
+    min_time: minimum time during transit before lines start smearing across resolution elements.
+
+    min_time_per_pixel: minimum time during transit before lines start smearing across a single pixel.
+    dphase_per_exp: the amount of phase (values from 0 to 1) taken up by a single exposure, given above constraints.
+    n_exp: number of exposures you can take during transit.
+
+    """
+
+    orbit = KeplerianOrbit(
+        m_star=mstar,  # solar masses!
+        r_star=rstar,  # solar radii!
+        #     m_planet_units = u.jupiterMass,
+        t0=0,  # reference transit at 0!
+        period=period,
+        ecc=e,
+        b=b,
+        omega=np.radians(peri),  # periastron, radians
+        Omega=np.radians(peri),  # periastron, radians
+        m_planet=mplanet * 0.0009543,
+    )
+    t = np.linspace(0, period * 1.1, 1000)  # days
+    z_vel = orbit.get_relative_velocity(t)[2].eval()
+
+    z_vel *= 695700 / 86400  # km / s
+
+    phases = t / period
+    if plot:
+        plt.plot(phases, z_vel)
+
+        plt.axvline(0.5, color="gray", label="Eclipse")
+        plt.axvline(0.25, label="Quadrature (should be maximized)", color="teal")
+
+        plt.axvline(0.75, label="Quadrature (should be maximized)", color="teal")
+
+        plt.legend()
+        plt.xlabel("Time (days)")
+        plt.ylabel("Radial velocity (km/s)")
+    acceleration = (
+        np.diff(z_vel) / np.diff((t * u.day).to(u.s).si.value) * u.km / (u.s**2)
+    )
+
+    if plot:
+        plt.plot(phases[:-1], acceleration)
+
+        plt.axvline(0.5, color="gray", label="Eclipse")
+        plt.axvline(0.25, label="Quadrature (should be minimized)", color="teal")
+
+        plt.axvline(0.75, label="Quadrature (should be minimized)", color="teal")
+        plt.xlabel("Orbital phase")
+        plt.ylabel("Radial acceleration (km/s^2)")
+
+        # cool, this is the acceleration.
+
+        # now want the pixel crossing time.
+
+        plt.legend()
+
+    # R = c / delta v
+    # delta v = c / R
+
+    delta_v = (const.c / R).to(u.km / u.s)
+    res_crossing_time = abs(delta_v / acceleration).to(u.s)
+    if plot:
+        plt.figure()
+        plt.plot(phases[:-1], res_crossing_time)
+        plt.axvline(0.5, color="gray", label="Eclipse")
+        plt.axvline(0.25, label="Quadrature (should be maximized)", color="teal")
+
+        plt.axvline(0.75, label="Quadrature (should be maximized)", color="teal")
+
+        plt.legend()
+        plt.xlabel("Orbital phase")
+        plt.ylabel("Resolution crossing time (s)")
+        plt.yscale("log")
+
+        plt.figure()
+        plt.plot(phases[:-1], res_crossing_time)
+        plt.legend()
+
+        plt.ylim(820, 900)
+        plt.xlabel("Orbital phase")
+        plt.ylabel("Resolution crossing time (s)")
+
+    within_phases = (phases[:-1] > phase_start) & (phases[:-1] < phase_end)
+
+    res_crossing_time_phases = res_crossing_time[within_phases]
+
+    max_time = np.min(res_crossing_time_phases)
+
+    max_time_per_pixel = max_time / pix_per_res
+    period = period * u.day
+    dphase_per_exp = (np.min(res_crossing_time_phases) / period).si
+
+    phasedur = phase_end - phase_start  # degrees. todo: calculate.
+    n_exp = phasedur / dphase_per_exp
+    # todo: actually get phase_start and phase_end in there
+
+    # then query https://igrins-jj.firebaseapp.com/etc/simple:
+    # this gives us, for the given exposure time, what the SNR is going to be.
+    # well that's more the maximum time
+    # so that's the maximum time, but we want more than that many exposures.
+    # don't have to worry about pixel-crossing time.
+    return max_time, max_time_per_pixel, dphase_per_exp, n_exp

src/scope/input.txt

@@ -26,7 +26,12 @@ snr_path               /home/asavel/scratch/scope/src/scope/data/snr_IGRINS_LTT-
 
 # Astrophysical Parameters. Can specify DATABASE to pull based on name for parameters with [DB]
 Rp                     DATABASE          # planetary radius in Jupiter radii. [DB]
+Mp                     DATABASE          # planetary mass in Jupiter masses. [DB]
+e                      DATABASE          # eccentricity of the planet's orbit. [DB]
+omega                  DATABASE          # argument of periastron of the planet's orbit, in degrees. [DB]
+b                      DATABASE          # impact parameter of the planet's orbit. [DB]
 Rstar                  DATABASE          # stellar radius, in solar radii. [DB]
+Mstar                  DATABASE          # stellar mass, in solar masses. [DB]
 kp                     128.1              # expected planetary orbital velocity assuming circular orbit, in km/s. input NULL if you'd like to calculate from orbital parameters.
 v_sys                  0.0          # systemic velocity, in km/s.  [DB]
 v_rot                  NULL              # equatorial rotational velocity. input NULL if you'd like to calculate from orbital parameters.
@@ -53,7 +58,7 @@ vary_throughput        True          # whether to include throughput variations.
 observation            transmission      # type of observation to perform. supported observations are ``emission`` and ``transmission``.
 phase_start            0.986222200994206           # phase of the beginning of the observations. 0 is center of transit, 0.5 is secondary eclipse. If DATABASE, just the transit duration. [DB]
 phase_end              1.0137777990057941           # phase of the end of the observations. 0 is center of transit, 0.5 is secondary eclipse. If DATABASE, just the transit duration. [DB]
-n_exposures            17             # number of exposures to simulate. sets the phases of the exposures.
+n_exposures            0             # number of exposures to simulate. sets the phases of the exposures. if set to 0, the minimum number of exposures that prevent pixel crossing for the provided instrument is used.
 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                    250           # signal-to-noise ratio of the observations, per pixel. I.e., what sets the photon counts at the detector.

src/scope/input_output.py

@@ -25,14 +25,18 @@ def __init__(self, message="scope input file error:"):
 # Mapping between input file parameters and database columns
 parameter_mapping = {
     "Rp": "pl_radj",
+    "Mp": "pl_massj",
     "Rstar": "st_rad",
+    "Mstar": "st_mass",
     "v_sys": "system_velocity",
     "a": "pl_orbsmax",
     "P_rot": "pl_orbper",
     "v_sys": "st_radv",
     "planet_name": "pl_name",
     "Rp_solar": "planet_radius_solar",
     "lambda_misalign": "pl_projobliq",
+    "e": "pl_orbeccen",
+    "peri": "pl_orblper",
 }
 
 
@@ -198,6 +202,9 @@ def parse_input_file(
         "a",
         "u1",
         "u2",
+        "Mstar",
+        "Mp",
+        "peri",
     ]
 
     # Convert values to appropriate types
@@ -302,7 +309,6 @@ def write_input_file(data, output_file_path="input.txt"):
 
 
 def parse_arguments():
-
     parser = argparse.ArgumentParser(description="Simulate observation")
 
     # Required parameters
@@ -413,5 +419,4 @@ def parse_arguments():
         "--input_file", type=str, default="input.txt", help="Input file with parameters"
     )
 
-
     return parser.parse_args()

src/scope/tests/test_calc_quantities.py

@@ -0,0 +1,58 @@
+import unittest
+from scope.calc_quantities import calc_crossing_time
+
+
+class TestCrossingTime(unittest.TestCase):
+    """
+    tests the exposure time calculation.
+    """
+
+    period = 1.80988198
+    mstar = 1.458
+    e = 0.000
+    inc = 89.623
+    mplanet = 0.894
+    rstar = 1.756
+    peri = 0
+    b = 0.027
+    R = 45000
+    pix_per_res = 3.3
+    plot = False
+
+    def test_crossing_time_moremassive_star(self):
+        """
+        in a greater potential well, the crossing time is shorter.
+        """
+        low_mass_time, _, _, _ = calc_crossing_time()
+
+        high_mass_time, _, _, _ = calc_crossing_time(mstar=10 * self.mstar)
+        assert low_mass_time > high_mass_time
+
+    def test_crossing_time_moremassive_star(self):
+        """
+        on a longer period for the same star, the crossing time is longer.
+        """
+        short_period_time, _, _, _ = calc_crossing_time()
+
+        long_period_time, _, _, _ = calc_crossing_time(period=10 * self.period)
+        assert short_period_time < long_period_time
+
+    def test_crossing_time_pixel_per_res(self):
+        """
+        more pixels per resolution element, shorter time to cross the pixels.
+        """
+        _, few_pixels_per_element, _, _ = calc_crossing_time()
+
+        _, many_pixels_per_element, _, _ = calc_crossing_time(
+            pix_per_res=10 * self.pix_per_res
+        )
+        assert many_pixels_per_element < few_pixels_per_element
+
+    def test_crossing_time_resolution(self):
+        """
+        higher resolution, shorter crossing time.
+        """
+        low_resolution_time, _, _, _ = calc_crossing_time()
+
+        high_resolution_time, _, _, _ = calc_crossing_time(R=10 * self.R)
+        assert high_resolution_time < low_resolution_time