Quickstart Dysmalpy example

This code can be pasted within a python session to quickly test Dysmalpy model creation.

Alternatively, the script can be run with:

python dysmalpy_quickstart_example.py

The output 3D cube will be saved as a FITS file (named dpy_test_model_3D.fits) in your current working directory. The 1D and 2D plot files will be shown interactively (or otherwise also saved to the current directory, if interactive plotting is not available).

# Example `Dysmalpy` model: saves output to current directory

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

from dysmalpy import galaxy, observation, models, instrument, aperture_classes, \
                     parameters, plotting

import numpy as np
import astropy.units as u
import astropy.io.fits as fits

import matplotlib as mpl
import matplotlib.pyplot as plt

if mpl.get_backend() == 'agg':
    havedisplay=False
else:
    havedisplay=True

# **Set function to tie scale height relative to effective radius**
def tie_sigz_reff(model_set):
    reff = model_set.components['disk+bulge'].r_eff_disk.value
    invq = model_set.components['disk+bulge'].invq_disk
    sigz = 2.0*reff/invq/2.35482
    return sigz


# ---------------------------------------------------------------

## Initialize galaxy, model set, observations, and instruments
gal = galaxy.Galaxy(z=2., name='galaxy')
mod_set = models.ModelSet()
obs = observation.Observation(name='OBS', tracer='LINE')
inst = instrument.Instrument()

## Set observation options:
obs.mod_options.oversample = 3
obs.fit_options.fit_flux = True   # Also plot the flux.
                                  # Velocity/dispersion included by default.

### 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 = 12.0
conc = 5.0

halo_fixed = {'mvirial': False,
              'conc': True}

halo_bounds = {'mvirial': (10, 13),
               'conc': (1, 20)}

halo = models.NFW(mvirial=mvirial, conc=conc, z=gal.z,
                  fixed=halo_fixed, bounds=halo_bounds, name='halo')

halo.mvirial.prior = parameters.BoundedGaussianPrior(center=11.5, stddev=0.5)


### Dispersion profile -- note the tracer must match the observation's tracer:
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

### Geometry -- note the "obs_name" attribute must match the observation's name:
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='OBS')


### 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)


### Set kinematic options for calculating velocity profile
mod_set.kinematic_options.adiabatic_contract = False
mod_set.kinematic_options.pressure_support = True


### Set up the 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.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

inst.moment = False                      # Use Gaussian fitting to extract to 1D/2D

# Set the beam kernel so it doesn't have to be calculated every step
inst.set_beam_kernel()
inst.set_lsf_kernel()


### Add the instrument to the observation:
obs.instrument = inst

### Add the model set, observation to the Galaxy
gal.model = mod_set
gal.add_observation(obs)


## Create models
f_cube = 'dpy_test_model_3D.fits'
if havedisplay:
    fileout1D = fileout2D = None
else:
    fileout1D = "dpy_test_model_1D.pdf"
    fileout2D = "dpy_test_model_2D.pdf"

### 3D model
gal.observations['OBS'].instrument.ndim = 3  # Set ndim of model
gal.create_model_data()
gal.observations['OBS'].model_cube.data.write(f_cube, overwrite=True)


### 2D model
gal.observations['OBS'].instrument.ndim = 2  # Set ndim of model
gal.create_model_data()
plotting.plot_model_2D(gal, inst_corr=True, fileout_base=fileout2D)

### 1D model
gal.observations['OBS'].instrument.ndim = 1  # Set ndim of model
# Define apertures for 1D profile: slit PA=kin PA; slit width=PSF FWHM
aper_arr = np.linspace(-(inst.fov[0]-1)/2., (inst.fov[0]-1)/2.,
                           num=inst.fov[0])*inst.pixscale.value
apertures = aperture_classes.setup_aperture_types(obs=gal.observations['OBS'],
                                                  profile1d_type='circ_ap_cube',
                                                  aper_centers=aper_arr,
                                                  slit_width=0.55, slit_pa=142.)
gal.observations['OBS'].instrument.apertures = apertures

gal.create_model_data()
plotting.plot_model_1D(gal, inst_corr=True, best_dispersion=sigma0, fileout_base=fileout1D)