# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
sbpy Activity Core Module
Core module functions and classes, especially for handling coma
geometries.
created on June 23, 2017
"""
__all__ = [
'Aperture',
'CircularAperture',
'AnnularAperture',
'RectangularAperture',
'GaussianAperture',
]
from abc import ABC, abstractmethod
import numpy as np
import astropy.units as u
from .. import data as sbd
from .. import units as sbu
[docs]class Aperture(ABC):
"""
Abstract base class for photometric apertures.
Notes
-----
The shape of the aperture must be passed as the first argument of
`__init__`, or else `as_length` and `as_angle` must be overridden.
"""
def __init__(self, dim):
if not dim.unit.is_equivalent((u.radian, u.meter)):
raise u.UnitTypeError(
'aperture must be defined with angles or lengths.')
self.dim = dim
def __str__(self):
"""Description of the aperture."""
# assumes preferred format for __repr__
return repr(self)[1:-1].replace('Aperture:', ' aperture,')
@abstractmethod
def __repr__(self):
"""Preferred format <ShapedAperture: size>"""
[docs] @sbd.dataclass_input(eph=sbd.Ephem)
@sbd.quantity_to_dataclass(eph=(sbd.Ephem, 'delta'))
def as_angle(self, eph):
"""This aperture in units of angle.
Parameters
----------
eph : dictionary-like, `~sbpy.data.Ephem`, or `~astropy.units.Quantity`
The observer-target distance (``delta``).
Returns
-------
aper
"""
dim = self.dim.to('arcsec', sbu.projected_size(eph))
return type(self)(dim)
[docs] @sbd.dataclass_input(eph=sbd.Ephem)
@sbd.quantity_to_dataclass(eph=(sbd.Ephem, 'delta'))
def as_length(self, eph):
"""This aperture in units of length.
Parameters
----------
eph : dictionary-like, `~sbpy.data.Ephem`, or `~astropy.units.Quantity`
The observer-target distance (``delta``).
Returns
-------
aper
"""
dim = self.dim.to('km', sbu.projected_size(eph))
return type(self)(dim)
[docs] @abstractmethod
def coma_equivalent_radius(self):
"""Circular aperture radius that yields same flux for a 1/ρ coma.
Returns
-------
rap : `~astropy.units.Quantity`
"""
[docs]class CircularAperture(Aperture):
"""Circular aperture projected at the distance of the target.
Parameters
----------
radius : `~astropy.units.Quantity`
Angular or projected linear radius for the aperture.
"""
def __init__(self, radius):
super().__init__(radius)
def __repr__(self):
return '<CircularAperture: radius {}>'.format(self.dim)
@property
def radius(self):
"""Aperture radius."""
return self.dim
[docs] def coma_equivalent_radius(self):
return self.radius
coma_equivalent_radius.__doc__ = Aperture.coma_equivalent_radius.__doc__
[docs]class AnnularAperture(Aperture):
"""Annular aperture projected at the distance of the target.
Parameters
----------
shape : `~astropy.units.Quantity`
A two-element `~astropy.units.Quantity` of angular or
projected linear size for the inner and outer radius of the
aperture.
"""
def __init__(self, shape):
if len(shape) != 2:
raise ValueError('shape must be 2-elements')
super().__init__(shape)
def __repr__(self):
return ('<AnnularAperture: radii {0[0].value:}–{0[1]:}>'
.format(self.dim))
@property
def shape(self):
"""Annulus inner and outer radii."""
return self.dim
[docs] def coma_equivalent_radius(self):
return max(self.dim) - min(self.dim)
coma_equivalent_radius.__doc__ = Aperture.coma_equivalent_radius.__doc__
[docs]class RectangularAperture(Aperture):
"""Rectangular aperture projected at the distance of the target.
Parameters
----------
shape : `~astropy.units.Quantity`
A two-element `~astropy.units.Quantity` of angular or
projected linear size for the width and height of the
aperture. The order is not significant.
"""
def __init__(self, shape):
if len(shape) != 2:
raise ValueError('shape must be 2-elements')
super().__init__(shape)
def __repr__(self):
return ("<RectangularAperture: dimensions {0[0].value:}×{0[1]:}>"
.format(self.dim))
@property
def shape(self):
"""Rectangle dimensions."""
return self.dim
[docs] def coma_equivalent_radius(self):
# Int_0^θ Int_0^r 1/r * r * dr dθ
# limits on r depend on θ; eq. of a line: x * sec(θ)
# --> Int_0^θ Int_0^x*sec(θ) dr dθ
# --> Int_0^θ x*sec(θ) dθ
# --> x * log(tan(θ) + sec(θ))
# First, integrate the 1/rho distribution over the first
# "octant" of the rectangle in polar coordinates. The
# azimuthal limits are 0 to arctan(y / x). Also, x, and y are
# the full rectangle dimensions, so they must be halved.
# th = np.arctan(y / x)
# I = (x / 2) * log(tan(th) + sec(th))
# sec(th) = cos(th)**-1
# cos(arctan(y / x)) = 1 / sqrt(1 + (y / x)**2)
# sec(th) = sqrt(1 + (y / x)**2)
# I1 = x / 2 * np.log(y / x + np.sqrt(1 + (y / x)**2))
# Then, integrate the second "octant": th = 0 to arctan(x / y)
# I2 = y / 2 * np.log(x / y + np.sqrt(1 + (x / y)**2))
# The two octants correspond to 1/4 the full area:
# I = 4 * (I1 + I2)
# For the circular aperture, the integral is 2 pi rho.
# rho = I / 2 / np.pi
# implement the above, moving constants around
x, y = self.shape
I1 = x * np.log(y / x + np.sqrt(1 + (y / x)**2))
I2 = y * np.log(x / y + np.sqrt(1 + (x / y)**2))
return (I1 + I2) / np.pi
coma_equivalent_radius.__doc__ = Aperture.coma_equivalent_radius.__doc__
[docs]class GaussianAperture(Aperture):
"""Gaussian-shaped aperture, e.g., for radio observations.
The aperture is normalized to 1.0 at the center.
Parameters
----------
sigma : `~astropy.units.Quantity`, optional
The width of the Gaussian beam (square-root of the variance)
as an angular or projected size.
fwhm : `~astropy.units.Quantity`, optional
The full-width at half-maximum of the Gaussian beam as an
angular or projected size.
Notes
-----
One of `sigma` or `fwhm` is required.
"""
def __init__(self, sigma=None, fwhm=None):
if (sigma is None) and (fwhm is None):
raise ValueError('One of `sigma` or `fwhm` must be defined')
if sigma is not None:
super().__init__(sigma)
else:
super().__init__(fwhm / 2.3548200450309493)
def __repr__(self):
return "<GaussianAperture: 1-σ width {}>".format(self.dim)
@property
def sigma(self):
"""Beam Gaussian width."""
return self.dim
@property
def fwhm(self):
"""Beam full-width at half-maximum."""
return self.dim * 2.3548200450309493
[docs] @sbd.dataclass_input
@sbd.quantity_to_dataclass(eph=(sbd.Ephem, 'delta'))
def __call__(self, rho, eph: sbd.Ephem=None):
"""Evaluate the aperture.
Parameters
----------
rho : `~astropy.units.Quantity`
Position to evaluate, in units of length or angle.
eph : dictionary-like, `~sbpy.data.Ephem`, or `~astropy.units.Quantity`, optional
The observer-target distance (``delta``). Use ``eph`` to
convert between angles and lengths, as needed.
"""
if eph is not None:
equiv = sbu.projected_size(eph)
else:
equiv = []
x = rho.to(self.dim.unit, equiv)
# normalize to 1.0 at the center
return np.exp(-x**2 / self.sigma**2 / 2)
[docs] def coma_equivalent_radius(self):
# This beam is normalized to 1.0 at the center.
return np.sqrt(np.pi / 2) * self.sigma
coma_equivalent_radius.__doc__ = Aperture.coma_equivalent_radius.__doc__