Source code for dysmalpy.models.baryons

# coding=utf8
# Copyright (c) MPE/IR-Submm Group. See LICENSE.rst for license information. 
#
# Baryon mass models for DysmalPy

from __future__ import (absolute_import, division, print_function,
                        unicode_literals)

# Standard library
import os, copy
import logging
import glob

# Third party imports
import numpy as np
import scipy.special as scp_spec
import scipy.integrate as scp_int
import scipy.optimize as scp_opt
import scipy.interpolate as scp_interp
import scipy.optimize as scp_opt
import scipy.integrate as scp_int
import astropy.constants as apy_con


from astropy.table import Table

# Local imports
from .base import MassModel, _LightMassModel, v_circular, \
                  sersic_mr, _I0_gaussring
# from .base import menc_from_vcirc
from dysmalpy.parameters import DysmalParameter

__all__ = ['Sersic', 'DiskBulge', 'LinearDiskBulge', 'ExpDisk', 'BlackHole',
           'GaussianRing',
           'surf_dens_exp_disk', 'menc_exp_disk', 'vcirc_exp_disk', 'sersic_menc_2D_proj',
           'mass_comp_conditional_ring',
           'NoordFlat', 'InfThinMassiveGaussianRing']

# NEW, RECALCULATED NOORDERMEER DIRECTORY: INCLUDES dlnrho/dlnr; MASS PROFILES
path = os.path.abspath(__file__)
dir_path = os.path.dirname(path)
# located one up:
dir_path = os.sep.join(dir_path.split(os.sep)[:-1])
_dir_deprojected_sersic_models = os.sep.join([dir_path, "data",
                                "deprojected_sersic_models_tables", ""])



# MASSIVE RING DIRECTORIES:
_dir_gaussian_ring_tables = os.getenv('GAUSSIAN_RING_PROFILE_DATADIR', None)

# CONSTANTS
G = apy_con.G
Msun = apy_con.M_sun
pc = apy_con.pc

# # +++++++++++++++++++++++++++++
# # TEMP:
# G = 6.67e-11 * u.m**3 / u.kg / (u.s**2)  #(unit='m3 / (kg s2)')
# Msun = 2e30 * u.kg
# pc = 3e16 * u.m
# # +++++++++++++++++++++++++++++

# LOGGER SETTINGS
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger('DysmalPy')
logger.setLevel(logging.INFO)


import warnings
warnings.filterwarnings("ignore")



[docs] def surf_dens_exp_disk(r, mass, rd): """ Radial surface density function for an infinitely thin exponential disk Parameters ---------- r : float or array Radius or radii at which to calculate the surface density mass : float Total mass of the disk rd : float Disk scale length. Returns ------- Sigr : float or array Surface density of a thin exponential disk at `r` """ Sig0 = mass / (2. * np.pi * rd**2) Sigr = Sig0 * np.exp(-r/rd) return Sigr
[docs] def menc_exp_disk(r, mass, rd): """ Enclosed mass function for exponential disk Parameters ---------- r : float or array Radius or radii at which to calculate the surface density mass : float Total mass of the disk rd : float Disk scale length. Returns ------- menc : float or array Enclosed mass of an exponential disk for the given `r` """ Sig0 = mass / (2. * np.pi * rd**2) menc = 2. * np.pi * Sig0 * rd**2 * ( 1 - np.exp(-r/rd)*(1.+r/rd) ) return menc
[docs] def vcirc_exp_disk(r, mass, rd): """ Rotation curve function for exponential disk Parameters ---------- r : float or array Radius or radii at which to calculate the surface density mass : float Total mass of the disk rd : float Disk scale length. Returns ------- vc : float or array Circular velocity of an exponential disk as a function of `r` """ #b1 = 1.6783469900166612 # scp_spec.gammaincinv(2.*n, 0.5), n=1 #rd = r_eff / b1 Sig0 = mass / (2. * np.pi * rd**2) y = r / (2.*rd) expdisk = y**2 * ( scp_spec.i0(y) * scp_spec.k0(y) - scp_spec.i1(y)*scp_spec.k1(y) ) VCsq = 4 * np.pi * G.cgs.value*Msun.cgs.value / (1000.*pc.cgs.value) * Sig0 * rd * expdisk VCsq[r==0] = 0. return np.sqrt(VCsq) / 1.e5
[docs] def sersic_menc_2D_proj(r, mass, n, r_eff): """ Enclosed mass as a function of r for a generic sersic model Parameters ---------- r : float or array Radius or radii at which to calculate the surface mass density mass : float Total mass of the Sersic component n : float Sersic index r_eff : float Effective radius Returns ------- menc : float or array Enclosed mass as a function of `r` Notes ----- This function is only valid in the case of an infinite cylinder """ bn = scp_spec.gammaincinv(2. * n, 0.5) integ = scp_spec.gammainc(2 * n, bn * (r / r_eff) ** (1. / n)) norm = mass menc = norm*integ return menc
def mass_comp_conditional_ring(param, modelset): """ Basic conditional prior on the mass of the other component(s) (i.e., a bulge or halo) when fitting both a massive ring and one or more other mass components, returning True/False This conditional could apply to the bulge mass, or to the halo mass (probably best to not do both). Intended to be passed as f_bounds function when setting, e.g., sersic_comp.total_mass.prior = ConditionalEmpiricalUniformPrior(f_cond=mass_comp_conditional_bounds_ring) """ # Double check param + model values are same: if param.value != modelset.components[param._model._name].__getattribute__(param._name).value: raise ValueError # Test radius array: rarr = np.arange(0., 10.1, 0.1) return np.all(np.isfinite(modelset.circular_velocity(rarr))) ##########################
[docs] class NoordFlat(object): """ Class to handle circular velocities / enclosed mass profiles for a thick Sersic component. Lookup tables are numerically calculated from the derivations provided in `Noordermeer 2008 <https://ui.adsabs.harvard.edu/abs/2008MNRAS.385.1359N/abstract>`_ which properly accounted for the thickness of the mass component. The lookup table provides rotation curves for Sersic components with `n` = 0.5 - 8 at steps of 0.1 and `invq` = [1, 2, 3, 4, 5, 6, 7, 8, 10, 20, 100]. If the given `n` and/or `invq` are not one of these values then the nearest ones are used. References ---------- `Noordermeer 2008 <https://ui.adsabs.harvard.edu/abs/2008MNRAS.385.1359N/abstract>`_ Parameters ---------- n : float Sersic index invq: float Sersic index """ def __init__(self, n=None, invq=None): self._n = n self._invq = invq self._n_current = None self._invq_current = None self.rho_interp_func = None self.dlnrhodlnr_interp_func = None self._rho_interp_type = None self._dlnrhodlnr_interp_type = None self._reset_interps() @property def n(self): return self._n @n.setter def n(self, value): if value < 0: raise ValueError("Sersic index can't be negative!") self._n = value # Reset vcirc interp: self._reset_interps() @property def invq(self): return self._invq @invq.setter def invq(self, value): if value < 0: raise ValueError("Invq can't be negative!") self._invq = value # Reset vcirc interp: self._reset_interps() def _reset_interps(self): self._set_vcirc_interp() self._set_menc_interp() if self.rho_interp_func is not None: self._set_rho_interp(interp_type=self._rho_interp_type) if self.dlnrhodlnr_interp_func is not None: self._set_dlnrhodlnr_interp(interp_type=self._dlnrhodlnr_interp_type)
[docs] def read_deprojected_sersic_table(self): # Use the "typical" collection of table values: table_n = np.arange(0.5, 8.1, 0.1) # Sersic indices table_invq = np.array([1., 2., 3., 4., 5., 6., 7., 8., 10., 20., 100., 1.11, 1.43, 1.67, 3.33, 0.5, 0.67]) # 1:1, 1:2, 1:3, ... flattening [also prolate 2:1, 1.5:1] nearest_n = table_n[ np.argmin( np.abs(table_n - self.n) ) ] nearest_invq = table_invq[ np.argmin( np.abs( table_invq - self.invq) ) ] file_sersic = _dir_deprojected_sersic_models + 'deproj_sersic_model_n{:0.1f}_invq{:0.2f}.fits'.format(nearest_n, nearest_invq) try: t = Table.read(file_sersic) except: # REMOVE BACKWARDS COMPATIBILITY! raise ValueError("File {} not found. _dir_deprojected_sersic_models={}.".format(file_sersic)) return t[0]
def _set_vcirc_interp(self): # SHOULD BE EXACTLY, w/in numerical limitations, EQUIV TO OLD CALCULATION table = self.read_deprojected_sersic_table() N2008_vcirc = table['vcirc'] N2008_rad = table['R'] self.N2008_Re = table['Reff'] self.N2008_mass = table['total_mass'] self.vcirc_interp = scp_interp.interp1d(N2008_rad, N2008_vcirc, fill_value="extrapolate") # vcirc = (v_interp(r / r_eff * N2008_Re) * np.sqrt( # mass / N2008_mass) * np.sqrt(N2008_Re / r_eff)) # return vcirc def _set_menc_interp(self): table = self.read_deprojected_sersic_table() table_Rad = table['R'] table_menc = table['menc3D_sph'] # Clean up values inside rmin: Add the value at r=0: menc=0 if table['R'][0] > 0.: table_Rad = np.append(0., table_Rad) table_menc = np.append(0., table_menc) self.menc_interp = scp_interp.interp1d(table_Rad, table_menc, fill_value="extrapolate")
[docs] def circular_velocity(self, r, r_eff, mass): """ Calculate circular velocity for a thick Sersic component, by interpolating Parameters ---------- r : float or array Radius or radii at which to calculate the circular velocity in kpc r_eff : float Effective radius of the Sersic component in kpc mass : float Total mass of the Sersic component Returns ------- vcirc : float or array Circular velocity at each given `r`, in km/s Notes ----- This function determines the circular velocity as a function of radius for a Sersic component with a total mass, `mass`, Sersic index, `n`, and an effective radius to scale height ratio, `invq`. This uses lookup tables numerically calculated from the derivations provided in `Noordermeer 2008 <https://ui.adsabs.harvard.edu/abs/2008MNRAS.385.1359N/abstract>`_ which properly account for the thickness of the mass component. References ---------- `Noordermeer 2008 <https://ui.adsabs.harvard.edu/abs/2008MNRAS.385.1359N/abstract>`_ """ vcirc = (self.vcirc_interp(r / r_eff * self.N2008_Re) * np.sqrt( mass / self.N2008_mass) * np.sqrt(self.N2008_Re / r_eff)) return vcirc
[docs] def enclosed_mass(self, r, r_eff, mass): """ Calculate enclosed mass for a thick Sersic component, by interpolating Parameters ---------- r : float or array Radius or radii at which to calculate the circular velocity in kpc r_eff : float Effective radius of the Sersic component in kpc mass : float Total mass of the Sersic component Returns ------- menc : float or array Enclosed mass (in a sphere) at each given `r`, in solar masses Notes ----- This function determines the spherical enclosed mass as a function of radius for a Sersic component with a total mass, `mass`, Sersic index, `n`, and an effective radius to scale height ratio, `invq`. This uses lookup tables numerically calculated from the derivations provided in `Noordermeer 2008 <https://ui.adsabs.harvard.edu/abs/2008MNRAS.385.1359N/abstract>`_ (as extended in `Price et al. 2022 <https://ui.adsabs.harvard.edu/abs/2022A%26A...665A.159P/abstract>`_) which properly account for the thickness of the mass component. References ---------- `Noordermeer 2008 <https://ui.adsabs.harvard.edu/abs/2008MNRAS.385.1359N/abstract>`_ `Price et al. 2022 <https://ui.adsabs.harvard.edu/abs/2022A%26A...665A.159P/abstract>`_ """ menc = self.menc_interp(r / r_eff * self.N2008_Re) * (mass / self.N2008_mass) # # TEST: # print("USING OLD vcirc from menc! (baryons.py)") # menc = menc_from_vcirc(self.circular_velocity(r, r_eff, mass), r) return menc
[docs] def rho(self, r, Reff, total_mass, interp_type='linear'): if (self._rho_interp_type != interp_type) | (self.rho_interp_func is None): # Update rho funcs: self._set_rho_interp(interp_type=interp_type) # Ensure it's an array: if isinstance(r*1., float): rarr = np.array([r]) else: rarr = np.array(r) # Ensure all radii are 0. or positive: rarr = np.abs(rarr) scale_fac = (total_mass / self.table_mass) * (self.table_Reff / Reff)**3 if interp_type.lower().strip() == 'cubic': rho_interp = np.zeros(len(rarr)) wh_in = np.where((rarr <= self.table_rad_rho.max()) & (rarr >= self.table_rad_rho.min()))[0] wh_extrap = np.where((rarr > self.table_rad_rho.max()) | (rarr < self.table_rad_rho.min()))[0] rho_interp[wh_in] = (self.rho_interp_func(rarr[wh_in] / Reff * self.table_Reff) * scale_fac ) rho_interp[wh_extrap] = (self.rho_interp_extrap_func(rarr[wh_extrap] / Reff * self.table_Reff) * scale_fac) elif interp_type.lower().strip() == 'linear': rho_interp = (self.rho_interp_func(rarr / Reff * self.table_Reff) * scale_fac ) else: raise ValueError("interp type '{}' unknown!".format(interp_type)) if (len(rarr) > 1): return rho_interp else: if isinstance(r*1., float): # Float input return rho_interp[0] else: # Length 1 array input return rho_interp
[docs] def dlnrho_dlnr(self, r, Reff, interp_type='linear'): """ Calculate log mass density gradient for a thick Sersic component, by interpolating. Can be used to determine an alternative pressure support correction. References ---------- `Noordermeer 2008 <https://ui.adsabs.harvard.edu/abs/2008MNRAS.385.1359N/abstract>`_ `Price et al. 2022 <https://ui.adsabs.harvard.edu/abs/2022A%26A...665A.159P/abstract>`_ """ if (self._dlnrhodlnr_interp_type != interp_type) | (self.dlnrhodlnr_interp_func is None): # Update rho funcs: self._set_dlnrhodlnr_interp(interp_type=interp_type) # Ensure it's an array: if isinstance(r*1., float): rarr = np.array([r]) else: rarr = np.array(r) # Ensure all radii are 0. or positive: rarr = np.abs(rarr) if interp_type.lower().strip() == 'cubic': dlnrho_dlnr_interp = np.zeros(len(rarr)) wh_in = np.where((rarr <= self.table_rad_dlnrhodlnr.max()) & \ (rarr >= self.table_rad_dlnrhodlnr.min()))[0] wh_extrap = np.where((rarr > self.table_rad_dlnrhodlnr.max()) | \ (rarr < self.table_rad_dlnrhodlnr.min()))[0] dlnrho_dlnr_interp[wh_in] = (self.dlnrhodlnr_interp_func(rarr[wh_in] / Reff * self.table_Reff) ) dlnrho_dlnr_interp[wh_extrap] = (self.dlnrhodlnr_interp_func_extrap(rarr[wh_extrap] / Reff * self.table_Reff)) elif interp_type.lower().strip() == 'linear': dlnrho_dlnr_interp = (self.dlnrhodlnr_interp_func(rarr / Reff * self.table_Reff) ) else: raise ValueError("interp type '{}' unknown!".format(interp_type)) if (len(rarr) > 1): return dlnrho_dlnr_interp else: if isinstance(r*1., float): # Float input return dlnrho_dlnr_interp[0] else: # Length 1 array input return dlnrho_dlnr_interp return dlnrho_dlnr_interp
def _set_rho_interp(self, interp_type='linear'): self._rho_interp_type = interp_type table = self.read_deprojected_sersic_table() table_rho = table['rho'] table_rad = table['R'] self.table_Reff = table['Reff'] self.table_mass = table['total_mass'] # Drop nonfinite parts: whfin = np.where(np.isfinite(table_rho))[0] table_rho = table_rho[whfin] table_rad = table_rad[whfin] self.table_rad_rho = table_rad if interp_type.lower().strip() == 'cubic': self.rho_interp_func = scp_interp.interp1d(table_rad, table_rho, fill_value=np.NaN, bounds_error=False, kind='cubic') self.rho_interp_extrap_func = scp_interp.interp1d(table_rad, table_rho, fill_value='extrapolate', kind='linear') elif interp_type.lower().strip() == 'linear': self.rho_interp_func = scp_interp.interp1d(table_rad, table_rho, fill_value='extrapolate', bounds_error=False, kind='linear') else: raise ValueError("interp type '{}' unknown!".format(interp_type)) def _set_dlnrhodlnr_interp(self, interp_type='linear'): self._dlnrhodlnr_interp_type = interp_type table = self.read_deprojected_sersic_table() table_dlnrho_dlnr = table['dlnrho_dlnR'] table_rad = table['R'] self.table_Reff = table['Reff'] self.table_mass = table['total_mass'] # Drop nonfinite parts: whfin = np.where(np.isfinite(table_dlnrho_dlnr))[0] table_dlnrho_dlnr = table_dlnrho_dlnr[whfin] table_rad = table_rad[whfin] self.table_rad_dlnrhodlnr = table_rad if interp_type.lower().strip() == 'cubic': self.dlnrhodlnr_interp_func = scp_interp.interp1d(table_rad, table_dlnrho_dlnr, fill_value=np.NaN, bounds_error=False, kind='cubic') self.dlnrhodlnr_interp_func_extrap = scp_interp.interp1d(table_rad, table_dlnrho_dlnr, fill_value='extrapolate', kind='linear') elif interp_type.lower().strip() == 'linear': self.dlnrhodlnr_interp_func = scp_interp.interp1d(table_rad, table_dlnrho_dlnr, fill_value='extrapolate', bounds_error=False, kind='linear') else: raise ValueError("interp type '{}' unknown!".format(interp_type))
##########################
[docs] class InfThinMassiveGaussianRing(object): """ Class to handle circular velocities / enclosed mass profiles for an infinitely thin, massive Gaussian ring. Lookup tables are numerically calculated, following B&T and Bovy online galaxies textbook. The lookup table provides rotation curves for Gaussian rings with `invh` = 0.5 - 8 at steps of 0.1 and `invq` = [1, 2, 3, 4, 5, 6, 7, 8, 10, 20, 100], where :math:`\mathrm{invh} = R_{\mathrm{peak}} / \mathrm{FWHM}_{\mathrm{ring}}` If the given `n` and/or `invq` are not one of these values then the nearest ones are used. Parameters ---------- invh : float Ratio of `R_peak` / `ring_FWHM` """ def __init__(self, invh=None): self._invh = invh self._invh_current = None self._reset_interps() @property def invh(self): return self._invh @invh.setter def invh(self, value): if value < 0: #raise ValueError("Invh can't be negative!") logger.warning('Invh is negative -- undefined, so interps will be all NaN!!') self._invh = value # Reset vcirc interp: self._reset_interps() def _reset_interps(self): # self._set_vcirc_interp() self._set_potential_gradient_interp() self._set_menc_interp()
[docs] def read_ring_table(self): # Use the "typical" collection of table values: #-------------------------------- # Glob values from path: table_invh = [] name_base = 'gauss_ring_profile_invh' fnames_glob = glob.glob(_dir_gaussian_ring_tables+name_base+'*.fits') for fn in fnames_glob: invh_str = fn.split(name_base)[-1].split('.fits')[0] table_invh.append(float(invh_str)) table_invh = np.array(table_invh) table_invh.sort() # #-------------------------------- # table_invh = np.array([0.01, 0.05, # 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., # 0.25, 0.75, # 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2., # 2.25, 2.5, 2.75, 3., 3.5, # 3.33, 6.67, # 4., 4.5, 5., 5.5, 6., 6.5, 7., 7.5, # 8., 8.5, 9., 9.5, 10., 12.5, # 15., 20., 25., 50.]) # #-------------------------------- nearest_invh = table_invh[ np.argmin( np.abs( table_invh - self.invh) ) ] if self.invh != nearest_invh: logger.warning('Using table for massive gausring with NON-EXACT invh!!' 'Model invh = {} ; table invh = {}'.format(self.invh, nearest_invh)) file_gausring = _dir_gaussian_ring_tables + 'gauss_ring_profile_invh{:0.2f}.fits'.format(nearest_invh) try: t = Table.read(file_gausring) except: raise ValueError("File {} not found. _dir_gaussian_ring_tables={}. Check that system var ${} is set correctly.".format(file_gausring, _dir_gaussian_ring_tables, 'GAUSSIAN_RING_PROFILE_DATADIR')) return t[0]
# def _set_vcirc_interp(self): # table = self.read_ring_table() # # tab_rad = table['R'] # tab_vcirc = table['vcirc'] # self.tab_invh = table['invh'] # self.tab_R_peak = table['R_peak'] # self.tab_ring_FWHM = table['ring_FWHM'] # self.tab_mass = table['total_mass'] # # self.vcirc_interp = scp_interp.interp1d(tab_rad, tab_vcirc, # fill_value="extrapolate") # # # scale_fac = np.sqrt(total_mass / table_mass) * np.sqrt(table_Rpeak / R_peak) # # vcirc_interp = v_interp(Rarr / R_peak * table_Rpeak) * scale_fac def _set_potential_gradient_interp(self): if self.invh < 0: self.potl_grad_interp = None else: table = self.read_ring_table() tab_rad = table['R'] tab_potl_grad = table['potential_gradient'] self.tab_invh = table['invh'] self.tab_R_peak = table['R_peak'] self.tab_ring_FWHM = table['ring_FWHM'] self.tab_mass = table['total_mass'] self.potl_grad_interp = scp_interp.interp1d(tab_rad, tab_potl_grad, fill_value="extrapolate") # scale_fac = (total_mass / table_mass) * (table_Rpeak / R_peak)**2 # potential_gradient_interp = potl_grad_interp(Rarr / R_peak * table_Rpeak) * scale_fac def _set_menc_interp(self): if self.invh < 0: self.menc_interp = None else: table = self.read_ring_table() tab_rad = table['R'] tab_menc = table['menc'] self.tab_invh = table['invh'] self.tab_R_peak = table['R_peak'] self.tab_ring_FWHM = table['ring_FWHM'] self.tab_mass = table['total_mass'] self.menc_interp = scp_interp.interp1d(tab_rad, tab_menc, fill_value="extrapolate") # scale_fac = (total_mass / table_mass) # menc_interp = m_interp(Rarr / R_peak * table_Rpeak) * scale_fac # def vcirc(self, R, R_peak, total_mass): # """ # Calculate circular velocity for a inf thin massive gaussian Ring # # Parameters # ---------- # R : float or array # Radius or radii at which to calculate the circular velocity in kpc # # R_peak : float # Peak of Gaussian ring in kpc # # total_mass : float # Total mass of the Gaussian ring component # # Returns # ------- # vcirc : float or array # Circular velocity at each given `R` # # Notes # ----- # This function determines the circular velocity as a function of radius for # a massive Gaussian ring component with a total mass, `total_mass`, # and a ring peak radius to ring FWHM ratio, `invh`. # This uses numerically calculated lookup tables. # # """ # scale_fac = np.sqrt(total_mass / self.tab_mass) * np.sqrt(self.tab_R_peak / R_peak) # vcirc = self.vcirc_interp(R / R_peak * self.tab_R_peak) * scale_fac # # # scale_fac = np.sqrt(total_mass / table_mass) * np.sqrt(table_Rpeak / R_peak) # # vcirc_interp = v_interp(Rarr / R_peak * table_Rpeak) * scale_fac # # return vcirc
[docs] def potential_gradient(self, R, R_peak, total_mass): """ Calculate potential gradient for a inf thin massive gaussian Ring Parameters ---------- R : float or array Radius or radii at which to calculate the potential gradient in kpc R_peak : float Peak of Gaussian ring in kpc total_mass : float Total mass of the Gaussian ring component Returns ------- potl_grad : float or array Potential gradient at each given `R` Notes ----- This function determines the potential gradient as a function of radius for a massive Gaussian ring component with a total mass, `total_mass`, and a ring peak radius to ring FWHM ratio, `invh`. This uses numerically calculated lookup tables. """ if self.potl_grad_interp is not None: scale_fac = total_mass / self.tab_mass * (self.tab_R_peak / R_peak)**2 potential_gradient_interp = self.potl_grad_interp(R / R_peak * self.tab_R_peak) * scale_fac return potential_gradient_interp else: return R*np.NaN
[docs] def enclosed_mass(self, R, R_peak, total_mass): """ Calculate enclosed mass for a inf thin massive gaussian Ring Parameters ---------- R : float or array Radius or radii at which to calculate the enclosed in kpc R_peak : float Peak of Gaussian ring in kpc total_mass : float Total mass of the Gaussian ring component Returns ------- vcirc : float or array Circular velocity at each given `R` Notes ----- This function determines the enclosed as a function of radius for a massive Gaussian ring component with a total mass, `total_mass`, and a ring peak radius to ring FWHM ratio, `invh`. This uses numerically calculated lookup tables. """ if self.menc_interp is not None: scale_fac = total_mass / self.tab_mass menc = self.menc_interp(R / R_peak * self.tab_R_peak) * scale_fac return menc else: return R*np.NaN
[docs] class BlackHole(MassModel): """ Central black hole. Treated as a point source at r = 0. Parameters ---------- BH_mass : float Log10 of the mass in solar units """ BH_mass = DysmalParameter(default=1, bounds=(0., 12.)) _subtype = 'baryonic' baryon_type = 'blackhole' def __init__(self, **kwargs): super(BlackHole, self).__init__(**kwargs)
[docs] @staticmethod def evaluate(r, BH_mass): """ Mass surface density of a BH (treat like delta function) """ # Ensure it's an array: if isinstance(r*1., float): rarr = np.array([r]) else: rarr = np.array(r) # Ensure all radii are 0. or positive: rarr = np.abs(rarr) mr = r * 0. wh0 = np.where((rarr == 0.))[0] mr[wh0] = BH_mass if (len(rarr) > 1): return mr else: if isinstance(r*1., float): # Float input return mr[0] else: # Length 1 array input return mr
[docs] def enclosed_mass(self, r): """ Central black hole enclosed mass (treat as step function) Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- menc : float or array Enclosed mass profile (unit: Msun) """ menc = r*0. + np.power(10.,self.BH_mass) return menc
[docs] def projected_enclosed_mass(self, r): # Point source: 2D is same as 3D return self.enclosed_mass(r)
[docs] def circular_velocity(self, r): """ Circular velocity as a function of radius Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- vcirc : float or array Circular velocity in km/s """ return super(BlackHole, self).circular_velocity(r)
[docs] def light_profile(self, r): """ Conversion from mass to light as a function of radius Assuming NO LIGHT emitted by central BH (eg, ignoring any emission in surrounding medium, eg flares) Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- light : float or array Relative line flux as a function of radius """ return r * 0.
[docs] class ExpDisk(MassModel, _LightMassModel): """ Infinitely thin exponential disk (i.e. Freeman disk) Parameters ---------- total_mass : float Log of total mass of the disk in solar units r_eff : float Effective radius in kpc baryon_type : {'gas+stars', 'stars', 'gas'} What type of baryons are included. Used for dlnrhogas/dlnr """ total_mass = DysmalParameter(default=1, bounds=(5, 14)) r_eff = DysmalParameter(default=1, bounds=(0, 50)) mass_to_light = DysmalParameter(default=1, fixed=True) _subtype = 'baryonic' tracer = 'mass' def __init__(self, baryon_type='gas+stars', **kwargs): self.baryon_type = baryon_type super(ExpDisk, self).__init__(**kwargs)
[docs] @staticmethod def evaluate(r, total_mass, r_eff, mass_to_light): """ Mass surface density of a thin exponential disk """ return surf_dens_exp_disk(r, 10.**total_mass, r_eff / 1.6783469900166612)
@property def rd(self): #b1 = 1.6783469900166612 # scp_spec.gammaincinv(2.*n, 0.5), n=1 return self.r_eff / 1.6783469900166612
[docs] def enclosed_mass(self, r): """ Enclosed mass as a function of radius Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- menc : float or array 1D enclosed mass profile """ return menc_exp_disk(r, 10**self.total_mass, self.rd)
[docs] def circular_velocity(self, r): """ Circular velocity as a function of radius Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- vcirc : float or array Circular velocity in km/s """ vcirc = vcirc_exp_disk(r, 10**self.total_mass, self.rd) return vcirc
[docs] def light_profile(self, r): """ Conversion from mass to light as a function of radius Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- light : float or array Relative line flux as a function of radius """ #light = surf_dens_exp_disk(r, 1.0, self.rd) light = surf_dens_exp_disk(r, (1./self.mass_to_light) * 10**self.total_mass, self.rd) return light
[docs] def rhogas(self, r): """ Mass surface density as a function of radius Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- surf_dens : float or array Mass surface density at `r` in units of Msun/kpc^2 """ return surf_dens_exp_disk(r, 10.**self.total_mass, self.rd)
[docs] def dlnrhogas_dlnr(self, r): """ Exponential disk asymmetric drift term Parameters ---------- r : float or array Radius in kpc Returns ------- log_drhodr : float or array Log surface density derivative as a function or radius Notes ----- See [1]_ for derivation and specificall Equations 3-11 References ---------- .. [1] https://ui.adsabs.harvard.edu/abs/2010ApJ...725.2324B/abstract """ # Shortcut for the exponential disk asymmetric drift term, from Burkert+10 eq 11: return -2. * (r / self.rd)
[docs] class Sersic(MassModel, _LightMassModel): """ Mass distribution following a Sersic profile Parameters ---------- total_mass : float Log10 of the total mass in solar units r_eff : float Effective (half-light) radius in kpc n : float Sersic index invq : float Ratio of the effective radius to the effective radius in the z-direction noord_flat : bool If True, use circular velocity profiles derived in Noordermeer 2008. If False, circular velocity is derived through `v_circular` baryon_type : {'gas+stars', 'stars', 'gas'} What type of baryons are included. Used for dlnrhogas/dlnr Notes ----- Model formula: .. math:: M(r) = M_e \exp \\left\{ -b_n \\left[ \\left( \\frac{r}{r_{\mathrm{eff}}} \\right)^{1/n} -1 \\right] \\right\} The constant :math:`b_n` is defined such that :math:`r_{\mathrm{eff}}` contains half the total mass, and can be solved for numerically. .. math:: \Gamma(2n) = 2\gamma (b_n,2n) Examples -------- .. plot:: :include-source: import numpy as np from dysmalpy.models import Sersic import matplotlib.pyplot as plt plt.figure() plt.subplot(111, xscale='log', yscale='log') s1 = Sersic(total_mass=10, r_eff=5, n=1) r=np.arange(0, 100, .01) for n in range(1, 10): s1.n = n plt.plot(r, s1(r), color=str(float(n) / 15)) plt.axis([1e-1, 30, 1e5, 1e10]) plt.xlabel('log Radius [kpc]') plt.ylabel('log Mass Surface Density [log Msun/kpc^2]') plt.text(.25, 8.e7, 'n=1') plt.text(.25, 3.e9, 'n=10') plt.show() """ total_mass = DysmalParameter(default=1, bounds=(5, 14)) r_eff = DysmalParameter(default=1, bounds=(0, 50)) n = DysmalParameter(default=1, bounds=(0, 8)) mass_to_light = DysmalParameter(default=1, fixed=True) _subtype = 'baryonic' tracer = 'mass' def __init__(self, invq=1.0, noord_flat=False, baryon_type='gas+stars', **kwargs): self.invq = invq self.baryon_type = baryon_type self._noord_flat = noord_flat super(Sersic, self).__init__(**kwargs) self._initialize_noord_flatteners() def __setstate__(self, state): super(Sersic, self).__setstate__(state) if 'baryon_type' in state.keys(): pass else: self.baryon_type = 'gas+stars' if '_noord_flat' in state.keys(): pass else: self._noord_flat = state['noord_flat'] self._initialize_noord_flatteners() @property def noord_flat(self): return self._noord_flat @noord_flat.setter def noord_flat(self, value): if type(value) is not bool: raise ValueError("noord_flat must be True/False!") self._noord_flat = value self._initialize_noord_flatteners() def _initialize_noord_flatteners(self): if self.noord_flat: # Initialize NoordFlat object: self.noord_flattener = NoordFlat(n=self.n.value, invq=self.invq) def _update_noord_flatteners(self): if self.n.value != self.noord_flattener._n: self.noord_flattener.n = self.n.value if self.invq != self.noord_flattener._invq: self.noord_flattener.invq = self.invq
[docs] @staticmethod def evaluate(r, total_mass, r_eff, n, mass_to_light): """ Sersic mass surface density """ return sersic_mr(r, 10**total_mass, n, r_eff)
[docs] def enclosed_mass(self, r): """ Sersic enclosed mass (linear) Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- menc : float or array Enclosed mass profile """ if self.noord_flat: # Check menc, invq are right: self._update_noord_flatteners() # Correct flattened mass profile: return self.noord_flattener.enclosed_mass(r, self.r_eff, 10**self.total_mass) else: return sersic_menc_2D_proj(r, 10**self.total_mass, self.n, self.r_eff)
[docs] def projected_enclosed_mass(self, r): return sersic_menc_2D_proj(r, 10**self.total_mass, self.n, self.r_eff)
[docs] def circular_velocity(self, r): """ Circular velocity as a function of radius Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- vcirc : float or array Circular velocity in km/s """ if self.noord_flat: # Check v, invq are right: self._update_noord_flatteners() vcirc = self.noord_flattener.circular_velocity(r, self.r_eff, 10**self.total_mass) else: vcirc = super(Sersic, self).circular_velocity(r) return vcirc
[docs] def light_profile(self, r): """ Conversion from mass to light as a function of radius Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- light : float or array Relative line flux as a function of radius """ #return sersic_mr(r, 1.0, self.n, self.r_eff) return sersic_mr(r, (1./self.mass_to_light) * 10**self.total_mass, self.n, self.r_eff)
[docs] def rhogas(self, r): """ Mass density as a function of radius (if noord_flat; otherwise surface density) Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- dens : float or array Mass density at `r` in units of Msun/kpc^3 (if noord_flat; otherwise surface density) """ if 'gas' in self.baryon_type.lower().strip(): if self.noord_flat: #rhogas = sersic_curve_rho(r, self.r_eff, 10**self.total_mass, self.n, self.invq) # Check v, invq are right: self._update_noord_flatteners() rhogas = self.noord_flattener.rho(r, self.r_eff, 10**self.total_mass) else: rhogas = sersic_mr(r, 10**self.total_mass, self.n, self.r_eff) else: rhogas = r * 0. return rhogas
[docs] def dlnrhogas_dlnr(self, r): """ Sersic asymmetric drift term Parameters ---------- r : float or array Radius in kpc Returns ------- log_drhodr : float or array Log surface density derivative as a function or radius """ if 'gas' in self.baryon_type.lower().strip(): if self.noord_flat: #dlnrhogas_dlnr_arr = sersic_curve_dlnrho_dlnr(r, self.r_eff, self.n, self.invq) # Check v, invq are right: self._update_noord_flatteners() dlnrhogas_dlnr_arr = self.noord_flattener.dlnrho_dlnr(r, self.r_eff) else: bn = scp_spec.gammaincinv(2. * self.n, 0.5) dlnrhogas_dlnr_arr = -2. * (bn / self.n) * np.power(r/self.r_eff, 1./self.n) else: dlnrhogas_dlnr_arr = r * 0. return dlnrhogas_dlnr_arr
[docs] class DiskBulge(MassModel, _LightMassModel): """ Mass distribution with a disk and bulge Parameters ---------- total_mass : float Log10 of the combined disk and bulge in solar units r_eff_disk : float Effective radius of the disk in kpc n_disk : float Sersic index of the disk r_eff_bulge : float Effective radius of the bulge n_bulge : float Sersic index of the bulge bt : float Bulge-to-total mass ratio invq_disk : float Effective radius to effective height ratio for the disk invq_bulge : float Effective radius to effective height ratio for the bulge noord_flat : bool If True, use circular velocity profiles derived in Noordermeer 2008. If False, circular velocity is derived through `v_circular` light_component : {'disk', 'bulge', 'total'} Which component to use as the flux profile gas_component : {'disk', 'total'} Which component contributes to dlnrhogas/dlnr baryon_type : {'gas+stars', 'stars', 'gas'} What type of baryons are included. Used for dlnrhogas/dlnr Notes ----- This model is the combination of 2 components, a disk and bulge, each described by a `Sersic`. The model is parametrized such that the B/T is a free parameter rather than the individual masses of the disk and bulge. """ total_mass = DysmalParameter(default=10, bounds=(5, 14)) r_eff_disk = DysmalParameter(default=1, bounds=(0, 50)) n_disk = DysmalParameter(default=1, fixed=True, bounds=(0, 8)) r_eff_bulge = DysmalParameter(default=1, bounds=(0, 50)) n_bulge = DysmalParameter(default=4., fixed=True, bounds=(0, 8)) bt = DysmalParameter(default=0.2, bounds=(0, 1)) mass_to_light = DysmalParameter(default=1, fixed=True) _subtype = 'baryonic' tracer = 'mass' def __init__(self, invq_disk=5, invq_bulge=1, noord_flat=False, light_component='disk', gas_component='disk', baryon_type='gas+stars', **kwargs): self.invq_disk = invq_disk self.invq_bulge = invq_bulge self.light_component = light_component self.gas_component = gas_component self.baryon_type = baryon_type self._noord_flat = noord_flat super(DiskBulge, self).__init__(**kwargs) self._initialize_noord_flatteners() def __setstate__(self, state): state_mod = copy.deepcopy(state) if 'noord_flat' in state.keys(): del state_mod['noord_flat'] state_mod['_noord_flat'] = state['noord_flat'] super(DiskBulge, self).__setstate__(state_mod) if 'baryon_type' in state_mod.keys(): pass else: self.baryon_type = 'gas+stars' self.gas_component = 'disk' if 'noord_flat' in state.keys(): self._initialize_noord_flatteners() @property def noord_flat(self): return self._noord_flat @noord_flat.setter def noord_flat(self, value): if type(value) is not bool: raise ValueError("noord_flat must be True/False!") self._noord_flat = value self._initialize_noord_flatteners() def _initialize_noord_flatteners(self): # Initialize NoordFlat objects: self.noord_flattener_disk = NoordFlat(n=self.n_disk.value, invq=self.invq_disk) self.noord_flattener_bulge = NoordFlat(n=self.n_bulge.value, invq=self.invq_bulge) def _update_noord_flatteners(self): if self.n_disk.value != self.noord_flattener_disk._n: self.noord_flattener_disk.n = self.n_disk.value if self.n_bulge.value != self.noord_flattener_bulge._n: self.noord_flattener_bulge.n = self.n_bulge.value if self.invq_disk != self.noord_flattener_disk._invq: self.noord_flattener_disk.invq = self.invq_disk if self.invq_bulge != self.noord_flattener_bulge._invq: self.noord_flattener_bulge.invq = self.invq_bulge
[docs] @staticmethod def evaluate(r, total_mass, r_eff_disk, n_disk, r_eff_bulge, n_bulge, bt, mass_to_light): """Disk+Bulge mass surface density""" print("consider if Noord flat: this will be modified") mbulge_total = 10**total_mass*bt mdisk_total = 10**total_mass*(1 - bt) mr_bulge = sersic_mr(r, mbulge_total, n_bulge, r_eff_bulge) mr_disk = sersic_mr(r, mdisk_total, n_disk, r_eff_disk) return mr_bulge+mr_disk
[docs] def enclosed_mass(self, r): """ Disk+Bulge total enclosed mass Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- menc : float or array Enclosed mass profile """ mbulge_total = 10 ** self.total_mass * self.bt mdisk_total = 10 ** self.total_mass * (1 - self.bt) if self.noord_flat: # Check menc, invq are right: self._update_noord_flatteners() # Correct flattened mass profiles: menc_bulge = self.noord_flattener_bulge.enclosed_mass(r, self.r_eff_bulge, mbulge_total) menc_disk = self.noord_flattener_disk.enclosed_mass(r, self.r_eff_disk, mdisk_total) else: # 2D projected: menc_bulge = sersic_menc_2D_proj(r, mbulge_total, self.n_bulge, self.r_eff_bulge) menc_disk = sersic_menc_2D_proj(r, mdisk_total, self.n_disk, self.r_eff_disk) return menc_disk+menc_bulge
[docs] def enclosed_mass_disk(self, r): """ Enclosed mass of the disk component Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- menc : float or array Enclosed mass profile """ mdisk_total = 10 ** self.total_mass * (1 - self.bt) if self.noord_flat: # Check menc, invq are right: self._update_noord_flatteners() # Correct flattened mass profiles: menc_disk = self.noord_flattener_disk.enclosed_mass(r, self.r_eff_disk, mdisk_total) else: # 2D projected: menc_disk = sersic_menc_2D_proj(r, mdisk_total, self.n_disk, self.r_eff_disk) return menc_disk
[docs] def enclosed_mass_bulge(self, r): """ Enclosed mass of the bulge component Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- menc : float or array Enclosed mass profile """ mbulge_total = 10 ** self.total_mass * self.bt if self.noord_flat: # Check menc, invq are right: self._update_noord_flatteners() # Correct flattened mass profiles: menc_bulge = self.noord_flattener_bulge.enclosed_mass(r, self.r_eff_bulge, mbulge_total) else: # 2D projected: menc_bulge = sersic_menc_2D_proj(r, mbulge_total, self.n_bulge, self.r_eff_bulge) return menc_bulge
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[docs] def projected_enclosed_mass(self, r): menc_disk = self.projected_enclosed_mass_disk(r) menc_bulge = self.projected_enclosed_mass_bulge(r) return menc_disk + menc_bulge
[docs] def projected_enclosed_mass_disk(self, r): mdisk_total = 10 ** self.total_mass * (1 - self.bt) return sersic_menc_2D_proj(r, mdisk_total, self.n_disk, self.r_eff_disk)
[docs] def projected_enclosed_mass_bulge(self, r): mbulge_total = 10 ** self.total_mass * self.bt return sersic_menc_2D_proj(r, mbulge_total, self.n_bulge, self.r_eff_bulge)
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[docs] def circular_velocity_disk(self, r): """ Circular velocity of the disk as a function of radius Parameters ---------- r : float or array Radii at which to calculate the circular velocity Returns ------- vcirc : float or array Circular velocity in km/s """ if self.noord_flat: mdisk_total = 10**self.total_mass*(1-self.bt) # Check v, invq are right: self._update_noord_flatteners() vcirc = self.noord_flattener_disk.circular_velocity(r, self.r_eff_disk, mdisk_total) else: mass_enc = self.enclosed_mass_disk(r) vcirc = v_circular(mass_enc, r) return vcirc
[docs] def circular_velocity_bulge(self, r): """ Circular velocity of the bulge as a function of radius Parameters ---------- r : float or array Radii at which to calculate the circular velocity Returns ------- vcirc : float or array Circular velocity in km/s """ if self.noord_flat: mbulge_total = 10**self.total_mass*self.bt # Check v, invq are right: self._update_noord_flatteners() vcirc = self.noord_flattener_bulge.circular_velocity(r, self.r_eff_bulge, mbulge_total) else: mass_enc = self.enclosed_mass_bulge(r) vcirc = v_circular(mass_enc, r) return vcirc
[docs] def circular_velocity(self, r): """ Total circular velocity as a function of radius Parameters ---------- r : float or array Radii at which to calculate the circular velocity in kpc Returns ------- vcirc : float or array Circular velocity in km/s """ vbulge = self.circular_velocity_bulge(r) vdisk = self.circular_velocity_disk(r) vcirc = np.sqrt(vbulge**2 + vdisk**2) return vcirc
[docs] def velocity_profile(self, r, modelset): """ Total rotational velocity due to the disk+bulge Parameters ---------- r : float or array Radii at which to calculate the circular velocity in kpc modelset : `ModelSet` Full ModelSet this component belongs to Returns ------- vrot : float or array Rotational velocity in km/s Notes ----- This method requires a `ModelSet` input to be able to apply the pressure support correction due to the gas turbulence. """ vcirc_sq = self.vcirc_sq(r) vrot_sq = modelset.kinematic_options.apply_pressure_support(r, modelset, vcirc_sq) vrot = np.sqrt(vrot_sq) return vrot
[docs] def velocity_profile_disk(self, r, modelset): """ Rotational velocity due to the disk Parameters ---------- r : float or array Radii at which to calculate the circular velocity in kpc modelset : `ModelSet` Full ModelSet this component belongs to Returns ------- vrot : float or array Rotational velocity in km/s Notes ----- This method requires a `ModelSet` input to be able to apply the pressure support correction due to the gas turbulence. """ vcirc_sq = self.circular_velocity_disk(r) ** 2 vrot_sq = modelset.kinematic_options.apply_pressure_support(r, modelset, vcirc_sq) vrot = np.sqrt(vrot_sq) return vrot
[docs] def velocity_profile_bulge(self, r, modelset): """ Rotational velocity due to the bulge Parameters ---------- r : float or array Radii at which to calculate the circular velocity in kpc modelset : `ModelSet` Full ModelSet this component belongs to Returns ------- vrot : float or array Rotational velocity in km/s Notes ----- This method requires a `ModelSet` input to be able to apply the pressure support correction due to the gas turbulence. """ vcirc_sq = self.circular_velocity_bulge(r) ** 2 vrot_sq = modelset.kinematic_options.apply_pressure_support(r, modelset, vcirc_sq) vrot = np.sqrt(vrot_sq) return vrot
[docs] def light_profile(self, r): """ Conversion from mass to light as a function of radius Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- light : float or array Relative line flux as a function of radius Notes ----- The resulting light profile depends on what `DiskBulge.light_component` is set to. If 'disk' or 'bulge' then only the mass associated with the disk or bulge will be converted into light. If 'total', then both components will be used. """ if self.light_component == 'disk': #flux = sersic_mr(r, 1.0, self.n_disk, self.r_eff_disk) flux = sersic_mr(r, (1./self.mass_to_light) * (1.0-self.bt) * 10**self.total_mass, self.n_disk, self.r_eff_disk) elif self.light_component == 'bulge': #flux = sersic_mr(r, 1.0, self.n_bulge, self.r_eff_bulge) flux = sersic_mr(r, (1./self.mass_to_light) * self.bt * 10**self.total_mass, self.n_bulge, self.r_eff_bulge) elif self.light_component == 'total': # flux_disk = sersic_mr(r, 1.0-self.bt, # self.n_disk, self.r_eff_disk) # flux_bulge = sersic_mr(r, self.bt, # self.n_bulge, self.r_eff_bulge) flux_disk = sersic_mr(r, (1./self.mass_to_light) * (1.0-self.bt) * 10**self.total_mass, self.n_disk, self.r_eff_disk) flux_bulge = sersic_mr(r, (1./self.mass_to_light) * self.bt * 10**self.total_mass, self.n_bulge, self.r_eff_bulge) flux = flux_disk + flux_bulge else: raise ValueError("light_component can only be 'disk', 'bulge', " "or 'total.'") return flux
[docs] def rhogas_disk(self, r): """ Mass density of the disk as a function of radius (if noord_flat; otherwise surface density) Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- dens : float or array Mass density at `r` in units of Msun/kpc^3 (if noord_flat; otherwise surface density) """ if 'gas' in self.baryon_type.lower().strip(): if self.gas_component in ['total', 'disk']: if self.noord_flat: mdisk_total = 10**self.total_mass*(1 - self.bt) # rhogas = sersic_curve_rho(r, self.r_eff_disk, mdisk_total, # self.n_disk, self.invq_disk) # Check v, invq are right: self._update_noord_flatteners() rhogas = self.noord_flattener_disk.rho(r, self.r_eff_disk, mdisk_total) else: mdisk_total = 10**self.total_mass*(1 - self.bt) # Just use the surface density as "rho", as this is the razor-thin case rhogas = sersic_mr(r, mdisk_total, self.n_disk, self.r_eff_disk) else: rhogas = r * 0. else: rhogas = r * 0. return rhogas
[docs] def rhogas_bulge(self, r): """ Mass density of the bulge as a function of radius (if noord_flat; otherwise surface density) Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- dens : float or array Mass density at `r` in units of Msun/kpc^3 (if noord_flat; otherwise surface density) """ if 'gas' in self.baryon_type.lower().strip(): # Only include bas in bulge if gas_component is 'total': if self.gas_component in ['total']: if self.noord_flat: mbulge_total = 10**self.total_mass*self.bt # rhogas = sersic_curve_rho(r, self.r_eff_bulge, mbulge_total, # self.n_bulge, self.invq_bulge) # Check v, invq are right: self._update_noord_flatteners() rhogas = self.noord_flattener_bulge.rho(r, self.r_eff_bulge, mbulge_total) else: mbulge_total = 10**self.total_mass*self.bt # Just use the surface density as "rho", as this is the razor-thin case rhogas = sersic_mr(r, mbulge_total, self.n_bulge, self.r_eff_bulge) else: rhogas = r * 0. else: rhogas = r * 0. return rhogas
[docs] def rhogas(self, r): """ Mass density as a function of radius (if noord_flat; otherwise surface density) Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- dens : float or array Mass density at `r` in units of Msun/kpc^3 (if noord_flat; otherwise surface density) """ # All cases handled internally in rhogas_disk, rhogas_bulge rhogas = self.rhogas_disk(r) + self.rhogas_bulge(r) return rhogas
[docs] def dlnrhogas_dlnr_disk(self, r): if self.noord_flat: #dlnrhogas_dlnr_arr = sersic_curve_dlnrho_dlnr(r, self.r_eff_disk, self.n_disk, self.invq_disk) # Check v, invq are right: self._update_noord_flatteners() dlnrhogas_dlnr_arr = self.noord_flattener_disk.dlnrho_dlnr(r, self.r_eff_disk) return dlnrhogas_dlnr_arr else: bn = scp_spec.gammaincinv(2. * self.n_disk, 0.5) return -2. * (bn / self.n_disk) * np.power(r/self.r_eff_disk, 1./self.n_disk)
[docs] def dlnrhogas_dlnr_bulge(self, r): if self.noord_flat: #dlnrhogas_dlnr_arr = sersic_curve_dlnrho_dlnr(r, self.r_eff_bulge, self.n_bulge, self.invq_bulge) # Check v, invq are right: self._update_noord_flatteners() dlnrhogas_dlnr_arr = self.noord_flattener_bulge.dlnrho_dlnr(r, self.r_eff_bulge) return dlnrhogas_dlnr_arr else: bn = scp_spec.gammaincinv(2. * self.n_bulge, 0.5) return -2. * (bn / self.n_bulge) * np.power(r/self.r_eff_bulge, 1./self.n_bulge)
[docs] def dlnrhogas_dlnr(self, r): """ Asymmetric drift term for the combined disk and bulge Parameters ---------- r : float or array Radius in kpc Returns ------- log_drhodr : float or array Log surface density derivative as a function or radius """ if 'gas' in self.baryon_type.lower().strip(): if self.gas_component == 'total': rhogasD = self.rhogas_disk(r) rhogasB = self.rhogas_bulge(r) dlnrhogas_dlnr_tot = (1./(rhogasD + rhogasB)) * \ (rhogasD*self.dlnrhogas_dlnr_disk(r) + rhogasB*self.dlnrhogas_dlnr_bulge(r)) elif self.gas_component == 'disk': dlnrhogas_dlnr_tot = self.dlnrhogas_dlnr_disk(r) else: dlnrhogas_dlnr_tot = r * 0. return dlnrhogas_dlnr_tot
[docs] class LinearDiskBulge(MassModel, _LightMassModel): """ Mass distribution with a disk and bulge Parameters ---------- total_mass : float Combined disk and bulge mass in solar units r_eff_disk : float Effective radius of the disk in kpc n_disk : float Sersic index of the disk r_eff_bulge : float Effective radius of the bulge n_bulge : float Sersic index of the bulge bt : float Bulge-to-total mass ratio invq_disk : float Effective radius to effective height ratio for the disk invq_bulge : float Effective radius to effective height ratio for the bulge noord_flat : bool If True, use circular velocity profiles derived in Noordermeer 2008. If False, circular velocity is derived through `v_circular` light_component : {'disk', 'bulge', 'total'} Which component to use as the flux profile baryon_type : {'gas+stars', 'stars', 'gas'} What type of baryons are included. Used for dlnrhogas/dlnr Notes ----- This model is the exactly the same as `DiskBulge` except that `total_mass` is in linear units instead of log. """ total_mass = DysmalParameter(default=10, bounds=(5, 14)) r_eff_disk = DysmalParameter(default=1, bounds=(0, 50)) n_disk = DysmalParameter(default=1, fixed=True, bounds=(0, 8)) r_eff_bulge = DysmalParameter(default=1, bounds=(0, 50)) n_bulge = DysmalParameter(default=4., fixed=True, bounds=(0, 8)) bt = DysmalParameter(default=0.2, bounds=(0, 1)) mass_to_light = DysmalParameter(default=1, fixed=True) _subtype = 'baryonic' tracer = 'mass' def __init__(self, invq_disk=5, invq_bulge=1, noord_flat=False, light_component='disk', baryon_type='gas+stars', **kwargs): self.invq_disk = invq_disk self.invq_bulge = invq_bulge self.light_component = light_component self.baryon_type = baryon_type self._noord_flat = noord_flat super(LinearDiskBulge, self).__init__(**kwargs) self._initialize_noord_flatteners() def __setstate__(self, state): super(LinearDiskBulge, self).__setstate__(state) if 'baryon_type' in state.keys(): pass else: self.baryon_type = 'gas+stars' if '_noord_flat' in state.keys(): pass else: self._noord_flat = state['noord_flat'] self._initialize_noord_flatteners() @property def noord_flat(self): return self._noord_flat @noord_flat.setter def noord_flat(self, value): if type(value) is not bool: raise ValueError("noord_flat must be True/False!") self._noord_flat = value self._initialize_noord_flatteners() def _initialize_noord_flatteners(self): # Initialize NoordFlat objects: self.noord_flattener_disk = NoordFlat(n=self.n_disk.value, invq=self.invq_disk) self.noord_flattener_bulge = NoordFlat(n=self.n_bulge.value, invq=self.invq_bulge) def _update_noord_flatteners(self): if self.n_disk.value != self.noord_flattener_disk._n: self.noord_flattener_disk.n = self.n_disk.value if self.n_bulge.value != self.noord_flattener_bulge._n: self.noord_flattener_bulge.n = self.n_bulge.value if self.invq_disk != self.noord_flattener_disk._invq: self.noord_flattener_disk.invq = self.invq_disk if self.invq_bulge != self.noord_flattener_bulge._invq: self.noord_flattener_bulge.invq = self.invq_bulge
[docs] @staticmethod def evaluate(r, total_mass, r_eff_disk, n_disk, r_eff_bulge, n_bulge, bt, mass_to_light): """Disk+Bulge mass surface density""" print("consider if Noord flat: this will be modified") mbulge_total = total_mass*bt mdisk_total = total_mass*(1 - bt) mr_bulge = sersic_mr(r, mbulge_total, n_bulge, r_eff_bulge) mr_disk = sersic_mr(r, mdisk_total, n_disk, r_eff_disk) return mr_bulge+mr_disk
[docs] def enclosed_mass(self, r): """ Disk+Bulge total enclosed mass Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- menc : float or array Enclosed mass profile """ mbulge_total = self.total_mass * self.bt mdisk_total = self.total_mass * (1 - self.bt) if self.noord_flat: # Check menc, invq are right: self._update_noord_flatteners() # Correct flattened mass profiles: menc_bulge = self.noord_flattener_bulge.enclosed_mass(r, self.r_eff_bulge, mbulge_total) menc_disk = self.noord_flattener_disk.enclosed_mass(r, self.r_eff_disk, mdisk_total) else: # 2D projected: menc_bulge = sersic_menc_2D_proj(r, mbulge_total, self.n_bulge, self.r_eff_bulge) menc_disk = sersic_menc_2D_proj(r, mdisk_total, self.n_disk, self.r_eff_disk) return menc_disk+menc_bulge
[docs] def enclosed_mass_disk(self, r): """ Enclosed mass of the disk component Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- menc : float or array Enclosed mass profile """ mdisk_total = self.total_mass * (1 - self.bt) if self.noord_flat: # Check menc, invq are right: self._update_noord_flatteners() # Correct flattened mass profiles: menc_disk = self.noord_flattener_disk.enclosed_mass(r, self.r_eff_disk, mdisk_total) else: # 2D projected: menc_disk = sersic_menc_2D_proj(r, mdisk_total, self.n_disk, self.r_eff_disk) return menc_disk
[docs] def enclosed_mass_bulge(self, r): """ Enclosed mass of the bulge component Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- menc : float or array Enclosed mass profile """ mbulge_total = self.total_mass * self.bt if self.noord_flat: # Check menc, invq are right: self._update_noord_flatteners() # Correct flattened mass profiles: menc_bulge = self.noord_flattener_bulge.enclosed_mass(r, self.r_eff_bulge, mbulge_total) else: # 2D projected: menc_bulge = sersic_menc_2D_proj(r, mbulge_total, self.n_bulge, self.r_eff_bulge) return menc_bulge
[docs] def circular_velocity_disk(self, r): """ Circular velocity of the disk as a function of radius Parameters ---------- r : float or array Radii at which to calculate the circular velocity Returns ------- vcirc : float or array Circular velocity in km/s """ if self.noord_flat: mdisk_total = self.total_mass*(1-self.bt) # Check v, invq are right: self._update_noord_flatteners() vcirc = self.noord_flattener_disk.circular_velocity(r, self.r_eff_disk, mdisk_total) else: mass_enc = self.enclosed_mass_disk(r) vcirc = v_circular(mass_enc, r) return vcirc
[docs] def circular_velocity_bulge(self, r): """ Circular velocity of the bulge as a function of radius Parameters ---------- r : float or array Radii at which to calculate the circular velocity Returns ------- vcirc : float or array Circular velocity in km/s """ if self.noord_flat: mbulge_total = self.total_mass*self.bt # Check v, invq are right: self._update_noord_flatteners() vcirc = self.noord_flattener_bulge.circular_velocity(r, self.r_eff_bulge, mbulge_total) else: mass_enc = self.enclosed_mass_bulge(r) vcirc = v_circular(mass_enc, r) return vcirc
[docs] def circular_velocity(self, r): """ Total Circular velocity as a function of radius Parameters ---------- r : float or array Radii at which to calculate the circular velocity Returns ------- vcirc : float or array Circular velocity in km/s """ vbulge = self.circular_velocity_bulge(r) vdisk = self.circular_velocity_disk(r) vcirc = np.sqrt(vbulge**2 + vdisk**2) return vcirc
[docs] def velocity_profile(self, r, modelset): """ Total rotational velocity due to the disk+bulge Parameters ---------- r : float or array Radii at which to calculate the circular velocity in kpc modelset : `ModelSet` Full ModelSet this component belongs to Returns ------- vrot : float or array Rotational velocity in km/s Notes ----- This method requires a `ModelSet` input to be able to apply the pressure support correction due to the gas turbulence. """ vcirc_sq = self.vcirc_sq(r) vrot_sq = modelset.kinematic_options.apply_pressure_support(r, modelset, vcirc_sq) vrot = np.sqrt(vrot_sq) return vrot
[docs] def velocity_profile_disk(self, r, modelset): """ Rotational velocity due to the disk Parameters ---------- r : float or array Radii at which to calculate the circular velocity in kpc modelset : `ModelSet` Full ModelSet this component belongs to Returns ------- vrot : float or array Rotational velocity in km/s Notes ----- This method requires a `ModelSet` input to be able to apply the pressure support correction due to the gas turbulence. """ vcirc_sq = self.circular_velocity_disk(r) ** 2 vrot_sq = modelset.kinematic_options.apply_pressure_support(r, modelset, vcirc_sq) vrot = np.sqrt(vrot_sq) return vrot
[docs] def velocity_profile_bulge(self, r, modelset): """ Rotational velocity due to the bulge Parameters ---------- r : float or array Radii at which to calculate the circular velocity in kpc modelset : `ModelSet` Full ModelSet this component belongs to Returns ------- vrot : float or array Rotational velocity in km/s Notes ----- This method requires a `ModelSet` input to be able to apply the pressure support correction due to the gas turbulence. """ vcirc_sq = self.circular_velocity_bulge(r) ** 2 vrot_sq = modelset.kinematic_options.apply_pressure_support(r, modelset, vcirc_sq) vrot = np.sqrt(vrot_sq) return vrot
[docs] def light_profile(self, r): """ Conversion from mass to light as a function of radius Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- light : float or array Relative line flux as a function of radius Notes ----- The resulting light profile depends on what `DiskBulge.light_component` is set to. If 'disk' or 'bulge' then only the mass associated with the disk or bulge will be converted into light. If 'total', then both components will be used. """ if self.light_component == 'disk': #flux = sersic_mr(r, 1.0, self.n_disk, self.r_eff_disk) flux = sersic_mr(r, (1./self.mass_to_light) * (1.0-self.bt) * self.total_mass, self.n_disk, self.r_eff_disk) elif self.light_component == 'bulge': #flux = sersic_mr(r, 1.0, self.n_bulge, self.r_eff_bulge) flux = sersic_mr(r, (1./self.mass_to_light) * self.bt * self.total_mass, self.n_bulge, self.r_eff_bulge) elif self.light_component == 'total': # flux_disk = sersic_mr(r, 1.0-self.bt, # self.n_disk, self.r_eff_disk) # flux_bulge = sersic_mr(r, self.bt, # self.n_bulge, self.r_eff_bulge) flux_disk = sersic_mr(r, (1./self.mass_to_light) * (1.0-self.bt) * self.total_mass, self.n_disk, self.r_eff_disk) flux_bulge = sersic_mr(r, (1./self.mass_to_light) * self.bt * self.total_mass, self.n_bulge, self.r_eff_bulge) flux = flux_disk + flux_bulge else: raise ValueError("light_component can only be 'disk', 'bulge', " "or 'total.'") return flux
[docs] class GaussianRing(MassModel, _LightMassModel): r""" Mass distribution following an infinitely thin Gaussian ring profile. Parameters ---------- total_mass : float Log10 of the total mass in solar units R_peak : float Peak of gaussian (radius) in kpc FWHM: float FWHM of gaussian, in kpc baryon_type : string What type of baryons are included. Used for dlnrhogas/dlnr. Options: {'gas+stars', 'stars', 'gas'} Notes ----- Model formula: .. math:: M(r)&=M_0\exp\left(\frac{(r-r_{\rm peak})^2}{2\sigma_R^2}\right) \sigma_R &= \mathrm{FWHM}/(2\sqrt{2\ln 2}) """ total_mass = DysmalParameter(default=1, bounds=(5, 14)) R_peak = DysmalParameter(default=1, bounds=(0, 50)) FWHM = DysmalParameter(default=1, bounds=(0, 50)) mass_to_light = DysmalParameter(default=1, fixed=True) _subtype = 'baryonic' tracer = 'mass' def __init__(self, baryon_type='gas+stars', **kwargs): self.baryon_type = baryon_type super(GaussianRing, self).__init__(**kwargs) self._initialize_ring_table() def __setstate__(self, state): super(GaussianRing, self).__setstate__(state) if 'baryon_type' in state.keys(): pass else: self.baryon_type = 'gas+stars' if 'ring_table' in state.keys(): pass else: self._initialize_ring_table() def _initialize_ring_table(self): # Initialize InfThinMassiveGaussianRing object: self.ring_table = InfThinMassiveGaussianRing(invh=self.ring_invh()) def _update_ring_table(self): if self.ring_invh() != self.ring_table._invh: self.ring_table.invh = self.ring_invh()
[docs] def sigma_R(self): return self.FWHM.value / (2.*np.sqrt(2.*np.log(2.)))
[docs] def ring_invh(self): return self.R_peak.value / self.FWHM.value
[docs] def ring_reff(self): # Find the effective radius explicitly by definition, using erf function # and in the dimensionless parameter x = r/Rpeak; with subbing u = 4*ln(2) h^2 (x-1)^2 # Solving for: # int_0^xeff t * exp{-4*ln(2)*invh**2*(t-1)**2) = 0.5 * int_0^inf t * exp{-4*ln(2)*invh**2*(t-1)**2) try: A = 4*np.log(2) func = lambda ueff: scp_spec.erf(np.sqrt(ueff)) - (np.pi * A * self.ring_invh() ** 2) ** -0.5 * np.exp(-ueff) - \ 0.5 * (1 - scp_spec.erf(np.sqrt(A * self.ring_invh() ** 2)) - (np.pi * A * self.ring_invh() ** 2) ** -0.5 * np.exp(-A * self.ring_invh() ** 2)) u_eff = scp_opt.brentq(lambda t: func(t), a=0., b=1.) x_eff = 1 + np.sqrt(u_eff / (A * self.ring_invh() ** 2)) return x_eff * self.R_peak.value except: logger.warning("Could not find ring_reff. Assuming reff=Rpeak instead...") return self.R_peak.value
[docs] def r_eff(self): return self.ring_reff()
[docs] @staticmethod def evaluate(r, total_mass, R_peak, FWHM, mass_to_light): """ Gaussian ring mass surface density """ sigma_R = FWHM/ (2.*np.sqrt(2.*np.log(2.))) I0 = _I0_gaussring(R_peak, sigma_R, 10**total_mass) return I0*np.exp(-(r-R_peak)**2/(2.*sigma_R**2))
[docs] def surface_density(self, r): """ Gaussian ring mass surface density """ return self.evaluate(r, self.total_mass.value, self.R_peak.value, self.FWHM.value, 1.)
[docs] def enclosed_mass(self, r): """ Sersic enclosed mass Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- menc : float or array Enclosed mass profile """ # Check invh is correct self._update_ring_table() return self.ring_table.enclosed_mass(r, self.R_peak.value, 10**self.total_mass.value)
[docs] def projected_enclosed_mass(self, r): """ Same as enclosed mass as this is infinitely thin gaussian ring """ return self.enclosed_mass(r)
[docs] def circular_velocity(self, r): """ Circular velocity as a function of radius Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- vcirc : float or array Circular velocity in km/s """ return np.sqrt(self.vcirc_sq(r))
[docs] def vcirc_sq(self, r): """ Square of circular velocity as a function of radius Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- vcirc_sq : float or array Square of circular velocity in km^2/s^2 Notes ----- Calculated as :math:`v_{\mathrm{circ}}^2(R) = R * \partial \Phi / \partial R` from the gradient of the potential, as the potential gradient has negative values. """ return r * self.potential_gradient(r)
[docs] def potential_gradient(self, r): r""" Method to evaluate the gradient of the potential, :math:`\Delta\Phi(r)/\Delta r`. Parameters ---------- r : float or array Radius or radii at which to calculate circular velocity in kpc Returns ------- dPhidr : float or array Gradient of the potential at `r` """ # Check invh is correct self._update_ring_table() return self.ring_table.potential_gradient(r, self.R_peak.value, 10**self.total_mass.value)
[docs] def light_profile(self, r): """ Conversion from mass to light as a function of radius Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- light : float or array Relative line flux as a function of radius """ #return 1.*np.exp(-(r-self.R_peak.value)**2/(2.*self.sigma_R()**2)) return self.surface_density(r) * (1./self.mass_to_light)
[docs] def rhogas(self, r): """ Mass density as a function of radius (if noord_flat; otherwise surface density) Parameters ---------- r : float or array Radii at which to calculate the enclosed mass Returns ------- dens : float or array Mass density at `r` in units of Msun/kpc^3 (if noord_flat; otherwise surface density) """ if 'gas' in self.baryon_type.lower().strip(): # Inf thin: rho(R,z=0) = Sigma(R) rhogas = self.surface_density(r) else: rhogas = r * 0. return rhogas
[docs] def dlnrhogas_dlnr(self, r): """ Sersic asymmetric drift term Parameters ---------- r : float or array Radius in kpc Returns ------- log_drhodr : float or array Log surface density derivative as a function or radius """ if 'gas' in self.baryon_type.lower().strip(): # Inf thin: rho(R,z=0) = Sigma(R) dlnrhogas_dlnr_arr = - r * (r-self.R_peak.value) / (self.sigma_R()**2) else: dlnrhogas_dlnr_arr = r * 0. return dlnrhogas_dlnr_arr