"""
Utility functions for photometry.
"""
import math
import warnings
import matplotlib.pyplot as plt
import spectres
import numpy as np
from species.core import box
from species.read import read_model, read_calibration, read_filter, read_planck
# from species.read import read_model, read_calibration, read_filter, read_planck, read_radtrans
[docs]def multi_photometry(datatype,
spectrum,
filters,
parameters):
"""
Parameters
----------
datatype : str
Data type ('model' or 'calibration').
spectrum : str
Spectrum name (e.g., 'drift-phoenix').
filters : tuple(str, )
Filter IDs.
parameters : dict
Parameters and values for the spectrum
Returns
-------
species.core.box.SynphotBox
Box with synthetic photometry.
"""
print('Calculating synthetic photometry...', end='', flush=True)
flux = {}
if datatype == 'model':
for item in filters:
if spectrum == 'planck':
readmodel = read_planck.ReadPlanck(filter_name=item)
else:
readmodel = read_model.ReadModel(spectrum, filter_name=item)
try:
flux[item] = readmodel.get_flux(parameters)[0]
except IndexError:
flux[item] = np.nan
warnings.warn(f'The wavelength range of the {item} filter does not match with '
f'the wavelength coverage of {spectrum}. The flux is set to NaN.')
elif datatype == 'calibration':
for item in filters:
readcalib = read_calibration.ReadCalibration(spectrum, filter_name=item)
flux[item] = readcalib.get_flux(parameters)[0]
print(' [DONE]')
return box.create_box('synphot', name='synphot', flux=flux)
[docs]def apparent_to_absolute(app_mag,
distance):
"""
Function for converting an apparent magnitude into an absolute magnitude. The uncertainty on
the distance is propagated into the uncertainty on the absolute magnitude.
Parameters
----------
app_mag : tuple(float, float), tuple(numpy.ndarray, numpy.ndarray)
Apparent magnitude and uncertainty (mag). The returned error on the absolute magnitude
is set to None if the error on the apparent magnitude is set to None, for example
``app_mag=(15., None)``.
distance : tuple(float, float), tuple(numpy.ndarray, numpy.ndarray)
Distance and uncertainty (pc). The error is not propagated into the error on the absolute
magnitude if set to None, for example ``distance=(20., None)``.
Returns
-------
float, numpy.ndarray
Absolute magnitude (mag).
float, numpy.ndarray, None
Uncertainty (mag).
"""
abs_mag = app_mag[0] - 5.*np.log10(distance[0]) + 5.
if app_mag[1] is not None and distance[1] is not None:
dist_err = distance[1] * (5./(distance[0]*math.log(10.)))
abs_err = math.sqrt(app_mag[1]**2 + dist_err**2)
elif app_mag[1] is not None and distance[1] is None:
abs_err = app_mag[1]
else:
abs_err = None
return abs_mag, abs_err
[docs]def absolute_to_apparent(abs_mag,
distance):
"""
Function for converting an absolute magnitude into an apparent magnitude.
Parameters
----------
abs_mag : tuple(float, float), tuple(numpy.ndarray, numpy.ndarray)
Absolute magnitude and uncertainty (mag). The same uncertainty is used for the
apparent magnitude.
distance : tuple(float, float), tuple(numpy.ndarray, numpy.ndarray)
Distance and uncertainty (pc).
Returns
-------
float, numpy.ndarray
Apparent magnitude (mag).
float, numpy.ndarray, None
Uncertainty (mag).
"""
app_mag = abs_mag[0] + 5.*np.log10(distance[0]) - 5.
return app_mag, abs_mag[1]
[docs]def get_residuals(datatype,
spectrum,
parameters,
filters,
objectbox,
inc_phot=True,
inc_spec=False,
**kwargs_radtrans):
"""
Parameters
----------
datatype : str
Data type ('model' or 'calibration').
spectrum : str
Name of the atmospheric model or calibration spectrum.
parameters : dict
Parameters and values for the spectrum
filters : tuple(str, )
Filter IDs. All available photometry of the object is used if set to None.
objectbox : species.core.box.ObjectBox
Box with the photometry and/or spectra of an object. A scaling and/or error inflation of
the spectra should be applied with :func:`~species.util.read_util.update_spectra`
beforehand.
inc_phot : bool
Include photometry.
inc_spec : bool
Include spectrum.
Keyword arguments
-----------------
kwargs_radtrans : dict
Dictionary with the keyword arguments for the ``ReadRadtrans`` object, containing
``line_species``, ``cloud_species``, and ``scattering``.
Returns
-------
species.core.box.ResidualsBox
Box with the photometry and/or spectrum residuals.
"""
if filters is None:
filters = objectbox.filters
if inc_phot:
model_phot = multi_photometry(datatype=datatype,
spectrum=spectrum,
filters=filters,
parameters=parameters)
res_phot = np.zeros((2, len(objectbox.flux)))
for i, item in enumerate(filters):
transmission = read_filter.ReadFilter(item)
res_phot[0, i] = transmission.mean_wavelength()
res_phot[1, i] = (objectbox.flux[item][0]-model_phot.flux[item])/objectbox.flux[item][1]
else:
res_phot = None
if inc_spec:
res_spec = {}
readmodel = None
for key in objectbox.spectrum:
wavel_range = (0.9*objectbox.spectrum[key][0][0, 0],
1.1*objectbox.spectrum[key][0][-1, 0])
wl_new = objectbox.spectrum[key][0][:, 0]
spec_res = objectbox.spectrum[key][3]
if spectrum == 'planck':
readmodel = read_planck.ReadPlanck(wavel_range=wavel_range)
model = readmodel.get_spectrum(model_param=parameters, spec_res=1000.)
flux_new = spectres.spectres(wl_new, model.wavelength, model.flux)
else:
if spectrum == 'petitradtrans':
# TODO change back
pass
# radtrans = read_radtrans.ReadRadtrans(line_species=kwargs_radtrans['line_species'],
# cloud_species=kwargs_radtrans['cloud_species'],
# scattering=kwargs_radtrans['scattering'],
# wavel_range=wavel_range)
#
# model = radtrans.get_model(parameters, spec_res=None)
#
# # separate resampling to the new wavelength points
#
# flux_new = spectres.spectres(wl_new, model.wavelength, model.flux)
else:
readmodel = read_model.ReadModel(spectrum, wavel_range=wavel_range)
# resampling to the new wavelength points is done in teh get_model function
model = readmodel.get_model(parameters,
spec_res=spec_res,
wavel_resample=wl_new,
smooth=True)
flux_new = model.flux
res_tmp = (objectbox.spectrum[key][0][:, 1]-flux_new)/objectbox.spectrum[key][0][:, 2]
res_spec[key] = np.column_stack([wl_new, res_tmp])
else:
res_spec = None
print('Calculating residuals... [DONE]')
print('Residuals (sigma):')
if res_phot is not None:
for i, item in enumerate(filters):
print(f' - {item}: {res_phot[1, i]:.2f}')
if res_spec is not None:
for key in objectbox.spectrum:
print(f' - {key}: min: {np.amin(res_spec[key]):.2f}, '
f'max: {np.amax(res_spec[key]):.2f}')
return box.create_box(boxtype='residuals',
name=objectbox.name,
photometry=res_phot,
spectrum=res_spec)