Merge pull request #85 from arjunsavel/final_cleans [skip ci]

clean up some language and make final funcs public

docs/pages/inputfile.rst

@@ -84,7 +84,7 @@ These parameters set the underlying simulated physical objects.
    * - u2
      - | Second quadratic limb darkening coefficient
        | Not used if LD=False or in emission mode
-   * - include_rm
+   * - include_rm (experimental)
      - | Enable/disable the Rossiter-McLaughlin effect
        | Only affects transmission observations
    * - v_rot_star
@@ -158,9 +158,9 @@ The only analysis currently implemented is Principal Components Analysis (PCA).
 
    * - n_princ_comp
      - Number of principal components to remove before cross-correlation.
-   * - divide_out_of_transit
+   * - divide_out_of_transit (experimental)
      - | Enable/disable division by median out-of-transit data
        | Only used in transmission mode
-   * - out_of_transit_dur
+   * - out_of_transit_dur (experimental)
      - | Duration of out-of-transit data in units of transit duration
        | Only used if divide_out_of_transit=True in transmission mode

docs/pages/quickstart.rst

@@ -3,3 +3,23 @@ Quickstart
 The aim of `scope` is to simulate observations of exoplanet atmospheres with ground-based, high-resolution spectroscopy.
 Simulating the atmosphere *itself* — performing the radiative transfer calculations to generate a spectrum
 that depends on atmospheric properties — is out of scope for `scope`.
+
+Once the package is installed and data are downloaded, the simulation is as simple as running, from the command line:
+
+.. code-block:: python
+
+    python run_simulation.py
+
+
+
+`scope` can also be run within a Python interpreter. The following code snippet demonstrates how to run a simulation
+(note that these are not real filepaths!):
+
+.. code-block:: python
+
+    from scope import run_simulation
+
+    simulate_observation(
+    planet_spectrum_path="./planet_spectrum.pkl",
+    star_spectrum_path="./stellar_spectrum.pkl",
+    data_cube_path="./data_cube.pkl",

pyproject.toml

@@ -3,7 +3,7 @@ requires = ["setuptools>=64", "setuptools_scm[toml]>=8", "setuptools-git-version
 build-backend = "setuptools.build_meta"
 
 [project]
-name = "scope"
+name = "scope-astr"
 authors = [
     {name = "Arjun Savel", email = "asavel@umd.edu"},
     {name= "Megan Bedell"},
@@ -61,18 +61,27 @@ docs = ["nbsphinx",
 
 
 [tool.setuptools_scm]
-write_to = "src/scope/__version.py"
+write_to = "src/scope/_version.py"
+version_scheme = "release-branch-semver"
+local_scheme = "no-local-version"
+git_describe_command = "git describe --tags --exact-match"
 
 [tool.setuptools-git-versioning]
 enabled = true
 template = "{tag}"
-dev_template = "{tag}.post{ccount}+git.{sha}"
-dirty_template = "{tag}.post{ccount}+git.{sha}.dirty"
+dev_template = "{tag}"
+dirty_template = "{tag}"
+
+[tool.setuptools]
+include-package-data = true
+
+[tool.setuptools.packages.find]
+where = ["src"]  # if using src layout
 
 [project.urls]
 Homepage = "https://github.com/arjunsavel/scope"
 Issues = "https://github.com/arjunsavel/scope/issues"
-Documentation = "scope-astr.readthedocs.io"
+Documentation = "https://scope-astr.readthedocs.io/en/latest/"
 
 [tool.black]
 target_version = ['py310', 'py311']

src/scope/broadening.py

@@ -1,5 +1,5 @@
 """
-file to calculate the intensity map and broadening said map. and the velocity profile.
+File to calculate the intensity map and broadening said map. and the velocity profile.
 """
 
 

src/scope/calc_quantities.py

@@ -1,3 +1,7 @@
+"""
+This module contains functions to calculate derived quantities from the input data.
+"""
+
 import astropy.units as u
 import numpy as np
 

src/scope/ccf.py

@@ -1,3 +1,7 @@
+"""
+Calculates the cross-correlation function (and log likelihood function from the Brogi & Line 2019 mapping)
+"""
+
 from functools import wraps
 
 import jax

src/scope/emcee_fit_hires.py

@@ -1,4 +1,7 @@
-import pickle
+"""
+Module to fit (simulated) HRCCS data with emcee, using MPI (multi-core).
+The fit is with respect to the velocity parameters (Kp, Vsys) and the scale factor.
+"""
 import sys
 
 import emcee

src/scope/grid.py

@@ -1,3 +1,7 @@
+"""
+An example parameter grid.
+"""
+
 import numpy as np
 from sklearn.model_selection import ParameterGrid
 

src/scope/input.txt

@@ -33,7 +33,7 @@ scale                  1.0               # scaling factor for the model spectrum
 LD                     True              # whether to include limb darkening in the simulation or not. only matters if observation is set to transmission.
 u1                     0.1               # first quadratic limb darkening coefficient. not used if limb_darkening is set to False or if observation is set to emission.
 u2                     0.1               # second quadratic limb darkening coefficient. not used if limb_darkening is set to False or if observation is set to emission.
-include_rm             False             # whether to include Rossiter-McLaughlin effect in the simulation or not. only matters if observation is set to transmission.
+include_rm             False             # *experimental*. Whether to include Rossiter-McLaughlin effect in the simulation or not. only matters if observation is set to transmission.
 v_rot_star             3                 # equatorial rotational velocity of the star in km/s. only matters if include_rm is set to True. and observation is set to transmission.
 lambda_misalign        0                 # misalignment angle of the planet's orbit with respect to the star's rotation axis, in degrees. only matters if include_rm is set to True. and observation is set to transmission. [DB]
 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.
@@ -58,5 +58,5 @@ time_dep_tell          False         # whether the tellurics are time-dependent
 
 # Analysis Parameters
 n_princ_comp           4             # number of principal components to remove from the simulated data before cross-correlation.
-divide_out_of_transit  False         # whether to divide the in-transit data by median out-of-transit or not. only used if observation is set to transmission.
-out_of_transit_dur     1.          # duration of the out-of-transit data, in units of transit duration. only used if observation is set to transmission and divide_out_of_transit is set to True.
+divide_out_of_transit  False         # *experimental*. Whether to divide the in-transit data by median out-of-transit or not. only used if observation is set to transmission.
+out_of_transit_dur     1.          # *experimental*. Duration of the out-of-transit data, in units of transit duration. only used if observation is set to transmission and divide_out_of_transit is set to True.

src/scope/input_output.py

@@ -1,3 +1,7 @@
+"""
+Module to handle the input files.
+"""
+
 import io
 import os
 from datetime import datetime

src/scope/rm_effect.py

@@ -1,3 +1,7 @@
+"""
+Module for simulating the Rossiter-McLaughlin effect. This is experimental for now.
+"""
+
 import matplotlib.pyplot as plt
 import numpy as np
 from astropy import constants as const

src/scope/run_simulation.py

@@ -1,10 +1,11 @@
-# import sys,os
-# sys.path.append(os.path.realpath('..'))
-# sys.path.append(os.path.realpath('.'))
+"""
+Main module for running a simulation of the HRCCS data.
+This module contains the main functions for simulating the data and
+calculating the log likelihood and cross-correlation function of the data given the model parameters.
+"""
+
 import os
-import sys
 
-import jax
 import numpy as np
 from tqdm import tqdm
 
@@ -346,10 +347,9 @@ def calc_log_likelihood(
     CCF = 0.0
     logL = 0.0
     for order in range(n_order):
-        wlgrid_order = np.copy(wlgrid[order,])  # Cropped wavelengths
-        model_flux_cube = np.zeros(
-            (n_exposure, n_pixel)
-        )  # "shifted" model spectra array at each phase
+        # grab the wavelengths from each order
+        wlgrid_order = np.copy(wlgrid[order,])
+        model_flux_cube = np.zeros((n_exposure, n_pixel))
 
         model_flux_cube = doppler_shift_planet_star(
             model_flux_cube,
@@ -374,19 +374,15 @@ def calc_log_likelihood(
             reprocessing=True,
         )
 
-        # ok now do the PCA. where does it fall apart?
         if do_pca:
-            # process the model same as the "data"!
+            # this is the "model reprocessing" step.
             model_flux_cube *= A_noplanet[order]
             model_flux_cube, _ = perform_pca(model_flux_cube, n_princ_comp, False)
-        # I = np.ones(n_pixel)
 
         logl, ccf = calc_ccf(model_flux_cube, flux_cube[order], n_pixel)
         CCF += ccf
         logL += logl
 
-        # # todo: airmass detrending reprocessing
-
     return logL, CCF  # returning CCF and logL values
 
 
@@ -489,6 +485,7 @@ def simulate_observation(
             star_flux, star_wave, v_rot_star, phases, Rstar, inc, lambda_misalign, a, Rp
         )
 
+    # todo: swap out
     Fstar_conv = get_star_spline(
         star_wave, star_flux, wl_model, instrument_kernel, smooth=False
     )

src/scope/utils.py

@@ -1,3 +1,7 @@
+"""
+Utility functions for simulating HRCCS data.
+"""
+
 import os
 import pickle
 
@@ -6,7 +10,8 @@
 import numpy as np
 from exoplanet.orbits.keplerian import KeplerianOrbit
 from numba import njit
-from scipy.interpolate import splev, splrep
+from scipy import signal
+from scipy.interpolate import interp1d
 from scipy.ndimage import gaussian_filter1d
 from tqdm import tqdm
 
@@ -57,8 +62,11 @@ def doppler_shift_planet_star(
             ) + 1.0
             if not reprocessing:
                 model_flux_cube[exposure,] *= flux_star * Rstar**2
-
-        else:  # in transmission, after we "divide out" (with PCA) the star and tellurics, we're left with Fp.
+        elif observation == "emission":
+            model_flux_cube[exposure,] = flux_planet * (Rp_solar * rsun) ** 2
+        elif (
+            observation == "transmission"
+        ):  # in transmission, after we "divide out" (with PCA) the star and tellurics, we're left with Fp.
             I = calc_limb_darkening(u1, u2, a, b, Rstar, phases[exposure], LD)
             model_flux_cube[exposure,] = 1.0 - flux_planet * I
             if not reprocessing:
@@ -84,12 +92,38 @@ def make_outdir(outdir):
         print("Directory already exists. Continuing!")
 
 
-def get_instrument_kernel(resolution_ratio=250000 / 45000):
-    xker = np.arange(41) - 20
-    sigma = resolution_ratio / (2.0 * np.sqrt(2.0 * np.log(2.0)))  # nominal
-    yker = np.exp(-0.5 * (xker / sigma) ** 2.0)
-    yker /= yker.sum()
-    return yker
+def get_instrument_kernel(resolution_ratio=250000 / 45000, kernel_size=41):
+    """
+    Creates a Gaussian kernel for instrument response using an alternative implementation.
+
+    Parameters
+    ----------
+    resolution_ratio : float
+        Ratio of resolutions (default: 250000/45000)
+    kernel_size : int
+        Size of the kernel (must be odd, default: 41)
+
+    Returns
+    -------
+    numpy.ndarray
+        Normalized Gaussian kernel
+    """
+    # Ensure kernel size is odd
+    if kernel_size % 2 == 0:
+        raise ValueError("Kernel size must be odd")
+
+    # Convert resolution ratio to standard deviation
+    std = resolution_ratio / (2.0 * np.sqrt(2.0 * np.log(2.0)))
+
+    # Create the Gaussian window
+    gaussian_window = signal.windows.gaussian(
+        kernel_size, std=std, sym=True  # Ensure symmetry
+    )
+
+    # Normalize to sum to 1
+    normalized_kernel = gaussian_window / gaussian_window.sum()
+
+    return normalized_kernel
 
 
 def save_data(outdir, run_name, flux_cube, flux_cube_nopca, A_noplanet, just_tellurics):
@@ -387,7 +421,8 @@ def calc_rvs(v_sys, v_sys_measured, Kp, Kstar, phases):
 
 def get_star_spline(star_wave, star_flux, planet_wave, yker, smooth=True):
     """
-    calculates the stellar spline. accounting for convolution and such.
+    Calculates the stellar spline using an alternative implementation with linear interpolation
+    instead of B-splines.
 
     Inputs
     ------
@@ -405,20 +440,31 @@ def get_star_spline(star_wave, star_flux, planet_wave, yker, smooth=True):
     Returns
     -------
     star_flux: array
-        Fluxes of the star. Measured in W/m^2/micron. todo: check.
+        Interpolated and processed stellar fluxes. Measured in W/m^2/micron.
     """
-    star_spline = splrep(star_wave, star_flux, s=0.0)
+    # Create interpolation function using scipy's interp1d
+    # Using cubic interpolation for smoother results
+    interpolator = interp1d(
+        star_wave,
+        star_flux,
+        kind="cubic",
+        bounds_error=False,
+        fill_value=(star_flux[0], star_flux[-1]),  # Extrapolate with edge values
+    )
 
-    star_flux = splev(
-        planet_wave, star_spline, der=0
-    )  # now on same wavelength grid as planet wave
+    # Interpolate onto planet wavelength grid
+    interpolated_flux = interpolator(planet_wave)
 
-    star_flux = np.convolve(star_flux, yker, mode="same")  # convolving star too
+    # Perform convolution
+    convolved_flux = np.convolve(interpolated_flux, yker, mode="same")
 
+    # Apply Gaussian smoothing if requested
     if smooth:
-        star_flux = gaussian_filter1d(star_flux, 200)
+        final_flux = gaussian_filter1d(convolved_flux, sigma=200)
+    else:
+        final_flux = convolved_flux
 
-    return star_flux
+    return final_flux
 
 
 def change_wavelength_solution(wl_cube_model, flux_cube_model, doppler_shifts):