API

krithika has several classes and functions that are used for different purposes. They are listed below:

To perform aperture photometry, the ApPhoto class can be used. It has several methods to perform aperture photometry on the data:

class krithika.ApPhoto(times, frames, errors, badpix, aprad=None, sky_rad1=None, sky_rad2=None, brightpix=False, nos_brightest=12, nos_faintest=None, minmax=None)[source]

Helper for aperture photometry on image time-series.

This class encapsulates common operations required for aperture photometry on a data cube of frames (shape (N_frames, Ny, Nx)), including centroid estimation, aperture and sky-mask construction, simple aperture flux extraction, and pixel-level decorrelation (PLD) basis construction.

Parameters:

times (array_like) – 1D array of time stamps corresponding to each frame.
frames (ndarray) – 3D data cube with shape (N_frames, Ny, Nx) containing image pixel values for each frame.
errors (ndarray) – Per-pixel uncertainties with the same shape as frames.
badpix (ndarray) – 2D bad-pixel mask for a single frame (shape (Ny, Nx)). A value of zero indicates a bad/ignored pixel; non-zero indicates a usable pixel.
aprad (float or None, optional) – Aperture radius in pixels for circular aperture photometry. Used by aperture_mask() and simple_aperture_photometry(). Ignored when brightpix is True. Default is None.
sky_rad1 (float or None, optional) – Inner radius (pixels) of the sky annulus. Used by sky_mask() and simple_aperture_photometry() for background estimation. If omitted and brightpix is False, an empty (zero) mask is returned. Default is None.
sky_rad2 (float or None, optional) – Outer radius (pixels) of the sky annulus. Used by sky_mask() and simple_aperture_photometry() for background estimation. Default is None.
brightpix (bool, optional) – If True, select nos_brightest brightest pixels to form an aperture instead of circular aperture. Used by aperture_mask() and sky_mask(). Default False.
nos_brightest (int, optional) – Number of brightest pixels to include in aperture when brightpix=True. Default 12.
nos_faintest (int, optional) – Number of faintest pixels to include in background aperture when it is not None and brightpix=True. When None, all pixels that are not in aperture are included. Default is None.
minmax (dict or None, optional) – Optional spatial bounds for brightest pixel selection with keys ‘rmin’, ‘rmax’, ‘cmin’, ‘cmax’. Default is None.

times, frames, errors, badpix

Type:: as above

aprad, sky_rad1, sky_rad2, brightpix, nos_brightest, minmax

Type:: as above

cen_r, cen_c

Per-frame row/column centroid positions (populated by find_center()).

Type:: ndarray

ap_bool_mask_pld

2D aperture boolean mask used by PLD routines.

Type:: ndarray

Psum, Phat, V, eigenvalues, PCA

Intermediate PLD quantities (pixel-sum, normalized fractions, PCA eigenvectors/values and projected time-series) populated by pixel_level_decorrelation().

Type:: ndarray

Main methods

------------

identify_crays(clip, niters)[source]: Flag cosmic-ray affected pixels and update badpix.

replace_nan(max_iter)[source]: Iteratively fill NaNs in frames from neighboring pixels.

find_center(rmin, rmax, cmin, cmax)[source]: Compute center-of-flux centroids per frame.

simple_aperture_photometry(method, plot, ...)[source]: Extract aperture photometry using manual, photutils, or median methods.

pixel_level_decorrelation(removeNan)[source]: Build normalized pixel level light curves and compute PCA basis for PLD.

Notes

The class stores both raw inputs and intermediate results used by subsequent calls; many methods modify the instance in-place and also return results for convenience.

aperture_mask(plot=False)[source]

Build a binary aperture mask (and masked frame/error arrays).

By default a circular aperture is created on the median image using the median centroid position computed from self.cen_r and self.cen_c. Alternatively, when self.brightpix=True, the mask is formed by selecting the self.nos_brightest brightest pixels from the median image (optionally limited by self.minmax bounds).

Parameters:

plot (bool, optional) – If True, produce a plot of median image with aperture plotted on the top (default False).

Returns:

aper_fl_mask (ndarray) – self.frames multiplied by the 2D aperture boolean mask (shape (N_frames, Ny, Nx)).
aper_err_mask (ndarray) – self.errors multiplied by the 2D aperture boolean mask.
ap_bool_msk (ndarray) – 2D binary (0/1) mask describing the aperture on a single frame (shape (Ny, Nx)).
fig (matplotlib.figure.Figure or None) – Figure of median image and aperture mask, if plot=True, otherwise None.
im (matplotlib.image.AxesImage or None) – Image object from the median image plot if plot=True, otherwise None.
axs (ndarray or None) – Axes of the figure if produced, otherwise None.

Notes

Uses instance attributes self.aprad, self.brightpix, self.nos_brightest, and self.minmax to control aperture construction.

find_center(rmin=None, rmax=None, cmin=None, cmax=None, plot=False)[source]

Compute center-of-flux centroids for each frame in the 3D data cube.

The centroid for each frame is computed as the center of flux (i.e. sum(x*I)/sum(I)) over either the full image or an optional subimage defined by the row/column limits rmin,rmax,cmin,cmax. Results are stored on the instance as self.cen_r and self.cen_c and also returned.

Parameters:

rmin (int or None, optional) – Optional integer bounds (rows and columns) defining a subimage to use for the centroid calculation. If omitted, the full frame is used.
rmax (int or None, optional) – Optional integer bounds (rows and columns) defining a subimage to use for the centroid calculation. If omitted, the full frame is used.
cmin (int or None, optional) – Optional integer bounds (rows and columns) defining a subimage to use for the centroid calculation. If omitted, the full frame is used.
cmax (int or None, optional) – Optional integer bounds (rows and columns) defining a subimage to use for the centroid calculation. If omitted, the full frame is used.
plot (bool, optional) – If True, produce a plot of centroids as a function of time (default False).

Returns:

cen_r, cen_c (ndarray) – 1D arrays (length N_frames) containing the row and column centroid positions (in pixel coordinates) for each frame.
fig (matplotlib.figure.Figure or None) – Figure of centroids vs time, if plot=True, otherwise None.
axs (ndarray or None) – Axes of the figure if produced, otherwise None.

growth_function(rmin, rmax, noise='pipe', plot=False, **kwargs)[source]

Compute the growth curve (flux vs aperture radius) and noise.

Parameters:

rmin (int) – Minimum aperture radius (inclusive) in pixels.
rmax (int) – Maximum aperture radius (exclusive) in pixels.
noise ({'pipe','std','rms','astropy'}, optional) – Method used to estimate scatter for each aperture radius. Default is ‘pipe’ (uses utils.pipe_mad).
plot (bool, optional) – If True, produce a plot of the growth function (default False).
**kwargs (dict) – Additional keyword arguments forwarded to simple_aperture_photometry (for example sky annulus parameters).

Returns:

ap_flux_radii (ndarray) – Median aperture flux for each radius in the tested range.
noise_radii (ndarray) – Estimated noise metric (in ppm) for each radius.
fig (matplotlib.figure.Figure or None) – Figure of the growth function, if plot=True, otherwise None.
axs (ndarray or None) – Axes of the figure if produced, otherwise None.

identify_crays(clip=5, niters=5)[source]

Identify cosmic-ray affected pixels and update the bad-pixel map: by comparing each data frame with the median frame.

Parameters:

clip (float, optional) – Sigma threshold used to flag cosmic-ray candidates (default 5).
niters (int, optional) – Number of iterative passes to update self.badpix (default 5).

Returns:

badpix – Updated 2D bad-pixel mask (shape (Ny, Nx)) where flagged pixels are set to zero and good pixels are one.

Return type:

ndarray

noise_with_pcs(pc_max, noise='pipe', plot=False)[source]

Estimate noise (MAD or other) as a function of included PCs.

Parameters:

pc_max (int) – Maximum number of principal components to test (the method will evaluate values from 1 to pc_max-1).
noise ({'pipe','std','rms','astropy'}, optional) – Noise estimator to use (default ‘pipe’ which uses utils.pipe_mad).
plot (bool, optional) – If True, display a plot of noise metric vs number of PCs.

Returns:

pcs_to_include (ndarray) – Array of numbers of PCs tested (1..pc_max-1).
mad_with_pcs (ndarray) – Noise metric computed for each tested number of PCs.
fig (matplotlib.figure.Figure or None) – Figure object if plot=True, else None.
axs (matplotlib.axes.Axes or None) – Axis object for the plot if produced, else None.

pixel_level_decorrelation(removeNan=False)[source]

Prepare pixel-level decorrelation (PLD) basis functions.

The method builds the normalized pixel-level light curves (Phat) inside the chosen aperture and performs a PCA (principal componene analysis) on them to produce a set of orthogonal basis vectors (self.PCA) with singular values and eigenvectors stored on the instance.

Parameters:

removeNan (bool, optional) – If True, remove frames that contain NaNs in the basis before performing PCA.

Returns:

V (ndarray) – Matrix of eigenvectors from the PCA (shape: n_components, n_pixels).
eigenvalues (ndarray) – Eigenvalues corresponding to each PCA component.
PCA (ndarray) – Projected PCA time-series (shape: n_components, n_frames).

pld_correction(pc_max=5, extra_vectors=None, plot=False)[source]

Apply PLD correction using the first pc_max principal components.

Parameters:

pc_max (int, optional) – Number of PCA components to include in the PLD fit (default 5).
extra_vectors (ndarray or None, optional) – Additional vectors to include in the PLD fit (shape: either 1D array of length n_frames or 2D array of shape (n_extra, n_frames)). If provided, these will be appended to the PCA basis before fitting. Default is None (no extra vectors).
plot (bool, optional) – If True, produce a diagnostic plot comparing the raw and PLD-predicted flux (default False).

Returns:

Psum (ndarray) – Sum of pixel fluxes (SAP flux) before correction.
prediction (ndarray) – PLD model prediction for the flux (same length as Psum).
fig (matplotlib.figure.Figure or None) – Diagnostic figure if plot=True, otherwise None.
axs (ndarray or None) – Axes of the diagnostic figure if produced, otherwise None.

plot_correlation_matrices()[source]

Plot correlation matrices before and after PCA on pixel-level light curves.

Returns:

fig (matplotlib.figure.Figure) – Figure containing the two correlation matrix subplots.
axs (ndarray) – Array of two Axes objects for the before/after PCA plots.
im_list (list) – List with the two AxesImage objects (useful for colorbars).

plot_eigenvectors()[source]

Visualize the first 10 PCA eigenvectors as spatial maps.

Returns:

fig (matplotlib.figure.Figure) – Figure containing the eigenvector images.
axs (ndarray) – Array of Axes for each eigenvector subplot.

plot_median_image()[source]

This function plots the median image frame

Returns:

fig (matplotlib.figure.Figure or None) – Figure of median image.
im (matplotlib.image.AxesImage or None) – Image object from the median image plot.
axs (ndarray or None) – Axes of the figure.

plot_pcs(nmax=10)[source]

Plot the PC time-series (principal components) up to nmax.

Parameters:

nmax (int, optional) – Number of principal components to plot (default 10).

Returns:

fig (matplotlib.figure.Figure) – Figure containing the stacked PC plots.
axs (ndarray) – Array of Axes objects, one per plotted PC.

replace_nan(max_iter=50)[source]

Fill NaN entries in the data array by the mean of neighbouring pixels.

On each iteration the array is rolled along each axis to build a collection of nearest-neighbour shifts; the mean of these shifted arrays is used to replace NaNs. Iteration stops when there are no remaining NaNs or when max_iter is reached.

Parameters:: max_iter (int, optional) – Maximum number of fill iterations to perform (default 50).
Returns:: frames – The data array with NaNs replaced (also modified in-place and returned).
Return type:: ndarray

simple_aperture_photometry(method='manual', plot=False, bkg_corr='mean', robust_err=False, **kwargs)[source]

Compute aperture photometry for each frame in the data cube.

The function supports three methods: - manual: sums pixels inside a binary aperture mask and subtracts

the estimated sky background per pixel computed from a sky annulus.

photutils: uses photutils aperture classes to compute aperture sums and background statistics per frame (requires photutils to be installed). Additional keyword arguments are forwarded to the photutils helpers.
median: returns the per-frame median of the aperture pixels (useful as a robust estimator or for quick diagnostics).

Parameters:

method ({'manual','photutils','median'}, optional) – Photometry method to use. Default is 'manual'.
plot (bool, optional) – If True, produce a plot of aperture photometry, flux vs time (default False).
bkg_corr (str, optional) – The method to compute sky background flux, can be either mean or median. Default is mean.
robust_err (bool, optional) – Computes robust errors on aperture photometry by including scatter in sky background flux. Default is False.
**kwargs (dict) – Additional keyword arguments forwarded to photutils routines when method='photutils' (for example, custom aperture or statistic options).

Returns:

ap_flux (ndarray) – 1D array (length N_frames) with the background-subtracted aperture flux for each frame (same units as frames).
ap_flux_err (ndarray) – 1D array (length N_frames) with the estimated uncertainty on the aperture flux for each frame.
tot_sky_bkg (ndarray) – 1D array (length N_frames) with the total sky background contribution (summed over aperture pixels) subtracted from each frame. For methods that do not estimate a sky background this will be zeros.
fig (matplotlib.figure.Figure or None) – Figure of aperture photometry vs time, if plot=True, otherwise None.
axs (ndarray or None) – Axes of the figure if produced, otherwise None.

Notes

The function expects centroids to be available as self.cen_r and self.cen_c (one per frame) before being called.
Uses instance attributes self.aprad, self.sky_rad1, self.sky_rad2, self.brightpix, self.nos_brightest, and self.minmax to control aperture and background extraction.
When using photutils, this routine constructs CircularAperture/CircularAnnulus objects per frame and uses ApertureStats and aperture_photometry; errors are estimated combining aperture sum errors and annulus statistics following standard propagation.

sky_mask(plot=False)[source]

Create a sky-annulus mask (and masked frame/error arrays).

When self.brightpix is False this builds a circular annulus defined by self.sky_rad1 (inner radius) and self.sky_rad2 (outer radius) around the median centroid position. When self.brightpix is True the function instead selects the self.nos_brightest brightest pixels (optionally limited by self.minmax) to define the aperture, and rest of the pixels as background.

Parameters:

plot (bool, optional) – If True, produce a plot of median image with sky mask plotted on the top (default False).

Returns:

sky_fl_mask (ndarray) – self.frames multiplied by the 2D sky boolean mask (shape (N_frames, Ny, Nx)).
sky_err_mask (ndarray) – self.errors multiplied by the 2D sky boolean mask.
sky_bool_msk (ndarray) – 2D binary (0/1) mask describing the sky annulus on a single frame (shape (Ny, Nx)).
fig (matplotlib.figure.Figure or None) – Figure of median image and sky mask, if plot=True, otherwise None.
im (matplotlib.image.AxesImage or None) – Image object from the median image plot if plot=True, otherwise None.
axs (ndarray or None) – Axes of the figure if produced, otherwise None.

Notes

Uses instance attributes self.sky_rad1, self.sky_rad2, self.brightpix, self.nos_brightest, and self.minmax to control sky mask construction.

The ApPhoto class is wrapped around several other classes, CHEOPSData, TESSData, and KeplerData, that allows the users to download CHEOPS, TESS, and Kepler data directly from command lines and use ApPhoto class to perform aperture photometry on these datasets. These classes have the following methods:

class krithika.CHEOPSData(object_name=None, pipe_filename=None)[source]

Manipulate CHEOPS time-series data from simple reading to extraction.

This class provides tools to read, process, and extract photometry from CHEOPS (CHaracterising ExOPlanet Satellite) observations. It supports reading PIPE-produced light curves, downloading data reduction pipeline (DRP) products, handling sub-array data, and performing aperture photometry on subarrays.

Parameters:

object_name (str, optional) – Name or identifier of the target object. Used for querying the CHEOPS database when downloading data products. Default is None.
pipe_filename (str, optional) – Path to a pipeline-produced FITS file for reading with pipe_data(). Default is None.

pipe_fname

Stored pipeline FITS filename.

Type:: str or None

object_name

Stored target object name.

Type:: str or None

drp_data

Dictionary containing DRP light curve data populated by get_drp_lightcurves() with keys for each instrument (e.g., ‘TG000001’), each containing ‘TIME’, ‘FLUX’, ‘FLUX_ERR’, ‘ROLL’, ‘BG’, ‘CONTAM’, ‘SMEAR’, ‘XC’, ‘YC’.

Type:: dict

time_sub, flux_sub, flux_err_sub, badpix_sub

Dictionaries containing sub-array data populated by get_subarrays().

Type:: dict

subarray_data

Dictionary that may be populated with processed subarray data.

Type:: dict, optional

Examples

Reading PIPE data:

>>> cheops = CHEOPSData(pipe_filename='path/to/file.fits')
>>> data = cheops.pipe_data(bgmin=300)

Downloading DRP light curves:

>>> cheops = CHEOPSData(object_name='WASP-189b')
>>> drp_data = cheops.get_drp_lightcurves(pout='/path/to/output')

Notes

Requires the dace_query library to download CHEOPS data products from the DACE (Data and Analysis Center for Exoplanets) database.

get_drp_lightcurves(pout='/home/docs/checkouts/readthedocs.org/user_builds/krithika/checkouts/latest/docs', load=False, save=False, aperture='default', visit_nos=None, filekey=None, bkg_clip=20)[source]

Download or load Data Reduction Pipeline light curves from CHEOPS.

This method queries the CHEOPS DACE database for light curve products (DRP level) for the target specified in self.object_name. It can download a specific visit (by visit_nos or filekey) or all available visits. Downloaded data are filtered to remove flagged points, high-background measurements, and cosmic rays (via EVENT/STATUS flags).

Flux is normalized by its median value per instrument.

Parameters:

pout (str, optional) – Output directory for downloaded/saved products. Created if it does not exist. Default is the current working directory.
load (bool, optional) – If True, skip downloading and instead load pre-downloaded files from pout. Useful for re-reading data. Default is False.
save (bool, optional) – If True, retain downloaded FITS files in pout after extraction. If False (default), files are deleted after reading.
aperture (str, optional) – Aperture type to download (‘default’, ‘optimal’, ‘rsup’, and ‘rinf). Default is 'default'.
visit_nos (int or None, optional) – Specific visit number to download (e.g., 701, 101, etc.). If provided, only that visit’s data are retrieved. Overridden by filekey. Default is None.
filekey (str or None, optional) – Specific file key (e.g., ‘CH_PR100009_TG000001_…’) for precise product selection. Takes precedence over visit_nos. Default is None.
bkg_clip (float, optional) – Sigma threshold for background clipping. Points with background above median(BG) + bkg_clip * MAD(BG) are removed. Default is 20.

Returns:

drp_data –

Nested dictionary with structure:

{instrument_name: {
    'TIME': ndarray,
    'FLUX': ndarray,
    'FLUX_ERR': ndarray,
    'ROLL': ndarray,
    'BG': ndarray,
    'CONTAM': ndarray,  # Contamination from nearby sources
    'SMEAR': ndarray,   # Smearing correction
    'XC', 'YC': ndarray # Centroids
}}

Also stored in self.drp_data.

Return type:

dict

Notes

Data filtering steps: 1. Remove non-finite flux values 2. Remove points with EVENT != 0 or STATUS != 0 3. Remove high-background points (sigma clipping)

The method creates a tar archive internally when downloading multiple files, which is extracted to pout then deleted.

Examples

Download all DRP light curves for a target:

>>> cheops = CHEOPSData(object_name='WASP-189b')
>>> drp_data = cheops.get_drp_lightcurves(pout='./cheops_data')

Download a specific visit:

>>> drp_data = cheops.get_drp_lightcurves(visit_nos=1, pout='./data')

Load previously downloaded data:

>>> drp_data = cheops.get_drp_lightcurves(load=True, pout='./data')

get_subarrays(pout='/home/docs/checkouts/readthedocs.org/user_builds/krithika/checkouts/latest/docs', download_all=False, load=False, visit_nos=None, filekey=None)[source]

Download or load CHEOPS sub-array image time-series data.

This method queries and downloads sub-array FITS files from CHEOPS data products, which contain the flux measurements as 2D image cubes for each frame in the time series. These are useful for custom photometry and pixel-level analysis.

Sub-array data are stored in self.time_sub, self.flux_sub, self.flux_err_sub, and self.badpix_sub.

Parameters:

pout (str, optional) – Output directory for downloaded/saved products. Created if needed. Default is the current working directory.
download_all (bool, optional) – If True, download all available file types (not just sub-arrays). If False, filter for sub-array products only. Default is False.
load (bool, optional) – If True, skip downloading and load pre-downloaded files from pout. Default is False.
visit_nos (int or None, optional) – Specific visit number to download. If provided with filekey=None, only that visit’s sub-arrays are retrieved. Default is None.
filekey (str or None, optional) – Specific file key for precise product selection. Takes precedence over visit_nos. Default is None.

Returns:

Results are stored in instance dictionaries:

self.time_sub[instrument]ndarray
Time stamps (shape: N_frames).
self.flux_sub[instrument]ndarray
Flux images (shape: N_frames, Ny, Nx).
self.flux_err_sub[instrument]ndarray
Flux error (shape: N_frames, Ny, Nx)
self.badpix_sub[instrument]ndarray
Bad-pixel mask (shape: Ny, Nx) where 1 = good, 0 = bad.

Return type:

None

Notes

The method attempts to load a bad-pixel map from the CHEOPS data products. If not found, a default all-good mask is created.

Sub-array file names are identified by the key ‘COR_SubArray’ in the filename.

Examples

Download all sub-arrays for a target:

>>> cheops = CHEOPSData(object_name='WASP-189b')
>>> cheops.get_subarrays(pout='./cheops_data')

Load previously downloaded sub-array data:

>>> cheops.get_subarrays(load=True, pout='./data')

pipe_data(bgmin=300, not_flg_arr=None)[source]

Extract and process light curve data from a PIPE FITS file.

This method reads a CHEOPS pipeline-produced FITS file (typically from self.pipe_fname) and extracts time series photometry along with auxiliary data such as centroids, background, and temperature. Data points are filtered by quality flags and background levels.

The flux is normalized by its median value.

Parameters:

bgmin (float, optional) – Minimum background level threshold in counts. Data points with background below this value are excluded. Default is 300.
not_flg_arr (array_like of int, optional) – Array of FLAG values to include in the output (i.e., non-zero flags are normally rejected, but values in this array are kept). If None, only FLAG==0 data are retained. Default is None.

Returns:

data – Dictionary containing extracted and processed data with keys:

’TIME’ndarray
Barycentric Julian Date time stamps.
’FLUX’ndarray
Normalized flux (normalized by median).
’FLUX_ERR’ndarray
Flux uncertainties (normalized by median).
’ROLL’ndarray
Roll angle in degrees.
’XC’, ‘YC’ndarray
Centroid positions in pixels.
’BG’ndarray
Background level in counts.
’TF2’ndarray
thermFront2 temperature in Kelvin.
’U1’ … ‘U16’ndarray, optional
Principal components from PCA.

Return type:

dict

Notes

Quality flags (FLAG column): - 0 : Good data - >0 : Data with issues (normally rejected)

The method performs two masking steps: 1. Removes flagged data and optionally includes specified flags 2. Removes high-background points (BG < bgmin)

Examples

>>> cheops = CHEOPSData(pipe_filename='file.fits')
>>> data = cheops.pipe_data(bgmin=300)
>>> print(data['TIME'], data['FLUX'])

class krithika.TESSData(object_name)[source]

Helper to query and fetch TESS products from MAST and juliet.

This lightweight utility class provides two convenience methods: get_lightcurves which delegates to juliet.utils.get_all_TESS_data (when available) to gather per-instrument light curves, and get_tpfs which queries MAST for Target Pixel Files (TPFs), downloads the products and populates simple dictionaries on the instance with the extracted time and pixel-cube data.

Parameters:: object_name (str) – Target identifier accepted by MAST (for example a TIC or common object name).

object_name

The provided target name.

Type:: str

tim_lc, fl_lc, fle_lc

Time/flux/flux-error collections populated by get_lightcurves when that method is called.

Type:: dict or None

out_dict_lc

Optional output dictionary returned by the juliet helper.

Type:: object or None

tim_tpf, fl_tpf, fle_tpf, badpix, quality

Dictionaries populated by get_tpfs() containing per-sector time arrays, pixel flux cubes, pixel flux error cubes, bad-pixel masks and quality arrays respectively.

Type:: dict

get_lightcurves(pdc=True, save=False, pout='/home/docs/checkouts/readthedocs.org/user_builds/krithika/checkouts/latest/docs', **kwargs)[source]

Collect TESS light curves for self.object_name.

This method attempts to use juliet.utils.get_all_TESS_data to retrieve light curves. The juliet helper may return either a simple tuple (tim, fl, fle) or an (out_dict, tim, fl, fle) tuple; both cases are handled and stored on the instance.

Parameters:

pdc (bool, optional) – If True request PDC-corrected fluxes when available. Default is True.
save (bool, optional) – If True the method will write ASCII files named LC_<object>_<inst>.dat to pout. Default False.
pout (str, optional) – Output directory used when save=True (defaults to the current working directory).
**kwargs (dict) – Forwarded to juliet.utils.get_all_TESS_data.

Returns:

tim_lc, fl_lc, fle_lc, out_dict_lc – Stored time/flux/fluxerr containers and an optional output dictionary (out_dict_lc may be None).

Return type:

tuple

get_tpfs(pout='/home/docs/checkouts/readthedocs.org/user_builds/krithika/checkouts/latest/docs', load=False, save=False, savefits=False)[source]

Download or load TESS Target Pixel Files (TPFs) for the target.

If load is False the method queries MAST for TESS time-series products for self.object_name, filters suitable products (TPFs), downloads them, reads the FITS tables and populates the instance dictionaries tim_tpf, fl_tpf, fle_tpf, quality and badpix. If save is True a small pickled dictionary per sector is written to pout. load=True simply reads these previously saved dictionaries.

Parameters:

pout (str, optional) – Output directory used when save=True. Default is the current working directory.
load (bool, optional) – If True, attempt to load previously saved TPF pickles instead of querying/downloading from MAST. (Loading logic is not implemented in this helper; calling with load=True will bypass the MAST query branch.)
save (bool, optional) – If True save extracted per-sector dictionaries as TPF_<object>_<sector>.pkl in pout. Default False.
savefits (bool, optional) – If True, save the downloaded fits files in pout folder. Default is False.

Returns:

tim_tpf, fl_tpf, fle_tpf, badpix, quality – The dictionaries populated on the instance containing the per-sector time arrays, pixel flux cubes, flux error cubes, quality arrays and a placeholder bad-pixel mask.

Return type:

tuple

class krithika.KeplerData(object_name)[source]

get_lightcurves(pdc=True, long_cadence=True, verbose=True, save=False, pout='/home/docs/checkouts/readthedocs.org/user_builds/krithika/checkouts/latest/docs')[source]

Collect Kepler/K2 light curves for self.object_name.

This method uses astroquery to retrieve light curves. The function returns a simple tuple (tim, fl, fle).

Parameters:

pdc (bool, optional) – If True request PDC-corrected fluxes when available. Default is True.
long_cadence (bool, optional) – If True the long cadence data will be downloaded. Default if False.
verbose (bool) – Boolean on whether to print updates. Default is True
save (bool, optional) – If True the method will write ASCII files named LC_<object>_<inst>.dat to pout. Default False.
pout (str, optional) – Output directory used when save=True (defaults to the current working directory).

get_tpfs(long_cadence=True, load=False, verbose=True, save=False, pout='/home/docs/checkouts/readthedocs.org/user_builds/krithika/checkouts/latest/docs', savefits=False)[source]

Collect Kepler/K2 target pixel files for self.object_name.

This method uses astroquery to retrieve target pixel files.

Parameters:

load (bool, optional) – If True request PDC-corrected fluxes when available. Default is True.
long_cadence (bool, optional) – If True the long cadence data will be downloaded. Default if False.
verbose (bool) – Boolean on whether to print updates. Default is True
save (bool, optional) – If True the method will write ASCII files named LC_<object>_<inst>.dat to pout. Default False.
pout (str, optional) – Output directory used when save=True (defaults to the current working directory).

If you used juliet for fitting the light curves, you can use the julietPlots class to plot the light curves and the fits. It has the following methods:

class krithika.julietPlots(input_folder=None, dataset=None, res=None, N=1000, **kwargs)[source]

Plotting helper for results produced by a juliet analysis.

This class loads a juliet input/output folder, computes models from the fitted posterior and provides a collection of convenience plotting and data-processing methods for light-curve visualisation. Typical usage is to instantiate the class with the input_folder used for a juliet fit and then call plotting helpers such as full_model_lc(), plot_gp() or detrend_data().

input_folder

Path to the juliet input/output folder provided at construction. The user need to provide either input_folder or dataset and res

Type:: str, optional

dataset

The object returned by juliet.load for the given input folder.

Type:: juliet.load object, optioanl

res

The result of calling dataset.fit(...).

Type:: juliet.fit, optional

N

Number of posterior samples drawn when evaluating sampled models.

Type:: int

\*\*kwargs

Any additional keywords provided to juliet.fit object.

Type:: dict, optional

Main methods

------------

full_model_lc(instruments=None, save=False, nrandom=50, quantile_models=True)[source]: Plot observed light curves with the full fitted model and residuals.

phase_folded_lc(phmin=0.8, instruments=None, highres=False, ...)[source]: Plot phase folder light curve with best-fitted planetary model for one or more instruments.

plot_gp(instruments=None, highres=False, one_plot=False, ...)[source]: Plot GP components and binned residuals for one or more instruments.

plot_fake_allan_deviation(instruments=None, binmax=10, method='pipe', timeunit=None)[source]: Plot “allan deviation” plot, which is noise as a function of binning, for one or more instruments

plot_corner(planet_only=False, save=True)[source]: Plot corner plots for fitted posterior samples, can use planet-only or all parameters.

compute_gp_model(instrument=None, highres=None, times=None, pred_time=None, resids=None, errors=None, parameter_vector=None)[source]: This function computes the GP models. It accesses the GP object from the juliet files, and then do GP.compute and GP.predict.

detrend_data(phmin, instruments, planet='p1')[source]: This function will generate detrended data from the raw data. It will also generate phases.

full_model_lc(instruments=None, save=False, nrandom=50, quantile_models=True)[source]

Plot the full fitted light-curve model for one or more instruments.

This method creates a figure per instrument showing the observed light curve with the fitted (full) model overplotted and a residuals panel beneath.

Parameters:

instruments (list or None) – List of instrument names to plot. If None, all instruments in the current juliet dataset plotted.
save (bool, optional) – If True, each figure is saved to <input_folder>/full_model_<instrument>.png. Default False.
nrandom (int, optional) – Number of random posterior-sample models to draw when quantile_models is False. Default is 50.
quantile_models (bool, optional) – If True, plot the 68%% credible interval as a filled band; if False, plot nrandom random posterior models instead.

Returns:

all_fig (list) – List of matplotlib Figure objects created (one per instrument).
all_axs1 (list) – List of top-panel Axes (data + model) for each figure.
all_axs2 (list) – List of bottom-panel Axes (residuals) for each figure.

Notes

The routine expects that self.models_all_ins was populated by a prior call (see the class constructor) and that the dataset arrays (times, data, errors) are available in self.dataset.

highres_planet_models_workaround(instrument, times_highres, GP_reg_highres, lin_reg_highres)[source]: The usual way of generating high time-resolution models is not working with juliet when light_travel_delay=True for occultations and phase curve modelling. This is a workaround for this.

phase_folded_lc(phmin=0.8, instruments=None, planet='p1', highres=False, nrandom=50, quantile_models=True, one_plot=None, figsize=(10.666666666666666, 6.0), pycheops_binning=False, nos_bin_tra=20, nos_bin_pc=30)[source]

Plot phase-folded light curves and models for the fitted dataset.

This method computes detrended data and planetary models by calling detrend_data and detrend_model internally, detects whether the fit contains transits, eclipses or phase-curves, and then produces matplotlib figures showing the data, model, and residuals. It can produce either one figure per instrument or a single combined figure with all instruments plotted together.

Parameters:

phmin (float, optional) – Minimum phase value used for generating orbital phase (default 0.8).
instruments (list or None, optional) – List of instrument names to plot. If None, all instruments in the juliet dataset (self.dataset.inames_lc) are used. One plot is produced if all instruments are provided; otherwise, one plot per instrument is created.
planet (str, optional) – Planet identifier used in the juliet fit (e.g., ‘p1’, ‘p2’, etc.) to select the appropriate planet model. Default is ‘p1’.
highres (bool, optional) – If True, compute and plot high time-resolution planet-only models. Default False.
nrandom (int, optional) – Number of random posterior-sample models to draw when quantile_models is False. Default 50.
quantile_models (bool, optional) – If True, display the 68%% credible interval as a filled band; if False, overplot nrandom random posterior samples. Default True.
one_plot (bool or None, optional) – If True, produce a single combined plot for all instruments (regardless of the number of instruments provided). If False, produce one plot per instrument (unless all instruments are provided). If None, the behaviour is as just described. Default None.
figsize (tuple, optional) – Figure size passed to matplotlib when creating figures. Default is (16/1.5, 9/1.5).
pycheops_binning (bool, optional) – If True, use binning as produced by pycheops the binned datapoints. If False, it will use default juliet binning. Default False.
nos_bin_tra (int, optional) – Number of total binned data points in transit plot. Default is 20.
nos_bin_pc (int, optional) – Number of total binned data points in phase curve plot. Default is 30.

Returns:

If a single combined plot is produced the function returns the
combined fig and the primary axes used. If multiple figures are
generated, the function returns lists (``figs_all, axs1_all,)
axs2_all, axs3_all, axs4_all`` corresponding to created figures and
their axes panels.

plot_corner(planet_only=False, param_list=None, save=True)[source]

Create a corner plot of selected posterior parameters.

Parameters:

planet_only (bool, optional) – If True, include only planetary parameters in the corner plot. If False (default), include instrumental and noise parameters as well when available.
param_list (list of str or None, optional) – If provided, a list of parameter names to include in the corner plot. If None (default), all parameters are included (subject to the planet_only filter).
save (bool, optional) – If True (default), save the generated figure to <input_folder>/corner_plot.png.

Returns:

fig – The figure produced by utils.corner_plot containing the marginal and pairwise plots for the selected parameters.

Return type:

matplotlib.figure.Figure

plot_fake_allan_deviation(instruments=None, binmax=10, method='pipe', timeunit=None)[source]

Compute and plot per-instrument noise-vs-bin-size curves.

This is a thin wrapper around utils.fake_allan_deviation that computes residuals for each instrument (data minus the full model) and returns the figures, axes and numeric results for further use.

Parameters:

instruments (list or None) – List of instrument names to include. If None, all instruments in the juliet dataset (self.dataset.inames_lc) are processed.
binmax (int, optional) – Passed to utils.fake_allan_deviation to control the maximum number of bins (default 10).
method ({'pipe','std','rms','astropy'}, optional) – Method used to estimate the scatter of the binned residuals: - ‘pipe’ : use pipe_mad (median absolute differences estimator) - ‘std’ : use np.nanstd - ‘rms’ : use the root-mean-square (rms helper) - ‘astropy’ : use astropy.stats.mad_std (robust MAD-based std)
timeunit (str or None, optional) – If provided, forwarded to utils.fake_allan_deviation to force the secondary x-axis time unit; otherwise the helper chooses an appropriate unit automatically.

Returns:

figs_all (list) – List of matplotlib Figure objects (one per instrument).
axs_all (list) – List of matplotlib Axes objects corresponding to the figures.
binsize_all (list) – List of binsize arrays returned for each instrument.
noise_all (list) – List of noise arrays (ppm) computed for each instrument.
white_noise_all (list) – List of white-noise expectation arrays (ppm) for each instrument.

plot_gp(instruments=None, highres=False, one_plot=False, figsize=(10.666666666666666, 6.0), pycheops_binning=False, xlabel='GP regressor', robust_errors=False, plot_errorbars=False, parameter_vector=None, nos_bin=20)[source]

Plot the Gaussian Process component(s) from the fitted juliet model.

This routine computes and plots the GP model (and binned data) for one or more instruments present in the current juliet dataset.

Parameters:

instruments (list or None) – List of instrument names to plot. If None, all instruments for which GP regressors exist in self.dataset.GP_lc_arguments are plotted.
highres (bool, optional) – If True, compute a high-resolution GP prediction (calls compute_gp_model). Default is False.
one_plot (bool, optional) – If True, combine all instruments into a single figure and return a single (fig, axs). If False, return lists of figures and axes (one per instrument). Default False.
figsize (tuple, optional) – Figure size passed to matplotlib. Default (16/1.5, 9/1.5).
pycheops_binning (bool, optional) – Use pycheops produced binning when set to True. Default False.
xlabel (str, optional) – X-axis label for the GP regressor. Default 'GP regressor'.
robust_errors (bool, optional) – If True, compute GP uncertainties using compute_gp_model instead of using juliet precomputed errors.
plot_errorbars (bool, optional) – If True, plot error bars on the GP data points. Default False.
parameter_vector (array-like or None, optional) – If provided, passed to the GP parameter setter before computing predictions (useful for evaluating GP at different parameter values).
nos_bin (int, optional) – Approximate number of bins to use when computing binned data per instrument. Default 20.

Returns:

If one_plot is True, returns a single tuple (fig, axs) for the combined figure. Otherwise returns (fig_all, axs_all) where each is a list containing one Figure/Axis per instrument.

Return type:

(fig, axs) or (fig_all, axs_all)

The SpectroscopicLC class can be used to perform analysis of spectroscopic light curves from, e.g., JWST. It has the following methods:

class krithika.SpectroscopicLC(times, lc, lc_errs, wavelengths, gp_pars=None, lin_pars=None, priors=None, pout='/home/docs/checkouts/readthedocs.org/user_builds/krithika/checkouts/latest/docs')[source]

Analyze and fit multi-wavelength time-series spectroscopic data.

This class provides tools for handling spectroscopic light curves where flux is measured as a function of both time and wavelength. It supports binning spectral data into channels, fitting individual wavelength bins using juliet, parallel processing of multiple channels, and visualization of spectral transit/eclipse signatures and systematics.

Parameters:

times (ndarray) – 1D array of time stamps (in BJD or similar) with shape (N_frames,).
lc (ndarray) – 2D array of flux measurements with shape (N_frames, N_wavelengths). Each row is a time frame; each column is a wavelength bin.
lc_errs (ndarray) – 2D array of flux uncertainties with the same shape as lc.
wavelengths (ndarray) – 1D array of wavelength values (in microns or other units) corresponding to each column of lc. Shape is (N_wavelengths,).
gp_pars (dict or None, optional) – Dictionary mapping instrument names to arrays of Gaussian Process (GP) regressor values (e.g., time, roll angle). Used when fitting with GP systematics. Default is None.
lin_pars (dict or None, optional) – Dictionary mapping instrument names to 2D arrays of linear regressor values (shape: N_frames x N_regressors). Used for linear detrending. Default is None.
priors (callable or None, optional) – A callable that accepts a channel name (string) and returns a juliet-format priors for that channel. Required for fitting via analyse_lc_parallel(). Default is None.
pout (str, optional) – Output directory for saving results, intermediate data, and pickled light curves. Default is the current working directory.

times

Input time array.

Type:: ndarray

lc

Input 2D light-curve array.

Type:: ndarray

lc_errs

Input 2D error array.

Type:: ndarray

wavelengths

Input wavelength array.

Type:: ndarray

gp_pars, lin_pars

Input systematics parameters.

Type:: dict or None

priors

Input priors function.

Type:: callable or None

pout

Output directory.

Type:: str

col_st, col_end

Start and end column indices for each binned channel (populated by col_spec()).

Type:: ndarray

spectral_lcs_dict

Dictionary containing binned spectral light curves and metadata, populated by generating_lightcurves(). Keys include ‘lc’, ‘err’, ‘wave’, and ‘wave_bin’.

Type:: dict

wav, wav_bin

Wavelength array and half-widths for each bin (populated by plot_parameter_spectrum()).

Type:: ndarray

pars_med, pars_up, pars_lo

Median and upper/lower credible interval bounds on fitted parameters as a function of wavelength (populated by plot_parameter_spectrum()).

Type:: ndarray

qua_white

Quantiles (median, upper, lower) of the white-light parameter (populated by plot_parameter_spectrum()).

Type:: tuple

Main methods

------------

generating_lightcurves(ch_nos=None)[source]: Extract and optionally bin spectral light curves into channels.

bin_lightcurve(unbinned_lc, unbinned_lc_err)[source]: Inverse-variance weighted binning of light curves along wavelength direction.

white_light_lc[source]: Compute the white-light (all-wavelength) light curve.

analyse_lc_parallel(nthreads, ch_nos=None, juliet_fit_kwargs={})[source]: End-to-end parallel analysis: bin, organize, and fit all channels.

plot_parameter_spectrum(parameter, bins=None, post_size=10000, ...)[source]: Plot fitted parameter values as a function of wavelength.

plot2Ddata(cmap='plasma')[source]: 2D heatmap of spectral light curves vs time.

plot2D_data_model_resids(cmap='plasma', detrend=False, ...)[source]: Side-by-side 2D heatmaps of data, model, and residuals.

joint_fake_allan_deviation(binmax=10, method='pipe', timeunit=None)[source]: Combined noise vs binning plot across all channels.

Notes

All fitting is performed using the juliet package (Espinoza et al. 2019).
Parallel processing uses Python’s multiprocessing.Pool.
Results are saved as pickled dictionaries and ASCII text files to pout.
Methods check for existing output files and skip re-computation if files are already present (useful for resuming interrupted analyses).

Examples

Create and analyze a spectroscopic dataset:

>>> import numpy as np
>>> times = np.linspace(2459000, 2459001, 100)
>>> wavelengths = np.linspace(1.1, 1.7, 50)
>>> lc = np.random.normal(1.0, 0.01, (100, 50))
>>> lc_errs = np.full_like(lc, 0.005)

>>> def get_priors(ch_name):
...     return juliet.utils.priors()

>>> spec = SpectroscopicLC(times, lc, lc_errs, wavelengths,
...                        priors=get_priors, pout='./results')
>>> spec.analyse_lc_parallel(nthreads=4, ch_nos=10)
>>> fig, ax = spec.plot_parameter_spectrum('p_p1', bins=10, plot_white=True)

analyse_lc_parallel(nthreads, ch_nos=None, juliet_fit_kwargs={})[source]

End-to-end parallel analysis: bin, organize, and fit lightcurves from all channels.

This high-level wrapper orchestrates the full analysis pipeline: 1. Bins the spectral light curve into channels via generating_lightcurves() 2. Organizes the data into a format suitable for fitting 3. Calls multi_fit_lcs() to fit all channels in parallel

Parameters:

nthreads (int) – Number of parallel threads to use for fitting.
ch_nos (int or None, optional) – Number of wavelength channels to bin into. If None, use native wavelength resolution. Default is None.
juliet_fit_kwargs (dict, optional) – Keyword arguments to forward to juliet.dataset.fit(...), such as sampler='dynamic_dynesty', nthreads=8, etc. Default is an empty dictionary.

Returns:

Results are written to pout and its subdirectories.

Return type:

None

Notes

Requires that self.priors is a callable that returns a priors for each channel name.

bin_lightcurve(unbinned_lc, unbinned_lc_err)[source]

Perform inverse-variance weighted binning of a light curve.

Each output bin is the weighted average of input values, where weights are the inverse squares of the input uncertainties. Output uncertainties are computed by propagating input uncertainties through the weighted average formula.

Parameters:

unbinned_lc (ndarray) – 2D array of unbinned flux values (shape: N_frames x N_wavelengths).
unbinned_lc_err (ndarray) – 2D array of flux uncertainties (same shape as unbinned_lc).

Returns:

bin_lc (ndarray) – 1D array of binned flux values (shape: N_frames,).
bin_lc_err (ndarray) – 1D array of binned flux uncertainties (same shape as bin_lc).

col_spec(ch_nos)[source]: Given a number of total number of channels, this function gives an array containing start and end column for each channel

fit_lc(lightcurves, ch_name)[source]: Inteernal helper function for multiprocessing. This function fits lightcurve for one spectroscopic channel

generating_lightcurves(ch_nos=None)[source]

Extract and optionally bin spectral light curves into channels.

This method either bins the original light curve into ch_nos channels (using inverse-variance weighting) or uses the original wavelength resolution. Results are saved as a pickled dictionary and stored in self.spectral_lcs_dict.

If a pickle file for the requested binning already exists, it is loaded instead of recomputing.

Parameters:

ch_nos (int or None, optional) – Number of channels to bin the light curve into. If None, use native wavelength resolution (no binning). Default is None.

Returns:

Results are stored in self.spectral_lcs_dict with keys:

’lc’ndarray
Binned light curves (shape: N_frames x N_channels).
’err’ndarray
Binned flux errors (same shape as ‘lc’).
’wave’ndarray
Central wavelength of each bin.
’wave_bin’ndarray
Half-width (or full width) of each wavelength bin.

Return type:

None

Notes

The output is pickled to pout/spectroscopic_lc_ch_<ch_nos>.pkl for future loading.

joint_fake_allan_deviation(binmax=10, method='pipe', timeunit=None)[source]

Plot combined noise-vs-binning curves across all spectroscopic channels.

Computes the “Allan deviation” (actually the MAD-based noise as a function of binning) for each channel’s residuals and overlays all curves on a single plot. Also shows the white-noise expectation curve.

Parameters:

binmax (int, optional) – Maximum number of bins to test. Default is 10.
method ({'pipe', 'std', 'rms', 'astropy'}, optional) – Noise estimator to use (passed to utils.fake_allan_deviation()). Default is 'pipe' (MAD-based).
timeunit ({'d', 'hr', 'min'} or None, optional) – Time unit for the secondary x-axis. If None, chosen automatically based on data span. Default is None.

Returns:

fig (matplotlib.figure.Figure) – The created figure.
axs (matplotlib.axes.Axes) – The plot axes (with primary x-axis for bin size and secondary x-axis for time period).

plot2D_data_model_resids(cmap='plasma', detrend=False, binwidth=None, ppm=False, **kwargs)[source]

Create side-by-side 2D heatmaps of data, model, and residuals.

Loads the fitted models and residuals from all channels and creates three synchronized heatmaps showing the full spectro-temporal evolution of the data, fit quality, and systematics.

Parameters:

cmap (str, optional) – Matplotlib colormap name. Default is 'plasma'.
detrend (bool, optional) – If True, load and plot detrended data (subtracting GP and linear systematics). Requires additional juliet processing and may be slow. If False, plot raw data. Default is False.
binwidth (float or None, optional) – Time binwidth (in days) for additional temporal binning of the data/model/residuals. If None, no binning is performed. Default is None.
ppm (bool, optional) – If True, express data and model in ppm relative to unity. Residuals are always in ppm. Default is False.
**kwargs (dict) – Forwarded to the julietPlots constructor when detrend=True (e.g., N=1000 for the number of posterior samples).

Returns:

fig (matplotlib.figure.Figure) – The created figure with three subplots.
axs (ndarray) – Array of three Axes objects (data, model, residuals).
im_list (list) – List of three AxesImage objects (useful for colorbars).

plot2Ddata(cmap='plasma')[source]

Create a 2D heatmap of spectral light curves (flux vs time) vs wavelength.

Parameters:

cmap (str, optional) – Matplotlib colormap name. Default is 'plasma'.

Returns:

fig (matplotlib.figure.Figure) – The created figure.
ax1 (matplotlib.axes.Axes) – The plot axes.
im1 (matplotlib.image.AxesImage) – The image object.
cbar (matplotlib.colorbar.Colorbar) – The colorbar object.

plot_parameter_spectrum(parameter, bins=None, post_size=10000, bin_method='mean', ppm=False, plot_white=False)[source]

Plot a fitted parameter as a function of wavelength (a spectrum).

Loads posterior samples for a parameter from all fitted channels, optionally re-bins them in the wavelength direction, and plots median with credible intervals as a function of wavelength.

Parameters:

parameter (str) – Name of the parameter to plot (e.g., ‘p_p1’, ‘fp_p1’, ‘C1_p1’). Should match the parameter naming convention used by juliet.
bins (int or None, optional) – Number of wavelength bins to use for the output spectrum. If None or equal to the native number of channels, no re-binning is performed. Default is None.
post_size (int, optional) – Number of posterior samples to draw from each channel. Default is 10000.
bin_method ({'mean', 'median', 'weighted_average'}, optional) – Method for combining posteriors across wavelength bins. Default is 'mean'.
ppm (bool, optional) – If True, scale parameter values to ppm (multiply by 1e6). Useful for small parameters like occultation depths. Default is False.
plot_white (bool, optional) – If True, overplot the white-light (all-wavelength) parameter value as a horizontal band with credible interval. Default is False.

Returns:

fig (matplotlib.figure.Figure) – The created figure.
axs (matplotlib.axes.Axes) – The plot axes.

Notes

Posterior samples are loaded from files named like pout/CH<i>/*_posteriors.pkl.

white_light_lc()[source]

Compute the white-light (all-wavelength) light curve.

Bins the entire 2D light curve (all wavelengths combined) into a single 1D time-series using inverse-variance weighting.

Returns:

white_lc (ndarray) – 1D array of white-light flux values.
white_lc_err (ndarray) – 1D array of white-light flux uncertainties.

The InvertCowanAgolPC class can be used to invert the Cowan & Agol (2008) phase curve model to get the thermal properties such as the brightness temperature map and Bond albedo. It has the following methods:

class krithika.InvertCowanAgolPC(E, C1, D1, C2, D2, rprs, bandpass=None, method='Nelder-Mead', teff_star=None, stellar_spec=None, nthreads=2, ntheta=49, nphi=100, pout='/home/docs/checkouts/readthedocs.org/user_builds/krithika/checkouts/latest/docs')[source]

Invert phase curve observations to physical properties using Cowan & Agol (2008) model.

This class implements tools to analyze phase-curve observations, computing thermal properties (temperature maps, day/night temperatures, Bond albedo, heat re-distribution efficiencies) and orbital properties (hotspot offset, phase offset) using the parameterized phase-curve model of Cowan & Agol (2008).

The class leverages the model’s Fourier expansion coefficients (A0, A1, B1, A2, B2) derived from eclipse depths (E) and phase-curve coefficients (C1, D1, C2, D2) to construct 2D thermal maps and compute various physical properties across posterior samples for uncertainty quantification.

Parameters:

E (ndarray) – Array of eclipse depths (F_p/F_*) from posterior samples.
C1 (ndarray) – Phase-curve cosine coefficient (first harmonic) from posterior samples.
D1 (ndarray) – Phase-curve sine coefficient (first harmonic) from posterior samples.
C2 (ndarray) – Phase-curve cosine coefficient (second harmonic) from posterior samples.
D2 (ndarray) – Phase-curve sine coefficient (second harmonic) from posterior samples.
rprs (float or ndarray) – Planet-to-star radius ratio. If scalar, same value used for all samples; if array, sampled from for each posterior sample.
bandpass (dict, optional) – Dictionary containing instrument bandpass data with keys ‘WAVE’ and ‘RESPONSE’. Required for computing brightness temperatures. Values must be Astropy Quantities. Default is None.
method (str, optional) – Optimization method for fitting temperatures to flux (passed to scipy.optimize.minimize). Default is 'Nelder-Mead'.
teff_star (float (astropy.units.Quantity) or None, optional) – Stellar effective temperature in Kelvin. Required for computing Bond albedo, heat redistribution efficiency, and average day/night temperatures. If not provided, a warning is printed. Default is None.
stellar_spec (dict or None, optional) – Dictionary containing stellar spectrum with keys ‘WAVE’ and ‘FLUX’. Used to interpolate stellar flux at bandpass wavelengths. If provided, overrides teff_star-based blackbody spectrum. Values must be Astropy Quantities. Default is None.
nthreads (int, optional) – Number of CPU threads to use for parallel processing of temperature maps and multiprocessing operations. Default is the total number of available CPU cores (multiprocessing.cpu_count()).
ntheta (int, optional) – Number of grid points across latitude (from -π/2 to π/2). Default is 49.
nphi (int, optional) – Number of grid points across longitude (from -π to π). Default is 100.
pout (str, optional) – Output directory for saving intermediate results (temperature maps, brightness temperatures, etc.). Default is the current working directory.

E, C1, D1, C2, D2

Stored input parameters (eclipse depths and phase-curve coefficients).

Type:: ndarray

rprs

Stored planet-to-star radius ratio.

Type:: float or ndarray

bandpass, method, teff_star, stellar_spec

Stored input parameters for spectral calculations.

Type:: as above

nthreads, ntheta, nphi, pout

Stored configuration parameters.

Type:: as above

A0, A1, B1, A2, B2

Fourier coefficients derived from input parameters according to Cowan & Agol (2008) model equations.

Type:: ndarray

wav_star, fl_star

Stellar wavelength array and flux spectrum (blackbody or interpolated).

Type:: Astropy Quantity

wav_instrument, response_instrument, trans_fun

Instrument bandpass wavelengths, response curve, and normalized transmission function evaluated on stellar wavelengths.

Type:: Astropy Quantity or ndarray

phi_ang, theta_ang

1D arrays of longitude and latitude grid points.

Type:: ndarray

temp_map

3D array of temperature maps (shape: N_samples, N_phi, N_theta) computed when temperature_map_distribution() is called.

Type:: ndarray

Main methods

------------

TdayTnight[source]: Compute dayside and nightside brightness temperatures.

calculate_0D_Ab_eps(ar)[source]: Compute 0D Bond albedo and heat redistribution efficiency from brightness temperatures.

phase_offsets(method='root')[source]: Compute hotspot offset and phase offset from phase curve.

temperature_map_distribution(nsamples=2000)[source]: Compute 2D temperature map distribution across posterior samples.

median_temperature_map(plot=False, cmap='plasma')[source]: Compute and optionally plot median temperature map.

equatorial_temp_map(plot=False, nsamples=2000)[source]: Compute and optionally plot temperature along equator.

forward_phase_curve_model(plot=False)[source]: Compute forward-model phase curve from 2D flux distribution.

albedo_eps_from_temp_map(a_by_Rst, uniform_nightside_temp=True, nsamples=2000)[source]: Compute Bond albedo and heat redistribution from temperature map distribution.

albedo_eps_from_fpfs_map(a_by_Rst)[source]: Compute Bond albedo and heat redistribution from flux map distribution.

albedo_from_kempton23(a_by_Rst, teq)[source]: Compute Bond albedo using Kempton et al. (2023) method.

Notes

All temperature computations use the bandpass-weighted brightness temperature assuming blackbody emission.
The class uses multiprocessing for computationally intensive operations (e.g., temperature map fitting).
Intermediate results are cached as .npy files in pout to avoid re-computation.
The Cowan & Agol (2008) model assumes a parameterized intensity map I(φ, θ) with latitudinal and longitudinal dependencies.

References

Cowan, N. B., & Agol, E. 2008. Inverting Phase Functions to Map Exoplanets. The Astrophysical Journal Letters, 678(2), L129.
Kempton, E. M. R., et al. 2023. A reflective, metal-rich atmosphere for GJ 1214b from its JWST phase curve, Nature, 620(7972), 67
Keating, D., Cowan, N. B., & Dang, L. 2019. Uniformly hot nightside temperatures on short-period gas giants. Nature Astronomy, 3, 1092.

Examples

Analyze eclipse depths and phase curve to extract thermal properties:

>>> # Load posterior samples
>>> E = posterior_samples['eclipse_depth']  # Array of eclipse depths
>>> C1 = posterior_samples['phase_c1']
>>> D1 = posterior_samples['phase_d1']
>>> C2 = posterior_samples['phase_c2']
>>> D2 = posterior_samples['phase_d2']
>>> rprs = posterior_samples['rprs']

>>> # Define bandpass and stellar spectrum
>>> bandpass = load_bandpass('TESS')  # Must have 'WAVE' and 'RESPONSE' keys
>>> teff_star = 5800 * u.K  # Kelvin

>>> # Create the inversion object
>>> invert = InvertCowanAgolPC(E=E, C1=C1, D1=D1, C2=C2, D2=D2, rprs=rprs,
...                             bandpass=bandpass, teff_star=teff_star)

>>> # Compute temperature maps
>>> temp_maps = invert.temperature_map_distribution(nsamples=1000)

>>> # Compute albedos and heat redistribution
>>> a_by_rst = 10.0  # Semi-major axis in units of stellar radius
>>> A_B, eps_kelp, T_day, T_night, eps_keating = invert.albedo_eps_from_temp_map(a_by_rst)

TdayTnight()[source]

Compute 0D dayside and nightside brightness temperatures.

Calculates integrated brightness temperatures for the dayside and nightside using the eclipse depths and phase-curve parameters. Results are cached to disk to avoid re-computation. Uses the BrightnessTemperatureCalculator helper class for computation.

The dayside brightness temperature is derived directly from the eclipse depth, while the nightside is estimated from E - 2*C1.

Returns:

bt_day (ndarray) – Dayside brightness temperatures (in Kelvin) for each posterior sample. Shape: (N_samples,).
bt_night (ndarray) – Nightside brightness temperatures (in Kelvin) for each posterior sample. Shape: (N_samples,).

Notes

Results are cached to {pout}/Brightness_temp_day.npy and {pout}/Brightness_temp_night.npy.
Requires that self.bandpass, self.teff_star (or self.stellar_spec), are properly set during initialization.
This function should be called before calculate_0D_Ab_eps().

Examples

>>> invert = InvertCowanAgolPC(E=E, C1=C1, D1=D1, C2=C2, D2=D2,
...                             rprs=0.1, bandpass=bp, teff_star=5800 * u.K)
>>> bt_day, bt_night = invert.TdayTnight()
>>> print(f"Dayside brightness temperature: {np.median(bt_day)} K")

albedo_eps_from_fpfs_map(a_by_Rst)[source]

Compute Bond albedo and heat redistribution from flux map distribution.

Calculates Bond albedo and heat redistribution efficiency directly from the 2D distribution of F_p/F_* maps (without intermediate temperature conversion). Provides an alternative approach to albedo_eps_from_temp_map() that does not require brightness temperature estimation.

Parameters:

a_by_Rst (float or ndarray) – Semi-major axis in units of stellar radius (a/R_*). If scalar, same value used for all samples; if array, sampled for each.

Returns:

A_Bond_all (ndarray) – Bond albedo for each posterior sample. Shape: (N_samples,).
eps_kelp_all (ndarray) – Heat redistribution efficiency (Kelp method). Shape: (N_samples,).
eps_keating_all (ndarray) – Heat redistribution efficiency (Keating et al. 2019 method). Shape: (N_samples,).

Notes

Does not require calls to temperature-mapping routines.
Computes F_p/F_* maps directly from Fourier coefficients.
Prints posterior median and 68% credible intervals for all outputs.
Uses _albedo_eps_fpfs() for computation per sample.

Examples

>>> a_by_rst = 10.0
>>> A_B, eps_k, eps_keat = invert.albedo_eps_from_fpfs_map(a_by_rst)

albedo_eps_from_temp_map(a_by_Rst, uniform_nightside_temp=False, nsamples=2000)[source]

Compute Bond albedo and heat redistribution from temperature map distribution.

Calculates Bond albedo (A_Bond), heat redistribution efficiency (ε) using two methods (Kelp and Keating et al. 2019), and dayside/nightside brightness temperatures from the distribution of 2D temperature maps.

Parameters:

a_by_Rst (float or ndarray) – Semi-major axis in units of stellar radius (a/R_*). If scalar, same value used for all samples; if array, sampled for each.
uniform_nightside_temp (bool, optional) – If True, replace zero nightside temperatures (from zero grid points) with the median of measured nightside temperatures for self-consistency. Default is False.
nsamples (int, optional) – Number of posterior samples to process when computing temperature maps. Default is 2000.

Returns:

A_Bond_all (ndarray) – Bond albedo for each posterior sample. Shape: (N_samples,).
eps_kelp_all (ndarray) – Heat redistribution efficiency (Kelp method) for each sample. Shape: (N_samples,). Computed as nightside/dayside flux ratio.
Tday_all (ndarray) – Dayside brightness temperature for each sample. Shape: (N_samples,).
Tnight_all (ndarray) – Nightside brightness temperature for each sample. Shape: (N_samples,).
eps_keating_all (ndarray) – Heat redistribution efficiency (Keating et al. 2019 method). Shape: (N_samples,). Computed as nightside flux / total flux.

Notes

Calls temperature_map_distribution() to compute temperature maps if not already cached.
Prints posterior median and 68% credible intervals for all output parameters.
Uses _albedo_eps_tdaynight() for computation per sample.=

References

Keating, D., Cowan, N. B., & Dang, L. 2019. Uniformly hot nightside temperatures on short-period gas giants. Nature Astronomy, 3, 1092.

Examples

>>> a_by_rst = 10.0
>>> A_B, eps_k, T_d, T_n, eps_keat = invert.albedo_eps_from_temp_map(a_by_rst)

albedo_from_kempton23(a_by_Rst, rst, teq)[source]

Compute Bond albedo using Kempton et al. (2023) method across posteriors.

Calculates Bond albedo for each posterior sample using the energy-balance method of Kempton et al. (2023). This approach directly relates observed eclipse depths to Bond albedo accounting for bandpass-dependent radiation.

Parameters:

a_by_Rst (float or ndarray) – Semi-major axis in units of stellar radius (a/R_*). If scalar, same value used for all samples; if array, sampled for each.
rst (float or ndarray (must be astropy.units.Quantity)) – Stellar radius in units of solar radius (R_*). If scalar, same value used for all samples; if array, sampled for each.
teq (float or ndarray (must be astropy.units.Quantity)) – Equilibrium temperature of the planet in Kelvin. If scalar, same value used for all samples; if array, sampled for each.

Returns:

A_Bond_all – Bond albedo for each posterior sample. Shape: (N_samples,). Printed to console as median ± 68% credible interval.

Return type:

ndarray

Notes

Uses _albedo_kempton23() for computation per sample.
Requires rst (stellar radius) to be provided as an argument.
Prints posterior median and 68% credible intervals.

References

Kempton, E. M. R., et al. 2023. A reflective, metal-rich atmosphere for GJ 1214b from its JWST phase curve, Nature, 620(7972), 67

Examples

>>> a_by_rst = 10.0
>>> rst = 1.0 * u.R_sun  # Solar radius
>>> teq = 1500 * u.K  # Kelvin
>>> A_Bond = invert.albedo_from_kempton23(a_by_rst, rst, teq)
>>> print(f"Bond albedo: {np.median(A_Bond):.3f}")

calculate_0D_Ab_eps(ar)[source]

Compute Bond albedo and heat redistribution efficiency from 0D brightness temperatures.

Calculates the 0D (calculated from spatially averaged temperatures) Bond albedo and heat redistribution efficiency using the dayside and nightside brightness temperatures computed from TdayTnight(). This implementation uses the relationship between the planet’s temperature, stellar radiation, and thermal emission to infer physical properties.

Parameters:

ar (float, or ndarray) – Semi-major axis in units of stellar radius (a/R_*).

Returns:

A_B (ndarray) – Bond albedo for each posterior sample (unitless). Shape: (N_samples,). Values typically range from 0 to 1.
eps (ndarray) – Heat redistribution efficiency for each posterior sample (unitless). Shape: (N_samples,). Values typically range from 0 to 1.

Notes

Requires TdayTnight() to be called first (or this function will call it).
Assumes the 0D approximation where the planet is treated as having uniform dayside and nightside temperatures.

Examples

>>> bt_day, bt_night = invert.TdayTnight()
>>> A_Bond, eps = invert.calculate_0D_Ab_eps(ar=10.0)
>>> print(f"Bond albedo: {np.median(A_Bond):.3f}")
>>> print(f"Heat redistribution: {np.median(eps):.3f}")

equatorial_temp_map(plot=False, nsamples=2000)[source]

Compute and optionally plot temperature profile along the equator.

Calculates a distribution of 1D temperature profiles along the equator (θ = 0) across posterior samples. Useful for understanding longitudinal variations in brightness temperature at the substellar latitude.

Parameters:

plot (bool, optional) – If True, generate a plot showing the median equatorial temperature profile with individual posterior samples overplotted. Default is False.
nsamples (int, optional) – Number of posterior samples to process. Default is 2000.

Returns:

temp_map_equatorial (ndarray) – 2D array of brightness temperatures along the equator with shape (nsamples, N_phi). Also saved to {pout}/Equatorial_Temperature_map.npy.
fig (matplotlib.figure.Figure or None) – Figure object if plot=True, otherwise None.
axs (matplotlib.axes.Axes or None) – Axes object if plot=True, otherwise None.

Notes

Uses multiprocessing to parallelize temperature calculations.
Results are cached and reloaded on subsequent calls.
Posterior samples are randomly selected if the total exceeds nsamples.
When plot=True, the plot shows 100 random posterior samples as light curves and the median as a heavy line.

Examples

>>> Tmap_equator = invert.equatorial_temp_map(plot=True, nsamples=1000)
>>> print(Tmap_equator.shape)
(1000, 100)

forward_phase_curve_model(plot=False)[source]

Compute forward-model phase curve from 2D flux distribution.

Converts the median 2D flux map into a forward-modeled phase curve for comparison with the fitted phase-curve parameterization. Provides a self-consistency check between the temperature-map representation and the fitted phase-curve model.

Parameters:

plot (bool, optional) – If True, generate a plot overlaying the forward model on the fitted parameterized phase curve. Default is False.

Returns:

FpFs_forward_model (ndarray) – Forward-model phase curve values (F_p/F_*) sampled at 1000 orbital phase angles from -π to π. Also returned when plot=True.
fig (matplotlib.figure.Figure or None) – Figure object if plot=True, otherwise None.
axs (matplotlib.axes.Axes or None) – Axes object if plot=True, otherwise None.

Notes

Uses the median of the flux map distribution across all posterior samples.
Fitted phase curve is computed from posterior median values of C1, D1, C2, D2.
Plots both curves for visual comparison of model consistency.

Examples

>>> FpFs_fwd = invert.forward_phase_curve_model(plot=True)
>>> print(f"Max forward model flux: {np.max(FpFs_fwd):.4f}")

median_temperature_map(plot=False, cmap='plasma')[source]

Compute and optionally plot the median temperature map across posteriors.

Computes the median of the flux distribution across all posterior samples, then converts it to a single 2D temperature map. This provides a robust point estimate of the planet’s thermal structure.

Parameters:

plot (bool, optional) – If True, generate a Mollweide projection plot of the median temperature map. Default is False.
cmap (str, optional) – Colormap name (passed to matplotlib) for the plot. Default is 'plasma'.

Returns:

temp_map_median (ndarray) – 2D array of median brightness temperatures with shape (N_phi, N_theta). Also saved to {pout}/Median_Temperature_map.npy.
fig (matplotlib.figure.Figure or None) – Figure object if plot=True, otherwise None.
ax (matplotlib.axes.Axes or None) – Axes object (Mollweide projection) if plot=True, otherwise None.
cax (matplotlib.collections.QuadMesh or None) – Mesh container from pcolormesh if plot=True, otherwise None.

Notes

Results are cached to {pout}/Median_Temperature_map.npy and reloaded on subsequent calls.
If plot=True, the temperature color scale is set to median ± 5σ (using MAD-based standard deviation) to highlight structure.
Uses a Mollweide projection to show the full planetary surface.

Examples

>>> Tmap_median = invert.median_temperature_map(plot=True, cmap='viridis')
>>> print(f"Max dayside temperature: {np.nanmax(Tmap_median)} K")

phase_offsets(method='root')[source]

Compute hotspot offset and phase offset from phase-curve parameters.

Calculates two complementary measures of the offset of thermal features from expected locations:

Hotspot offset: The longitude (in degrees) where the brightness temperature maximum occurs, relative to the substellar point.
Phase offset: The orbital phase (in degrees relative to secondary eclipse) where the observed phase-curve flux peaks.

These quantities help constrain heat transport and atmospheric dynamics by revealing where thermal features are located relative to the tidally-locked dayside.

Parameters:

method (str, optional) – Root-finding method to use for identifying extrema. Options are 'root' (default, uses scipy.optimize.root) or 'fsolve' (uses scipy.optimize.fsolve). Default is 'root'.

Returns:

hotspot_off (ndarray) – Hotspot offset in degrees for each posterior sample. Shape: (N_samples,). Must be computed in the range [-180, 180] degrees.
phase_off (ndarray) – Phase offset in degrees for each posterior sample. Shape: (N_samples,). Positive values indicate the phase-curve maximum occurs before the secondary eclipse; negative values indicate it occurs after.

Notes

Computation uses _find_hotspot_off() and _find_phase_off() for each posterior sample.
This method uses progress bars (tqdm) to show computation progress.
Results are not cached; re-computing will recalculate offsets from scratch.

Examples

>>> hotspot_off, phase_off = invert.phase_offsets(method='root')
>>> print(f"Median hotspot offset: {np.median(hotspot_off):.1f} degrees")
>>> print(f"Median phase offset: {np.median(phase_off):.1f} degrees")

temperature_map_distribution(nsamples=2000)[source]

Compute 2D temperature map distribution across posterior samples.

Builds a distribution of 2D temperature maps (shape: N_samples, N_phi, N_theta) by converting the posterior distribution of flux maps into temperature maps via brightness temperature fitting. Results are cached and can be re-used by other methods.

This method randomly samples from the posterior distribution if more posterior samples exist than requested, improving computational efficiency while maintaining statistical representation.

Parameters:: nsamples (int, optional) – Number of posterior samples to process (randomly selected from all available samples). Default is 2000.
Returns:: temp_map – 3D array of brightness temperature maps with shape (nsamples, N_phi, N_theta). Also stored as self.temp_map.
Return type:: ndarray

Notes

Uses multiprocessing with self.nthreads processes for parallel computation of temperature maps.
Results are cached to {pout}/Temperature_map.npy and reloaded on subsequent calls.
Uses progress bars (tqdm) to show computation progress.
This is a computationally intensive operation and typically requires several minutes to hours for large grids.

Examples

>>> temp_maps = invert.temperature_map_distribution(nsamples=1000)
>>> print(temp_maps.shape)
(1000, 100, 49)

>>> # Get median map
>>> temp_median = np.median(temp_maps, axis=0)

The SelectLinDetrend class can be used to perform linear detrending and model selection using log-evidence or scatter in the residuals. It has the following methods:

class krithika.SelectLinDetrend(time, flux, flux_err, priors, linear_regressors, gp_regressors=None, roll_degree=None, roll=None, noise_method='astropy', juliet_fit_kwargs={}, nthreads=2, pout='/home/docs/checkouts/readthedocs.org/user_builds/krithika/checkouts/latest/docs')[source]

Select optimal linear detrending vectors using model selection with Bayesian evidence.

This class provides tools to systematically identify the best combination of linear detrending regressors (e.g., time polynomials, roll angle modulation, background variations) for multi-instrument light-curve fitting. It employs model selection based on log-evidence improvements, starting from a base model with no linear regressors and iteratively adding vectors that maximize evidence gain. A regressor is accepted only if it improves log-evidence by at least 2 units; the search terminates when no remaining regressor meets this threshold.

The class interfaces with the juliet package for Bayesian model fitting and supports parallel processing of model comparisons through concurrent fitting. It automatically handles multi-instrument datasets where different instruments may have different numbers and types of regressors available.

Linear regressors can be passed directly or be generated from roll angle values via Fourier series expansion. All data (time, flux, errors, regressors) are expected as per-instrument dictionaries, enabling straightforward multi-instrument analysis.

Parameters:

time (dict) – Time stamps for each instrument. Keys are instrument names; values are 1-D numpy arrays of time points (in BJD or similar). All arrays should have the same length as corresponding flux and flux_err arrays.
flux (dict) – Flux (or normalized flux) measurements for each instrument. Keys are instrument names; values are 1-D numpy arrays of flux measurements.
flux_err (dict) – Flux uncertainties for each instrument. Keys are instrument names; values are 1-D numpy arrays of flux uncertainties (same length as flux).
priors (callable) –
A function mapping an instrument name (keyword argument instrument) to a tuple of three lists:
- param_names : list of str Names of juliet-compatible model parameters (e.g., ‘t0_p1’, ‘p_p1’, ‘P_p1’).
- distributions : list of str juliet-compatible distribution names (e.g., ‘normal’, ‘uniform’, ‘loguniform’).
- hyperparams : list Hyperparameter lists for each distribution. For example, a ‘normal’ distribution takes [mean, std], and a ‘uniform’ takes [lower, upper].
The function is called once per instrument during initialization and results are cached internally.
linear_regressors (dict) –
Linear regressors to consider for detrending. Keys are instrument names; values are dicts where each key is a regressor name (str) and value is a 1-D numpy array (same length as time for that instrument). Examples include ‘time’, ‘time^2’, ‘background’, ‘systematics’, etc.

If roll_degree and roll are both provided, roll-angle Fourier series regressors (‘ROLL_SIN_1’, ‘ROLL_COS_1’, etc.) are automatically added to linear_regressors during initialization.
gp_regressors (dict or None, optional) – Gaussian Process (GP) regressors to use for systematics fitting. Keys are instrument names; values are 1-D numpy arrays of regressor values (e.g., roll angle, time). If provided, these are sorted and passed to juliet for GP handling. Default is None.
roll_degree (int or None, optional) – Maximum degree of roll-angle Fourier series to include as linear regressors. If provided (and roll is also provided), Fourier series regressors up to degree roll_degree are automatically generated and added to linear_regressors. Must be paired with roll parameter. Default is None.
roll (dict or None, optional) – Roll angle values (in degrees) for each instrument. Keys are instrument names; values are 1-D numpy arrays of roll angles. Required only if roll_degree is specified. Default is None.
noise_method (str, optional) –
Method to compute noise/scatter in residuals. Options are:
- ’astropy’ : Uses median absolute deviation (astropy.stats.mad_std).
- ’std’ : Uses standard deviation (numpy.std).
- ’pipe’ : Uses the PIPE photometric uncertainty (pipe_mad function).
This method is applied to residuals after each model fit to estimate per-instrument noise levels. Default is 'astropy'.
juliet_fit_kwargs (dict, optional) – Additional keyword arguments to pass to juliet.fit() during model fits. For example, {'sampler': 'dynamic_dynesty', 'nlive': 500}. Default is an empty dict {}.
nthreads (int, optional) – Total number of CPU threads available for fitting. Shared among parallel workers if n_parallel > 1 is used in select_optimal_parameters(). Default is the total number of available CPU cores.
pout (str, optional) – Output directory for saving juliet fit results. Each model fit creates a subdirectory named after the included regressors (e.g., ‘linsel_base’, ‘linsel_time’, ‘linsel_time_background’). Default is the current working directory.

time, flux, flux_err

Stored input light-curve data (per-instrument).

Type:: dict

instruments

Extracted list of instrument names from the keys of time dict.

Type:: list of str

priors

Stored prior function.

Type:: callable

linear_regressors

Stored linear regressors (or regressors with added roll terms if roll_degree was provided).

Type:: dict

GP_regressors

Stored GP regressors.

Type:: dict or None

roll_degree, roll

Stored roll parameters.

Type:: int or None, dict or None

noise_method

Stored noise method name.

Type:: str

noise_func

Cached noise function (e.g., mad_std, np.std, rms, or pipe_mad).

Type:: callable

juliet_fit_kwargs

Stored juliet fitting kwargs.

Type:: dict

nthreads

Stored thread count.

Type:: int

pout

Stored output directory.

Type:: str

_sort_idx, _tim_s, _fl_s, _fle_s

Per-instrument sort indices and sorted data arrays (cached by _precompute()).

Type:: dict

_gp_pars_s, _gp_pars_s_juliet

Cached GP regressor arrays (sorted, or None if not applicable).

Type:: dict or None

_cached_priors

Per-instrument prior tuples cached from the priors callable.

Type:: dict

Main Methods

----------

select_optimal_parameters(n_parallel=1)[source]: Execute forward selection to identify optimal linear detrending vectors based on log-evidence maximization.

Notes

All fitting is performed using the juliet Bayesian fitting package.
Log-evidence improvements are computed via juliet.fit() and returned as posteriors['lnZ'].
The model selection algorithm is greedy: it selects the regressor with the maximum evidence gain at each step, then moves to the next round.
A regressor is accepted only if it improves log-evidence by >= delta_lnZ_threshold (default 2.0). This threshold reflects a Bayes factor of ~exp(2) ≈ 7.4, corresponding to “strong” evidence in Bayesian hypothesis testing.
When n_parallel > 1 in select_optimal_parameters(), concurrent fitting is performed using concurrent.futures.ProcessPoolExecutor. To keep total thread usage under control, the per-worker nthreads budget is automatically adjusted as nthreads // n_parallel.
Roll-angle Fourier series regressors (if generated) follow the naming convention ‘ROLL_SIN_k’ and ‘ROLL_COS_k’ for degree k = 1, 2, …, roll_degree.
Regressors that do not exist for a particular instrument (e.g., a regressor available only for one camera in a multi-camera setup) are automatically skipped for that instrument during model construction.
All per-instrument data (time, flux, errors, regressors) are length-checked during juliet loading but not explicitly in __init__. Ensure consistency before instantiation.

Examples

Basic forward selection with two instruments:

>>> import numpy as np
>>> time = {
...     'camera_1': np.linspace(2459000, 2459001, 200),
...     'camera_2': np.linspace(2459000, 2459001, 150),
... }
>>> flux = {
...     'camera_1': 1.0 + 0.01 * np.random.randn(200),
...     'camera_2': 1.0 + 0.005 * np.random.randn(150),
... }
>>> flux_err = {
...     'camera_1': np.full(200, 0.005),
...     'camera_2': np.full(150, 0.003),
... }
>>> linear_regressors = {
...     'camera_1': {
...         'time': time['camera_1'],
...         'time^2': time['camera_1']**2,
...     },
...     'camera_2': {
...         'time': time['camera_2'],
...         'time^2': time['camera_2']**2,
...     },
... }
>>> def priors_func(instrument):
...     return (['t0_p1', 'p_p1', 'P_p1'],
...             ['normal', 'normal', 'uniform'],
...             [[2459000.5, 0.1], [3.5, 0.01], [0.01, 0.1]])
>>>
>>> selector = SelectLinDetrend(
...     time=time,
...     flux=flux,
...     flux_err=flux_err,
...     priors=priors_func,
...     linear_regressors=linear_regressors,
...     pout='/path/to/results'
... )
>>> selected_regs, lnZ_hist, scatter_hist = selector.select_optimal_parameters(n_parallel=2)

With roll-angle detrending:

>>> roll_angles = {
...     'camera_1': np.random.uniform(0, 360, 200),
...     'camera_2': np.random.uniform(0, 360, 150),
... }
>>> selector = SelectLinDetrend(
...     time=time,
...     flux=flux,
...     flux_err=flux_err,
...     priors=priors_func,
...     linear_regressors=linear_regressors,
...     roll_degree=3,
...     roll=roll_angles,
...     pout='/path/to/results'
... )
>>> # Now linear_regressors includes 'ROLL_SIN_1', 'ROLL_COS_1', ..., 'ROLL_SIN_3', 'ROLL_COS_3'
>>> selected_regs, lnZ_hist, scatter_hist = selector.select_optimal_parameters()

It is possible to use scatter in the residuals (instead of log-evidence) for model selection by setting selection_method=’scatter’ in select_optimal_parameters. In this case, the regressor that minimizes the scatter (as computed by noise_method) is selected at each step.

>>> # Using scatter-based selection instead of log-evidence
>>> selected_regs, lnZ_hist, scatter_hist = selector.select_optimal_parameters(selection_method='scatter')

select_optimal_parameters(n_parallel=1, delta_lnZ_threshold=2.0, selection_method='lnZ', scatter_threshold=5.0)[source]

Select optimal linear detrending vectors using forward selection based on Bayesian evidence or light-curve scatter.

The algorithm starts with a base model that has no linear regressors and iteratively adds the vector that yields the greatest improvement according to the chosen selection metric. The search stops when no remaining vector clears the acceptance threshold.

Unique regressor names are collected across all instruments. When a regressor is selected it is included for every instrument that carries it; instruments that do not have that regressor are left unchanged.

Parameters:

n_parallel (int, optional) –
Number of model fits to run concurrently within each selection round. Default is 1 (sequential).

When n_parallel > 1 the fits in the inner loop are dispatched to a ProcessPoolExecutor for true CPU parallelism (no GIL contention). The juliet nthreads budget is automatically divided by n_parallel so that the total thread count across all workers stays within self.nthreads. The module-level _lindetrend_worker function is used as the callable so that no part of this instance (including the potentially un-picklable priors callable) needs to be serialised.
delta_lnZ_threshold (float, optional) – Minimum improvement in log-evidence required to accept a new regressor when selection_method=’lnZ’. A candidate is added only if its lnZ exceeds the current base model’s lnZ by at least this value. Default is 2.0, which corresponds to “strong evidence” on the Jeffreys scale. Ignored when selection_method=’scatter’.
selection_method ({'lnZ', 'scatter'}, optional) –
Metric used to rank and accept candidate regressors.
- 'lnZ' (default) — select the regressor that maximises the Bayesian log-evidence improvement. Works for any number of instruments.
- 'scatter' — select the regressor that minimises the MAD scatter of the residuals for the single instrument. If more than one instrument is provided this option is not supported and the method automatically falls back to 'lnZ' with a warning.
scatter_threshold (float, optional) – Minimum reduction in scatter (in ppm) required to accept a new regressor when selection_method=’scatter’. A candidate is accepted only if its residual scatter is at least this many ppm lower than the current base model’s scatter. Default is 5.0 ppm. Ignored when selection_method=’lnZ’.

Returns:

selected_regressors (dict) – Dict with instrument names as keys. Each value is either a dict of {regressor_name: array} for the selected regressors, or None if no regressor was selected for that instrument.
lnZ_history (list of float) – Log-evidence at each accepted step (index 0 is the base model).
scatter_history (list) – Scatter (MAD) values at each accepted step (index 0 is the base model).

krithika has a simple light weight class to view the 2D, 3D, or 4D data cubes, with functionalities similar to DS9. The NDImageViewer class has the following methods:

class krithika.NDImageViewer(data, cmap='magma')[source]

Interactive viewer for 2D, 3D, and 4D image data with real-time analysis tools.

This class provides an interactive matplotlib-based interface for viewing and analyzing multi-dimensional image data. It supports 2D images, 3D image cubes (e.g., time series), and 4D data (e.g., time series of grouped observations). Features include dynamic scaling modes, cut profile extraction, and draggable cut tools for pixel-level analysis.

Parameters:

data (array-like) – Input image data with shape (ny, nx) for 2D, (nt, ny, nx) for 3D, or (nt, ng, ny, nx) for 4D, where nt=time, ng=group, ny/nx=spatial dimensions.
cmap (str, optional) – Matplotlib colormap name for image display. Default is 'magma'.

data

The input image data.

Type:: ndarray

ndim

Number of dimensions in the data (2, 3, or 4).

Type:: int

i_image

Current image/frame index for 3D/4D data (default 0).

Type:: int

i_group

Current group index for 4D data (default 0).

Type:: int

scale_mode

Current scaling mode (‘linear’, ‘log’, ‘asinh’, or ‘zscale’).

Type:: str

cuts

List of extracted cut profiles, each containing ‘p1’, ‘p2’, ‘line’, ‘color’, and ‘visible’ keys.

Type:: list of dict

active_cut

Cut currently being drawn (awaiting two clicks).

Type:: dict or None

fig

The interactive figure object.

Type:: matplotlib.figure.Figure

ax_img

Main image display axes.

Type:: matplotlib.axes.Axes

ax_prof

Profile plot axes for cut data.

Type:: matplotlib.axes.Axes

im

The displayed image object.

Type:: matplotlib.image.AxesImage

cbar

Colorbar for the image display.

Type:: matplotlib.colorbar.Colorbar

Notes

Data must be 2D, 3D, or 4D; other dimensions will raise a ValueError.
Finite data values are automatically detected for scaling limits.
The viewer is interactive: sliders update in real-time, and cuts can be drawn and modified by clicking/dragging on the image.

Examples

View a 2D image:

>>> import numpy as np
>>> data_2d = np.random.rand(256, 256)
>>> viewer = NDImageViewer(data=data_2d)
>>> viewer.show()

View a 3D time series:

>>> data_3d = np.random.rand(100, 256, 256)  # 100 frames
>>> viewer = NDImageViewer(data=data_3d, cmap='viridis')
>>> viewer.show()

View 4D data and extract cuts:

>>> data_4d = np.random.rand(10, 8, 256, 256)  # 10 times, 8 groups
>>> viewer = NDImageViewer(data=data_4d)
>>> viewer.show()
>>> # Click "Add cut" button, then click two points on the image to define a cut

show()[source]

Display the interactive viewer window.

Calls plt.show() to render the figure and start the interactive loop.

Apart from these classes, there are several functions and classes that can be used for different purposes. They are listed below:

krithika.utils.t14(per, ar, rprs, b, ecc=0, omega=90, transit=True)[source]

To compute transit/eclipse duration from Period, a/R*, Rp/R*, b, eccentricity, and omega

Parameters:

perfloat, or numpy.ndarray: Orbital period (in days) of the planet
aRfloat, or numpy.ndarray: Scaled semi-major axis, a/R*
rprsfloat, or numpy.ndarray: Planet-to-star radius ratio, Rp/R*
bbfloat, or numpy.ndarray: Impact parameter
ecc: float, or numpy.ndarray: Eccentricity of the orbit
omega: float, or numpy.ndarray: Argument of periastron (deg)
transit: bool: Whether to compute transit or occultation duration

krithika.utils.b_to_inc(b, ar, ecc=0.0, omega=90.0, transit=True)[source]: Computing inclination from impact parameter

krithika.utils.inc_to_b(inc, ar, ecc=0.0, omega=90.0, transit=True)[source]: To compute impact parameter from inclination

krithika.utils.rho_to_ar(rho, per)[source]: To compute a/R* from rho

krithika.utils.pipe_mad(x, axis=0)[source]

krithika.utils.rms(x)[source]

krithika.utils.fake_allan_deviation(times, residuals, binmax=10, method='pipe', timeunit=None, plot=True)[source]

Estimate and optionally plot a noise-vs-bin-size curve (a proxy for Allan deviation).

The function computes the scatter of binned residuals for a sequence of integer bin sizes and returns both the computed noise (in ppm) and a white-noise expectation for comparison. When plot=True a log-log figure is produced with a secondary x-axis that maps bin size to a time-unit (days/hours/minutes) according to the total span of times.

Parameters:

times (array-like) – Time stamps of the observations (assumed in days).
residuals (array-like) – Residual fluxes (same length as times).
binmax (int, optional) – Parameter controlling the maximum number of bins. The function will evaluate binsizes from 1..int(len(times)/binmax). Default 10.
method ({'pipe', 'std', 'rms', 'astropy'}, optional) – Method used to estimate the scatter of the binned residuals: - ‘pipe’ : use pipe_mad (median absolute differences estimator) - ‘std’ : use np.nanstd - ‘rms’ : use the root-mean-square (rms helper) - ‘astropy’ : use astropy.stats.mad_std (robust MAD-based std)
timeunit ({'d', 'hr', 'min'}, optional) – If provided, forces the time unit used in the secondary x-axis to days (‘d’), hours (‘hr’), or minutes (‘min’). By default the function selects the most appropriate unit based on the total span of times.
plot (bool, optional) – If True (default), create and display a matplotlib figure with the computed noise curve and the white-noise expectation.

Returns:

fig (matplotlib.figure.Figure or None) – The created figure if plot=True, otherwise None.
axs (matplotlib.axes.Axes or None) – The axes object for the figure if plot=True, otherwise None.
binsize (ndarray) – Integer array of evaluated bin sizes.
noise_ppm (ndarray) – Computed scatter of binned residuals converted to parts-per-million (ppm), corresponding to binsize.
white_noise_expec (ndarray) – White-noise expectation (ppm) computed as noise_func(residuals) / sqrt(binsize) (with a small-sample correction applied). Values for binsize==0 will be inf or undefined.

krithika.utils.corner_plot(samples, labels, **kwargs)[source]

Generate a corner plot from MCMC samples.

Parameters:

samples (ndarray) – 2D array of shape (n_samples, n_parameters) containing the MCMC samples.
labels (list of str) – List of parameter names for labeling the axes.
figsize (tuple, optional) – Size of the figure to create. Default is (9, 9).

Returns:

fig – The generated corner plot figure.

Return type:

matplotlib.figure.Figure

krithika.utils.make_psd(times, flux, nos_freq_points=100000, plot=True, plot_max_freq=True, timeunit=None)[source]

Compute a Lomb-Scargle power spectral density (PSD) for a light curve.

The routine converts times (assumed in days) to seconds and the input flux to parts-per-million (ppm), computes a dense frequency grid, evaluates a Lomb-Scargle PSD on that grid and (optionally) produces a log-log plot with a secondary axis showing period in an appropriate time unit.

Parameters:

times (array-like) – Time stamps of the observations (in days).
flux (array-like) – Relative flux measurements (unitless, e.g. normalised to 1.0). The function converts these to ppm internally.
nos_freq_points (int, optional) – Number of frequency points to evaluate the PSD on. Default is 100,000.
plot (bool, optional) – If True (default) create and return a matplotlib figure and axes containing the PSD plot. If False, no plot is created and fig and axs in the return tuple will be None.
plot_max_freq (bool, optional) – If True (default) mark the frequency with maximum power on the plot with a vertical dashed line.
timeunit ({'d','hr','min'} or None, optional) – Force the secondary x-axis unit for the period conversion. If None (default) the function selects days/hours/minutes automatically based on the frequency grid span.

Returns:

freq_grid (astropy.units.Quantity (Hz)) – Frequency grid used to evaluate the PSD.
psd1 (astropy.units.Quantity) – Lomb-Scargle power spectral density evaluated on freq_grid (units: ppm^2 Hz^-1).
fig (matplotlib.figure.Figure or None) – Figure with the PSD plot when plot=True, otherwise None.
axs (matplotlib.axes.Axes or None) – Axes for the PSD plot when plot=True, otherwise None.
per_max_pow (astropy.units.Quantity) – Period corresponding to the maximum power in the PSD (converted to days).

Notes

The function uses astropy.timeseries.LombScargle for the PSD and multiplies the input flux by 1e6 to report power in ppm^2 Hz^-1.

krithika.utils.generate_times_with_gaps(times, efficiency)[source]

This function takes the times array with a regular cadence and generate another times array with data gaps in it with a given efficiency to mimick CHEOPS observations. Efficiency here is the fraction of the actual time spent on the star.

Parameters:

times (numpy.ndarray) – Regularly spaced time array
efficiency (float) – Efficiency of the observations in per cent.

Returns:

tim (numpy.ndarray) – Time array with gaps in it.
roll (numpy.ndarray) – The corresponding roll angle array.

class krithika.BrightnessTemperatureCalculator(fp, rprs, bandpass, nthreads=2, method='Nelder-Mead', teff_star=None, stellar_spec=None, pout='/home/docs/checkouts/readthedocs.org/user_builds/krithika/checkouts/latest/docs')[source]

Compute brightness temperature(s) from eclipse depths and a bandpass.

Parameters:

fp (float or array-like) – Eclipse depths (planet/star flux ratio). If array-like, temperature distribution will be computed (multiprocessing is used in this case).
rprs (float or array-like) – Planet-to-star radius ratio (Rp/R*). If scalar and fp is array-like, the same rprs is applied to all elements.
bandpass (dict) – Bandpass dict with keys ‘WAVE’ (astropy Quantity) and ‘RESPONSE’ (array).
nthreads (int, optional) – Number of worker processes for multiprocessing. Default is cpu count.
method (str, optional) – Minimizer method passed to scipy.optimize.minimize. Default ‘Nelder-Mead’.
teff_star (float or None) – Stellar effective temperature in K (used if stellar_spec not provided).
stellar_spec (dict or None) – Stellar spectrum dict with keys ‘WAVE’ and ‘FLUX’ (astropy Quantities).
pout (str, optional) – Output directory used to cache results. Default current working dir.

compute()[source]: Run the calculation and return temperature(s) in Kelvin (float or 1D ndarray).

compute()[source]

Compute and return brightness temperature(s) in Kelvin.

Returns:: temp_pl – Brightness temperature (K) for scalar input or 1D ndarray for array input.
Return type:: float or ndarray