Source code for autogalaxy.profiles.mass_profiles.stellar_mass_profiles

import copy
import numpy as np
from scipy.special import wofz
from scipy.integrate import quad
from typing import List, Tuple

import autoarray as aa

from autogalaxy.profiles.mass_profiles import MassProfile
from autogalaxy.profiles.mass_profiles.mass_profiles import (
    MassProfileMGE,
    MassProfileCSE,
)

from autogalaxy.profiles.mass_profiles.mass_profiles import psi_from


class StellarProfile:

    pass


[docs]class EllGaussian(MassProfile, StellarProfile):
[docs] def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), elliptical_comps: Tuple[float, float] = (0.0, 0.0), intensity: float = 0.1, sigma: float = 1.0, mass_to_light_ratio: float = 1.0, ): """ The elliptical Gaussian light profile. Parameters ---------- centre The (y,x) arc-second coordinates of the profile centre. elliptical_comps The first and second ellipticity components of the elliptical coordinate system, (see the module `autogalaxy -> convert.py` for the convention). intensity Overall intensity normalisation of the light profile (units are dimensionless and derived from the data the light profile's image is compared too, which is expected to be electrons per second). sigma The sigma value of the Gaussian. """ super(EllGaussian, self).__init__( centre=centre, elliptical_comps=elliptical_comps ) super(MassProfile, self).__init__( centre=centre, elliptical_comps=elliptical_comps ) self.mass_to_light_ratio = mass_to_light_ratio self.intensity = intensity self.sigma = sigma
def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike): """ Calculate the deflection angles at a given set of arc-second gridded coordinates. Parameters ---------- grid The grid of (y,x) arc-second coordinates the deflection angles are computed on. """ return self.deflections_2d_via_analytic_from(grid=grid) @aa.grid_dec.grid_2d_to_structure @aa.grid_dec.transform @aa.grid_dec.relocate_to_radial_minimum def deflections_2d_via_analytic_from(self, grid: aa.type.Grid2DLike): """ Calculate the deflection angles at a given set of arc-second gridded coordinates. Parameters ---------- grid The grid of (y,x) arc-second coordinates the deflection angles are computed on. """ deflections = ( self.mass_to_light_ratio * self.intensity * self.sigma * np.sqrt((2 * np.pi) / (1.0 - self.axis_ratio ** 2.0)) * self.zeta_from(grid=grid) ) return self.rotate_grid_from_reference_frame( np.multiply( 1.0, np.vstack((-1.0 * np.imag(deflections), np.real(deflections))).T ) ) @aa.grid_dec.grid_2d_to_structure @aa.grid_dec.transform @aa.grid_dec.relocate_to_radial_minimum def deflections_2d_via_integral_from(self, grid: aa.type.Grid2DLike): """ Calculate the deflection angles at a given set of arc-second gridded coordinates. Parameters ---------- grid The grid of (y,x) arc-second coordinates the deflection angles are computed on. Note: sigma is divided by sqrt(q) here. """ def calculate_deflection_component(npow, index): deflection_grid = self.axis_ratio * grid[:, index] for i in range(grid.shape[0]): deflection_grid[i] *= ( self.intensity * self.mass_to_light_ratio * quad( self.deflection_func, a=0.0, b=1.0, args=( grid[i, 0], grid[i, 1], npow, self.axis_ratio, self.sigma / np.sqrt(self.axis_ratio), ), )[0] ) return deflection_grid deflection_y = calculate_deflection_component(1.0, 0) deflection_x = calculate_deflection_component(0.0, 1) return self.rotate_grid_from_reference_frame( np.multiply(1.0, np.vstack((deflection_y, deflection_x)).T) ) @staticmethod def deflection_func(u, y, x, npow, axis_ratio, sigma): eta_u = np.sqrt(axis_ratio) * np.sqrt( (u * ((x ** 2) + (y ** 2 / (1 - (1 - axis_ratio ** 2) * u)))) ) return np.exp(-0.5 * np.square(np.divide(eta_u, sigma))) / ( (1 - (1 - axis_ratio ** 2) * u) ** (npow + 0.5) ) @aa.grid_dec.grid_2d_to_structure @aa.grid_dec.transform @aa.grid_dec.relocate_to_radial_minimum def convergence_2d_from(self, grid: aa.type.Grid2DLike): """Calculate the projected convergence at a given set of arc-second gridded coordinates. Parameters ---------- grid The grid of (y,x) arc-second coordinates the convergence is computed on. """ return self.convergence_func(self.grid_to_eccentric_radii(grid)) def convergence_func(self, grid_radius: float) -> float: return self.mass_to_light_ratio * self.image_2d_via_radii_from(grid_radius) @aa.grid_dec.grid_2d_to_structure def potential_2d_from(self, grid: aa.type.Grid2DLike): return np.zeros(shape=grid.shape[0]) def image_2d_via_radii_from(self, grid_radii: np.ndarray): """Calculate the intensity of the Gaussian light profile on a grid of radial coordinates. Parameters ---------- grid_radii The radial distance from the centre of the profile. for each coordinate on the grid. Note: sigma is divided by sqrt(q) here. """ return np.multiply( self.intensity, np.exp( -0.5 * np.square( np.divide(grid_radii, self.sigma / np.sqrt(self.axis_ratio)) ) ), ) @property def axis_ratio(self): axis_ratio = super().axis_ratio return axis_ratio if axis_ratio < 0.9999 else 0.9999 def zeta_from(self, grid: aa.type.Grid2DLike): q2 = self.axis_ratio ** 2.0 ind_pos_y = grid[:, 0] >= 0 shape_grid = np.shape(grid) output_grid = np.zeros((shape_grid[0]), dtype=np.complex128) scale_factor = self.axis_ratio / (self.sigma * np.sqrt(2.0 * (1.0 - q2))) xs_0 = grid[:, 1][ind_pos_y] * scale_factor ys_0 = grid[:, 0][ind_pos_y] * scale_factor xs_1 = grid[:, 1][~ind_pos_y] * scale_factor ys_1 = -grid[:, 0][~ind_pos_y] * scale_factor output_grid[ind_pos_y] = -1j * ( wofz(xs_0 + 1j * ys_0) - np.exp(-(xs_0 ** 2.0) * (1.0 - q2) - ys_0 * ys_0 * (1.0 / q2 - 1.0)) * wofz(self.axis_ratio * xs_0 + 1j * ys_0 / self.axis_ratio) ) output_grid[~ind_pos_y] = np.conj( -1j * ( wofz(xs_1 + 1j * ys_1) - np.exp(-(xs_1 ** 2.0) * (1.0 - q2) - ys_1 * ys_1 * (1.0 / q2 - 1.0)) * wofz(self.axis_ratio * xs_1 + 1j * ys_1 / self.axis_ratio) ) ) return output_grid def with_new_normalization(self, normalization): mass_profile = copy.copy(self) mass_profile.mass_to_light_ratio = normalization return mass_profile
# noinspection PyAbstractClass class AbstractEllSersic(MassProfile, MassProfileMGE, MassProfileCSE, StellarProfile): def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), elliptical_comps: Tuple[float, float] = (0.0, 0.0), intensity: float = 0.1, effective_radius: float = 0.6, sersic_index: float = 0.6, mass_to_light_ratio: float = 1.0, ): """ The Sersic mass profile, the mass profiles of the light profiles that are used to fit and subtract the lens \ model_galaxy's light. Parameters ---------- centre The (y,x) arc-second coordinates of the profile centre. elliptical_comps The first and second ellipticity components of the elliptical coordinate system, (see the module `autogalaxy -> convert.py` for the convention). intensity Overall flux intensity normalisation in the light profiles (electrons per second). effective_radius The radius containing half the light of this profile. sersic_index Controls the concentration of the profile (lower -> less concentrated, higher -> more concentrated). mass_to_light_ratio The mass-to-light ratio of the light profiles """ super(AbstractEllSersic, self).__init__( centre=centre, elliptical_comps=elliptical_comps ) super(MassProfile, self).__init__( centre=centre, elliptical_comps=elliptical_comps ) super(MassProfileMGE, self).__init__() super(MassProfileCSE, self).__init__() self.mass_to_light_ratio = mass_to_light_ratio self.intensity = intensity self.effective_radius = effective_radius self.sersic_index = sersic_index def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike): return self.deflections_2d_via_cse_from(grid=grid) @aa.grid_dec.grid_2d_to_structure @aa.grid_dec.transform @aa.grid_dec.relocate_to_radial_minimum def deflections_2d_via_mge_from(self, grid: aa.type.Grid2DLike): """ Calculate the projected 2D deflection angles from a grid of (y,x) arc second coordinates, by computing and summing the convergence of each individual cse used to decompose the mass profile. The cored steep elliptical (cse) decomposition of a the elliptical NFW mass profile (e.g. `decompose_convergence_via_cse`) is using equation (12) of Oguri 2021 (https://arxiv.org/abs/2106.11464). Parameters ---------- grid The grid of (y,x) arc-second coordinates the convergence is computed on. """ return self._deflections_2d_via_mge_from( grid=grid, sigmas_factor=np.sqrt(self.axis_ratio) ) @aa.grid_dec.grid_2d_to_structure @aa.grid_dec.transform @aa.grid_dec.relocate_to_radial_minimum def deflections_2d_via_cse_from(self, grid: aa.type.Grid2DLike): """ Calculate the projected 2D deflection angles from a grid of (y,x) arc second coordinates, by computing and summing the convergence of each individual cse used to decompose the mass profile. The cored steep elliptical (cse) decomposition of a the elliptical NFW mass profile (e.g. `decompose_convergence_via_cse`) is using equation (12) of Oguri 2021 (https://arxiv.org/abs/2106.11464). Parameters ---------- grid The grid of (y,x) arc-second coordinates the convergence is computed on. """ return self._deflections_2d_via_cse_from(grid=grid) @aa.grid_dec.grid_2d_to_structure @aa.grid_dec.transform @aa.grid_dec.relocate_to_radial_minimum def convergence_2d_from(self, grid: aa.type.Grid2DLike): """Calculate the projected convergence at a given set of arc-second gridded coordinates. Parameters ---------- grid The grid of (y,x) arc-second coordinates the convergence is computed on. """ return self.convergence_func(self.grid_to_eccentric_radii(grid)) @aa.grid_dec.grid_2d_to_structure @aa.grid_dec.transform @aa.grid_dec.relocate_to_radial_minimum def convergence_2d_via_mge_from(self, grid: aa.type.Grid2DLike): """ Calculate the projected convergence at a given set of arc-second gridded coordinates. Parameters ---------- grid The grid of (y,x) arc-second coordinates the convergence is computed on. """ eccentric_radii = self.grid_to_eccentric_radii(grid=grid) return self._convergence_2d_via_mge_from(grid_radii=eccentric_radii) @aa.grid_dec.grid_2d_to_structure @aa.grid_dec.transform @aa.grid_dec.relocate_to_radial_minimum def convergence_2d_via_cse_from(self, grid: aa.type.Grid2DLike): """ Calculate the projected 2D convergence from a grid of (y,x) arc second coordinates, by computing and summing the convergence of each individual cse used to decompose the mass profile. The cored steep elliptical (cse) decomposition of a the elliptical NFW mass profile (e.g. `decompose_convergence_via_cse`) is using equation (12) of Oguri 2021 (https://arxiv.org/abs/2106.11464). Parameters ---------- grid The grid of (y,x) arc-second coordinates the convergence is computed on. """ elliptical_radii = self.grid_to_elliptical_radii(grid=grid) return self._convergence_2d_via_cse_from(grid_radii=elliptical_radii) def convergence_func(self, grid_radius: float) -> float: return self.mass_to_light_ratio * self.image_2d_via_radii_from(grid_radius) @aa.grid_dec.grid_2d_to_structure def potential_2d_from(self, grid: aa.type.Grid2DLike): return np.zeros(shape=grid.shape[0]) def image_2d_via_radii_from(self, radius: np.ndarray): """ Returns the intensity of the profile at a given radius. Parameters ---------- radius The distance from the centre of the profile. """ return self.intensity * np.exp( -self.sersic_constant * (((radius / self.effective_radius) ** (1.0 / self.sersic_index)) - 1) ) def decompose_convergence_via_mge(self) -> Tuple[List, List]: radii_min = self.effective_radius / 100.0 radii_max = self.effective_radius * 20.0 def sersic_2d(r): return ( self.mass_to_light_ratio * self.intensity * np.exp( -self.sersic_constant * (((r / self.effective_radius) ** (1.0 / self.sersic_index)) - 1.0) ) ) return self._decompose_convergence_via_mge( func=sersic_2d, radii_min=radii_min, radii_max=radii_max ) def decompose_convergence_via_cse(self,) -> Tuple[List, List]: """ Decompose the convergence of the Sersic profile into cored steep elliptical (cse) profiles. This decomposition uses the standard 2d profile of a Sersic mass profile. Parameters ---------- func The function representing the profile that is decomposed into CSEs. radii_min: The minimum radius to fit radii_max: The maximum radius to fit total_cses The number of CSEs used to approximate the input func. sample_points: int (should be larger than 'total_cses') The number of data points to fit Returns ------- Tuple[List, List] A list of amplitudes and core radii of every cored steep elliptical (cse) the mass profile is decomposed into. """ upper_dex, lower_dex, total_cses, sample_points = cse_settings_from( effective_radius=self.effective_radius, sersic_index=self.sersic_index, sersic_constant=self.sersic_constant, mass_to_light_gradient=0.0, ) scaled_effective_radius = self.effective_radius / np.sqrt(self.axis_ratio) radii_min = scaled_effective_radius / 10.0 ** lower_dex radii_max = scaled_effective_radius * 10.0 ** upper_dex def sersic_2d(r): return ( self.mass_to_light_ratio * self.intensity * np.exp( -self.sersic_constant * ( ((r / scaled_effective_radius) ** (1.0 / self.sersic_index)) - 1.0 ) ) ) return self._decompose_convergence_via_cse_from( func=sersic_2d, radii_min=radii_min, radii_max=radii_max, total_cses=total_cses, sample_points=sample_points, ) @property def sersic_constant(self): """A parameter derived from Sersic index which ensures that effective radius contains 50% of the profile's total integrated light. """ return ( (2 * self.sersic_index) - (1.0 / 3.0) + (4.0 / (405.0 * self.sersic_index)) + (46.0 / (25515.0 * self.sersic_index ** 2)) + (131.0 / (1148175.0 * self.sersic_index ** 3)) - (2194697.0 / (30690717750.0 * self.sersic_index ** 4)) ) @property def ellipticity_rescale(self): return 1.0 - ((1.0 - self.axis_ratio) / 2.0) @property def elliptical_effective_radius(self): """ The effective_radius of a Sersic light profile is defined as the circular effective radius. This is the \ radius within which a circular aperture contains half the profiles's total integrated light. For elliptical \ systems, this won't robustly capture the light profile's elliptical shape. The elliptical effective radius instead describes the major-axis radius of the ellipse containing \ half the light, and may be more appropriate for highly flattened systems like disk galaxies. """ return self.effective_radius / np.sqrt(self.axis_ratio) def with_new_normalization(self, normalization): mass_profile = copy.copy(self) mass_profile.mass_to_light_ratio = normalization return mass_profile
[docs]class EllSersic(AbstractEllSersic, MassProfileMGE, MassProfileCSE): @aa.grid_dec.grid_2d_to_structure @aa.grid_dec.transform @aa.grid_dec.relocate_to_radial_minimum def deflections_2d_via_integral_from(self, grid: aa.type.Grid2DLike): """ Calculate the deflection angles at a given set of arc-second gridded coordinates. Parameters ---------- grid The grid of (y,x) arc-second coordinates the deflection angles are computed on. """ def calculate_deflection_component(npow, index): sersic_constant = self.sersic_constant deflection_grid = self.axis_ratio * grid[:, index] for i in range(grid.shape[0]): deflection_grid[i] *= ( self.intensity * self.mass_to_light_ratio * quad( self.deflection_func, a=0.0, b=1.0, args=( grid[i, 0], grid[i, 1], npow, self.axis_ratio, self.sersic_index, self.effective_radius, sersic_constant, ), )[0] ) return deflection_grid deflection_y = calculate_deflection_component(1.0, 0) deflection_x = calculate_deflection_component(0.0, 1) return self.rotate_grid_from_reference_frame( np.multiply(1.0, np.vstack((deflection_y, deflection_x)).T) ) @staticmethod def deflection_func( u, y, x, npow, axis_ratio, sersic_index, effective_radius, sersic_constant ): eta_u = np.sqrt(axis_ratio) * np.sqrt( (u * ((x ** 2) + (y ** 2 / (1 - (1 - axis_ratio ** 2) * u)))) ) return np.exp( -sersic_constant * (((eta_u / effective_radius) ** (1.0 / sersic_index)) - 1) ) / ((1 - (1 - axis_ratio ** 2) * u) ** (npow + 0.5))
[docs]class SphSersic(EllSersic):
[docs] def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), intensity: float = 0.1, effective_radius: float = 0.6, sersic_index: float = 0.6, mass_to_light_ratio: float = 1.0, ): """ The Sersic mass profile, the mass profiles of the light profiles that are used to fit and subtract the lens model_galaxy's light. Parameters ---------- centre The (y,x) arc-second coordinates of the profile centre intensity Overall flux intensity normalisation in the light profiles (electrons per second) effective_radius The circular radius containing half the light of this profile. sersic_index Controls the concentration of the profile (lower -> less concentrated, higher -> more concentrated). mass_to_light_ratio The mass-to-light ratio of the light profile. """ super().__init__( centre=centre, elliptical_comps=(0.0, 0.0), intensity=intensity, effective_radius=effective_radius, sersic_index=sersic_index, mass_to_light_ratio=mass_to_light_ratio, )
[docs]class EllExponential(EllSersic):
[docs] def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), elliptical_comps: Tuple[float, float] = (0.0, 0.0), intensity: float = 0.1, effective_radius: float = 0.6, mass_to_light_ratio: float = 1.0, ): """ The EllExponential mass profile, the mass profiles of the light profiles that are used to fit and subtract the lens model_galaxy's light. Parameters ---------- centre The (y,x) arc-second coordinates of the profile centre. elliptical_comps The first and second ellipticity components of the elliptical coordinate system, (see the module `autogalaxy -> convert.py` for the convention). intensity Overall flux intensity normalisation in the light profiles (electrons per second). effective_radius The circular radius containing half the light of this profile. mass_to_light_ratio The mass-to-light ratio of the light profiles """ super().__init__( centre=centre, elliptical_comps=elliptical_comps, intensity=intensity, effective_radius=effective_radius, sersic_index=1.0, mass_to_light_ratio=mass_to_light_ratio, )
[docs]class SphExponential(EllExponential):
[docs] def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), intensity: float = 0.1, effective_radius: float = 0.6, mass_to_light_ratio: float = 1.0, ): """ The Exponential mass profile, the mass profiles of the light profiles that are used to fit and subtract the lens model_galaxy's light. Parameters ---------- centre The (y,x) arc-second coordinates of the profile centre. intensity Overall flux intensity normalisation in the light profiles (electrons per second). effective_radius The circular radius containing half the light of this profile. mass_to_light_ratio The mass-to-light ratio of the light profiles. """ super().__init__( centre=centre, elliptical_comps=(0.0, 0.0), intensity=intensity, effective_radius=effective_radius, mass_to_light_ratio=mass_to_light_ratio, )
[docs]class EllDevVaucouleurs(EllSersic):
[docs] def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), elliptical_comps: Tuple[float, float] = (0.0, 0.0), intensity: float = 0.1, effective_radius: float = 0.6, mass_to_light_ratio: float = 1.0, ): """ The EllDevVaucouleurs mass profile, the mass profiles of the light profiles that are used to fit and subtract the lens model_galaxy's light. Parameters ---------- centre The (y,x) arc-second coordinates of the profile centre. elliptical_comps The first and second ellipticity components of the elliptical coordinate system, (see the module `autogalaxy -> convert.py` for the convention). intensity Overall flux intensity normalisation in the light profiles (electrons per second). effective_radius The radius containing half the light of this profile. mass_to_light_ratio The mass-to-light ratio of the light profile. """ super().__init__( centre=centre, elliptical_comps=elliptical_comps, intensity=intensity, effective_radius=effective_radius, sersic_index=4.0, mass_to_light_ratio=mass_to_light_ratio, )
[docs]class SphDevVaucouleurs(EllDevVaucouleurs):
[docs] def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), intensity: float = 0.1, effective_radius: float = 0.6, mass_to_light_ratio: float = 1.0, ): """ The DevVaucouleurs mass profile, the mass profiles of the light profiles that are used to fit and subtract the lens model_galaxy's light. Parameters ---------- centre The (y,x) arc-second coordinates of the profile centre. intensity Overall flux intensity normalisation in the light profiles (electrons per second). effective_radius The circular radius containing half the light of this profile. mass_to_light_ratio The mass-to-light ratio of the light profiles. """ super().__init__( centre=centre, elliptical_comps=(0.0, 0.0), intensity=intensity, effective_radius=effective_radius, mass_to_light_ratio=mass_to_light_ratio, )
[docs]class EllSersicRadialGradient(AbstractEllSersic):
[docs] def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), elliptical_comps: Tuple[float, float] = (0.0, 0.0), intensity: float = 0.1, effective_radius: float = 0.6, sersic_index: float = 0.6, mass_to_light_ratio: float = 1.0, mass_to_light_gradient: float = 0.0, ): """ Setup a Sersic mass and light profiles. Parameters ---------- centre The (y,x) arc-second coordinates of the profile centre. elliptical_comps The first and second ellipticity components of the elliptical coordinate system, (see the module `autogalaxy -> convert.py` for the convention). intensity Overall flux intensity normalisation in the light profiles (electrons per second). effective_radius The circular radius containing half the light of this profile. sersic_index Controls the concentration of the profile (lower -> less concentrated, higher -> more concentrated). mass_to_light_ratio The mass-to-light ratio of the light profile. mass_to_light_gradient The mass-to-light radial gradient. """ super().__init__( centre=centre, elliptical_comps=elliptical_comps, intensity=intensity, effective_radius=effective_radius, sersic_index=sersic_index, mass_to_light_ratio=mass_to_light_ratio, ) self.mass_to_light_gradient = mass_to_light_gradient
@aa.grid_dec.grid_2d_to_structure @aa.grid_dec.transform @aa.grid_dec.relocate_to_radial_minimum def deflections_2d_via_integral_from(self, grid: aa.type.Grid2DLike): """ Calculate the deflection angles at a given set of arc-second gridded coordinates. Parameters ---------- grid The grid of (y,x) arc-second coordinates the deflection angles are computed on. """ def calculate_deflection_component(npow, index): sersic_constant = self.sersic_constant deflection_grid = self.axis_ratio * grid[:, index] for i in range(grid.shape[0]): deflection_grid[i] *= ( self.intensity * self.mass_to_light_ratio * quad( self.deflection_func, a=0.0, b=1.0, args=( grid[i, 0], grid[i, 1], npow, self.axis_ratio, self.sersic_index, self.effective_radius, self.mass_to_light_gradient, sersic_constant, ), )[0] ) return deflection_grid deflection_y = calculate_deflection_component(1.0, 0) deflection_x = calculate_deflection_component(0.0, 1) return self.rotate_grid_from_reference_frame( np.multiply(1.0, np.vstack((deflection_y, deflection_x)).T) ) @staticmethod def deflection_func( u, y, x, npow, axis_ratio, sersic_index, effective_radius, mass_to_light_gradient, sersic_constant, ): eta_u = np.sqrt(axis_ratio) * np.sqrt( (u * ((x ** 2) + (y ** 2 / (1 - (1 - axis_ratio ** 2) * u)))) ) return ( (((axis_ratio * eta_u) / effective_radius) ** -mass_to_light_gradient) * np.exp( -sersic_constant * (((eta_u / effective_radius) ** (1.0 / sersic_index)) - 1) ) / ((1 - (1 - axis_ratio ** 2) * u) ** (npow + 0.5)) ) @aa.grid_dec.grid_2d_to_structure @aa.grid_dec.transform @aa.grid_dec.relocate_to_radial_minimum def convergence_2d_from(self, grid: aa.type.Grid2DLike): """Calculate the projected convergence at a given set of arc-second gridded coordinates. Parameters ---------- grid The grid of (y,x) arc-second coordinates the convergence is computed on. """ return self.convergence_func(self.grid_to_eccentric_radii(grid)) def convergence_func(self, grid_radius: float) -> float: return ( self.mass_to_light_ratio * ( ((self.axis_ratio * grid_radius) / self.effective_radius) ** -self.mass_to_light_gradient ) * self.image_2d_via_radii_from(grid_radius) ) def decompose_convergence_via_mge(self): radii_min = self.effective_radius / 100.0 radii_max = self.effective_radius * 20.0 def sersic_radial_gradient_2D(r): return ( self.mass_to_light_ratio * self.intensity * ( ((self.axis_ratio * r) / self.effective_radius) ** -self.mass_to_light_gradient ) * np.exp( -self.sersic_constant * (((r / self.effective_radius) ** (1.0 / self.sersic_index)) - 1.0) ) ) return self._decompose_convergence_via_mge( func=sersic_radial_gradient_2D, radii_min=radii_min, radii_max=radii_max ) def decompose_convergence_via_cse(self) -> Tuple[List, List]: """ Decompose the convergence of the Sersic profile into singular isothermal elliptical (sie) profiles. This decomposition uses the standard 2d profile of a Sersic mass profile. Parameters ---------- func The function representing the profile that is decomposed into CSEs. radii_min: The minimum radius to fit radii_max: The maximum radius to fit total_sies The number of SIEs used to approximate the input func. sample_points: int (should be larger than 'total_sies') The number of data points to fit Returns ------- Tuple[List, List] A list of amplitudes and core radii of every singular isothernal ellipsoids (sie) the mass profile is decomposed into. """ upper_dex, lower_dex, total_cses, sample_points = cse_settings_from( effective_radius=self.effective_radius, sersic_index=self.sersic_index, sersic_constant=self.sersic_constant, mass_to_light_gradient=self.mass_to_light_gradient, ) scaled_effective_radius = self.effective_radius / np.sqrt(self.axis_ratio) radii_min = scaled_effective_radius / 10.0 ** lower_dex radii_max = scaled_effective_radius * 10.0 ** upper_dex def sersic_radial_gradient_2D(r): return ( self.mass_to_light_ratio * self.intensity * ( ((self.axis_ratio * r) / scaled_effective_radius) ** -self.mass_to_light_gradient ) * np.exp( -self.sersic_constant * ( ((r / scaled_effective_radius) ** (1.0 / self.sersic_index)) - 1.0 ) ) ) return self._decompose_convergence_via_cse_from( func=sersic_radial_gradient_2D, radii_min=radii_min, radii_max=radii_max, total_cses=total_cses, sample_points=sample_points, )
[docs]class SphSersicRadialGradient(EllSersicRadialGradient):
[docs] def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), intensity: float = 0.1, effective_radius: float = 0.6, sersic_index: float = 0.6, mass_to_light_ratio: float = 1.0, mass_to_light_gradient: float = 0.0, ): """ Setup a Sersic mass and light profiles. Parameters ---------- centre The (y,x) arc-second coordinates of the profile centre. intensity Overall flux intensity normalisation in the light profiles (electrons per second). effective_radius The circular radius containing half the light of this profile. sersic_index Controls the concentration of the profile (lower -> less concentrated, higher -> more concentrated). mass_to_light_ratio The mass-to-light ratio of the light profile. mass_to_light_gradient The mass-to-light radial gradient. """ super().__init__( centre=centre, elliptical_comps=(0.0, 0.0), intensity=intensity, effective_radius=effective_radius, sersic_index=sersic_index, mass_to_light_ratio=mass_to_light_ratio, mass_to_light_gradient=mass_to_light_gradient, )
class EllSersicCore(EllSersic): def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), elliptical_comps: Tuple[float, float] = (0.0, 0.0), effective_radius: float = 0.6, sersic_index: float = 4.0, radius_break: float = 0.01, intensity_break: float = 0.05, gamma: float = 0.25, alpha: float = 3.0, mass_to_light_ratio: float = 1.0, ): """ The elliptical cored-Sersic light profile. Parameters ---------- centre The (y,x) arc-second coordinates of the profile centre. elliptical_comps The first and second ellipticity components of the elliptical coordinate system, (see the module `autogalaxy -> convert.py` for the convention). intensity Overall intensity normalisation of the light profile (units are dimensionless and derived from the data the light profile's image is compared too, which is expected to be electrons per second). effective_radius The circular radius containing half the light of this profile. sersic_index Controls the concentration of the profile (lower -> less concentrated, higher -> more concentrated). radius_break The break radius separating the inner power-law (with logarithmic slope gamma) and outer Sersic function. intensity_break The intensity at the break radius. gamma The logarithmic power-law slope of the inner core profiles alpha : Controls the sharpness of the transition between the inner core / outer Sersic profiles. """ super().__init__( centre=centre, elliptical_comps=elliptical_comps, intensity=intensity_break, effective_radius=effective_radius, sersic_index=sersic_index, mass_to_light_ratio=mass_to_light_ratio, ) self.radius_break = radius_break self.intensity_break = intensity_break self.alpha = alpha self.gamma = gamma def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike): return self.deflections_2d_via_mge_from(grid=grid) def image_2d_via_radii_from(self, grid_radii: np.ndarray): """ Calculate the intensity of the cored-Sersic light profile on a grid of radial coordinates. Parameters ---------- grid_radii The radial distance from the centre of the profile. for each coordinate on the grid. """ return np.multiply( np.multiply( self.intensity_prime, np.power( np.add( 1, np.power(np.divide(self.radius_break, grid_radii), self.alpha), ), (self.gamma / self.alpha), ), ), np.exp( np.multiply( -self.sersic_constant, ( np.power( np.divide( np.add( np.power(grid_radii, self.alpha), (self.radius_break ** self.alpha), ), (self.effective_radius ** self.alpha), ), (1.0 / (self.alpha * self.sersic_index)), ) ), ) ), ) def decompose_convergence_via_mge(self): radii_min = self.effective_radius / 50.0 radii_max = self.effective_radius * 20.0 def core_sersic_2D(r): return ( self.mass_to_light_ratio * self.intensity_prime * (1.0 + (self.radius_break / r) ** self.alpha) ** (self.gamma / self.alpha) * np.exp( -self.sersic_constant * ( (r ** self.alpha + self.radius_break ** self.alpha) / self.effective_radius ** self.alpha ) ** (1.0 / (self.sersic_index * self.alpha)) ) ) return self._decompose_convergence_via_mge( func=core_sersic_2D, radii_min=radii_min, radii_max=radii_max ) @property def intensity_prime(self): """Overall intensity normalisation in the rescaled Core-Sersic light profiles (electrons per second)""" return ( self.intensity_break * (2.0 ** (-self.gamma / self.alpha)) * np.exp( self.sersic_constant * ( ((2.0 ** (1.0 / self.alpha)) * self.radius_break) / self.effective_radius ) ** (1.0 / self.sersic_index) ) ) class SphSersicCore(EllSersicCore): def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), effective_radius: float = 0.6, sersic_index: float = 4.0, radius_break: float = 0.01, intensity_break: float = 0.05, gamma: float = 0.25, alpha: float = 3.0, ): """ The elliptical cored-Sersic light profile. Parameters ---------- centre The (y,x) arc-second coordinates of the profile centre. intensity Overall intensity normalisation of the light profile (units are dimensionless and derived from the data the light profile's image is compared too, which is expected to be electrons per second). effective_radius The circular radius containing half the light of this profile. sersic_index Controls the concentration of the profile (lower -> less concentrated, higher -> more concentrated). radius_break The break radius separating the inner power-law (with logarithmic slope gamma) and outer Sersic function. intensity_break The intensity at the break radius. gamma The logarithmic power-law slope of the inner core profiles alpha : Controls the sharpness of the transition between the inner core / outer Sersic profiles. """ super().__init__( centre=centre, elliptical_comps=(0.0, 0.0), effective_radius=effective_radius, sersic_index=sersic_index, radius_break=radius_break, intensity_break=intensity_break, gamma=gamma, alpha=alpha, ) self.radius_break = radius_break self.intensity_break = intensity_break self.alpha = alpha self.gamma = gamma
[docs]class EllChameleon(MassProfile, StellarProfile):
[docs] def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), elliptical_comps: Tuple[float, float] = (0.0, 0.0), intensity: float = 0.1, core_radius_0: float = 0.01, core_radius_1: float = 0.02, mass_to_light_ratio: float = 1.0, ): """ The elliptical Chamelon mass profile. Parameters ---------- centre The (y,x) arc-second coordinates of the profile centre. elliptical_comps The first and second ellipticity components of the elliptical coordinate system, (see the module `autogalaxy -> convert.py` for the convention). intensity Overall intensity normalisation of the light profile (units are dimensionless and derived from the data the light profile's image is compared too, which is expected to be electrons per second). core_radius_0 : the core size of the first elliptical cored Isothermal profile. core_radius_1 : core_radius_0 + core_radius_1 is the core size of the second elliptical cored Isothermal profile. We use core_radius_1 here is to avoid negative values. Profile form: mass_to_light_ratio * intensity *\ (1.0 / Sqrt(x^2 + (y/q)^2 + core_radius_0^2) - 1.0 / Sqrt(x^2 + (y/q)^2 + (core_radius_0 + core_radius_1)**2.0)) """ super(EllChameleon, self).__init__( centre=centre, elliptical_comps=elliptical_comps ) super(MassProfile, self).__init__( centre=centre, elliptical_comps=elliptical_comps ) self.mass_to_light_ratio = mass_to_light_ratio self.intensity = intensity self.core_radius_0 = core_radius_0 self.core_radius_1 = core_radius_1
def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike): return self.deflections_2d_via_analytic_from(grid=grid) @aa.grid_dec.grid_2d_to_structure @aa.grid_dec.transform @aa.grid_dec.relocate_to_radial_minimum def deflections_2d_via_analytic_from(self, grid: aa.type.Grid2DLike): """ Calculate the deflection angles at a given set of arc-second gridded coordinates. Following Eq. (15) and (16), but the parameters are slightly different. Parameters ---------- grid The grid of (y,x) arc-second coordinates the deflection angles are computed on. """ factor = ( 2.0 * self.mass_to_light_ratio * self.intensity / (1 + self.axis_ratio) * self.axis_ratio / np.sqrt(1.0 - self.axis_ratio ** 2.0) ) core_radius_0 = np.sqrt( (4.0 * self.core_radius_0 ** 2.0) / (1.0 + self.axis_ratio) ** 2 ) core_radius_1 = np.sqrt( (4.0 * self.core_radius_1 ** 2.0) / (1.0 + self.axis_ratio) ** 2 ) psi0 = psi_from( grid=grid, axis_ratio=self.axis_ratio, core_radius=core_radius_0 ) psi1 = psi_from( grid=grid, axis_ratio=self.axis_ratio, core_radius=core_radius_1 ) deflection_y0 = np.arctanh( np.divide( np.multiply(np.sqrt(1.0 - self.axis_ratio ** 2.0), grid[:, 0]), np.add(psi0, self.axis_ratio ** 2.0 * core_radius_0), ) ) deflection_x0 = np.arctan( np.divide( np.multiply(np.sqrt(1.0 - self.axis_ratio ** 2.0), grid[:, 1]), np.add(psi0, core_radius_0), ) ) deflection_y1 = np.arctanh( np.divide( np.multiply(np.sqrt(1.0 - self.axis_ratio ** 2.0), grid[:, 0]), np.add(psi1, self.axis_ratio ** 2.0 * core_radius_1), ) ) deflection_x1 = np.arctan( np.divide( np.multiply(np.sqrt(1.0 - self.axis_ratio ** 2.0), grid[:, 1]), np.add(psi1, core_radius_1), ) ) deflection_y = np.subtract(deflection_y0, deflection_y1) deflection_x = np.subtract(deflection_x0, deflection_x1) return self.rotate_grid_from_reference_frame( np.multiply(factor, np.vstack((deflection_y, deflection_x)).T) ) @aa.grid_dec.grid_2d_to_structure @aa.grid_dec.transform @aa.grid_dec.relocate_to_radial_minimum def convergence_2d_from(self, grid: aa.type.Grid2DLike): """Calculate the projected convergence at a given set of arc-second gridded coordinates. Parameters ---------- grid The grid of (y,x) arc-second coordinates the convergence is computed on. """ return self.convergence_func(self.grid_to_elliptical_radii(grid)) def convergence_func(self, grid_radius: float) -> float: return self.mass_to_light_ratio * self.image_2d_via_radii_from(grid_radius) @aa.grid_dec.grid_2d_to_structure def potential_2d_from(self, grid: aa.type.Grid2DLike): return np.zeros(shape=grid.shape[0]) def image_2d_via_radii_from(self, grid_radii: np.ndarray): """Calculate the intensity of the Chamelon light profile on a grid of radial coordinates. Parameters ---------- grid_radii The radial distance from the centre of the profile. for each coordinate on the grid. """ axis_ratio_factor = (1.0 + self.axis_ratio) ** 2.0 return np.multiply( self.intensity / (1 + self.axis_ratio), np.add( np.divide( 1.0, np.sqrt( np.add( np.square(grid_radii), (4.0 * self.core_radius_0 ** 2.0) / axis_ratio_factor, ) ), ), -np.divide( 1.0, np.sqrt( np.add( np.square(grid_radii), (4.0 * self.core_radius_1 ** 2.0) / axis_ratio_factor, ) ), ), ), ) @property def axis_ratio(self): axis_ratio = super().axis_ratio return axis_ratio if axis_ratio < 0.99999 else 0.99999 def with_new_normalization(self, normalization): mass_profile = copy.copy(self) mass_profile.mass_to_light_ratio = normalization return mass_profile
[docs]class SphChameleon(EllChameleon):
[docs] def __init__( self, centre: Tuple[float, float] = (0.0, 0.0), intensity: float = 0.1, core_radius_0: float = 0.01, core_radius_1: float = 0.02, mass_to_light_ratio: float = 1.0, ): """ The spherica; Chameleon mass profile. Profile form: mass_to_light_ratio * intensity *\ (1.0 / Sqrt(x^2 + (y/q)^2 + core_radius_0^2) - 1.0 / Sqrt(x^2 + (y/q)^2 + (core_radius_0 + core_radius_1)**2.0)) Parameters ---------- centre The (y,x) arc-second coordinates of the profile centre. elliptical_comps The first and second ellipticity components of the elliptical coordinate system, (see the module `autogalaxy -> convert.py` for the convention). intensity Overall intensity normalisation of the light profile (units are dimensionless and derived from the data the light profile's image is compared too, which is expected to be electrons per second). core_radius_0 : the core size of the first elliptical cored Isothermal profile. core_radius_1 : core_radius_0 + core_radius_1 is the core size of the second elliptical cored Isothermal profile. We use core_radius_1 here is to avoid negative values. """ super().__init__( centre=centre, elliptical_comps=(0.0, 0.0), intensity=intensity, core_radius_0=core_radius_0, core_radius_1=core_radius_1, mass_to_light_ratio=mass_to_light_ratio, )
def cse_settings_from( effective_radius, sersic_index, sersic_constant, mass_to_light_gradient ): if mass_to_light_gradient > 0.5: if effective_radius > 0.2: lower_dex = 6.0 upper_dex = np.min( [np.log10((18.0 / sersic_constant) ** sersic_index), 1.1] ) if sersic_index <= 1.2: total_cses = 50 sample_points = 80 elif sersic_index > 3.8: total_cses = 40 sample_points = 50 lower_dex = 6.5 else: total_cses = 30 sample_points = 50 else: if sersic_index <= 1.2: upper_dex = 1.0 total_cses = 50 sample_points = 80 lower_dex = 4.5 elif sersic_index > 3.8: total_cses = 40 sample_points = 50 lower_dex = 6.0 upper_dex = 1.5 else: upper_dex = 1.1 lower_dex = 6.0 total_cses = 30 sample_points = 50 else: upper_dex = np.min( [ np.log10((23.0 / sersic_constant) ** sersic_index), 0.85 - np.log10(effective_radius), ] ) if (sersic_index <= 0.9) and (sersic_index > 0.8): total_cses = 50 sample_points = 80 upper_dex = np.log10((18.0 / sersic_constant) ** sersic_index) lower_dex = 4.3 + np.log10(effective_radius) elif sersic_index <= 0.8: total_cses = 50 sample_points = 80 upper_dex = np.log10((16.0 / sersic_constant) ** sersic_index) lower_dex = 4.0 + np.log10(effective_radius) elif sersic_index > 3.8: total_cses = 40 sample_points = 50 lower_dex = 4.5 + np.log10(effective_radius) else: lower_dex = 3.5 + np.log10(effective_radius) total_cses = 30 sample_points = 50 return upper_dex, lower_dex, total_cses, sample_points