Skip to content

Commit

Permalink
Merge pull request #22 from tomasstolker/photon_detector
Browse files Browse the repository at this point in the history 
Support for filters from photon-counting detectors
  • Loading branch information
Tomas Stolker committed Jul 30, 2020
2 parents 163b0df + 31c1df6 commit 532bcd3
Show file tree
Hide file tree
Showing 6 changed files with 197 additions and 51 deletions.
62 changes: 38 additions & 24 deletions species/analysis/photometry.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,9 @@

class SyntheticPhotometry:
"""
Class for calculating synthetic photometry from a spectrum.
Class for calculating synthetic photometry from a spectrum and also for conversion between
magnitudes and fluxes. Note that depending on the detector type (energy- or photon-counting)
the integral for the filter-weighted flux contains an additional wavelength factor.
"""

@typechecked
Expand All @@ -32,7 +34,7 @@ def __init__(self,
----------
filter_name : str
Filter name as listed in the database. Filters from the SVO Filter Profile Service are
downloaded and added to the database.
automatically downloaded and added to the database.
Returns
-------
Expand All @@ -53,6 +55,9 @@ def __init__(self,

self.database = config['species']['database']

read_filt = read_filter.ReadFilter(self.filter_name)
self.det_type = read_filt.detector_type()

@typechecked
def zero_point(self) -> np.float64:
"""
Expand Down Expand Up @@ -113,13 +118,13 @@ def spectrum_to_flux(self,
Wavelength points (um).
flux : np.ndarray
Flux (W m-2 um-1).
error : np.ndarray
Uncertainty (W m-2 um-1). Not used if set to None.
error : np.ndarray, None
Uncertainty (W m-2 um-1). Not used if set to ``None``.
threshold : float, None
Transmission threshold (value between 0 and 1). If the minimum transmission value is
larger than the threshold, a NaN is returned. This will happen if the input spectrum
does not cover the full wavelength range of the filter profile. Not used if set to
None.
``None``.
Returns
-------
Expand All @@ -129,6 +134,11 @@ def spectrum_to_flux(self,
Uncertainty (W m-2 um-1).
"""

if error is not None:
# The error calculation requires the original spectrum because spectrum_to_flux is used
wavel_error = wavelength.copy()
flux_error = flux.copy()

if self.filter_interp is None:
transmission = read_filter.ReadFilter(self.filter_name)
self.filter_interp = transmission.interpolate_filter()
Expand All @@ -150,7 +160,6 @@ def spectrum_to_flux(self,
'point. Photometry is set to NaN.')

else:

if threshold is None and (wavelength[0] > self.wavel_range[0] or
wavelength[-1] < self.wavel_range[1]):

Expand All @@ -166,9 +175,6 @@ def spectrum_to_flux(self,
wavelength = wavelength[indices]
flux = flux[indices]

if error is not None:
error = error[indices]

transmission = self.filter_interp(wavelength)

if threshold is not None and \
Expand All @@ -189,8 +195,15 @@ def spectrum_to_flux(self,
indices = np.isnan(transmission)
indices = np.logical_not(indices)

integrand1 = transmission[indices]*flux[indices]
integrand2 = transmission[indices]
if self.det_type == 'energy':
# Energy counting detector
integrand1 = transmission[indices]*flux[indices]
integrand2 = transmission[indices]

if self.det_type == 'photon':
# Photon counting detector
integrand1 = wavelength[indices]*transmission[indices]*flux[indices]
integrand2 = wavelength[indices]*transmission[indices]

integral1 = np.trapz(integrand1, wavelength[indices])
integral2 = np.trapz(integrand2, wavelength[indices])
Expand All @@ -201,11 +214,12 @@ def spectrum_to_flux(self,
phot_random = np.zeros(200)

for i in range(200):
spec_random = flux + np.random.normal(loc=0.,
scale=1.,
size=wavelength.shape[0])*error
# Use the original spectrum size (i.e. wavel_error and flux_error)
spec_random = flux_error + np.random.normal(loc=0.,
scale=1.,
size=wavel_error.shape[0])*error

phot_random[i] = self.spectrum_to_flux(wavelength,
phot_random[i] = self.spectrum_to_flux(wavel_error,
spec_random,
error=None,
threshold=threshold)[0]
Expand Down Expand Up @@ -240,20 +254,20 @@ def spectrum_to_magnitude(self,
error : np.ndarray, list(np.ndarray), None
Uncertainty (W m-2 um-1).
distance : tuple(float, float), None
Distance and uncertainty (pc). No absolute magnitude is calculated if set to None.
No error on the absolute magnitude is calculated if the uncertainty is set to None.
Distance and uncertainty (pc). No absolute magnitude is calculated if set to ``None``.
No error on the absolute magnitude is calculated if the uncertainty is set to ``None``.
threshold : float, None
Transmission threshold (value between 0 and 1). If the minimum transmission value is
larger than the threshold, a NaN is returned. This will happen if the input spectrum
does not cover the full wavelength range of the filter profile. Not used if set to
None.
``None``.
Returns
-------
tuple(float, float)
Apparent magnitude and uncertainty (mag).
Apparent magnitude and uncertainty.
tuple(float, float)
Absolute magnitude and uncertainty (mag).
Absolute magnitude and uncertainty.
"""

zp_flux = self.zero_point()
Expand Down Expand Up @@ -312,9 +326,9 @@ def magnitude_to_flux(self,
Parameters
----------
magnitude : float
Magnitude (mag).
Magnitude.
error : float, None
Error (mag). Not used if set to ``None``.
Error on the magnitude. Not used if set to ``None``.
zp_flux : float, None
Zero-point flux (W m-2 um-1). The value is calculated if set to ``None``.
Expand Down Expand Up @@ -369,9 +383,9 @@ def flux_to_magnitude(self,
Returns
-------
tuple(float, float), tuple(np.ndarray, np.ndarray)
Apparent magnitude and uncertainty (mag).
Apparent magnitude and uncertainty.
tuple(float, float), tuple(np.ndarray, np.ndarray)
Absolute magnitude and uncertainty (mag).
Absolute magnitude and uncertainty.
"""

zp_flux = self.zero_point()
Expand Down
25 changes: 17 additions & 8 deletions species/data/database.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -209,19 +209,26 @@ def add_dust(self) -> None:
@typechecked
def add_filter(self,
filter_name: str,
filename: Optional[str] = None) -> None:
filename: Optional[str] = None,
detector_type: str = 'photon') -> None:
"""
Function for adding a filter profile to the database, either from the SVO Filter profile
Service or input file.
Service or from an input file.
Parameters
----------
filter_name : str
Filter name from the SVO Filter Profile Service (e.g., 'Paranal/NACO.Lp').
Filter name from the SVO Filter Profile Service (e.g., 'Paranal/NACO.Lp') or a
user-defined name if a ``filename`` is specified.
filename : str
Filename of the filter profile. The first column should contain the wavelength
(um) and the second column the transmission. The profile is downloaded from the SVO
Filter Profile Service if set to None.
(um) and the second column the fractional transmission. The profile is downloaded from
the SVO Filter Profile Service if the argument of ``filename`` is set to ``None``.
detector_type : str
The detector type ('photon' or 'energy'). The argument is only used if a ``filename``
is provided. Otherwise, for filters that are fetched from the SVO website, the detector
type is read from the SVO data. The detector type determines if a wavelength factor
is included in the integral for the synthetic photometry.
Returns
-------
Expand Down Expand Up @@ -250,11 +257,13 @@ def add_filter(self,
transmission = data[:, 1]

else:
wavelength, transmission = filters.download_filter(filter_name)
wavelength, transmission, detector_type = filters.download_filter(filter_name)

if wavelength is not None and transmission is not None:
h5_file.create_dataset(f'filters/{filter_name}',
data=np.column_stack((wavelength, transmission)))
dset = h5_file.create_dataset(f'filters/{filter_name}',
data=np.column_stack((wavelength, transmission)))

dset.attrs['det_type'] = str(detector_type)

h5_file.close()

Expand Down
29 changes: 16 additions & 13 deletions species/data/filters.py
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
"""
Module for downloading filter data from the SVO website.
Module for downloading filter data from the website of the SVO Filter Profile Service.
"""

import os
Expand All @@ -15,21 +15,26 @@


@typechecked
def download_filter(filter_id: str) -> Tuple[Optional[np.ndarray], Optional[np.ndarray]]:
def download_filter(filter_id: str) -> Tuple[Optional[np.ndarray],
Optional[np.ndarray],
Optional[str]]:
"""
Function for downloading filter transmission data from the SVO Filter Profile Service.
Function for downloading filter profile data from the SVO Filter Profile Service.
Parameters
----------
filter_id : str
Filter name as listed on the SVO website.
Filter name as listed on the website of the SVO Filter Profile Service
(see http://svo2.cab.inta-csic.es/svo/theory/fps/).
Returns
-------
np.ndarray
Wavelength (um).
np.ndarray
Transmission.
Fractional transmission.
str
Detector type ('energy' or 'photon').
"""

if filter_id == 'LCO/VisAO.Ys':
Expand All @@ -41,6 +46,9 @@ def download_filter(filter_id: str) -> Tuple[Optional[np.ndarray], Optional[np.n
wavelength = wavelength[:-7]
transmission = transmission[:-7]

# Not sure if energy- or photon-counting detector
det_type = 'energy'

os.remove('VisAO_Ys_filter_curve.dat')

else:
Expand All @@ -56,6 +64,7 @@ def download_filter(filter_id: str) -> Tuple[Optional[np.ndarray], Optional[np.n
except IndexError:
wavelength = None
transmission = None
det_type = None

warnings.warn(f'Filter \'{filter_id}\' is not available on the SVO Filter Profile '
f'Service.')
Expand All @@ -77,15 +86,9 @@ def download_filter(filter_id: str) -> Tuple[Optional[np.ndarray], Optional[np.n
det_type = 'photon'

except KeyError:
# Energy-counting detector if the DetectorType key is not present
det_type = 'energy'

if det_type == 'photon':
raise ValueError(f'The detector of the {filter_id} filter is a photon counter, '
f'therefore the transmission profile has to be multiplied with '
f'the wavelength when calculating average fluxes. This is '
f'currently not implemented in species. Please open an issue on '
f'Github if needed.')

wavelength *= 1e-4 # (um)

os.remove('filter.xml')
Expand Down Expand Up @@ -120,4 +123,4 @@ def download_filter(filter_id: str) -> Tuple[Optional[np.ndarray], Optional[np.n
wavelength = wavelength[indices]
transmission = transmission[indices]

return wavelength, transmission
return wavelength, transmission, det_type
33 changes: 27 additions & 6 deletions species/read/read_filter.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,16 @@ def __init__(self,

self.database = config['species']['database']

h5_file = h5py.File(self.database, 'r')

if 'filters' not in h5_file or self.filter_name not in h5_file['filters']:
h5_file.close()
species_db = database.Database()
species_db.add_filter(self.filter_name)
h5_file = h5py.File(self.database, 'r')

h5_file.close()

@typechecked
def get_filter(self) -> np.ndarray:
"""
Expand All @@ -59,12 +69,6 @@ def get_filter(self) -> np.ndarray:

h5_file = h5py.File(self.database, 'r')

if 'filters' not in h5_file or self.filter_name not in h5_file['filters']:
h5_file.close()
species_db = database.Database()
species_db.add_filter(self.filter_name)
h5_file = h5py.File(self.database, 'r')

data = np.asarray(h5_file[f'filters/{self.filter_name}'])

if data.shape[0] == 2 and data.shape[1] > data.shape[0]:
Expand Down Expand Up @@ -149,3 +153,20 @@ def filter_fwhm(self) -> float:
root2 = np.amin(diff[diff > 0.])

return root2 - root1

@typechecked
def detector_type(self) -> str:
"""
Function for returning the detector type.
Returns
-------
str
Detector type ('energy' or 'photon').
"""

with h5py.File(self.database, 'r') as h5_file:
dset = h5_file[f'filters/{self.filter_name}']
det_type = dset.attrs['det_type']

return det_type
Loading

0 comments on commit 532bcd3

Please sign in to comment.