"""Export coronagraph performance curves to external formats.
Provides functions to export coronagraph performance data in EXOSIMS FITS
format and AYO's CSV format.
All functions operate on a Coronagraph instance passed as the first argument.
"""
from __future__ import annotations
import csv
import json
from pathlib import Path
from typing import TYPE_CHECKING
import astropy.io.fits as pyfits
import astropy.units as u
import numpy as np
from lod_unit import lod
from .logger import logger
from .performance import compute_radial_average
if TYPE_CHECKING:
from .coronagraph import Coronagraph
[docs]
def export_exosims(
coro: Coronagraph,
aperture_radius_lod: float = 0.7,
fit_gaussian_for_core_area: bool = False,
use_phot_aperture_as_min: bool = False,
units: str = "LAMBDA/D",
) -> dict:
"""Export performance curves in EXOSIMS format.
Writes individual FITS files and a ``specs.json`` to
``<yip_path>/exosims/``.
Returns the EXOSIMS specs dictionary.
"""
# Validate that interpolators exist
required = [
"throughput_interp",
"raw_contrast_interp",
"occ_trans_interp",
"core_area_interp",
"core_intensity_interp",
]
missing = [a for a in required if not hasattr(coro, a)]
if missing:
raise ValueError(
f"Performance curves not computed. Missing: {missing}. "
"Call compute_all_performance_curves() first."
)
# Build separations
x_offsets = np.array(coro.offax.x_offsets)
separations = np.abs(x_offsets)
if hasattr(coro.offax, "max_offset_in_image"):
max_sep = coro.offax.max_offset_in_image.to(lod).value
separations = separations[separations <= max_sep]
separations = np.sort(np.unique(separations))
# Evaluate interpolators
throughput = coro.throughput_interp(separations)
raw_contrast = coro.raw_contrast_interp(separations)
core_area = coro.core_area_interp(separations)
# Clip small negatives from spline artifacts
if np.any(throughput < 0):
n_neg = np.sum(throughput < 0)
min_val = np.min(throughput)
if min_val < -0.01:
raise ValueError(
f"Found {n_neg} negative throughput values. Min: {min_val:.6f}"
)
logger.debug(
f"Clipping {n_neg} small negative throughput values (min: {min_val:.6f})"
)
throughput = np.clip(throughput, 0, None)
if np.any(throughput > 1):
raise ValueError(f"Throughput > 1: max = {np.max(throughput):.3f}")
# Occulter transmission
sky_trans_data = coro.sky_trans()
sep_occ_trans, occ_trans = compute_radial_average(
sky_trans_data, coro.pixel_scale_arcsec.value
)
# Core mean intensity
core_intensities = (
coro.core_intensity_dict
if hasattr(coro, "core_intensity_dict")
else coro.core_mean_intensity_curve(plot=False)[1]
)
first_diam = next(iter(core_intensities.keys()))
stellar_psf = coro.stellar_intens(first_diam)
dims = stellar_psf.shape
nbins = int(np.floor(np.max(dims) / 2))
center = [coro.stellar_intens.center_x, coro.stellar_intens.center_y]
y, x = np.indices(stellar_psf.shape)
r = np.sqrt((x - center[0]) ** 2 + (y - center[1]) ** 2)
max_radius = np.max(r)
bins = np.linspace(0, max_radius, nbins + 1)
bin_centers = (bins[1:] + bins[:-1]) / 2
sep_core_intensity = bin_centers * coro.pixel_scale_arcsec.value
# Write FITS files
_save_to_exosims_format(
coro,
sep=separations,
throughput=throughput,
raw_contrast=raw_contrast,
core_area=core_area,
sep_occ_trans=sep_occ_trans,
occ_trans=occ_trans,
sep_core_intensity=sep_core_intensity,
core_intensities=core_intensities,
aperture_radius_lod=aperture_radius_lod,
fit_gaussian_for_core_area=fit_gaussian_for_core_area,
use_phot_aperture_as_min=use_phot_aperture_as_min,
units=units,
)
# IWA / OWA
IWA = coro.IWA
OWA = coro.OWA
to_arcsec = units.lower() == "arcsec"
if to_arcsec:
angunit = ((coro.header.lambda0) / (coro.header.diameter)).to(
u.arcsec, equivalencies=u.dimensionless_angles()
)
IWA_output = (IWA * angunit).to(u.arcsec).value
OWA_output = (OWA * angunit).to(u.arcsec).value
else:
IWA_output = IWA.to_value(lod)
OWA_output = OWA.to_value(lod)
# Core area filename or scalar
core_area_fname = (
"core_area.fits"
if fit_gaussian_for_core_area
else aperture_radius_lod**2 * np.pi
)
# deltaLam
deltaLam = None
if coro.header.maxlam is not None and coro.header.minlam is not None:
deltaLam = (coro.header.maxlam - coro.header.minlam).to(u.nm).value
# Specs dict
outdict = {
"pupilDiam": coro.header.diameter.to(u.m).value,
"obscurFac": coro.header.obscured,
"starlightSuppressionSystems": [
{
"name": coro.name,
"lam": coro.header.lambda0.to(u.nm).value,
"deltaLam": deltaLam,
"occ_trans": "occ_trans.fits",
"core_thruput": "core_thruput.fits",
"core_mean_intensity": "core_mean_intensity.fits",
"core_area": core_area_fname,
"IWA": IWA_output,
"OWA": OWA_output,
"input_angle_units": units,
}
],
}
exosims_dir = Path(coro.yip_path, "exosims")
specs_file = exosims_dir / "specs.json"
with open(specs_file, "w") as f:
json.dump(outdict, f, indent=2)
logger.info(f"EXOSIMS specs saved to {specs_file}")
return outdict
[docs]
def export_ayo_csv(
coro: Coronagraph,
output_path,
sep_min: float = 0.125,
sep_max: float = 32.0,
sep_step: float = 0.25,
contrast_floor: float = 1e-10,
ppf: float = 30.0,
) -> Path:
"""Export performance curves in AYO CSV format.
Returns the path to the saved CSV file.
"""
output_path = Path(output_path)
separations = np.arange(sep_min, sep_max + sep_step / 2, sep_step)
raw_contrast = np.abs(coro.raw_contrast(separations))
contrast = np.maximum(raw_contrast, contrast_floor)
noise_floor_val = contrast / ppf
throughput = coro.throughput(separations)
occ_trans = coro.occulter_transmission(separations)
with open(output_path, "w", newline="") as f:
writer = csv.writer(f)
writer.writerow(
[
"Sep (l/D)",
"Contrast (for a point source)",
"Noise floor (point source 1-sigma)",
"Core throughput",
"Skytrans",
]
)
for i, sep in enumerate(separations):
writer.writerow(
[
f"{sep:.6f}",
f"{contrast[i]:.6e}",
f"{noise_floor_val[i]:.6e}",
f"{throughput[i]:.6e}",
f"{occ_trans[i]:.6f}",
]
)
logger.info(f"AYO-format CSV saved to {output_path}")
return output_path
