Skip to content

First commit of the exo-Venus imaging tutorial #31

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ladsantos wants to merge 1 commit into
base: main
Choose a base branch
from
Open

First commit of the exo-Venus imaging tutorial #31

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
312 changes: 312 additions & 0 deletions tutorials/imaging_early-vs-modern_Venus.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
}
Loading