Source code for sbpy.activity.gas.core

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""activity.gas core

"""

__all__ = [
    "photo_lengthscale",
    "photo_timescale",
    "fluorescence_band_strength",
    "Haser",
    "VectorialModel",
]

from abc import ABC, abstractmethod

# distutils is deprecated in python 3.10 and will be removed in 3.12 (PEP 632).
# Migration from distutils.log -> logging
from dataclasses import dataclass
from typing import Callable, Tuple

import numpy as np
import astropy.units as u

try:
    import scipy
    from scipy import special
    from scipy.integrate import quad, dblquad, romberg
    from scipy.interpolate import CubicSpline, PPoly
except ImportError:
    scipy = None
    PPoly = None

from ... import bib
from ... import data as sbd
from ... import units as sbu
from ...utils.decorators import requires
from ..core import (
    Aperture,
    RectangularAperture,
    GaussianAperture,
    AnnularAperture,
    CircularAperture,
)


[docs] def photo_lengthscale(species, source=None): """Photodissociation lengthscale for a gas species. Parameters ---------- species : string The species to look up. source : string, optional Retrieve values from this source (case insensitive). See references for keys. Returns ------- gamma : `~astropy.units.Quantity` The lengthscale at 1 au. Examples -------- >>> from sbpy.activity import photo_lengthscale >>> gamma = photo_lengthscale('OH') References ---------- [CS93] H2O and OH from Table IV of Cochran & Schleicher 1993, Icarus 105, 235-253. Quoted for intermediate solar activity. """ from .data import photo_lengthscale as data default_sources = { "H2O": "CS93", "OH": "CS93", } if species not in data: summary = "" for k, v in sorted(data.items()): summary += "\n{} [{}]".format(k, ", ".join(v.keys())) raise ValueError( "Invalid species {}. Choose from:{}".format(species, summary)) gas = data[species] source = default_sources[species] if source is None else source if source not in gas: raise ValueError( "Source key {} not available for {}. Choose from: {}".format( source, species, ", ".join(gas.keys()) ) ) gamma, bibcode = gas[source] bib.register(photo_lengthscale, bibcode) return gamma
[docs] def photo_timescale(species, source=None): """Photodissociation timescale for a gas species. Parameters ---------- species : string Species to look up. source : string, optional Retrieve values from this source. See references for keys. Returns ------- tau : `~astropy.units.Quantity` The timescale at 1 au. May be a two-element array: (quiet Sun, active Sun). Examples -------- >>> from sbpy.activity import photo_timescale >>> tau = photo_timescale('OH') References ---------- [CS93] Table IV of Cochran & Schleicher 1993, Icarus 105, 235-253. Quoted for intermediate solar activity. [C94] Crovisier 1994, JGR 99, 3777-3781. [CE83] Crovisier & Encrenaz 1983, A&A 126, 170-182. [H92] Huebner et al. 1992, Astroph. & Space Sci. 195, 1-294. """ from .data import photo_timescale as data default_sources = { "H2O": "CS93", "OH": "CS93", "HCN": "C94", "CH3OH": "C94", "H2CO": "C94", "CO2": "CE83", "CO": "CE83", "CN": "H92", } if species not in data: summary = "" for k, v in sorted(data.items()): summary += "\n{} [{}]".format(k, ", ".join(v.keys())) raise ValueError( "Invalid species {}. Choose from:{}".format(species, summary)) gas = data[species] source = default_sources[species] if source is None else source if source not in gas: raise ValueError( "Source key {} not available for {}. Choose from: {}".format( source, species, ", ".join(gas.keys()) ) ) tau, bibcode = gas[source] bib.register(photo_timescale, bibcode) return tau
[docs] @sbd.dataclass_input(eph=sbd.Ephem) def fluorescence_band_strength(species, eph=None, source=None): """Fluorescence band strength. Parameters ---------- species : string Species to look up. eph : `~astropy.units.Quantity`, `~sbpy.data.Ephem` or `dict` optional The target ephemeris. The strength is scaled to the given heliocentric distance, if present. Some species require heliocentric radial velocity ('rdot'). source : string, optional Retrieve values from this source (case insensitive). See references for keys. Returns ------- LN : `~astropy.units.Quantity` Luminosity per molecule, scaled to rh, if provided. Examples -------- >>> import astropy.units as u >>> from sbpy.activity import fluorescence_band_strength >>> >>> eph = {'rh': 1 * u.au, 'rdot': -1 * u.km / u.s} >>> LN = fluorescence_band_strength('OH 0-0', eph, 'SA88') >>> print(LN) # doctest: +FLOAT_CMP [1.54e-15] erg / s """ from .data import fluorescence_band_strength as data default_sources = { "OH 0-0": "SA88", "OH 1-0": "SA88", "OH 1-1": "SA88", "OH 2-2": "SA88", "OH 0-1": "SA88", "OH 0-2": "SA88", "OH 2-0": "SA88", "OH 2-1": "SA88", } if species not in data: raise ValueError( "No data available for {}. Choose one of: {}".format( species, ", ".join(data.keys()) ) ) band = data[species] source = default_sources[species] if source is None else source if source not in band: raise ValueError( "No source {} for {}. Choose one of: {}".format( source, species, ", ".join(band.keys()) ) ) LN, note, bibcode = band[source] if bibcode is not None: bib.register(fluorescence_band_strength, bibcode) return LN(eph)
class GasComa(ABC): """Abstract base class for gas coma models. Parameters ---------- Q : `~astropy.units.Quantity` Production rate, number per time. v : `~astropy.units.Quantity` Radial outflow speed, distance per time. """ @u.quantity_input(Q=(u.s**-1, u.mol / u.s), v=u.m / u.s) def __init__(self, Q, v): self.Q = Q self.v = v @u.quantity_input(r=u.m) def volume_density(self, r): """Coma volume density. Parameters ---------- r : `~astropy.units.Quantity` Linear distance to the nucleus. Returns ------- n : `~astropy.units.Quantity` Local number density. """ return self._volume_density(r.to_value("m")) / u.m**3 @sbd.dataclass_input(eph=sbd.Ephem) @sbd.quantity_to_dataclass(eph=(sbd.Ephem, "delta")) def column_density(self, rho, eph=None): """Coma column density at a projected distance from nucleus. Parameters ---------- rho : `~astropy.units.Quantity` Projected distance to the region of interest on the plane of the sky in units of length or angle. eph : dictionary-like, `~sbpy.data.Ephem`, `~astropy.units.Quantity`, optional Target-observer distance, or ephemeris with ``'delta'`` field. Required to convert rho to a projected size. Returns ------- sigma : `~astropy.units.Quantity` Coma column density along the line of sight at a distance rho. """ equiv = [] if eph is not None: equiv = sbu.projected_size(eph) rho = rho.to_value("m", equiv) return self._column_density(rho) / u.m**2 @sbd.dataclass_input(eph=sbd.Ephem) @sbd.quantity_to_dataclass(eph=(sbd.Ephem, "delta")) def total_number(self, aper, eph=None): """Total number of molecules in aperture. Parameters ---------- aper : `~astropy.units.Quantity`, `~sbpy.activity.Aperture` Observation aperture. May be a circular aperture radius with units of length or angle. eph : dictionary-like, `~sbpy.data.Ephem`, `~astropy.units.Quantity` Target-observer distance, or ephemeris with `delta`. Required if the aperture has angular units. Returns ------- N : float Total number of molecules within the aperture. """ if not isinstance(aper, Aperture): aper = CircularAperture(aper) if eph is not None: aper = aper.as_length(eph) return self._total_number(aper) def _total_number(self, aper): """Total number of molecules in aperture. Sub-classes of ``GasComa`` may override this method, instead of ``total_number``, which avoids reusing the boiler plate aperture conversions. Parameters ---------- aper : `~sbpy.activity.Aperture` Observation aperture in units of length. """ return self._integrate_column_density(aper)[0] @abstractmethod def _volume_density(self, r): """Unitless volume density function. Parameters ---------- r : float Linear distance to the nucleus in meters. Returns ------- n : float Local number density in inverse cubic-meters. """ @abstractmethod def _column_density(self, rho): """Unitless column density function. Parameters ---------- rho : float Projected distance of the region of interest on the plane of the sky in units of meters. Returns ------- sigma : float Coma column density along the line of sight at a distance rho in units of inverse square-meters. """ @requires("scipy") def _integrate_volume_density(self, rho, epsabs=1.49e-8): """Integrate volume density along the line of sight. Parameters ---------- rho : float Projected distance of the region of interest on the plane of the sky in units of meters epsabs : float, int, optional Absolute and relative error tolerance for integrals. See `scipy.integrate.quad`. Returns ------- sigma : float Coma column density along ``rho`` in units of inverse square-meters. err : float Estimated integration error. """ def f(s, rho2): r = np.sqrt(rho2 + s**2) return self._volume_density(r) # quad diverges integrating to infinity, but 1e6 × rho is good # enough limit = 30 points = rho * np.logspace(-4, 4, limit // 2) sigma, err = quad( f, 0, 1e6 * rho, args=(rho**2,), limit=limit, points=points, epsabs=epsabs ) # spherical symmetry sigma *= 2 err *= 2 return sigma, err @requires("scipy") def _integrate_column_density(self, aper, epsabs=1.49e-8): """Integrate column density over an aperture. Parameters ---------- aper : `~sbpy.activity.Aperture` Aperture, in units of length. epsabs : float, int, optional Absolute and relative error tolerance for integrals. See `scipy.integrate.quad` (circular, annular, Gaussian) and `~scipy.integrate.dblquad` (rectangular) for details. Returns ------- N : float Total number. err : float Estimated integration error. """ if isinstance(aper, (CircularAperture, AnnularAperture)): if isinstance(aper, CircularAperture): limits = (0, aper.radius.to_value("m")) else: limits = aper.shape.to_value("m") # integrate in polar coordinates def f(rho): """Column density integration in polar coordinates. rho in m, column_density in m**-2 """ return rho * self._column_density(rho) N, err = quad(f, *limits, epsabs=epsabs) N *= 2 * np.pi err *= 2 * np.pi elif isinstance(aper, RectangularAperture): shape = aper.shape.to_value("m") def f(rho, th): """Column density integration in polar coordinates. rho in m, column_density in m**-2 th is ignored (azimuthal symmetry) """ return rho * self._column_density(rho) # first "octant"; rho1 and rho2 are the limits of the # integration def rho1(th): "Lower limit" return 0 def rho2(th): "Upper limit (a line)" return shape[0] / 2 / np.cos(th) th = np.arctan(shape[1] / shape[0]) N1, err1 = dblquad(f, 0, th, rho1, rho2, epsabs=epsabs) # second "octant" def rho2(th): "Upper limit (a line)" return shape[1] / 2 / np.cos(th) th = np.arctan(shape[0] / shape[1]) N2, err2 = dblquad(f, 0, th, rho1, rho2, epsabs=epsabs) # N1 + N2 constitute 1/4th of the rectangle N = 4 * (N1 + N2) err = 4 * (err1 + err2) elif isinstance(aper, GaussianAperture): # integrate in polar coordinates def f(rho, sigma): """Column density integration in polar coordinates. rho and sigma in m, column_density in m**-2 """ return ( rho * np.exp(-(rho**2) / sigma**2 / 2) * self._column_density(rho) ) sigma = aper.sigma.to_value("m") N, err = quad(f, 0, np.inf, args=(sigma,), epsabs=epsabs) N *= 2 * np.pi err *= 2 * np.pi return N, err
[docs] class Haser(GasComa): """Haser coma model. Some functions require `scipy`. Parameters ---------- Q : `~astropy.units.Quantity` Production rate, per time. v : `~astropy.units.Quantity` Radial outflow speed, distance per time. parent : `~astropy.units.Quantity` Coma lengthscale of the parent species. daughter : `~astropy.units.Quantity`, optional Coma lengthscale of the daughter species. References ---------- Haser 1957, Bulletin de la Societe Royale des Sciences de Liege 43, 740. Newburn and Johnson 1978, Icarus 35, 360-368. """ @bib.cite({"model": "1957BSRSL..43..740H"}) @u.quantity_input(parent=u.m, daughter=u.m) def __init__(self, Q, v, parent, daughter=None): super().__init__(Q, v) self.parent = parent self.daughter = daughter def _volume_density(self, r): n = (self.Q / self.v).to_value("1/m") / r**2 / 4 / np.pi parent = self.parent.to_value("m") if self.daughter is None or self.daughter == 0: # parent only n *= np.exp(-r / parent) else: daughter = self.daughter.to_value("m") n *= ( daughter / (parent - daughter) * (np.exp(-r / parent) - np.exp(-r / daughter)) ) return n def _iK0(self, x): """Integral of the modified Bessel function of 2nd kind, 0th order.""" return special.iti0k0(x)[1] def _K1(self, x): """Modified Bessel function of 2nd kind, 1st order.""" return special.k1(x) @requires("scipy") @bib.cite({"model": "1978Icar...35..360N"}) def _column_density(self, rho): sigma = (self.Q / self.v).to_value("1/m") / rho / 2 / np.pi parent = self.parent.to_value("m") if self.daughter is None or self.daughter == 0: sigma *= np.pi / 2 - self._iK0(rho / parent) else: daughter = self.daughter.to_value("m") sigma *= ( daughter / (parent - daughter) * (self._iK0(rho / daughter) - self._iK0(rho / parent)) ) return sigma @requires("scipy") def _total_number(self, aper): # Inspect aper and handle as appropriate if isinstance(aper, (RectangularAperture, GaussianAperture)): return self._integrate_column_density(aper)[0] elif isinstance(aper, AnnularAperture): N0 = self._total_number(CircularAperture(aper.shape[0])) N1 = self._total_number(CircularAperture(aper.shape[1])) return N1 - N0 # Solution for the circular aperture of radius rho: bib.register(self.total_number, {"model": "1978Icar...35..360N"}) rho = aper.radius parent = self.parent.to(rho.unit) x = (rho / parent).to_value(u.dimensionless_unscaled) N = (self.Q * rho / self.v).to_value(u.dimensionless_unscaled) if self.daughter is None or self.daughter == 0: N *= 1 / x - self._K1(x) + np.pi / 2 - self._iK0(x) else: daughter = self.daughter.to(rho.unit) y = (rho / daughter).to_value("") N *= (daughter / (parent - daughter)).to_value("") * ( self._iK0(y) - self._iK0(x) + x**-1 - y**-1 + self._K1(y) - self._K1(x) ) return N
@dataclass class VMParent: """ Physical information about the parent necessary for the vectorial model Parameters ---------- v_outflow : float See VectorialModel documentation. tau_d : float See VectorialModel documentation. tau_T : float See VectorialModel documentation. sigma : float See VectorialModel documentation. """ v_outflow: float tau_d: float tau_T: float sigma: float @dataclass class VMFragment: """ Physical information about the fragment necessary for the vectorial model Parameters ---------- v_photo : float See VectorialModel documentation. tau_T : float See VectorialModel documentation. """ v_photo: float tau_T: float @dataclass class VMGridParams: """ Vectorial model gridding parameters to control how finely the space around the comet should be gridded, and how detailed the outflow sampling should be. Parameters ---------- radial_points : int See VectorialModel documentation. angular_points : int See VectorialModel documentation. radial_substeps : int See VectorialModel documentation. """ radial_points: int angular_points: int radial_substeps: int @dataclass class VMParams: """ Vectorial model parameters unrelated to physical inputs, dealing with the model's assumptions and detalis of the calculations. Parameters ---------- parent_destrution_level : float See VectorialModel documentation. fragment_destruction_level : float See VectorialModel documentation. max_fragment_lifetimes : float See VectorialModel documentation. """ parent_destruction_level: float fragment_destruction_level: float max_fragment_lifetimes: float @dataclass class VMFragmentSputterPolar: """ Describes the distribution of fragment volume density (``fragment_density[i][j]``) as function of (``rs[i]``, ``thetas[j]``) in a spherical coordinate system, given a column of parent molecules flowing outward along the azimuthal (z) axis. Parameters ---------- rs : `np.ndarray` List of radii (`astropy.units.Quantity`). thetas: `np.ndarray` List of polar angles theta (radians) in a spherical coordinate system. fragment_density: `np.ndarray` List of fragment densities at the corresponding ``rs[i]`` and ``thetas[j]``. """ rs: np.ndarray thetas: np.ndarray fragment_density: np.ndarray @dataclass class VMResult: """ Dataclass to hold a collection of vectorial model details and results in a language- and model-independent way. Parameters ---------- volume_density_grid : `np.ndarray` List of radii (`~astropy.units.Quantity`) that the model used for volume density calculations. volume_density : `np.ndarray` List of volume densities (`~astropy.units.Quantity`) computed at the corresponding radius: at radius = ``volume_density_grid[i]``, the volume density is ``volume_density[i]``. column_density_grid : `np.ndarray` List of radii (`~astropy.units.Quantity`) that the model used for column density calculations. column_density : `np.ndarray` List of column densities (`~astropy.units.Quantity`), computed at the corresponding radius: at radius = ``column_density_grid[i]``, the column density is ``column_density[i]``. fragment_sputter : `VMFragmentSputterPolar` Describes the distribution of fragment volume density as function of (r, theta) in a spherical coordinate system, given a column of parent molecules flowing outward along the azimuthal (z) axis. solid_angle_sputter : `VMFragmentSputterPolar` Similar to ``fragment_sputter``, but adjusted by a factor of ``sin(theta)``. volume_density_interpolation : `scipy.interpolate.PPoly` Function that takes radius as a float in meters (no astropy units) and returns the volume density in 1/m^3. Not reliable beyond ``max_grid_radius``, and questionable for radii less than ``collision_sphere_radius``. column_density_interpolation : `scipi.interpolate.PPoly` Function that takes radius as a float in meters (no astropy units) and returns the column density in 1/m^2. Not relaible beyond ``max_grid_radius``, and questionable for radii less than ``collision_sphere_radius``. collision_sphere_radius : `astropy.units.Quantity` The radius of the collision sphere as described and calculated in Festou (1981): the radius beyond which a molecule can expect to see one collision over its lifetime. max_grid_radius : `astropy.units.Quantity` Cutoff radius for the model calculations. Attempting to take results from the model beyond this radius will be wrong or unreliable at best. coma_radius : `astropy.units.Quantity` Cutoff radius beyond which we take the parents to be entirely dissociated. num_fragments_theory : `float` Total number of fragments that we expect based on a steady-state calculation. Unreliable when time dependence is strong. num_fragments_grid : `float` Total number of fragments based on the actual results of the model. Agreement with ``num_fragments_theory`` is generally very good in the steady-state case. t_perm_flow : `float` Time for the comet to reach a steady state/permanent flow regime, where the number of fragments produced is equal to the number of fragments lost to dissociation. In a time-dependent production context, this will also measure how long the effects of a single outburst can affect the density. """ volume_density_grid: np.ndarray = None volume_density: np.ndarray = None column_density_grid: np.ndarray = None column_density: np.ndarray = None fragment_sputter: VMFragmentSputterPolar = None solid_angle_sputter: VMFragmentSputterPolar = None volume_density_interpolation: PPoly = None column_density_interpolation: PPoly = None collision_sphere_radius: u.Quantity = None max_grid_radius: u.Quantity = None coma_radius: u.Quantity = None num_fragments_theory: float = None num_fragments_grid: float = None t_perm_flow: u.Quantity = None
[docs] class VectorialModel(GasComa): """ Vectorial model for fragments in a coma produced with a dissociative energy kick. Parameters ---------- base_q : `~astropy.units.Quantity` Base production rate, per time parent: `~sbpy.data.Phys` Object with the following physical property fields: * ``tau_T``: total lifetime of the parent molecule * ``tau_d``: photodissociative lifetime of the parent molecule * ``v_outflow``: outflow velocity of the parent molecule * ``sigma``: cross-sectional area of the parent molecule fragment: `~sbpy.data.Phys` Object with the following physical property fields: * ``tau_T``: total lifetime of the fragment molecule * ``v_photo``: velocity of fragment resulting from photodissociation of the parent q_t: callable, optional Calculates the parent production rate as a function of time: ``q_t(t)``. The argument ``t`` is the look-back time as a float in units of seconds. The return value is the production rate in units of inverse seconds. If provided, this value is added to ``base_q``. If no time-dependence function is given, the model will run with steady production at ``base_q`` stretching infinitely far into the past. radial_points: int, optional Number of radial grid points the model will use radial_substeps: int, optional Number of points along the outflow axis to integrate over angular_points: int, optional Number of angular grid points the model will use parent_destruction_level: float, optional Model will attempt to track parents until this percentage has dissociated fragment_destruction_level: float, optional Model will attempt to track fragments until this percentage has dissociated max_fragment_lifetimes: float, optional Fragments traveling through the coma will be ignored if they take longer than this to arrive and contribute to the density at any considered point. print_progress: bool, optional Print progress while calculating. References ---------- The density distribution of neutral compounds in cometary atmospheres. I - Models and equations, Festou, M. C. 1981, Astronomy and Astrophysics, vol. 95, no. 1, Feb. 1981, p. 69-79. """ @bib.cite({"model": "1981A&A....95...69F"}) @u.quantity_input(base_q=(u.s**-1, u.mol / u.s)) def __init__( self, base_q, parent, fragment, q_t=None, radial_points=50, radial_substeps=80, angular_points=30, parent_destruction_level=0.99, fragment_destruction_level=0.95, max_fragment_lifetimes=8.0, print_progress=False, ): super().__init__(base_q, parent["v_outflow"][0]) # Calculations are done internally in meters and seconds to match the # base GasComa class # Convert to unitless value of production per second self.base_q = base_q.to(1 / u.s).value # Copy time dependence or create a steady production function if q_t is None: self.q_t = self._make_steady_production() else: self.q_t = q_t # Copy parent info, stripping astropy units and converting to meters # and seconds self.parent = VMParent( tau_T=parent["tau_T"][0].to(u.s).value, tau_d=parent["tau_d"][0].to(u.s).value, v_outflow=parent["v_outflow"][0].to(u.m / u.s).value, sigma=parent["sigma"][0].to(u.m**2).value, ) # Same for the fragment info self.fragment = VMFragment( tau_T=fragment["tau_T"][0].to(u.s).value, v_photo=fragment["v_photo"][0].to(u.m / u.s).value, ) # Grid settings self.grid = VMGridParams( radial_points=radial_points, radial_substeps=radial_substeps, angular_points=angular_points, ) self.model_params = VMParams( # Helps define cutoff for radial grid at this percentage of parents # lost to decay parent_destruction_level=parent_destruction_level, # Default here is lower than parents because they are born farther # from nucleus, tracking them too long will stretch the radial grid # a bit too much fragment_destruction_level=fragment_destruction_level, # If a fragment has to travel longer than this many lifetimes to # contribute to the density at a point, ignore it max_fragment_lifetimes=max_fragment_lifetimes, ) self.print_progress = print_progress self.vmr = VMResult() # Calculate up a few things self._setup_calculations() # Build the radial grid self.fast_voldens_grid = self._make_radial_logspace_grid() self.vmr.volume_density_grid = self.fast_voldens_grid * (u.m) # Angular grid self.d_alpha = self.epsilon_max / self.grid.angular_points # Make array of angles adjusted up away from zero, to keep from # encroaching on the outflow axis self.angular_grid = np.linspace( 0, self.epsilon_max, num=self.grid.angular_points, endpoint=False ) self.angular_grid += self.d_alpha / 2 # Makes a 2d array full of zero values self.fragment_sputter = np.zeros( (self.grid.radial_points, self.grid.angular_points) ) # Do the main computation self._compute_fragment_density() # Turn our volume density into a function of r self._interpolate_volume_density() # Get the column density at our grid points and interpolate self._compute_column_density() # Count up the number of fragments in the grid versus theoretical value self.vmr.num_fragments_theory = self._calc_num_fragments_theory() self.vmr.num_fragments_grid = self._calc_num_fragments_grid() # Convert fragment sputters to group of one-dimensional arrays # of the form r_i, theta_i, fragment_density_i sputterlist = [] for (i, j), frag_dens in np.ndenumerate(self.fragment_sputter): sputterlist.append( [ self.vmr.volume_density_grid[i].to(u.m).value, self.angular_grid[j], frag_dens, ] ) sputter = np.array(sputterlist) # fill in the fragment sputter results self.vmr.fragment_sputter = VMFragmentSputterPolar( rs=sputter[:, 0] * u.m, thetas=sputter[:, 1], fragment_density=sputter[:, 2] / u.m**3, ) normsputterlist = [] for (i, j), norm_dens in np.ndenumerate(self.solid_angle_sputter): normsputterlist.append( [ self.vmr.volume_density_grid[i].to(u.m).value, self.angular_grid[j], norm_dens, ] ) normsputter = np.array(normsputterlist) self.vmr.solid_angle_sputter = VMFragmentSputterPolar( rs=normsputter[:, 0] * u.m, thetas=normsputter[:, 1], fragment_density=normsputter[:, 2] / u.m**3, ) if print_progress: print("Vectorial model calculations complete!")
[docs] @classmethod def binned_production(cls, qs, fragment, parent, ts, **kwargs): """Alternate constructor for vectorial model Parameters ---------- qs : `~astropy.units.Quantity` List of steady production rates, per time, with length equal to that of ``ts``. parent: `~sbpy.data.Phys` Same as __init__ fragment: `~sbpy.data.Phys` Same as __init__ ts : `~astropy.units.Quantity` List of times corresponding to when the production qs begin, with positive times indicating the past. kwargs: variable, optional Any additional parameters in kwargs are passed on to __init__, which are documented above and may be passed in here. Returns ------- VectorialModel Instance of the VectorialModel class Examples -------- This specifies that from 30 days ago to 7 days ago, the production was 1.e27, changes to 3.e27 between 7 and 5 days ago, then falls to 2.e27 from 5 days ago until now: >>> q_example = [1.e27, 3.e27, 1.e27] * (1/u.s) >>> t_example = [30, 7, 5] * u.day Notes ----- Preserves Festou's original fortran method of describing time dependence in the model - time bins of steady production at specified intervals. The base production of the model is taken from the first element in the production array, which assumes the arrays are time-ordered from oldest to most recent. The base production extends backward in time to infinity, so take care when using this method for time dependence if that is not what is intended. """ return cls( base_q=qs[0], q_t=VectorialModel._make_binned_production(qs, ts), fragment=fragment, parent=parent, **kwargs )
def _make_binned_production(qs, ts) -> Callable[[np.float64], np.float64]: """Produces a time dependence function out of lists given to binned_production constructor. Parameters ---------- qs : `astropy.units.Quantity` See binned_production for description ts : `astropy.units.Quantity` See binned_production for description Returns ------- q_t : function See __init__ for description Notes ----- We create a model-compatible function for time dependence out of the information specified in the arrays qs and ts The resulting function gives a steady specified production within the given time windows. """ base_q = qs[0] q_variations = qs - base_q q_invsecs = q_variations.to_value(1 / (u.s)) t_at_p_secs = ts.to_value(u.s) # extend the arrays to simplify our comparisons in binned_q_t # last bin stops at observation time, t = 0 t_at_p_secs = np.append(t_at_p_secs, [0]) # junk value for production because this is in the future q_invsecs = np.append(q_invsecs, [0]) # this function represents the deviation from base_q at any time t, so # we need to handle times outside of the range specified in 'ts' def q_t(t): # too long in the past? if t > t_at_p_secs[0]: return 0 # in the future? if t < 0: return 0 # find which bin the given time falls in, and return the # corresponding production for that interval for i in range(len(q_invsecs) - 1): if t < t_at_p_secs[i] and t > t_at_p_secs[i + 1]: return q_invsecs[i] return q_t def _make_steady_production(self) -> Callable[[np.float64], np.float64]: """Produces a time dependence function that contributes no extra parents at any time. Returns ------- q_t : function See __init__ for description Notes ----- If no q_t is given, we use this as our time dependence as the model needs a q_t to run """ # No additional production at any time def q_t(_): return 0.0 return q_t def _setup_calculations(self) -> None: """Miscellaneous calculations to inform the model later. Notes ----- Calculates the collision sphere radius, coma radius, time to permanent flow regime, the maximum radius our grid could possibly need to extend out to, and the maximum angle that a fragment's trajectory can deviate from its parent's trajectory (which is assumed to be radial). """ """ Calculate collision sphere radius based on the first production rate, Eq. (5) in Festou 1981 Note that this is only calculated with the base production rate, because it is assumed that the base production rate has had roughly enough time to reach a steady state before letting production vary with time. """ if self.print_progress: print("Performing setup calculations...") # This v_therm factor comes from molecular flux of ideal gas moving # through a surface, in our case the surface of the collision sphere # The gas is collisional inside this sphere and free beyond it, so we # can model this as effusion v_therm = self.parent.v_outflow * 0.25 q = self.base_q vp = self.parent.v_outflow vf = self.fragment.v_photo # Eq. 5 of Festou 1981 self.vmr.collision_sphere_radius = ( (self.parent.sigma * q * v_therm) / (vp * vp) ) * u.m # Calculates the radius of the coma (defined as the space the parents # occupy) # NOTE: Equation (16) of Festou 1981 where alpha is the percent # destruction of molecules parent_beta_r = - \ np.log(1.0 - self.model_params.parent_destruction_level) parent_r = parent_beta_r * vp * self.parent.tau_T self.vmr.coma_radius = parent_r * u.m fragment_beta_r = - \ np.log(1.0 - self.model_params.fragment_destruction_level) # Calculate the time needed to hit a steady, permanent production of # fragments perm_flow_radius = self.vmr.coma_radius.value + ( (vp + vf) * fragment_beta_r * self.fragment.tau_T ) t_secs = self.vmr.coma_radius.value / vp + ( perm_flow_radius - self.vmr.coma_radius.value ) / (vp + vf) self.vmr.t_perm_flow = (t_secs * u.s).to(u.day) # This is the total radial size that parents & fragments occupy, beyond # which we assume zero density self.vmr.max_grid_radius = perm_flow_radius * u.m # Two cases for angular range of ejection of fragment based on relative # velocities of parent and fragment species if vf < vp: self.epsilon_max = np.arcsin(vf / vp) else: self.epsilon_max = np.pi
[docs] def production_at_time(self, t: np.float64) -> np.float64: """Get production rate at time t. Parameters ---------- t : numpy.float64 Time in seconds, with positive values representing the past Returns ------- numpy.float64 Production rate, unitless, at the specified time """ return self.base_q + self.q_t(t)
def _make_radial_logspace_grid(self) -> np.ndarray: """Create an appropriate radial grid based on the model parameters. Returns ------- ndarray Logarithmically spaced samples of the radial space around the coma, out to a maximum distance. Notes ----- Creates a grid (in meters) with numpy's logspace function that covers the expected radial size, stretching from 2 times the collision sphere radius (near the nucleus be dragons) out to the calculated max. If we get too close to the nucleus things go very badly so don't do it, dear reader. """ start_power = np.log10(self.vmr.collision_sphere_radius.value * 2) end_power = np.log10(self.vmr.max_grid_radius.value) return np.logspace( start_power, end_power, num=self.grid.radial_points, endpoint=True ) def _outflow_axis_sampling( self, x: np.float64, y: np.float64, theta: np.float64 ) -> Tuple[np.ndarray, np.ndarray]: """Construct a list of points along the outflow axis, sampled to be more dense around the minimum distance to (x, y). Parameters ---------- x : numpy.float64 x-coordinate of the space where we are calculating the fragment density y : numpy.float64 y-coordinate of the space where we are calculating the fragment density theta: numpy.float64 polar spherical coordinate of (x, y) to avoid calculating it from (x, y) Returns ------- tuple(numpy.ndarray, numpy.ndarray) Tuple with list of ejection sites in the first index and their extents in the second index Notes ----- Returns list of radial points along outflow axis to sample for fragment ejection, with a list describing the size dr that each ejection point occupies along the outflow axis. """ outflow_axis_edge = self.vmr.coma_radius.value # calc angle to edge of the grid, and use it if it's less than # epsilon max. without this, when epsilon_max approaches pi we # get values of r outside the edge of the grid grid_edge_angle = np.arctan2(x, (y - outflow_axis_edge)) if grid_edge_angle < self.epsilon_max: max_subangle = grid_edge_angle else: max_subangle = self.epsilon_max # Trace angles from our grid point at x, y and find where dissociations # would have to occur along the outflow axis. This samples shorter # trajectories more finely along the axis, which is necessary because # the exponential decay along their transit means they contribute the # most to the fragment density ces = np.vectorize(lambda epsilon: y - x / np.tan(epsilon)) # array one element larger than radial_substeps because we will # use the midpoints defined by this array subangles = np.linspace( theta, max_subangle, num=self.grid.radial_substeps + 1, endpoint=True ) # Find where this set of angles intersects the outflow axis ejection_grid = ces(subangles) # Find dr for radial integration drs = np.diff(ejection_grid) # This puts the sampling of the ejection sites at the midpoints # of our ejection grid ejection_sites = (ejection_grid[1:] + ejection_grid[:-1]) / 2 return ejection_sites, drs def _fragment_sputter(self, r: np.float64, theta: np.float64) -> np.float64: """Compute the fragment density at (r, theta) in a spherical coordinate system where theta is the polar angle and the parents flow radially outward along the z axis. Parameters ---------- r : numpy.float64 radial coordinate in meters where we are calculating the fragment density. theta : numpy.float64 polar spherical coordinate where we are calculating the fragment density. Returns ------- np.float64 Fragment density in 1/m**3, no astropy units attached. """ sputter = 0.0 vp = self.parent.v_outflow vf = self.fragment.v_photo p_tau_T = self.parent.tau_T f_tau_T = self.fragment.tau_T x = r * np.sin(theta) y = r * np.cos(theta) ejection_sites, drs = self._outflow_axis_sampling(x, y, theta) # Loop over these ejection sites that contribute to x, y for slice_r, dr in zip(ejection_sites, drs): # Distance from dissociation site to our grid point sep_dist = np.sqrt(x**2 + (slice_r - y) ** 2) cos_eject = (y - slice_r) / sep_dist sin_eject = x / sep_dist # Parent extinction when traveling along to the # dissociation site p_extinction = np.e ** (-slice_r / (p_tau_T * vp)) # Calculate sqrt(vR^2 - u^2 sin^2 gamma) v_factor = np.sqrt(vf * vf - (vp * vp) * sin_eject**2) # The geometry of the problem can admit one or two solutions for # the velocity of the fragment v_one = vp * cos_eject + v_factor if vf > vp: velocities = [v_one] else: v_two = vp * cos_eject - v_factor velocities = [v_one, v_two] # TODO: this shouldn't be necessary if we only pick r, theta inside ejection if v_one == 0.0: continue for v in velocities: # Time taken to travel from the dissociation point at v, reject # if the time is too large and fragments have decayed beyond # self.fragment_destruction_level percent t_frag = sep_dist / v if t_frag > self.time_limit: continue # total time between parent emission from nucleus and fragment # arriving at our point of interest, which we then use to look # up Q at that time in the past t_total = (slice_r / vp) + t_frag # Division by parent velocity makes this production per # unit distance for radial integration q(r, epsilon) # given by eq. 32, Festou 1981 q = self.production_at_time(t_total) / vp q_r_eps = (v**2 * q) / (vf * np.abs(v - vp * cos_eject)) # Fragment extinction when traveling at speed v from # dissociation site to x, y f_extinction = np.e ** (-t_frag / f_tau_T) # differential addition to the density integrating along dr, # similar to eq. (36) Festou 1981 n_r = (p_extinction * f_extinction * q_r_eps) / (sep_dist**2 * v) sputter += n_r * dr return sputter def _compute_fragment_density(self) -> None: """Computes the density of fragments as a function of radius. Notes ----- Computes the density of fragments sputtered around the comet due to an outflow of parents along the z-axis, performing the integration in eq. (36), Festou 1981. The resulting radial fragment density will be in units of 1/(m^3) as we work in m, s, and m/s. We then interpolate the fragment density as a function of arbitrary radius. We use our results to calculate the total number of fragments in the coma for comparison to the theoretical number we expect, to provide the user with a rough idea of how well the chosen radial and angular grid sizes have captured the appropriate amount of particles. Note that some level of disagreement is expected because the parent_destruction_level and fragment_destruction_level parameters cut the grid off before all parents can dissociate, and thus some escape the model and come up missing in the fragment count based on the grid. """ if self.print_progress: print("Starting fragment density computations...") # Follow fragments until they have been totally destroyed self.time_limit = self.model_params.max_fragment_lifetimes * self.fragment.tau_T # More factors to fill out integral similar to eq. (36) Festou 1981 integration_factor = ( (1 / (4 * np.pi * self.parent.tau_d)) * self.d_alpha / (4.0 * np.pi) ) # vectorize and apply to each combination of (r, theta) sputter_func = np.vectorize(self._fragment_sputter) rs, thetas = np.meshgrid( self.fast_voldens_grid, self.angular_grid, indexing="ij" ) self.fragment_sputter = integration_factor * sputter_func(rs, thetas) # Make array to hold our data, no units self.fast_voldens = np.zeros(self.grid.radial_points) # Integration factors from angular part of integral, similar to # eq. (36) Festou 1981 # integrate over theta to produce radial fragment volume density # axis = 1 sums over second column, theta # Equivalent to summing over j for sin(theta[j]) * # fragment_sputter[i][j] with numpy magic self.solid_angle_sputter = np.sin(thetas) * self.fragment_sputter self.fast_voldens = 2.0 * np.pi * \ np.sum(self.solid_angle_sputter, axis=1) # Tag with proper units self.vmr.volume_density = self.fast_voldens / (u.m**3) @requires("scipy") def _interpolate_volume_density(self) -> None: """Interpolate the volume density as a function of radial distance from the nucleus. Takes our fragment density grid and constructs density as a function of arbitrary radius. """ if self.print_progress: print("Interpolating radial fragment density...") # Interpolate this radial density grid with a cubic spline for lookup # at non-grid radii, input in m, output in 1/m^3 self.vmr.volume_density_interpolation = CubicSpline( self.fast_voldens_grid, self.fast_voldens, bc_type="natural" ) def _column_density_at_rho(self, rho: np.float64) -> np.float64: """ Calculate the column density of fragments at a given impact parameter. Parameters ---------- rho : np.float64 Impact parameter of the column density integration, in meters Returns ------- np.float64 Column density at the given impact parameter in m^-2, no astropy units attached Notes ----- We return zero column density beyond a certain distance past the grid edge, which can lead to strange graphing results and sharp cutoffs. """ reach_factor = 2.0 # _volume_density() approximates beyond the edge of the grid, so allow # an arbirtary limit beyond which we just return zero column density r_max = self.vmr.max_grid_radius.value * reach_factor if rho > r_max: return 0 rhosq = rho**2 z_max = np.sqrt(r_max**2 - rhosq) def column_density_integrand(z): return self._volume_density(np.sqrt(z**2 + rhosq)) c_dens = 2 * romberg(column_density_integrand, 0, z_max, rtol=0.0001, divmax=50) # result is in 1/m^2 return c_dens @requires("scipy") def _compute_column_density(self) -> None: """Compute the column density on the grid and interpolate the results. Notes ----- The interpolator returns column density in m^-2, no astropy units attached. """ if self.print_progress: print("Computing column densities...") # make a grid for the column density values column_density_grid = self._make_radial_logspace_grid() # vectorize so we can hand the whole grid to our function cd_vectorized = np.vectorize(self._column_density_at_rho) # array now holds corresponding column density values column_densities = cd_vectorized(column_density_grid) self.fast_column_density_grid = column_density_grid self.vmr.column_density_grid = column_density_grid * u.m self.vmr.column_density = column_densities / (u.m**2) # Interpolation gives column density in m^-2 self.vmr.column_density_interpolation = CubicSpline( column_density_grid, column_densities, bc_type="natural" ) def _calc_num_fragments_theory(self) -> np.float64: """The total number of fragment species we expect in the coma. Returns ------- np.float64 Total number of fragment species we expect in the coma theoretically Notes ----- Outbursts/time dependent production in general will make this result poor due to the grid being sized to capture a certain fraction of parents/fragments at the oldest (first) production rate. The farther you get from this base production, the farther the model will deviate from capturing the requested percentage of particles. """ vp = self.parent.v_outflow vf = self.fragment.v_photo p_tau_T = self.parent.tau_T f_tau_T = self.fragment.tau_T p_tau_d = self.parent.tau_d t_perm = self.vmr.t_perm_flow.to(u.s).value alpha = f_tau_T * p_tau_T / p_tau_d max_r = self.vmr.max_grid_radius.value edge_adjust = np.pi * max_r * max_r * (vf + vp) * self.fast_voldens[-1] num_time_slices = 1000 # array starting at tperm (farthest back in time we expect something to # hang around), stepped down to zero (when the observation takes place) time_slices = np.linspace(t_perm, 0, num_time_slices, endpoint=True) # estimate based on discrete-time source and sink model theory_total = 0 for i, t in enumerate(time_slices[:-1]): extinction_one = t / p_tau_T extinction_two = time_slices[i + 1] / p_tau_T mult_factor = -np.e ** (-extinction_one) + \ np.e ** (-extinction_two) theory_total += self.production_at_time(t) * mult_factor return theory_total * alpha - edge_adjust def _calc_num_fragments_grid(self) -> np.float64: """Total number of fragments in the coma. Calculates the total number of fragments by integrating the density grid over its volume Returns ------- np.float64 Number of fragments in the coma based on our grid calculations Notes ----- Outbursts/time dependent production in general will make this result poor due to the grid being sized to capture a certain fraction of parents/fragments at the oldest (first) production rate. The farther you get from this base production, the farther the model will deviate from capturing the requested percentage of particles. """ max_r = self.vmr.max_grid_radius.value def vol_integrand(r, r_func): return r_func(r) * r**2 r_int = romberg( vol_integrand, 0, max_r, args=(self.vmr.volume_density_interpolation,), rtol=0.0001, divmax=20, ) return 4 * np.pi * r_int def _column_density(self, rho) -> np.float64: """Gives fragment column density at arbitrary impact parameter. Parameters ---------- rho : np.float64 Impact parameter, in meters, no astropy units attached. Returns ------- np.float64 Fragment column density at given impact parameter, in m^-2, no astropy units. """ if rho < self.fast_column_density_grid[0]: return self.vmr.column_density[0].value if rho > self.vmr.max_grid_radius.value: return 0 return self.vmr.column_density_interpolation(rho) def _volume_density(self, r) -> np.float64: """Gives fragment volume density at arbitrary radius. Parameters ---------- r : np.float64 Distance from nucleus, in meters, no astropy units attached. Returns ------- np.float64 Fragment volume density at specified radius, in m^-3, no astropy units. Notes ----- When asked for a value at a radius smaller than the first grid point, we have two choices: return a value based on the interpolation of the volume density, or return the value of the closest grid point to the nucleus. This function returns the closest grid point to the nucleus, because the model does not say anything about what happens inside the collision sphere - so we approximate by clamping the value as constant. We can't trust the interpolation outside of the grid because the density values can get very large or become negative. Outside the radius of the grid, the volume density is approximated with exponential decay based on the total lifetime of the fragment. """ if r < self.fast_voldens_grid[0]: return self.fast_voldens[0] if r > self.vmr.max_grid_radius.value: diff = r - self.vmr.max_grid_radius.value guess = self.fast_voldens[-1] * np.exp( -diff / (self.fragment.tau_T * self.fragment.v_photo) ) return guess return self.vmr.volume_density_interpolation(r)