Example Dysmalpy model
This notebook shows how to create a mock galaxy with user-specific settings and the visualization of the models
Includes the following components:
Disk + Bulge
NFW halo
Constant velocity dispersion
The structure of the notebook is the following:
Setup steps
Initialize galaxy and model set
Set up observations
Create models
1) Setup steps
Import modules
from __future__ import (absolute_import, division, print_function, unicode_literals)
import numpy as np
import astropy.units as u
import astropy.io.fits as fits
import copy
from dysmalpy import galaxy, models, instrument, parameters, \
plotting, observation, aperture_classes
INFO:numexpr.utils:Note: NumExpr detected 10 cores but "NUMEXPR_MAX_THREADS" not set, so enforcing safe limit of 8.
INFO:numexpr.utils:NumExpr defaulting to 8 threads.
Setup notebook
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline
mpl.rcParams['figure.dpi']= 300
mpl.rc("savefig", dpi=300)
Import functions to tie scale height relative to effective radius, and the halo mass to the \(\boldsymbol{f_\mathrm{DM}}\) value
from dysmalpy.fitting_wrappers.tied_functions import tie_sigz_reff, tie_lmvirial_to_fdm, \
tied_geom_lambda
2) Initialize galaxy & model set
gal = galaxy.Galaxy(z=2., name='galaxy')
mod_set = models.ModelSet()
Baryonic component: Combined Disk+Bulge
total_mass = 10.5 # M_sun
bt = 0.3 # Bulge-Total ratio
r_eff_disk = 4.0 # kpc
n_disk = 1.0
invq_disk = 5.0
r_eff_bulge = 1.0 # kpc
n_bulge = 4.0
invq_bulge = 1.0
noord_flat = True # Switch for applying Noordermeer flattening
# Fix components
bary_fixed = {'total_mass': False,
'r_eff_disk': False,
'n_disk': True,
'r_eff_bulge': True,
'n_bulge': True,
'bt': False}
# Set bounds
bary_bounds = {'total_mass': (10, 13),
'r_eff_disk': (1.0, 30.0),
'n_disk': (1, 8),
'r_eff_bulge': (1, 5),
'n_bulge': (1, 8),
'bt': (0, 1)}
bary = models.DiskBulge(total_mass=total_mass, bt=bt,
r_eff_disk=r_eff_disk, n_disk=n_disk,
invq_disk=invq_disk,
r_eff_bulge=r_eff_bulge, n_bulge=n_bulge,
invq_bulge=invq_bulge,
noord_flat=noord_flat,
name='disk+bulge',
fixed=bary_fixed, bounds=bary_bounds)
bary.r_eff_disk.prior = parameters.BoundedGaussianPrior(center=5.0, stddev=1.0)
Halo component
mvirial = -99
conc = 5.0
fdm = 0.5
halo_fixed = {'mvirial': False,
'conc': True,
'fdm': False}
halo_bounds = {'mvirial': (10, 13),
'conc': (1, 20),
'fdm': (0, 1)}
halo = models.NFW(mvirial=mvirial, conc=conc, fdm=fdm, z=gal.z,
fixed=halo_fixed, bounds=halo_bounds, name='halo')
halo.mvirial.tied = tie_lmvirial_to_fdm
# The tied component must have "fixed=False", and then a specified tied function
Dispersion profile
sigma0 = 39. # km/s
disp_fixed = {'sigma0': False}
disp_bounds = {'sigma0': (5, 300)}
disp_prof = models.DispersionConst(sigma0=sigma0, fixed=disp_fixed,
bounds=disp_bounds, name='dispprof', tracer='line')
z-height profile
sigmaz = 0.9 # kpc
zheight_fixed = {'sigmaz': False}
zheight_prof = models.ZHeightGauss(sigmaz=sigmaz, name='zheightgaus',
fixed=zheight_fixed)
zheight_prof.sigmaz.tied = tie_sigz_reff
# The tied component must have "fixed=False", and then a specified tied function
Geometry
Use the same settings for the 1D and 2D mocks:
inc = 62. # degrees
pa = 142. # degrees, blue-shifted side CCW from north
xshift = 0 # pixels from center
yshift = 0 # pixels from center
geom_fixed = {'inc': False,
'pa': True,
'xshift': True,
'yshift': True}
geom_bounds = {'inc': (0, 90),
'pa': (90, 180),
'xshift': (0, 4),
'yshift': (-10, -4)}
geom = models.Geometry(inc=inc, pa=pa, xshift=xshift, yshift=yshift,
fixed=geom_fixed, bounds=geom_bounds, name='geom',
obs_name='mock')
geom_1D = models.Geometry(inc=inc, pa=pa, xshift=xshift, yshift=yshift,
fixed=geom_fixed, bounds=geom_bounds, name='geom_1D',
obs_name='mock_1D')
# Tie the inc / PA to the 2D geom:
geom_1D.inc.tied = tied_geom_lambda('geom', 'inc')
geom_1D.pa.tied = tied_geom_lambda('geom', 'pa')
Add all model components to ModelSet
mod_set.add_component(bary, light=True)
mod_set.add_component(halo)
mod_set.add_component(disp_prof)
mod_set.add_component(zheight_prof)
mod_set.add_component(geom)
mod_set.add_component(geom_1D)
Set kinematic options for calculating velocity profile
mod_set.kinematic_options.adiabatic_contract = False
mod_set.kinematic_options.pressure_support = True
3) Set up observations
Set up the 1D + 2D observation & instruments:
obs = observation.Observation(name='mock', tracer='line')
inst = instrument.Instrument()
beamsize = 0.55*u.arcsec # FWHM of beam
sig_inst = 45*u.km/u.s # Instrumental spectral resolution
beam = instrument.GaussianBeam(major=beamsize)
lsf = instrument.LSF(sig_inst)
inst.ndim = 2
inst.beam = beam
inst.lsf = lsf
inst.pixscale = 0.125*u.arcsec # arcsec/pixel
inst.fov = [33, 33] # (nx, ny) pixels
inst.spec_type = 'velocity' # 'velocity' or 'wavelength'
inst.spec_step = 10*u.km/u.s # Spectral step
inst.spec_start = -1000*u.km/u.s # Starting value of spectrum
inst.nspec = 201 # Number of spectral pixels
# Set the beam kernel so it doesn't have to be calculated every step
inst.set_beam_kernel()
inst.set_lsf_kernel()
obs_1D = observation.Observation(name='mock_1D', tracer='line')
inst_1D = copy.deepcopy(inst)
inst_1D.ndim = 1
Add the model set, observations to the Galaxy
gal.model = mod_set
obs.instrument = inst
gal.add_observation(obs)
obs_1D.instrument = inst_1D
obs_1D_options = observation.ObsModOptions(oversample=3)
obs_1D.mod_options = obs_1D_options
# Define apertures:
aper_arr = np.linspace(-(inst_1D.fov[0]-1)/2., (inst_1D.fov[0]-1)/2.,
num=inst_1D.fov[0])*inst_1D.pixscale.value
apertures = aperture_classes.setup_aperture_types(obs=obs_1D, profile1d_type='circ_ap_cube',
slit_width = 0.55, aper_centers=aper_arr, slit_pa=142.)
obs_1D.instrument.apertures = apertures
gal.add_observation(obs_1D)
4) Create models
Creates models for all observations – including both the 1D and 2D cases.
gal.create_model_data()
# Alternatively: can create models on a per-obs basis, using the obs name:
# gal.create_model_data(obs_list=['mock'])
# gal.create_model_data(obs_list=['mock_1D'])
2D model
plotting.plot_model_2D(gal, inst_corr=True)
1D model
plotting.plot_model_1D(gal, inst_corr=True, best_dispersion=sigma0)
