from pathlib import Path
from dataclasses import dataclass
from typing import Any
import numpy as np
from astropy import units as u
from astropy.io import ascii as ioascii
from astropy.table import Table
from synphot import SpectralElement
from synphot.models import Empirical1D
from ..effects import ter_curves_utils as ter_utils
from .surface_utils import (make_emission_from_emissivity,
make_emission_from_array, extract_type_from_unit)
from ..utils import (get_meta_quantity, quantify, from_currsys,
convert_table_comments_to_dict, find_file, get_logger)
logger = get_logger(__name__)
[docs]@dataclass
class PoorMansSurface:
"""Solely used by SurfaceList."""
# FIXME: Use correct types instead of Any
emission: Any
throughput: Any
meta: Any
[docs]class SpectralSurface:
"""
Initialised by a file containing one or more of the following columns.
transmission, emissivity, reflection. The column wavelength must be given.
Alternatively kwargs for the above mentioned quantities can be passed as
arrays. If they are not Quantities, then a unit should also be passed with
the <array_name>_unit syntax (i.e. emission_unit or wavelength_unit)
If temperature is not given as a Quantity, it defaults to degrees Celsius.
"""
def __init__(self, filename=None, cmds=None, **kwargs):
filename = find_file(filename)
self.meta = {"filename": filename,
"temperature": -270 * u.deg_C, # deg C
"emission_unit": "",
"wavelength_unit": u.um}
self.cmds = cmds
self.table = Table()
if filename is not None and Path(filename).exists():
self.table = ioascii.read(filename)
tbl_meta = convert_table_comments_to_dict(self.table)
if isinstance(tbl_meta, dict):
self.meta.update(tbl_meta)
self.meta.update(kwargs)
@property
def area(self):
if "area" in self.meta:
the_area = self.from_meta("area", u.m ** 2)
elif "outer" in self.meta:
outer_diameter = self.from_meta("outer", u.m)
the_area = np.pi * (0.5 * outer_diameter) ** 2
if "inner" in self.meta:
inner_diameter = self.from_meta("inner", u.m)
the_area -= np.pi * (0.5 * inner_diameter) ** 2
else:
the_area = None
return the_area
@property
def mirror_angle(self):
if "angle" in self.meta:
mirr_angle = self.from_meta("angle", u.deg)
else:
mirr_angle = 0 * u.deg
return mirr_angle
@property
def wavelength(self):
return self._get_array("wavelength")
@property
def throughput(self):
return self._get_ter_property(self.meta.get("action", "transmission"))
@property
def transmission(self):
return self._get_ter_property("transmission")
@property
def emissivity(self):
return self._get_ter_property("emissivity")
@property
def reflection(self):
return self._get_ter_property("reflection")
@property
def emission(self):
"""
Look for an emission array in self.meta.
If it doesn't find this, it defaults to creating a blackbody and
multiplies this by the emissivity. Assumption is that
``self.meta["temperature"]`` is in ``deg_C``, unless it is a
``u.Quantity`` with temperature unit attached. Return units are in
``PHOTLAM arcsec^-2``, even though ``arcsec^-2`` is not given.
"""
flux = self._get_array("emission")
if flux is None and "temperature" not in self.meta:
return None
if flux is not None:
wave = self._get_array("wavelength")
flux = make_emission_from_array(flux, wave, meta=self.meta)
elif "temperature" in self.meta:
emiss = self.emissivity # SpectralElement [0..1]
temp = from_currsys(self.meta["temperature"], self.cmds)
if not isinstance(temp, u.Quantity):
temp = quantify(temp, u.deg_C)
temp = temp.to(u.Kelvin, equivalencies=u.temperature())
flux = make_emission_from_emissivity(temp, emiss)
has_solid_angle = extract_type_from_unit(flux.meta["solid_angle"],
"solid angle")[1] != u.Unit("")
if flux is not None and has_solid_angle:
conversion_factor = flux.meta["solid_angle"].to(u.arcsec ** -2)
flux = flux * conversion_factor
flux.meta["solid_angle"] = u.arcsec ** -2
flux.meta["history"].append(
f"Converted to arcsec-2: {conversion_factor}")
if flux is not None and "rescale_emission" in self.meta:
dic = from_currsys(self.meta["rescale_emission"], self.cmds)
amplitude = dic["value"] * u.Unit(dic["unit"])
filter_name = dic["filter_name"]
if "filename_format" in dic:
filter_name = dic["filename_format"].format(filter_name)
flux = ter_utils.scale_spectrum(flux, filter_name, amplitude)
return flux
def _get_ter_property(self, ter_property, fmt="synphot"):
"""
Look for arrays for transmission, emissivity, or reflection.
Parameters
----------
ter_property : str
``transmission``, ``emissivity``, ``reflection``
fmt : str
Returns
-------
response_curve : ``synphot.SpectralElement``
"""
compliment_names = ["transmission", "emissivity", "reflection"]
ii = np.where([ter_property == name for name in compliment_names])[0][0]
compliment_names.pop(ii)
wave = self._get_array("wavelength")
value_arr = self._get_array(ter_property)
if value_arr is None:
value_arr = self._compliment_array(*compliment_names)
if value_arr is not None and wave is not None and fmt == "synphot":
response_curve = SpectralElement(Empirical1D, points=wave,
lookup_table=value_arr)
elif fmt == "array":
response_curve = value_arr
else:
response_curve = None
logger.warning("Both wavelength and %s must be set", ter_property)
return response_curve
def _compliment_array(self, colname_a, colname_b):
"""
Return an complimentary array using: ``a + b + c = 1``.
E.g. ``Emissivity = 1 - (Transmission + Reflection)``
Parameters
----------
colname_a : str
Name of the first TER property to look for
colname_b
Name of the second TER property to look for
Returns
-------
actual : ``synphot.SpectralElement``
Complimentary spectrum to those given
"""
compliment_a = self._get_array(colname_a)
compliment_b = self._get_array(colname_b)
if compliment_a is not None and compliment_b is not None:
actual = 1 * compliment_a.unit - (compliment_a + compliment_b)
elif compliment_a is not None and compliment_b is None:
actual = 1 * compliment_a.unit - compliment_a
elif compliment_b is not None and compliment_a is None:
actual = 1 * compliment_b.unit - compliment_b
elif compliment_b is None and compliment_a is None:
actual = None
return actual
def _get_array(self, colname):
"""
Look for an array in either the self.meta or self.table attributes.
Order of search goes: 1. self.meta, 2. self.table
Parameters
----------
colname : str
Array column (or key) name
Returns
-------
val_out : array-like Quantity
"""
if colname in self.meta:
val = self.meta[colname]
elif colname in self.table.colnames:
val = self.table[colname].data
else:
logger.debug("%s not found in either '.meta' or '.table': [%s]",
colname, self.meta.get("name", self.meta["filename"]))
return None
col_units = colname + "_unit"
if isinstance(val, u.Quantity):
units = val.unit
elif col_units in self.meta:
units = u.Unit(self.meta[col_units])
else:
units = u.Unit("")
if isinstance(val, u.Quantity):
val_out = val.to(units)
elif isinstance(val, (list, tuple, np.ndarray)):
val_out = val * units
elif val is None:
val_out = None
else:
raise ValueError(f"{colname} must be of type: Quantity, array, "
f"list, tuple, but is {type(val)}")
return val_out
def __repr__(self):
msg = (f"{self.__class__.__name__}({self.meta['filename']}, "
f"**{self.meta!r})")
return msg
def __str__(self):
meta = self.meta
name = meta["name"] if "name" in meta else meta["filename"]
cols = "".join([col[0].upper() for col in self.table.colnames])
msg = f"SpectralSurface [{cols}] \"{name}\""
return msg