From c81899bddfd401951994f7d1039b35ec1f3afcb1 Mon Sep 17 00:00:00 2001 From: arjunsavel Date: Mon, 11 Nov 2024 13:55:16 -0500 Subject: [PATCH 1/6] clean up some language and make final funcs public --- pyproject.toml | 19 +++++++--- src/scope/run_simulation.py | 19 +++------- src/scope/utils.py | 71 ++++++++++++++++++++++++++++--------- 3 files changed, 74 insertions(+), 35 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 4ec7a11..c4ca0ad 100644 --- a/pyproject.toml +++ b/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'] diff --git a/src/scope/run_simulation.py b/src/scope/run_simulation.py index ad57c4b..d04f66a 100644 --- a/src/scope/run_simulation.py +++ b/src/scope/run_simulation.py @@ -1,10 +1,5 @@ -# import sys,os -# sys.path.append(os.path.realpath('..')) -# sys.path.append(os.path.realpath('.')) import os -import sys -import jax import numpy as np from tqdm import tqdm @@ -346,10 +341,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 +368,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 +479,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 ) diff --git a/src/scope/utils.py b/src/scope/utils.py index 4cf2200..ba7ae4a 100644 --- a/src/scope/utils.py +++ b/src/scope/utils.py @@ -6,7 +6,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 @@ -84,12 +85,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 +414,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 +433,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): From 3be9b6e0d252ca19331cbb13c1e10c86699c1f59 Mon Sep 17 00:00:00 2001 From: arjunsavel Date: Mon, 11 Nov 2024 14:10:11 -0500 Subject: [PATCH 2/6] little more verbiage to the quickstart --- docs/pages/quickstart.rst | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/docs/pages/quickstart.rst b/docs/pages/quickstart.rst index 78a58bd..1604696 100644 --- a/docs/pages/quickstart.rst +++ b/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", From 2c1a00098bb0260f6f0271323b9fa63122477639 Mon Sep 17 00:00:00 2001 From: arjunsavel Date: Tue, 12 Nov 2024 10:22:10 -0500 Subject: [PATCH 3/6] note what is experimental --- src/scope/input.txt | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/scope/input.txt b/src/scope/input.txt index 7bc6b88..818640f 100644 --- a/src/scope/input.txt +++ b/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. From 0dec86657114904ed078fd6469bf25717dfcca7d Mon Sep 17 00:00:00 2001 From: arjunsavel Date: Tue, 12 Nov 2024 10:23:05 -0500 Subject: [PATCH 4/6] add the inputfile experimental tags --- docs/pages/inputfile.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/pages/inputfile.rst b/docs/pages/inputfile.rst index ae5bc32..1cd1225 100644 --- a/docs/pages/inputfile.rst +++ b/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 From 9b3584f429f38d624dd8314ac5984b69d514bbc7 Mon Sep 17 00:00:00 2001 From: arjunsavel Date: Tue, 12 Nov 2024 10:29:01 -0500 Subject: [PATCH 5/6] complete module docstrings. thought i did this already lol --- src/scope/broadening.py | 2 +- src/scope/calc_quantities.py | 4 ++++ src/scope/ccf.py | 4 ++++ src/scope/emcee_fit_hires.py | 5 ++++- src/scope/grid.py | 4 ++++ src/scope/input_output.py | 4 ++++ src/scope/rm_effect.py | 4 ++++ src/scope/run_simulation.py | 6 ++++++ src/scope/utils.py | 4 ++++ 9 files changed, 35 insertions(+), 2 deletions(-) diff --git a/src/scope/broadening.py b/src/scope/broadening.py index 60b2588..dabb860 100644 --- a/src/scope/broadening.py +++ b/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. """ diff --git a/src/scope/calc_quantities.py b/src/scope/calc_quantities.py index 0aaa540..118b305 100644 --- a/src/scope/calc_quantities.py +++ b/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 diff --git a/src/scope/ccf.py b/src/scope/ccf.py index e4e6bb6..409cd25 100644 --- a/src/scope/ccf.py +++ b/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 diff --git a/src/scope/emcee_fit_hires.py b/src/scope/emcee_fit_hires.py index daaa3fc..5b0b22e 100644 --- a/src/scope/emcee_fit_hires.py +++ b/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 diff --git a/src/scope/grid.py b/src/scope/grid.py index d23ad2f..8c7012c 100644 --- a/src/scope/grid.py +++ b/src/scope/grid.py @@ -1,3 +1,7 @@ +""" +An example parameter grid. +""" + import numpy as np from sklearn.model_selection import ParameterGrid diff --git a/src/scope/input_output.py b/src/scope/input_output.py index 2267409..ea527fc 100644 --- a/src/scope/input_output.py +++ b/src/scope/input_output.py @@ -1,3 +1,7 @@ +""" +Module to handle the input files. +""" + import io import os from datetime import datetime diff --git a/src/scope/rm_effect.py b/src/scope/rm_effect.py index 9dc3034..f86760d 100644 --- a/src/scope/rm_effect.py +++ b/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 diff --git a/src/scope/run_simulation.py b/src/scope/run_simulation.py index d04f66a..6d4d512 100644 --- a/src/scope/run_simulation.py +++ b/src/scope/run_simulation.py @@ -1,3 +1,9 @@ +""" +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 numpy as np diff --git a/src/scope/utils.py b/src/scope/utils.py index ba7ae4a..00d506d 100644 --- a/src/scope/utils.py +++ b/src/scope/utils.py @@ -1,3 +1,7 @@ +""" +Utility functions for simulating HRCCS data. +""" + import os import pickle From 36e71c879fd8328e9dc4eba85311e79eda2a5795 Mon Sep 17 00:00:00 2001 From: arjunsavel Date: Tue, 12 Nov 2024 10:32:16 -0500 Subject: [PATCH 6/6] ah yes scale the planet correctly --- src/scope/utils.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/scope/utils.py b/src/scope/utils.py index 00d506d..4821e99 100644 --- a/src/scope/utils.py +++ b/src/scope/utils.py @@ -62,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: