# coding=utf-8
import numpy as np
from scipy import interpolate
try:
import bottleneck as bn
except ImportError:
import numpy as bn
from .hyperspectral_cube import HyperspectralCube as HyperCube
from .galaxy_parameters import GalaxyParameters
from .model_utilities import modelPlots, ModelExt
# CONSTANTS
G = 6.67384e-11 # Gravitational Constant (m^3.kg^-1.s^-2)
SOL_MASS = 2e30 # Solar Mass (kg)
PARSEC = 3e16 # Parsec (m)
# LOGGING CONFIGURATION
import logging
logging.basicConfig(level=logging.INFO)
[docs]class DiskModel(modelPlots, ModelExt):
"""
The first and default galpak `Model` (when unspecified).
It simulates a simple disk galaxy using `DiskParameters`.
flux_profile: 'exponential' | 'gaussian' | 'de_vaucouleurs' | 'sersicN2' | 'user'
The flux profile of the observed galaxy. The default is 'exponential'.
See http://en.wikipedia.org/wiki/Sersic%27s_law
If 'user' provided flux map, specify flux_image. Currently unsupported
thickness_profile: 'gaussian' [default] | 'exponential' | 'sech2' | 'none'
The perpendicular flux profile of the disk galaxy.
The default is 'gaussian'.
If none, the galaxy will be a cylinder
rotation_curve: 'arctan' | 'mass' | 'exp' | 'tanh' | 'isothermal'
The profile of the velocity v(r) can be in
- 'arctan' (default) : ~ Vmax arctan(r/rV), rV=turnover radius
- 'tanh' : Vmax tanh(r/rV), rV=turnover radius
- 'exp' : inverted exponential, 1-exp(-r/rV)
- 'mass' : a constant light-to-mass ratio v(r)=sqrt(G m(<r) / r)
- 'isothermal' : an isothermal rotation curve v(r)=Vmax * sqrt(1-rV/r * arctan(r/rV)
dispersion_profile: 'thick' [default] | 'thin'
The local disk dispersion from the rotation curve and disk thickness
from Binney & Tremaine 2008, see Genzel et al. 2008.
GalPak has 3 components for the dispersion:
- a component from the rotation curve arising from mixing velocities
of a disk with non-zero thickness.
- a component from the local disk dispersion specified by
`disk_dispersion`
- a spatially constant dispersion,
which is the output parameter `velocity_dispersion`.
line: None[default] | dict to fit doublets, use a dictionary with
line['wave']=[l, lref] eg. [3726.2, 3728.9] # Observed or Rest
line['ratio']=[0.8, 1.0] eg. [0.8, 1.0] # The primary line for redshifts is the reddest
redshift: float
The redshift of the observed galaxy, used in mass calculus.
Will override the redshift provided at init ; this is a convenience parameter.
If this is not set anywhere, GalPak will not try to compute the mass.
"""
# fixme: use integers instead of strings here and enjoy the performance gain
FLUX_VALID = ['gaussian', 'exponential', 'de_vaucouleurs', 'sersicN2', 'user']
THICKNESS_VALID = [ 'gaussian', 'exponential', 'sech2', 'none']
DISPERSION_VALID = [ 'thick', 'thin']
CURVE_VALID = ['arctan', 'mass', 'exponential', 'tanh', 'isothermal' ]
def __init__(self,
flux_profile='exponential',
thickness_profile='gaussian',
dispersion_profile='thick',
rotation_curve= 'arctan',
flux_image=None,
line=None,
redshift=None,
cosmology='planck15'):
self.flux_profile = flux_profile
self.thickness_profile = thickness_profile
self.dispersion_profile = dispersion_profile
self.rotation_curve = rotation_curve
self.line = line
self.redshift = redshift
self.cosmology= cosmology
self.logger = logging.getLogger('GalPaK: DiskModel:')
#if using input image
if self.flux_profile == 'user':
self.flux_map_user = flux_image
if self.redshift is not None and self.cosmology is not None:
self.set_cosmology()
self.Parameters = GalaxyParameters #just a copy of Parameter class
#@fixme to be removed
self.hz_profile = thickness_profile #for backwards compatibility
self.disk_dispersion = dispersion_profile #for backwards compatibility
self.GalaxyParameters = GalaxyParameters #for backwards compatibility
#From v2.0beta
def set_cosmology(self):
try:
from colossus.cosmology import cosmology as cosmo
from colossus import halo
from colossus.halo import mass_defs, mass_adv, mass_so, concentration
self.logger.info("setting cosmology with colossus package")
cosmo=cosmo.setCosmology(self.cosmology)
#put h in
h = cosmo.h
self.Ez = cosmo.Ez(self.redshift)
DA = cosmo.angularDiameterDistance(self.redshift) / h ####!!
self.DeltaVir = halo.mass_so.deltaVir(self.redshift)
self.halo = halo #Diemer Halo module for NFW halo models
except ImportError:
try:
from astLib import astCalc
#using astLib for cosmology
self.logger.info("Using AstLib for a 737 cosmology, h,OM, OL= (0.7, 0.3, 0.7)")
h = astCalc.H0 / 100.
self.Ez = astCalc.Ez(self.redshift)
DA = astCalc.da(self.redshift)#for h=0.7
self.DeltaVir = astCalc.DeltaVz(self.redshift)
self.halo = None
self.logger.warning('Best to use colossus package from Diemer& Kratsov for this module')
except ImportError:
#using astropy
from astropy import cosmology
self.logger.info("Using astropy.cosmology WMAP9 only for cosmology")
cosmo = cosmology.default_cosmology.get()
h = cosmo.h
DA = cosmo.angular_diameter_distance(self.redshift)
self.Ez = cosmo.H(self.redshift).value/cosmo.H0
self.DeltaVir = None### not available
self.halo = None
self.logger.warning('Best to use colossus package from Diemer& Kratsov for this module')
self.kpc = DA * 1e3 * np.radians(1./3600)# in unit kpc
#self.hkpc = self.kpc * h
def initial_parameters(self, runner):
"""
Returns an instance of a class extending `ModelParameters`.
You may omit parameters, they will be automatically set to the mean of
the max and min.
"""
# Default initial parameters
return GalaxyParameters(
radius=3.0,
turnover_radius=1.0,
inclination=30,
pa=np.random.uniform(0,360),
maximum_velocity=100.0,
velocity_dispersion=5.0,
flux=runner.flux_est
)
def min_boundaries(self, runner):
"""
Returns an instance of a class extending `ModelParameters`.
You MUST provide all the parameters.
"""
# Default boundaries
return GalaxyParameters(
x=runner.cube.shape[2] / 3.,
y=runner.cube.shape[1] / 3.,
z=runner.cube.shape[0] / 3.,
flux=0.,
radius=0.199,
inclination=-0.5,
pa=0,
turnover_radius=0.01,
maximum_velocity=0.,
velocity_dispersion=0.
)
def max_boundaries(self, runner):
"""
Returns an instance of a class extending `ModelParameters`.
You MUST provide all the parameters.
"""
# changed max radius to 3/8th of cube size.
return GalaxyParameters(
x=runner.cube.shape[2] * 2 / 3.,
y=runner.cube.shape[1] * 2 / 3.,
z=runner.cube.shape[0] * 2 / 3.,
flux=3. * runner.flux_est,
radius= 0.5 * (runner.cube.shape[1]+runner.cube.shape[2]) * 3 / 8.,
inclination=90.5,
pa=360,
turnover_radius=0.5 * (runner.cube.shape[1]+runner.cube.shape[2]) * 3 / 8. / 2.,
maximum_velocity=350.,
velocity_dispersion=180.
)
def sanitize_parameters(self, parameters):
"""
Mutates the parameters with custom logic.
This is called right after parameter jumping in the loop.
"""
inc = parameters['inclination']
if inc<0:
parameters.pa +=180
parameters.pa = parameters.pa % 360. #to deal with circularity
def sanitize_chain(self, chain):
"""
Mutates the `chain` with custom logic.
This is called at the end of the loop, before storing the chain.
"""
pa = chain['pa']
pa_corr = np.where(pa > 180, pa - 360, pa)
chain['pa'] = pa_corr
def setup_random_amplitude(self, amplitude):
"""
Mutates the random `amplitude` of the parameter jump with custom logic.
This is called during the setup of the MCMC runner.
The amplitude is already filled with :
sqrt((min_boundaries - max_boundaries) ** 2 / 12.) * p / v
where
p = number of parameters in the model
v = number of voxels in the cube
"""
amplitude['flux'] *= 0.5 # flux since we start at ~flux since flux is an integrated quantity
amplitude['pa'] *= 1.5
amplitude['maximum_velocity'] *= 1.5 #
amplitude['velocity_dispersion'] *= 1.5 # velocity dispersion
amplitude *= 3 # fixme: document this coefficient
#CUSTOM METHODS
def _set_flux_profile(self, galaxy, radius_cube):
"""
creates a flux_cube in 3D (x,y,z)
returns flux_cube
"""
rad = galaxy['radius']
if self.flux_profile == 'de_vaucouleurs':
flux_cube = self.flux_de_vaucouleurs(radius_cube, rad)
elif self.flux_profile == 'gaussian':
flux_cube = self.flux_gaussian(radius_cube, rad)
elif self.flux_profile == 'exponential':
flux_cube = self.flux_exponential(radius_cube, rad)
elif self.flux_profile == 'sersicN2':
flux_cube = self.flux_sersic(radius_cube, rad, 2.0)
elif self.flux_profile == 'user':
self.logger.error("Not implemented")
#should be face-on view...
flux_cube = self.flux_map_user[None,:,:] #@fixme need thickness?
else:
raise NotImplementedError("Flux profile not supported. Should be '%s' " %
self.model.FLUX_VALID)
return flux_cube
def _set_velocity_profile(self, galaxy, radius_cube):
"""
creates a velocity profile for circular orbits
return v_profile
"""
max_vel = galaxy['maximum_velocity']
rt = galaxy['turnover_radius']
if self.rotation_curve == 'arctan':
v_profile = np.arctan(radius_cube / rt) * max_vel * 2. / np.pi
# For Mass=Light profile
elif self.rotation_curve == 'mass':
v_profile = self._fv_newton(radius_cube, galaxy.radius, max_vel)
# For Feng's inv_exp profile
elif self.rotation_curve == 'exponential':
v_profile = max_vel * (1.0 - np.exp(-radius_cube / rt))
elif self.rotation_curve == 'tanh' :
v_profile = max_vel * np.tanh(radius_cube / rt)
elif self.rotation_curve == 'isothermal' :
v_profile = max_vel * np.sqrt(1 - rt/radius_cube*np.arctan(radius_cube/rt))
else:
raise NotImplementedError("Rotational Curve not supported. Should be '%s' " %
self.model.CURVE_VALID)
return v_profile
# DEFAULT METHODS
def _set_thickness_profile(self, nz, hz, flux_cube):
"""
add the thickness profile to the flux cube
returns flux_cube
"""
if self.thickness_profile == 'gaussian':
flux_cube *= np.exp(-nz**2. / 2. / hz**2.)
elif self.thickness_profile == 'exponential':
flux_cube *= np.exp(-np.sqrt(nz**2.) / hz) # exponential z-profile
elif self.thickness_profile == 'sech2':
flux_cube *= np.cosh(-nz / hz)**(-2.) # sech^2 z-profile
elif self.thickness_profile == 'none':
self.logger.warning("Using no disk thickness profile: "
"this galaxy will become a cylinder!")
flux_cube = flux_cube
else:
raise NotImplementedError("Context profile not supported. Should be '%s'" %
self.model.THICKNESS_VALID)
return flux_cube
def _set_dispersion_profile(self, hz, vtot, radius_cube):
"""
set disk dispersion profile
"""
### ISOTROPIC but this IS the 1-d dispersion in z-direction
# A) using R/h=0.5 (V/sigma)^2 for very thin disk
if self.dispersion_profile == 'thin':
# dimension of r [3D-spatial]
sig_ale_disk = np.sqrt(hz * vtot ** 2. / radius_cube)
# B) using R/h=V/Sigma for compact disk
elif self.dispersion_profile == 'thick':
# dimension of r [3D-spatial]
sig_ale_disk = hz * vtot / radius_cube
# C) Using Binney Merrifield
# elif self.dispersion_profile == 'infinitely_thin':
# dimension of r [3D-spatial]
# sig_ale_disk = np.sqrt(math.tau * Surf_density * G * hz_kpc)
else:
raise ValueError("Disk dispersion is not valid. Should be '%s' " %
self.DISPERSION_VALID)
return sig_ale_disk
###########################################
#HEART of model
###########################################
def _create_maps(self, galaxy, shape):
"""
Creates and returns:
- a cube sized like (`shape`) containing a clean simulation
of a galaxy for the `parameters`.
- the velocity map, which is an image.
- the dispersion (sigma) map, which is an image.
This is called on each loop, to compute the likelihood of each generated
set of parameters. Therefore, it needs to be fast.
parameters: DiskParameters
The parameters upon which the simulated disk galaxy will be built.
shape: Tuple of 3
The 3D shape of the resulting cube.
Eg: (21,21,21)
Returns a tuple.
"""
# Create radial velocities and radial dispersions
flux_cube, vz, vz_map, s_map, sig_map, sigz_disk_map, sig_intr = \
self._compute_galaxy_model(galaxy, shape)
# normalize to real flux
flux_map = flux_cube.sum(0)
flux_map = galaxy.flux * flux_map / flux_map.sum()
return flux_map, vz_map, s_map
def _create_cube(self, shape, f_map, v_map, s_map, z_step_kms, zo=None):
"""
Creates a 3D cube from Flux, Velocity and Dispersion maps.
galaxy: GalaxyParameters
shape: tuple
Defines the created cube's shape
f_map: ndarray
2D image of the flux
v_map: ndarray
2D image of the velocity
s_map: ndarray
2D image of the dispersion
z_step_kms: float
in km/s
zo: int
The centroid (in pixel) in the z-direction
If not set, will be the middle of the z-axis
"""
vmax = z_step_kms * (shape[0] - 1) / 2.
if zo is None:
vzero = vmax
else:
vzero = z_step_kms * zo # fixme: shouldn't this be with zo-0.5 ?
# Array of velocities : array of indices from to 0 to Vmax
velocities = np.linspace(0, 2 * vmax, shape[0]) - vzero
# Create a copy of Vmap in each z plane
tmpcube = np.outer(np.ones_like(velocities), v_map)
vcube = tmpcube.reshape(shape)
# Create a copy of Smap in each z plane
tmpcube = np.outer(np.ones_like(velocities), s_map)
scube = tmpcube.reshape(shape)
# Create a copy of velocities in each location
if z_step_kms > 0:
vbegin = np.min(velocities)
vfinal = np.max(velocities)
else:
vbegin = np.max(velocities)
vfinal = np.min(velocities)
# Create grid of velocities
vgrid, y, x = np.mgrid[
vbegin: vfinal + 0.9 * z_step_kms: z_step_kms,
0: np.shape(v_map)[0],
0: np.shape(v_map)[1]
]
# Apply Gaussian
finalcube = np.exp(-0.5 * (vgrid - vcube) ** 2 / (scube ** 2))
# For doublets, add the blue component
# deltaV is defined as (l2-l1)/l2 * c ie from the reddest line
if self.line is not None:
# Delta between emission peaks (in km/s ?)
# line['wave']=[l, lref] eg. [3726.2, 3728.9] # Observed or Rest;
# line['ratio']=[0.8, 1.0] eg. [0.8, 1.0] ; the primary line for redshifts is the reddest
delta = 3e5 * (self.line['wave'][1]-self.line['wave'][0])/(self.line['wave'][0]+self.line['wave'][1])*2
ratio = self.line['ratio'][0] / self.line['ratio'][1]
# Add the blue doublet
finalcube = finalcube + ratio * np.exp(-0.5 * (vgrid - (vcube - delta)) ** 2 / (scube ** 2))
# delta = 217. # km/s OII
# delta = -945 # km/s Ha 6564.614 NII 6585.27
# delta = 2872 # km/s OIII 4960.295 5008.240
# delta = 641 # km/s SII 6718.29 6732.67
# Normalize from amplitude ratios int f dv = int f dz
ampsum = bn.nansum(finalcube,0)
#amplitude = f_map / bn.nansum(finalcube, 0)
#amplitude[np.isnan(amplitude)] = 0
amplitude = f_map / ampsum
amplitude[np.isfinite(amplitude)==False] = 0
#amplitude = np.where(np.isfinite(amplitude),f_map / ampsum, 0)
# Create a cube of Amplitude in each z plane
amplitude_cube = np.resize(amplitude, shape)
finalcube = finalcube * amplitude_cube
finalcube[np.isfinite(finalcube)==False] = 0
cube = HyperCube(data=finalcube)
return cube
def _compute_galaxy_model(self, galaxy, shape):
"""
Create a disk model in a cube
galaxy: GalaxyParameters
The parameters of the galaxy
Returns
3D flux cube,
3D Vz cube,
2D Vmap (1st moment of Vz),
2D S_map,
2D (2nd moment of Vz),
2D disk_dispersion,
2D intrinsic
"""
# Return indices in source frame; vz in image frame
hz = galaxy.radius * 0.15
# trick to go faster, as we know the max z shape the resulting model will fit in
# z_shape = np.ceil(shape[1] / 2.5 * 0.15 * 6.) ## 6 sigmas
# z_shape = np.ceil(self.max_boundaries.radius / 2.5 * 0.15 * 6.) ## 6 sigmas
# or
# z_shape = np.ceil(self.max_boundaries.radius * 2 * 3 * 0.15) # 3 hz_max on each side
# this short cut introduced a bug in the flux_deconvolved maps! fixed after 1.6.0
z_shape = int(np.ceil ( (shape[1] + shape[2]) / 2.)) #mean of spatial dimensions
self.logger.debug(" using z_shape %d in compute_galaxy_model" % (z_shape) )
self.logger.debug(" cube is %d in z- " % (shape[0]) )
# Minimum shape for z direction -- collapsed, so big sizes unimportant
# z_shape = 6 hz(sigma)_max;
# hz_max = size_box/2.5 (proxy for radius max) * 0.15
# Create arrays of indices in image frame; this cube has spatial dimentions in x,y,z; so should be almost a perfect cube in shape.
ind_shape = tuple([z_shape, shape[1], shape[2]])
z, y, x = np.indices(ind_shape)
#adding 0.5 to be centered on pixel center
nx, ny, nz, vz_cube, vtot = self._map_indices(x, y, z, hz, galaxy)
radius_cube = np.hypot(nx,ny)
flux_cube = self._set_flux_profile(galaxy, radius_cube)
flux_cube = self._set_thickness_profile(nz, hz, flux_cube)
# Flux weighted
flux_map = bn.nansum(flux_cube, 0)
vz_map = bn.nansum(vz_cube * flux_cube / flux_map, 0)
v2_map = bn.nansum(vz_cube ** 2 * flux_cube / flux_map, 0)
Var_map = v2_map - vz_map ** 2
sig_map_disk = np.sqrt(Var_map)
### Add Sigma_Aleatoire_disk
sig_ale_disk = self._set_dispersion_profile(hz, vtot, radius_cube)
# Compute weighted mean for Sig_disk
# a) 2D map = flux-weighted map
#Sigz_disk_map = bn.nansum(Sig_disk*A,0)/Flux_map
# b) 2D map = flux-weighted map rms
#sigz_map_ale = np.sqrt(bn.nansum(sig_ale_disk**2*A,0)/Flux_map)
# c) Normalized with flux squared
sigz_map_ale = np.sqrt(bn.nansum(sig_ale_disk ** 2 * flux_cube ** 2, 0) / bn.nansum(flux_cube * flux_cube, 0))
# Intrinsic Dispersion
sig_intr = np.ones_like(sigz_map_ale) * galaxy.sigma0
# Total S-map
s_map = np.sqrt(sig_map_disk ** 2 + sigz_map_ale ** 2 + sig_intr ** 2)
# Normalize by flux
total = bn.nansum(flux_cube)
flux_cube = np.where(total>0, flux_cube / total, 0)
return flux_cube, vz_cube, vz_map, s_map, \
sig_map_disk, sigz_map_ale, sig_intr
def _map_indices(self, xx, yy, zz, hz, galaxy ):
"""
Takes arrays of indices, galaxy parameters, and rotates accordingly.
Returns xx, yy, zz indices of object in image frame
Returns z-component of velocity of object in image frame
"""
inclination = np.radians(galaxy.inclination)
pa = np.radians(galaxy.pa - 90)
zo = (np.size(zz, 0) - 1) / 2. # center along z
# Rotation inclination (around x)
rot_x = lambda i: np.array([[1, 0, 0], [0, np.cos(i), -np.sin(i)], [0, np.sin(i), np.cos(i)]])
rot_i = rot_x(inclination)
rot_mi = rot_x(-inclination)
# Rotation PA (around z) anti-clockwise
rot_z = lambda a: np.array([[np.cos(a), np.sin(a), 0], [-np.sin(a), np.cos(a), 0], [0, 0, 1]])
rot_pa = rot_z(pa)
rot_mpa = rot_z(-pa)
# Transformation 1 from sky coord to disk plane
xx = xx - galaxy.x
yy = yy - galaxy.y
zz = zz - zo
# # Rotation around z for PA
# nx = rot_pa[0, 0] * xx + rot_pa[0, 1] * yy
# ny = rot_pa[1, 0] * xx + rot_pa[1, 1] * yy
# nz = rot_pa[2, 2] * zz
# # Rotation around x for inclination
# x = rot_i[0, 0] * nx
# y = rot_i[1, 1] * ny + rot_i[1, 2] * nz
# z = rot_i[2, 1] * ny + rot_i[2, 2] * nz
coord = np.array([xx, yy, zz])
rot = np.linalg.multi_dot([rot_i, rot_pa, coord.reshape(3, -1)])
x, y, z = rot.reshape(coord.shape)
# Compute the radius cube
# radius_cube, _ = self._set_flux_profile(galaxy, x, y, z)
radius_cube = np.hypot(x, y)
# Circular orbits -- we're DIVIDING BY ZERO ; it's okay, need for speed
vx = y / radius_cube
vy = - x / radius_cube
vx[~np.isfinite(vx)] = 0
vy[~np.isfinite(vy)] = 0
v_profile = self._set_velocity_profile(galaxy, radius_cube)
# Remove NaNs where radius is 0
v_profile[radius_cube == 0] = 0.
# 3D velocity vectors
vx *= v_profile
vy *= v_profile
# # Rotation of -inclination around x
# vvx = rot_mi[0, 0] * vx
# vvy = rot_mi[1, 1] * vy
# vvz = rot_mi[2, 1] * vy
# # Rotation of -pa around z
# vx = rot_mpa[0, 0] * vvx + rot_mpa[0, 1] * vvy
# vy = rot_mpa[1, 0] * vvx + rot_mpa[1, 1] * vvy
# vz = rot_mpa[2, 2] * vvz
# v = np.sqrt(vx ** 2 + vy ** 2 + vz ** 2)
coord = np.array([vx, vy, vy]).reshape(3, -1)
coord *= rot_mi[[0,1,2], [0,1,1]][:, None]
rot = np.dot(rot_mpa, coord)
vz = rot[2].reshape(vx.shape)
v = np.sqrt(np.sum(rot**2, axis=0))
v = v.reshape(vx.shape)
return x, y, z, vz, v
####################################
def _fv_newton(self, r, radius, vmax):
"""
Returns velocity given the r-cube
The cumulative flux I(<r) profile is analytically calculated given the profile shape
"""
mass = 1 # enclosed_mass for Gaussian
if self.__dict__['flux_profile'] == 'de_vaucouleurs':
#(-4 (5040 + 5040 x^(1/4) + 2520 Sqrt[x] + 840 x^(3/4) + 210 x + 42 x^(5/4) + 7 x^(3/2) + x^(7/4)))/E^x^(1/4)
int_x_dev = lambda x: (-4. * (5040. + 5040. * x ** (1. / 4) + 2520. * np.sqrt(x) + 840. * x ** (3. / 4) + 210. * x + 42. * x ** (5. / 4) + 7. * x ** (3. / 2) + x ** (7. / 4))) * np.exp(-x ** (1. / 4))
mtot = mass * (int_x_dev(r / radius) - int_x_dev(0))
elif self.__dict__['flux_profile'] == 'gaussian':
mtot = mass * (1. - np.exp(-r ** 2 / 2. / (2 * radius / 2.35) ** 2))
elif self.__dict__['flux_profile'] == 'exponential':
#was wrong:
#mtot = mass * (1. - ( r / radius) * np.exp(- r / radius) - np.exp(- r / radius))
#should be:
mtot = mass * (1. - (1.68 * r / radius) * np.exp(-1.68 * r / radius) - np.exp(-1.68 * r / radius))
else:
raise NotImplementedError("Flux profile not supported. Should be '%s' " %
self.model.FLUX_VALID)
v = np.sqrt(mtot / (r + 1e-9)) # km/s
# Normalization set by Vmax (sets the mass)
if np.size(v) > 1:
v = v / bn.nanmax(v) * vmax
else:
v = v / np.max(v) * vmax
return v
def _compute_Mdyn_at_Rhalf(self, galaxy, radius=None): #, flux_profile, rotation_curve):
"""
From galaxy, computes V profile and dynamic_mass
radius: float
In pixels, the half-light radius
"""
rv, max_vel = galaxy['turnover_radius'], galaxy['maximum_velocity']
if radius is None:
radius = galaxy['radius']
try:
from astropy.cosmology import WMAP5
except ImportError:
raise ImportError('The package `astropy` is required.')
if self.redshift is None:
raise RuntimeError("Cannot compute mass if redshift is not provided.")
redshift = self.redshift
pixscale = self.pixscale
# Compute angular distance in Mpc
dang = WMAP5.angular_diameter_distance(redshift).value
# Convert it to kpc/arcsec
dang *= 1e3 * np.radians(1. / 3600)
rhalf_kpc = radius * pixscale * dang # in kpc
r_kpc = np.r_[0:200:100j] # in kpc
r = r_kpc / dang # in arcsec
r_pix = r / pixscale # in pixel
if self.__dict__['rotation_curve'] == 'arctan':
fv = lambda q: np.arctan(q / rv) * max_vel * 2. / np.pi # r in pix
elif self.__dict__['rotation_curve'] == 'mass':
s = r_pix
fv = lambda q: self._fv_newton(q, s, max_vel) # r in pix
elif self.__dict__['rotation_curve'] == 'exponential': # Feng inv_exp profile
fv = lambda q: max_vel * (1.0 - np.exp(-q / rv))
elif self.__dict__['rotation_curve'] == 'tanh' :
fv = lambda q: max_vel * np.tanh(q / rv)
elif self.__dict__['rotation_curve'] == 'isothermal':
fv = lambda q: max_vel * np.sqrt(1-rv/r*np.arctan(r/rv))
else:
raise ValueError("Unsupported Rotational Curve. Should be '%s'" % self.model.CURVE_VALID)
# v² = (GM/r) in km/s
compute_mass = lambda vol, rad: (1e3 * vol) ** 2 * (PARSEC * rad * 1e3) / G / SOL_MASS
v_profile = fv(r_pix)
v_rhalf = interpolate.spline(r_kpc, v_profile, rhalf_kpc)
dynamic_mass = compute_mass(v_rhalf, rhalf_kpc)#computed at R1/2.
#mass_profile = compute_mass(v_profile, r_kpc)
#v_2rhalf = interpolate.spline(r_kpc, v_profile, 2. * rhalf_kpc)
#rvir_kpc = 180 # compute_Rvir(Vmax)
#Vvir = interpolate.spline(r_kpc, v_profile, rvir_kpc)
#Mvir = 0.83 * (max_vel / 200.) ** 3 * (1 + context['redshift']) ** 1.5 * 1e12
return dynamic_mass
### STATIC METHODS #########################################################
@staticmethod
def flux_gaussian(nx, ny, nz, rhwhm):
"""
Disk Gaussian profile with Gaussian thickness
"""
radius = np.sqrt(nx ** 2 + ny ** 2)
energy = np.exp(-radius ** 2 / 2. / (2. * rhwhm / 2.35) ** 2) #* np.exp(-nz ** 2 / 2. / hz ** 2)
return radius, energy
@staticmethod
def flux_sersic(radius, size, index):
"""
Sersic profile with Gaussian thickness
See http://en.wikipedia.org/wiki/Sersic_profile
"""
beta = 1.9992 * index - 0.3271
energy = np.exp(-beta * (radius / size) ** (1. / index))
#energy = np.exp(-beta * (radius / size) ** (1. / index)) * np.exp(-nz ** 2 / 2. / hz ** 2)
#energy = np.exp(-beta * (radius / size) ** (1. / index)) * np.exp(-np.sqrt(nz ** 2) / hz ) #exponential z-profile
#energy = np.exp(-beta * (radius / size) ** (1. / index)) * np.cosh(-nz / hz )**(-2.) #sech^2 z-profile
return energy
@staticmethod
def flux_exponential(radius, size):
"""
Disk exponential profile with Gaussian thickness
This is merely a Sérsic profile with index = 1
"""
return DiskModel.flux_sersic(radius, size, 1.)
@staticmethod
def flux_de_vaucouleurs(radius, size):
"""
De Vaucouleurs profile with Gaussian thickness
See http://en.wikipedia.org/wiki/De_Vaucouleurs%27_law
This is merely a Sérsic profile with index = 4
"""
return DiskModel.flux_sersic(radius, size, 4.)
def __str__(self):
dic = vars(self)
dic['model'] = self.__name__()
return str(dic)
def __name__(self):
return self.__class__.__name__
