From 593a70a0b7e85921d96f7ecfeac134ea15d2f65f Mon Sep 17 00:00:00 2001 From: Leonardo dos Santos Date: Wed, 9 Sep 2026 11:12:07 -0400 Subject: [PATCH] First commit of the exo-Venus imaging tutorial --- tutorials/imaging_early-vs-modern_Venus.ipynb | 312 ++++++++++++++++++ 1 file changed, 312 insertions(+) create mode 100644 tutorials/imaging_early-vs-modern_Venus.ipynb diff --git a/tutorials/imaging_early-vs-modern_Venus.ipynb b/tutorials/imaging_early-vs-modern_Venus.ipynb new file mode 100644 index 0000000..0341197 --- /dev/null +++ b/tutorials/imaging_early-vs-modern_Venus.ipynb @@ -0,0 +1,312 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "# pyEDITH Tutorial: exo-Venus in Imaging Mode\n", + "\n", + "This tutorial will guide you through using pyEDITH in imaging mode to discern an early exo-Venus from a modern exo-Venus. We'll explore this science case in various filters.\n", + "\n", + "\n", + "## Before we start\n", + "\n", + "Make sure you follow the instructions on the [Installation](https://pyedith.readthedocs.io/en/latest/installation.html) page.\n", + "\n", + "## 1. Setup and Imports\n", + "\n", + "First, let's import the necessary modules and set up our environment:\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from astropy import units as u\n", + "from pyEDITH import parse_input, calculate_texp, calculate_snr, AstrophysicalScene, Observation, Observatory, set_verbosity,calculate_exposure_time_or_snr,Filter\n", + "from pyEDITH.units import *\n", + "\n", + "# Set verbosity to INFO, showing info, warnings and errors. Other options are \"warning\" (warnings and errors), \n", + "# \"quiet\" (only errors), and \"debug\" (all logs)\n", + "set_verbosity(level='info') \n", + "\n", + "# Set the necessary environment variables --> REPLACE WITH YOUR PATHS. \n", + "# You can also open your .bashrc (or .zshrc) and type:\n", + "# export SCI_ENG_DIR=\"/path/to/Sci-Eng-Interface/hwo_sci_eng\"\n", + "# export YIP_CORO_DIR=\"/path/to/yips\"\n", + "\n", + "# Loading HWO style package to make pretty plots\n", + "import hwostyle\n", + "hwostyle.use(\"light\")\n", + "colors = hwostyle.palette\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Since the exo-Venus theoretical spectra are in geometric albedo, we will need a function that converts that into contrast ratio as a function of phase angle, planet radius and separation." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Function to calculate contrast ratio as a function of\n", + "# geometric albedo, phase angle, planet radius and angular separation\n", + "# Astropy units will be needed for the planet radius and separation\n", + "def contrast_ratio(geometric_albedo, phase_angle, planet_radius, separation):\n", + " phase_function = (np.sin(phase_angle) + (np.pi - phase_angle) * np.cos(phase_angle)) / np.pi\n", + " contrast = geometric_albedo * phase_function * (planet_radius / separation).decompose().value ** 2\n", + " return contrast" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will also need a function that calculates the projected angular separation based on the planet's semi-major axis, distance, eccentricity, orbital inclination, argument of periastron and true anomaly. We define this function below." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def projected_angular_separation(\n", + " a_au: float,\n", + " d_pc: float,\n", + " ecc: float,\n", + " inc_deg: float,\n", + " omega_deg: float,\n", + " nu_deg: float,\n", + ") -> float:\n", + " \"\"\"Calculates the instantaneous projected angular separation of an exoplanet.\n", + "\n", + " Parameters:\n", + " a_au : Semi-major axis (Astronomical Units, AU)\n", + " d_pc : Distance to host star (Parsecs, pc)\n", + " ecc : Orbital eccentricity (0 <= ecc < 1)\n", + " inc_deg : Orbital inclination in degrees (0 = face-on, 90 = edge-on)\n", + " omega_deg : Argument of periastron in degrees\n", + " nu_deg : True anomaly (position in orbit) in degrees\n", + "\n", + " Returns:\n", + " Projected angular separation in arcseconds (arcsec)\n", + " \"\"\"\n", + " # Convert angles from degrees to radians\n", + " inc = np.radians(inc_deg)\n", + " omega = np.radians(omega_deg)\n", + " nu = np.radians(nu_deg)\n", + "\n", + " # 1. Calculate the instantaneous physical separation r(t) in AU\n", + " r_au = (a_au * (1 - ecc**2)) / (1 + ecc * np.cos(nu))\n", + "\n", + " # 2. Compute 3D position vector in orbital plane coordinates\n", + " # Primary axis aligned with line of nodes\n", + " u = omega + nu\n", + "\n", + " # 3. Apply projection geometry (face-on vs edge-on tilt)\n", + " # Projected physical separation in the sky plane (in AU)\n", + " r_proj_au = r_au * np.sqrt((np.cos(u)) ** 2 + (np.sin(u) * np.cos(inc)) ** 2)\n", + "\n", + " # 4. Convert projected physical distance (AU) at distance (pc) to arcseconds\n", + " theta_arcsec = r_proj_au / d_pc\n", + "\n", + " return theta_arcsec" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we read the Venus spectra from the HWO Tools GitHub repository and plot their corresponding contrast ratios." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Read spectra\n", + "BASE_URL = \"https://raw.githubusercontent.com/spacetelescope/hwo-tools/main/coron_model/planets/\"\n", + "\n", + "early_venus = np.loadtxt(BASE_URL + 'EarlyVenus_geo_albedo.txt')\n", + "early_venus_wl = early_venus[:, 0]\n", + "early_venus_alb = early_venus[:, 1]\n", + "modern_venus = np.loadtxt(BASE_URL + 'Venus_geo_albedo.txt')\n", + "modern_venus_wl = modern_venus[:, 0]\n", + "modern_venus_alb = modern_venus[:, 1]\n", + "\n", + "# Set Venus's parameters\n", + "ap = 0.7233 * u.au\n", + "rp = 0.95 * u.earthRad\n", + "dist = 10 * u.pc\n", + "\n", + "# Calculate contrast ratios\n", + "early_venus_fpfs = contrast_ratio(early_venus_alb, np.pi / 2, rp, ap)\n", + "modern_venus_fpfs = contrast_ratio(modern_venus_alb, np.pi / 2, rp, ap)\n", + "\n", + "# Define wavelength band of the filter, this will be used for calculations later\n", + "# but, for now, we only use it for plotting purposes\n", + "filter_name = 'UVIS/F750W'\n", + "wavelength_center = 0.75\n", + "bandwidth = 0.2\n", + "waveband = np.array([wavelength_center - bandwidth / 2, wavelength_center + bandwidth / 2])\n", + "\n", + "plt.plot(early_venus_wl, early_venus_fpfs, label='Early Venus')\n", + "plt.plot(modern_venus_wl, modern_venus_fpfs, label='Modern Venus')\n", + "plt.axvspan(xmin=waveband[0], xmax=waveband[1], color='k', alpha=0.1)\n", + "plt.xlabel(r'Wavelength ($\\mu$m)')\n", + "plt.ylabel('Fp/Fs')\n", + "plt.legend()\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 2. Defining Input Parameters\n", + "\n", + "Let's set up the parameters for a Venus-like planet around a Sun-like star. Our objective is to estimate the exposure time to differentiate between modern and early exo-Venus. To do that, first we need to estimate the required SNR." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Handy numpy array slicing to work with the model within the filter bandpass\n", + "sect_early = np.logical_and(early_venus_wl > waveband[0], early_venus_wl < waveband[1])\n", + "sect_modern = np.logical_and(modern_venus_wl > waveband[0], modern_venus_wl < waveband[1])\n", + "\n", + "# Average contrast ratios inside the wavelength band\n", + "Fp_Fs_band_earlyVenus = np.mean(early_venus_fpfs[sect_early])\n", + "Fp_Fs_band_modernVenus = np.mean(modern_venus_fpfs[sect_modern])\n", + "\n", + "# Ratio between modern and early Venus\n", + "ratio_modern_early = Fp_Fs_band_modernVenus / Fp_Fs_band_earlyVenus\n", + "print('Early and modern Venus differ by a factor of {:.2f} in Fp/Fs in the band {}-{} micron.'.format(ratio_modern_early, waveband[0], waveband[1]))\n", + "\n", + "# Required SNR for a 5-sigma differentiation between modern and early Venus\n", + "required_snr_5sigma = 5 / (ratio_modern_early - 1)\n", + "print('The required SNR to differ between Modern and Early Venus with 5-sigma confidence is thus roughly {:.2f}.'.format(required_snr_5sigma))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 1.3 Running the ETC\n", + "\n", + "We can now calculate the exposure time. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Calculate separation in arcsec\n", + "separation = projected_angular_separation(ap.value, dist.value, ecc=0.0, inc_deg=45.0, omega_deg=0.0, nu_deg=0.)\n", + "\n", + "imaging_params = {\n", + " 'wavelength': wavelength_center, # Wavelength in microns\n", + " 'snr': required_snr_5sigma, # Desired signal-to-noise ratio\n", + " 'CRb_multiplier': 2.0, # Count rate ratio multiplier (assuming differential imaging for PSF subtraction)\n", + " 'psf_trunc_ratio': 0.3, # PSF Truncation Ratio to calculate photometric aperture of solid angle Omega. \n", + " 'distance': dist.value, # Distance to star in parsecs\n", + " 'FstarV_10pc': 122.9279, # Stellar flux at 10 pc in the V band [ph/cm2/s/nm]\n", + " 'Fstar_10pc': 115.59984, # Stellar flux at 10 pc in the observed band [ph/cm2/s/nm]\n", + " 'Fp/Fs': Fp_Fs_band_earlyVenus, # Planet-to-star contrast\n", + " 'stellar_radius': 1, # Stellar radius in solar radii\n", + " 'nzodis': 3.0, # Number of zodiacal light disks\n", + " 'ra': 236.00757736823, # Right ascension of star [deg]\n", + " 'dec': 2.51516683165, # Declination of star [deg]\n", + " 'separation': separation, # Separation between star and planet in arcseconds\n", + " 'observatory_preset': 'EAC1', # Preset observatory configuration\n", + " 'observing_mode': 'IMAGER', # Observing mode\n", + " 'filter_list': Filter(filter_name, \n", + " center=wavelength_center, \n", + " bandwidth=bandwidth, \n", + " type=\"IMAGER\") # filter in which to perform the observation\n", + "}\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Make the parameters be the shape that the code desires \n", + "parsed_parameters= parse_input.parse_parameters(imaging_params)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Calculate Exposure time\n", + "texp, validation_output = calculate_texp(imaging_params)\n", + "for filter_name, result in texp.items():\n", + " print(f\"Filter: {filter_name}\")\n", + " print(f\" Wavelength: {result['wavelength']}\")\n", + " print(f\" Calculated exposure time: {result['exposure_time'].to(u.hr)}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`validation_output` contains some interesting quantities that you can use to double check or validate, or to plot additional quantities." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "validation_output" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.14.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}