from types import SimpleNamespace
from typing import Union
from pathlib import Path
import astropy.units as u
import numpy as np
from astropy.coordinates import SkyCoord
from astropy.table import Table
from astropy.time import Time
from astropy.units import Quantity
from sunpy.coordinates import HeliographicStonyhurst
from sunpy.time import TimeRange
from xrayvision.visibility import Visibilities, VisMeta
from stixpy.calibration.energy import get_elut
from stixpy.calibration.grid import get_grid_transmission
from stixpy.calibration.livetime import get_livetime_fraction
from stixpy.config.instrument import STIX_INSTRUMENT, _get_uv_points_data
from stixpy.coordinates.frames import STIXImaging
from stixpy.coordinates.transforms import get_hpc_info
from stixpy.product.sources import CompressedPixelData, RawPixelData, SummedCompressedPixelData
__all__ = [
"get_subcollimator_info",
"create_meta_pixels",
"create_visibility",
"get_uv_points_data",
"calibrate_visibility",
]
_PIXEL_SLICES = {
"top": slice(0, 4),
"bot": slice(4, 8),
"top+bot": slice(0, 8),
"all": slice(None),
"small": slice(8, None),
}
[docs]
@u.quantity_input
def get_uv_points_data(d_det: u.Quantity[u.mm] = 47.78 * u.mm, d_sep: u.Quantity[u.mm] = 545.30 * u.mm):
r"""
Return the STIX (u,v) points coordinates defined in [1], ordered with respect to the detector index.
It will return the precalculated cached data if the default values for `d_det` and `d_sep` are used.
Parameters
----------
d_det: `u.Quantity[u.mm]` optional
Distance between the rear grid and the detector plane (in mm). Default, 47.78 * u.mm
d_sep: `u.Quantity[u.mm]` optional
Distance between the front and the rear grid (in mm). Default, 545.30 * u.mm
Returns
-------
A dictionary containing sub-collimator indices, sub-collimator labels and coordinates of the STIX (u,v) points (defined in arcsec :sup:`-1`)
References
----------
[1] Massa et al., 2023, The STIX Imaging Concept, Solar Physics, 298,
https://doi.org/10.1007/s11207-023-02205-7
"""
if d_det == 47.78 * u.mm and d_sep == 545.30 * u.mm:
# If the default values are used, return the cached data
return STIX_INSTRUMENT.vis_info
return _get_uv_points_data(d_det=d_det, d_sep=d_sep)
[docs]
def get_subcollimator_info():
r"""
Return resolutions, orientations, and labels for sub-collimators.
Returns
-------
"""
grids = {
9: np.array([3, 20, 22]) - 1,
8: np.array([16, 14, 32]) - 1,
7: np.array([21, 26, 4]) - 1,
6: np.array([24, 8, 28]) - 1,
5: np.array([15, 27, 31]) - 1,
4: np.array([6, 30, 2]) - 1,
3: np.array([25, 5, 23]) - 1,
2: np.array([7, 29, 1]) - 1,
1: np.array([12, 19, 17]) - 1,
0: np.array([11, 13, 18]) - 1,
}
labels = {i: [str(i + 1) + "a", str(i) + "b", str(i) + "c"] for i in range(10)}
res32 = np.zeros(32)
res32[grids[9]] = 178.6
res32[grids[8]] = 124.9
res32[grids[7]] = 87.3
res32[grids[6]] = 61.0
res32[grids[5]] = 42.7
res32[grids[4]] = 29.8
res32[grids[3]] = 20.9
res32[grids[2]] = 14.6
res32[grids[1]] = 10.2
res32[grids[0]] = 7.1
res10 = [7.1, 10.2, 14.6, 20.9, 29.8, 42.7, 61.0, 87.3, 124.9, 178.6]
o32 = np.zeros(32)
o32[grids[9]] = [150, 90, 30]
o32[grids[8]] = [170, 110, 50]
o32[grids[7]] = [10, 130, 70]
o32[grids[6]] = [30, 150, 90]
o32[grids[5]] = [50, 170, 110]
o32[grids[4]] = [70, 10, 130]
o32[grids[3]] = [90, 30, 150]
o32[grids[2]] = [110, 50, 170]
o32[grids[1]] = [130, 70, 10]
o32[grids[0]] = [150, 90, 30]
# g03_10 = [g03, g04, g05, g06, g07, g08, g09, g10]
# g01_10 = [g01, g02, g03, g04, g05, g06, g07, g08, g09, g10]
# g_plot = [g10, g05, g09, g04, g08, g03, g07, g02, g06, g01]
# l_plot = [l10, l05, l09, l04, l08, l03, l07, l02, l06, l01]
res = SimpleNamespace(res32=res32, o32=o32, res10=res10, label=labels)
return res
def get_elut_correction(e_ind, pixel_data):
r"""
Get ELUT correction factors
Only correct the low and high energy edges as internal bins are contiguous.
Parameters
----------
e_ind : `np.ndarray`
Energy channel indices
pixel_data : `~stixpy.products.sources.CompressedPixelData`
Pixel data
Returns
-------
"""
energy_mask = pixel_data.energy_masks.energy_mask.astype(bool)
elut = get_elut(pixel_data.time_range.center)
ebin_edges_low = np.zeros((32, 12, 32), dtype=float)
ebin_edges_low[..., 1:] = elut.e_actual
ebin_edges_low = ebin_edges_low[..., energy_mask]
ebin_edges_high = np.zeros((32, 12, 32), dtype=float)
ebin_edges_high[..., 0:-1] = elut.e_actual
ebin_edges_high[..., -1] = np.nan
ebin_edges_high = ebin_edges_high[..., energy_mask]
ebin_widths = ebin_edges_high - ebin_edges_low
ebin_sci_edges_low = elut.e[..., 0:-1].value
ebin_sci_edges_low = ebin_sci_edges_low[..., energy_mask]
ebin_sci_edges_high = elut.e[..., 1:].value
ebin_sci_edges_high = ebin_sci_edges_high[..., energy_mask]
e_cor_low = (ebin_edges_high[..., e_ind[0]] - ebin_sci_edges_low[..., e_ind[0]]) / ebin_widths[..., e_ind[0]]
e_cor_high = (ebin_sci_edges_high[..., e_ind[-1]] - ebin_edges_low[..., e_ind[-1]]) / ebin_widths[..., e_ind[-1]]
return e_cor_high, e_cor_low
[docs]
def create_visibility(meta_pixels):
r"""
Create visibilities from meta-pixels
Parameters
----------
meta_pixels
Returns
-------
"""
abcd_rate_kev_cm = meta_pixels["abcd_rate_kev_cm"]
abcd_rate_error_kev_cm = meta_pixels["abcd_rate_error_kev_cm"]
real = abcd_rate_kev_cm[:, 2] - abcd_rate_kev_cm[:, 0]
imag = abcd_rate_kev_cm[:, 3] - abcd_rate_kev_cm[:, 1]
real_err = np.sqrt(abcd_rate_error_kev_cm[:, 2] ** 2 + abcd_rate_error_kev_cm[:, 0] ** 2).value * real.unit
imag_err = np.sqrt(abcd_rate_error_kev_cm[:, 3] ** 2 + abcd_rate_error_kev_cm[:, 1] ** 2).value * real.unit
vis_info = get_uv_points_data()
# Compute raw phases
phase = np.arctan2(imag, real)
# Compute amplitudes
observed_amplitude = np.sqrt(real**2 + imag**2)
# Compute error on amplitudes
amplitude_error = np.sqrt((real / observed_amplitude * real_err) ** 2 + (imag / observed_amplitude * imag_err) ** 2)
# Apply pixel correction
pix_set = meta_pixels.get("pixels", None)
if pix_set in {"top", "bot", "top+bot"}:
phase += 46.1 * u.deg # Center of large pixel in terms morie pattern phase
elif pix_set == "all":
phase += 45 * u.deg
elif pix_set == "small":
phase += 22.5 * u.deg
else:
raise ValueError(
f"Set of pixels used to make meta pixels not recognised, {pix_set} should be one of {list(_PIXEL_SLICES.keys())}"
)
obsvis = (np.cos(phase) + np.sin(phase) * 1j) * observed_amplitude
imaging_detector_indices = vis_info["isc"] - 1
vis_meta = VisMeta(
instrumet="STIX",
spectral_range=meta_pixels["energy_range"],
time_range=Time([meta_pixels["time_range"].start, meta_pixels["time_range"].end]),
vis_labels=vis_info["label"].tolist(),
isc=vis_info["isc"].value,
calibrated=False,
)
vis = Visibilities(
obsvis[imaging_detector_indices],
u=vis_info["u"],
v=vis_info["v"],
amplitude=observed_amplitude[imaging_detector_indices],
amplitude_uncertainty=amplitude_error[imaging_detector_indices],
meta=vis_meta,
)
return vis
[docs]
def calibrate_visibility(
vis: Visibilities, flare_location: SkyCoord = SkyCoord(0 * u.arcsec, 0 * u.arcsec, frame=STIXImaging)
):
"""
Calibrate visibility phase and amplitudes.
Parameters
----------
vis :
Uncalibrated Visibilities
flare_location
The location of the flare the visibility are shifted to have the phase center at this location if not given will
default to sun center `Helioprojective(0,0)`
Returns
-------
"""
modulation_efficiency = np.pi**3 / (8 * np.sqrt(2))
# Grid correction factors
grid_cal_file = Path(__file__).parent.parent / "config" / "data" / "grid" / "GridCorrection.csv"
grid_cal = Table.read(grid_cal_file, header_start=2, data_start=3)
grid_corr = grid_cal["Phase correction factor"][vis.meta["isc"] - 1] * u.deg
# Phase correction
phase_cal_file = Path(__file__).parent.parent / "config" / "data" / "grid" / "PhaseCorrFactors.csv"
phase_cal = Table.read(phase_cal_file, header_start=3, data_start=4)
phase_corr = phase_cal["Phase correction factor"][vis.meta["isc"] - 1] * u.deg
tr = TimeRange(vis.meta.time_range)
if not isinstance(flare_location, SkyCoord):
raise ValueError("'flare_location' is not a SkyCoord object")
if not isinstance(flare_location.frame, STIXImaging) or flare_location.obstime != tr.center:
roll, solo_heeq, stix_pointing = get_hpc_info(vis.meta.time_range[0], vis.meta.time_range[1])
solo_coord = HeliographicStonyhurst(solo_heeq, representation_type="cartesian", obstime=tr.center)
flare_location = flare_location.transform_to(STIXImaging(obstime=tr.center, observer=solo_coord))
phase_mapcenter_corr = -2 * np.pi * (flare_location.Tx * vis.u + flare_location.Ty * vis.v) * u.rad
ovis = vis.visibilities[:]
ovis_real = np.real(ovis)
ovis_imag = np.imag(ovis)
ovis_amp = np.sqrt(ovis_real**2 + ovis_imag**2)
ovis_phase = np.arctan2(ovis_imag, ovis_real).to("deg")
calibrated_amplitude = ovis_amp * modulation_efficiency
calibrated_phase = ovis_phase + grid_corr + phase_corr + phase_mapcenter_corr
# Compute error on amplitudes
systematic_error = 5 * u.percent # from G. Hurford 5% systematics
calibrated_amplitude_error = np.sqrt(
(vis.amplitude_uncertainty * modulation_efficiency) ** 2
+ (systematic_error * calibrated_amplitude).to(ovis_amp.unit) ** 2
)
calibrated_visibility = (np.cos(calibrated_phase) + np.sin(calibrated_phase) * 1j) * calibrated_amplitude
vis.meta["calibrated"] = True
vis.meta["offset"] = flare_location
cal_vis = Visibilities(
calibrated_visibility,
u=vis.u,
v=vis.v,
amplitude=calibrated_amplitude,
amplitude_uncertainty=calibrated_amplitude_error,
meta=vis.meta,
)
return cal_vis