Source code for galpak.lib.disk_model

# coding=utf-8

import numpy as np
from scipy import interpolate
from matplotlib import pyplot as plt

try:
    import bottleneck as bn
except ImportError:
    import numpy as bn

from .hyperspectral_cube import HyperspectralCube as HyperCube
from .galaxy_parameters import  GalaxyParameters

# 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)
logger = logging.getLogger('GalPaK: DiskModel:')


[docs]class DiskModel: """ 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'] 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', line=None, redshift=None): 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 #@fixme to be removed self.hz_profile = thickness_profile #for backwards compatibility self.disk_dispersion = dispersion_profile #for backwards compatibility 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, x, y, z): """ creates a flux_cube in 3D (x,y,z) returns flux_cube """ if self.__dict__['flux_profile'] == 'de_vaucouleurs': radius_cube, flux_cube = self.flux_de_vaucouleurs(x, y, z, galaxy.radius) elif self.__dict__['flux_profile'] == 'gaussian': radius_cube, flux_cube = self.flux_gaussian(x, y, z, galaxy.radius) elif self.__dict__['flux_profile'] == 'exponential': radius_cube, flux_cube = self.flux_exponential(x, y, z, galaxy.radius) elif self.__dict__['flux_profile'] == 'sersicN2': radius_cube, flux_cube = self.flux_sersic(x, y, z, galaxy.radius, 2.0) else: raise NotImplementedError("Flux profile not supported. Should be '%s' " % self.model.FLUX_VALID) return radius_cube, flux_cube def _set_velocity_profile(self, galaxy, radius_cube): """ creates a velocity profile for circular orbits return v_profile """ if self.__dict__['rotation_curve'] == 'arctan': v_profile = np.arctan(radius_cube / galaxy.rv) * galaxy.max_vel * 2. / np.pi # For Mass=Light profile elif self.__dict__['rotation_curve'] == 'mass': v_profile = self._fv_newton(radius_cube, galaxy.radius, galaxy.max_vel) # For Feng's inv_exp profile elif self.__dict__['rotation_curve'] == 'exponential': v_profile = galaxy.max_vel * (1.0 - np.exp(-radius_cube / galaxy.rv)) elif self.__dict__['rotation_curve'] == 'tanh' : v_profile = galaxy.max_vel * np.tanh(radius_cube / galaxy.rv) elif self.__dict__['rotation_curve'] == 'isothermal' : v_profile = galaxy.max_vel * np.sqrt(1-galaxy.rv/radius_cube*np.arctan(radius_cube/galaxy.rv)) 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.__dict__['thickness_profile'] == 'gaussian': flux_cube *= np.exp(-nz**2. / 2. / hz**2.) elif self.__dict__['thickness_profile'] == 'exponential': flux_cube *= np.exp(-np.sqrt(nz**2.) / hz) # exponential z-profile elif self.__dict__['thickness_profile'] == 'sech2': flux_cube *= np.cosh(-nz / hz)**(-2.) # sech^2 z-profile elif self.__dict__['thickness_profile'] == 'none': 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.__dict__['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.__dict__['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 logger.debug(" using z_shape %d in compute_galaxy_model" % (z_shape) ) 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, flux_cube = self._set_flux_profile(galaxy, nx, ny, nz) 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 # Compute the radius cube radius_cube, _ = self._set_flux_profile(galaxy, x, y, z) # Circular orbits -- we're DIVIDING BY ZERO ; it's okay, need for speed vx = np.where(np.isfinite(y / radius_cube), +y / radius_cube, 0) vy = np.where(np.isfinite(x / radius_cube), -x / radius_cube, 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 = vx * v_profile vy = 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) 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.diskdisk_model_dict_model['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, radius, rv, max_vel, pixscale): #, flux_profile, rotation_curve): """ From radius, rv, max_vel, computes V profile and dynamic_mass radius: float In pixels, the half-light 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.cube.xy_step # 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 ##################################### def plot_vprofile(self, parameters): max_vel = parameters.maximum_velocity from matplotlib import pyplot as plt x = np.r_[0.1:parameters.radius*5:0.1] v_x = self._set_velocity_profile(parameters, x) plt.plot(x,v_x,'k',label='V(r) ' + self.rotation_curve) sini = np.sin(np.radians(parameters.inclination)) plt.axhline(max_vel,label='Vmax') plt.ylabel('V(km/s) deproj.') plt.xlabel('x (pix)') plt.legend() plt.show() # def plot_AM(self, parameters): """ Compute and plot the Angular Momentum : cumulative specific (normalized and not), non cumulative specific (normalized and not) the normalization is (Vmax*Rh) for the cumulative one and (Vvir*Rvir) for the non cumulative one :param parameters: run parameters :return: AM spec and cumulative """ """ :param parameters: :return: """ max_vel = parameters.maximum_velocity rh = parameters.radius from matplotlib import pyplot as plt r = np.r_[0.1:rh*10:0.1] v_x = self._set_velocity_profile(parameters, r) j = v_x * r rr, flux = self._set_flux_profile(parameters,x=r/(2**.5),y=r/(2**.5),z=r/(2**.5)) # nx, ny = radius / sqrt(2) because r = sqrt(nx**2 + ny**2) in the function flux_sersic j_cum = np.cumsum(flux * v_x * r**2) / np.cumsum(flux * r) # AM specific j = v_x*r # AM #z = self.redshift #print z #M = flux #r_vir = 200*1e-3*(M/1e6)**(1./3.) * (10/(1+z)) # In kPc (M in Solar Mass)---------------------- ??? plt.figure(figsize=(15,10)) plt.title('Angular momentum') plt.subplot(2,2,1) plt.plot(r/rh,j_cum) plt.ylabel('$j_*(<r)$') plt.xlabel('$r/r_{1/2}$') plt.subplot(2,2,2) plt.plot(r/rh,j) plt.ylabel('$j_*(r)$') plt.xlabel('$r/r_{1/2}$') plt.subplot(2,2,3) plt.plot(r/rh,j_cum/(max_vel*rh)) plt.ylabel('$j_*(<r)\ /\ (V_{max}\ r_{1/2})$') plt.xlabel('$r/r_{1/2}$') plt.subplot(2,2,4) plt.plot(r/rh,j/(max_vel*rh)) plt.ylabel('$j_*(r)\ /\ (V_{max}\ r_{1/2})$') plt.xlabel('$r/r_{1/2}$') plt.subplots_adjust(hspace=.3) plt.tight_layout plt.show() ### 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(nx, ny, nz, size, index): """ Sersic profile with Gaussian thickness See http://en.wikipedia.org/wiki/Sersic_profile """ beta = 1.9992 * index - 0.3271 radius = np.sqrt(nx ** 2 + ny ** 2) 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 radius, energy @staticmethod def flux_exponential(nx, ny, nz, size): """ Disk exponential profile with Gaussian thickness This is merely a Sérsic profile with index = 1 """ return DiskModel.flux_sersic(nx, ny, nz, size, 1.) @staticmethod def flux_de_vaucouleurs(nx, ny, nz, 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(nx, ny, nz, size, 4.) def __str__(self): return str(vars(self)) def __dict__(self): return vars(self)