Example Dysmalpy 2D fitting, using fitting wrapper
In this example, we use dysmalpy to measure the kinematics of galaxy GS4_43501 at \(z=1.613\) in 2D, using a fitting wrapper which simplifies the implementation of the fitting algorithm to the user. In this specific case, the fittign method is \(\texttt{MPFIT}\), as specified at the bottom of the fitting_2D_mpfit.params file. The notebook shows how to find the best fit models for the two-dimensional velocity and velocity dispersion profiles (\(v(x,y)\) and \(\sigma(x,y)\)). These fits allow us to measure quantities such as the total mass (disk+bulge), the effective radius \(r_\mathrm{eff}\), dark matter fraction \(f_\mathrm{DM}\) and velocity dispersion \(\sigma_0\).
The fitting includes the following components:
Disk + Bulge
NFW halo
Constant velocity dispersion
The structure of the notebook is the following:
Setup steps (and load the params file)
Run
Dysmalpyfitting: 2D wrapper, with fit method= MPFITExamine results
1) Setup steps
Import modules
from __future__ import (absolute_import, division, print_function,
unicode_literals)
from dysmalpy.fitting_wrappers import dysmalpy_fit_single
from dysmalpy.fitting_wrappers import utils_io
from dysmalpy import fitting
import os
import numpy as np
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
# Setup plotting
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline
mpl.rcParams['figure.dpi']= 300
mpl.rc("savefig", dpi=300)
from IPython.core.display import Image
Set data, output paths
Note this will override the
datadirandoutdirspecified in the param file.(This is useful for the example here. When running from command line, it’s recommended to properly set the directories in the param file.)
# Data directory (datadir = /YOUR/DATA/PATH/)
filepath = os.path.abspath(fitting.__file__)
datadir = os.sep.join(os.path.dirname(filepath).split(os.sep)[:-1]+["tests", "test_data", ""])
# Load the parameters file from the examples directory
param_path = os.sep.join(os.path.dirname(filepath).split(os.sep)[:-1]+["examples", "examples_param_files", ""])
param_filename = param_path+'fitting_2D_mpfit.params'
print(param_filename)
# Where to save output files (output = /YOUR/OUTPUTS/PATH)
outdir = '/Users/sedona/data/dysmalpy_test_examples/JUPYTER_OUTPUT_2D_FITTING_WRAPPER/'
outdir = '/Users/jespejo/Dropbox/Postdoc/Data/dysmalpy_test_examples/JUPYTER_OUTPUT_2D_FITTING_WRAPPER/'
outdir_mpfit = outdir + 'MPFIT/'
/Users/jespejo/anaconda3/envs/test_dysmalpy/lib/python3.11/site-packages/dysmalpy/examples/examples_param_files/fitting_2D_mpfit.params
Settings in parameter file:
Note there are many commented out options / parameters. These given an more complete overview of the settings & parameters that can be specified with the fitting wrapper parameter files.
with open(param_filename, 'r') as f:
print(f.read())
# Example parameters file for fitting a single object with 1D data
# Note: DO NOT CHANGE THE NAMES IN THE 1ST COLUMN AND KEEP THE COMMAS!!
# See README for a description of each parameter and its available options.
# ******************************* OBJECT INFO **********************************
galID, GS4_43501 # Name of your object
z, 1.613 # Redshift
# ****************************** DATA INFO *************************************
datadir, None # Optional: Full path to data directory.
fdata_vel, GS4_43501_Ha_vm.fits # Full path to vel map. Alternatively, just the filename if 'datadir' is set.
fdata_verr, GS4_43501_Ha_vm_err.fits # Full path to vel. err map. Alternatively, just the filename if 'datadir' is set.
fdata_disp, GS4_43501_Ha_dm.fits # Full path to disp map. Alternatively, just the filename if 'datadir' is set.
fdata_derr, GS4_43501_Ha_dm_err.fits # Full path to disp. err map. Alternatively, just the filename if 'datadir' is set.
fdata_mask, GS4_43501_Ha_m.fits # Full path to mask
# -- strongly recommended to have a mask
# Alternatively, just the filename if 'datadir' is set.
data_inst_corr, True
symmetrize_data, False # Symmetrize data before fitting?
smoothing_type, median # Is the data median smoothed
# before extracting maps?
smoothing_npix, 3 # Number of pixels for smoothing aperture
## GALAXY CENTER:
# IMPORTANT: 0 indexed, so coordinates go from [0, nX-1] and [0, nY-1].
# So if using QFitsView, will need to subtract 1 from both coords (as QFitsView goes from [1,nX] and [1,nY])
xcenter, None # Center position in maps, x coord. Default: (nX-1)/2.
ycenter, None # Center position in maps, y coord. Default: (nY-1)/2.
# ***************************** OUTPUT *****************************************
outdir, GS4_43501_2D_out/ # Full path for output directory
# ***************************** OBSERVATION SETUP ******************************
# Instrument Setup
# ------------------
pixscale, 0.125 # Pixel scale in arcsec/pixel
fov_npix, 37 # Number of pixels on a side of model cube
spec_type, velocity # DON'T CHANGE!
spec_start, -1000. # Starting value for spectral axis // generally don't change
spec_step, 10. # Step size for spectral axis in km/s // generally don't change
nspec, 201 # Number of spectral steps // generally don't change
# LSF Setup
# ---------
use_lsf, True # True/False if using an LSF
sig_inst_res, 51.0 # Instrumental dispersion in km/s
# PSF Setup
# ---------
psf_type, Gaussian # Gaussian, Moffat, or DoubleGaussian
psf_fwhm, 0.55 # PSF FWHM in arcsecs
psf_beta, -99. # Beta parameter for a Moffat PSF
# ## ELLIPTICAL PSF:
# psf_type, Gaussian # Gaussian, Moffat, or DoubleGaussian
# psf_fwhm_major, 0.55 # PSF major axis FWHM in arcsecs
# psf_fwhm_minor, 0.25 # PSF minor axis FWHM in arcsecs
# psf_PA, 0. # PA of PSF major axis, in deg E of N. (0=N, 90=E)
# psf_beta, -99. # Beta parameter for a Moffat PSF
# # DoubleGaussian: settings instead of psf_fwhm
# psf_type, DoubleGaussian
# psf_fwhm1, 0.16 # FWHM of PSF component 1, in arcsecs. SINFONI AO: 0.16
# psf_fwhm2, 0.48 # FWHM of PSF component 1, in arcsecs. SINFONI AO: 0.48
# psf_scale1, 0.368 # Flux scaling (*not* peak height) of component 1. SINFONI AO: 0.368
# psf_scale2, 0.632 # Flux scaling (*not* peak height) of component 2. SINFONI AO: 0.632
# **************************** SETUP MODEL *************************************
# Model Settings
# -------------
# List of components to use: SEPARATE WITH SPACES
## MUST always keep: geometry zheight_gaus
## RECOMMENDED: always keep: disk+bulge const_disp_prof
components_list, disk+bulge const_disp_prof geometry zheight_gaus halo
# possible options:
# disk+bulge, sersic, blackhole
# const_disp_prof, geometry, zheight_gaus, halo,
# radial_flow, uniform_planar_radial_flow, uniform_bar_flow, uniform_wedge_flow,
# unresolved_outflow, biconical_outflow,
# CAUTION: azimuthal_planar_radial_flow, variable_bar_flow, spiral_flow
# List of components that emit light. SEPARATE WITH SPACES
## Current options: disk+bulge / bulge / disk [corresponding to the mass disk+bulge component],
## also: light_sersic, light_gaussian_ring
light_components_list, disk
# NOTE: if a separate light profile (eg light_sersic) is used,
# this MUST be changed to e.g., 'light_components_list, light_sersic'
adiabatic_contract, False # Apply adiabatic contraction?
pressure_support, True # Apply assymmetric drift correction?
noord_flat, True # Apply Noordermeer flattenning?
oversample, 1 # Spatial oversample factor
oversize, 1 # Oversize factor
moment_calc, False # If False, observed maps fit with GAUSSIANS
zcalc_truncate, True # Truncate in zgal direction when calculating or not
n_wholepix_z_min, 3 # Minimum number of whole pixels in zgal dir, if zcalc_truncate=True
# ********************************************************************************
# DISK + BULGE
# ------------
# Initial Values
total_mass, 11.0 # Total mass of disk and bulge log(Msun)
bt, 0.3 # Bulge-to-Total Ratio
r_eff_disk, 5.0 # Effective radius of disk in kpc
n_disk, 1.0 # Sersic index for disk
invq_disk, 5.0 # disk scale length to zheight ratio for disk
n_bulge, 4.0 # Sersic index for bulge
invq_bulge, 1.0 # disk scale length to zheight ratio for bulge
r_eff_bulge, 1.0 # Effective radius of bulge in kpc
# Fixed? True if its a fixed parameter, False otherwise
total_mass_fixed, False
r_eff_disk_fixed, False
bt_fixed, True
n_disk_fixed, True
r_eff_bulge_fixed, True
n_bulge_fixed, True
# Parameter bounds. Lower and upper bounds
total_mass_bounds, 10.0 13.0
bt_bounds, 0.0 1.0
r_eff_disk_bounds, 0.1 30.0
n_disk_bounds, 1.0 8.0
r_eff_bulge_bounds, 1.0 5.0
n_bulge_bounds, 1.0 8.0
# # ********************************************************************************
# # BLACK HOLE
# # ------------
#
# # Initial Values
# BH_mass, 11. # log(Msun)
#
# # Fixed? True if its a fixed parameter, False otherwise
# BH_mass_fixed, False
#
# # Parameter bounds. Lower and upper bounds
# BH_mass_bounds, 6. 18.
# # ********************************************************************************
# # Separate light profile: (Truncated) Sersic profile
# # ------------
# # Initial values
# L_tot_sersic, 1. # arbitrary units
# lr_eff, 4. # kpc
# lsersic_n, 1. # Sersic index of light profile
# lsersic_rinner, 0. # [kpc] Inner truncation radius of sersic profile. 0 = no truncation
# lsersic_router, inf # [kpc] Outer truncation radius of sersic profile. inf = no truncation
#
# # Fixed? True if its a fixed parameter, False otherwise
# L_tot_sersic_fixed, True
# lr_eff_fixed, False
# lsersic_n_fixed, True
# lsersic_rinner_fixed, True
# lsersic_router_fixed, True
#
# # Parameter bounds. Lower and upper bounds
# L_tot_sersic_bounds, 0. 2. # arbitrary units
# lr_eff_bounds, 0.5 15. # kpc
# lsersic_n_bounds, 0.5 8.
# lsersic_rinner_bounds, 0. 5. # kpc
# lsersic_router_bounds, 4. 20. # kpc
# # ********************************************************************************
# # Separate light profile: Gaussian ring
# # ------------
# # Initial values
# L_tot_gaus_ring, 1. # arbitrary units
# R_peak_gaus_ring, 6. # kpc
# FWHM_gaus_ring, 1. # kpc
#
# # Fixed? True if its a fixed parameter, False otherwise
# L_tot_gaus_ring_fixed, True
# R_peak_gaus_ring_fixed, True
# FWHM_gaus_ring_fixed, True
#
# # Parameter bounds. Lower and upper bounds
# L_tot_gaus_ring_bounds, 0. 2. # arbitrary units
# R_peak_gaus_ring_bounds, 0. 15. # kpc
# FWHM_gaus_ring_bounds, 0.1 10. # kpc
# ********************************************************************************
# ********************************************************************************
# ********************************************************************************
# DARK MATTER HALO
# ----------------
# Halo type: options: NFW / twopowerhalo / burkert / einasto / dekelzhao
halo_profile_type, NFW
# ** NOTE **: Uncomment the section below corresponding to the selected halo type.
# ********************************************************************************
# NFW halo
# Initial Values
mvirial, 11.0 # Halo virial mass in log(Msun)
halo_conc, 5.0 # Halo concentration parameter
fdm, 0.5 # Dark matter fraction at r_eff_disk
# Fixed? True if its a fixed parameter, False otherwise. Also set False if it will be tied (below)
mvirial_fixed, True
halo_conc_fixed, True
fdm_fixed, False
# Parameter bounds. Lower and upper bounds
mvirial_bounds, 10.0 13.0
halo_conc_bounds, 1.0 20.0
fdm_bounds, 0.0 1.0
# Tie the parameters?
fdm_tied, True # for NFW, fdm_tied=True determines fDM from Mvirial (+baryons)
mvirial_tied, False # for NFW, mvirial_tied=True determines Mvirial from fDM (+baryons)
# ********************************************************************************
# # ********************************************************************************
# # Two-power halo
#
# # Initial Values
# mvirial, 11.5 # Halo virial mass in log(Msun)
# halo_conc, 5.0 # Halo concentration parameter
# fdm, 0.5 # Dark matter fraction at r_eff_disk
# alpha, 1. # TPH: inner slope. NFW has alpha=1
# beta, 3. # TPH: outer slope. NFW has beta=3
#
# # Fixed? True if its a fixed parameter, False otherwise. Also set False if it will be tied (below)
# mvirial_fixed, False
# halo_conc_fixed, True
# fdm_fixed, False
# alpha_fixed, False
# beta_fixed, True
#
# # Parameter bounds. Lower and upper bounds
# mvirial_bounds, 10.0 13.0
# halo_conc_bounds, 1.0 20.0
# fdm_bounds, 0.0 1.0
# alpha_bounds, 0.0 3.0
# beta_bounds, 1.0 4.0
#
# # Tie the parameters?
# fdm_tied, True # for non-NFW, fdm_tied=True determines fDM from other halo params (+baryons)
# mvirial_tied, True # for non-NFW, mvirial_tied=True determines Mvirial from SMHM+fgas + baryon total_mass
# alpha_tied, False # for TPH, alpha_tied=True determines alpha from free fDM + other parameters.
#
# ### OTHER SETTINGS:
# mhalo_relation, Moster18 ## SMHM relation to use for tying Mvir to Mbar. options: Moster18 / Behroozi13
#
# fgas, 0.5 # Gas fraction for SMHM inference of Mvir if 'mvirial_tied=True'
# lmstar, -99. # Currently code uses fgas to infer lmstar
# # from fitting baryon total_mass for SMHM relation
# # ********************************************************************************
# # ********************************************************************************
# # Burkert halo
#
# # Initial Values
# mvirial, 11.5 # Halo virial mass in log(Msun)
# halo_conc, 5.0 # Halo concentration parameter
# fdm, 0.5 # Dark matter fraction at r_eff_disk
# rB, 10. # Burkert: Halo core radius, in kpc
#
# # Fixed? True if its a fixed parameter, False otherwise. Also set False if it will be tied (below)
# mvirial_fixed, False
# halo_conc_fixed, True
# fdm_fixed, False
# rB_fixed, False
#
# # Parameter bounds. Lower and upper bounds
# mvirial_bounds, 10.0 13.0
# halo_conc_bounds, 1.0 20.0
# fdm_bounds, 0.0 1.0
# rB_bounds, 1.0 20.0
#
# # Tie the parameters?
# fdm_tied, True # for non-NFW, fdm_tied=True determines fDM from other halo params (+baryons)
# mvirial_tied, True # for non-NFW, mvirial_tied=True determines Mvirial from SMHM+fgas + baryon total_mass
# rB_tied, False # for Burkert, rB_tied=True determines rB from free fDM + other parameters.
#
# ### OTHER SETTINGS:
# mhalo_relation, Moster18 ## SMHM relation to use for tying Mvir to Mbar. options: Moster18 / Behroozi13
#
# fgas, 0.5 # Gas fraction for SMHM inference of Mvir if 'mvirial_tied=True'
# lmstar, -99. # Currently code uses fgas to infer lmstar
# # from fitting baryon total_mass for SMHM relation
# # ********************************************************************************
# # ********************************************************************************
# # Einasto halo
# # Initial Values
# mvirial, 11.5 # Halo virial mass in log(Msun)
# halo_conc, 5.0 # Halo concentration parameter
# fdm, 0.5 # Dark matter fraction at r_eff_disk
# alphaEinasto, 1. # Einasto: Halo profile index
#
# # Fixed? True if its a fixed parameter, False otherwise. Also set False if it will be tied (below)
# mvirial_fixed, False
# halo_conc_fixed, True
# fdm_fixed, False
# alphaEinasto_fixed, False
#
# # Parameter bounds. Lower and upper bounds
# mvirial_bounds, 10.0 13.0
# halo_conc_bounds, 1.0 20.0
# fdm_bounds, 0.0 1.0
# alphaEinasto_bounds, 0.0 2.0
#
# # Tie the parameters?
# fdm_tied, True # for non-NFW, fdm_tied=True determines fDM from other halo params (+baryons)
# mvirial_tied, True # for non-NFW, mvirial_tied=True determines Mvirial from SMHM+fgas + baryon total_mass
# alphaEinasto_tied, False # for Einasto, alphaEinasto_tied=True determines alphaEinasto from free fDM + other params.
#
# ### OTHER SETTINGS:
# mhalo_relation, Moster18 ## SMHM relation to use for tying Mvir to Mbar. options: Moster18 / Behroozi13
#
# fgas, 0.5 # Gas fraction for SMHM inference of Mvir if 'mvirial_tied=True'
# lmstar, -99. # Currently code uses fgas to infer lmstar
# # from fitting baryon total_mass for SMHM relation
# # ********************************************************************************
# # ********************************************************************************
# # Dekel-Zhao halo
# # Initial Values
# mvirial, 12.0 # Halo virial mass in log(Msun)
# s1, 1.5 # Inner logarithmic slope (at resolution r1=0.01*Rvir)
# c2, 25.0 # Concentration parameter (defined relative to c, a)
# fdm, 0.5 # Dark matter fraction at r_eff_disk
#
# # Fixed? True if its a fixed parameter, False otherwise. Also set False if it will be tied (below)
# mvirial_fixed, False
# s1_fixed, False
# c2_fixed, False
# fdm_fixed, False
#
# # Parameter bounds. Lower and upper bounds
# mvirial_bounds, 10.0 13.0 # log(Msun)
# s1_bounds, 0.0 2.0
# c2_bounds, 0.0 40.0
# fdm_bounds, 0.0 1.0
#
# # Tie the parameters?
# mvirial_tied, True # mvirial_tied=True determines Mvirial from fDM, s1, c2.
# s1_tied, True # Tie the s1 to M*/Mvir using best-fit Freundlich+20 (Eqs 45, 47, 48, Table 1)
# c2_tied, True # Tie the c2 to M*/Mvir using best-fit Freundlich+20 (Eqs 47, 49, Table 1)
# fdm_tied, False # for non-NFW, fdm_tied=True determines fDM from other halo params (+baryons)
#
# ### OTHER SETTINGS:
# lmstar, 10.5 # Used to infer s1, c2 if s1_tied or c2_tied = True
#
# # ********************************************************************************
# ********************************************************************************
# INTRINSIC DISPERSION PROFILE
# ------------------
# Initial Values
sigma0, 39.0 # Constant intrinsic dispersion value
# Fixed? True if its a fixed parameter, False otherwise
sigma0_fixed, False
# Parameter bounds. Lower and upper bounds
sigma0_bounds, 5.0 300.0
# ********************************************************************************
# ********************************************************************************
# ********************************************************************************
# HIGHER ORDER COMPONENTS: INFLOW, OUTFLOW
# ----------------------------------------
# # ********************************************************************************
# # UNIFORM SPHERICAL RADIAL FLOW -- in rhat direction in spherical coordinates
# # radial_flow
# # -------------------
#
# vr, -90. # Radial flow [km/s]. Positive: Outflow. Negative: Inflow.
# # ********************************************************************************
# # UNIFORM PLANAR RADIAL FLOW -- in Rhat direction in cylindrical coordinates
# # (eg, radial in galaxy midplane)
# # uniform_planar_radial_flow
# # -------------------
#
# vr, -90. # Radial flow [km/s]. Positive: Outflow. Negative: Inflow.
# # ********************************************************************************
# # UNIFORM BAR FLOW -- in xhat direction along bar in cartesian coordinates,
# # with bar at an angle relative to galaxy major axis (blue)
# # uniform_bar_flow
# # -------------------
#
# vbar, -90. # Bar flow [km/s]. Positive: Outflow. Negative: Inflow.
# phi, 90. # Azimuthal angle of bar [degrees], counter-clockwise from blue major axis.
# # Default is 90 (eg, along galaxy minor axis)
# bar_width, 2 # Width of the bar perpendicular to bar direction.
# # Bar velocity only is nonzero between -bar_width/2 < ygal < bar_width/2.
# # ********************************************************************************
# # UNIFORM WEDGE FLOW -- in planar radial flow in cylindrical coordinates, restricted to pos, neg wedges
# # uniform_wedge_flow
# # -------------------
#
# vr, -90. # Radial flow [km/s]. Positive: Outflow. Negative: Inflow.
# theta, 60. # Opening angle of wedge [deg]. (the full angular span)
# phi, 90. # Angle offset relative to the galaxy angle, so the wedge center is at phi.
# # Default: 90 deg, so centered along minor axis
# # ********************************************************************************
# # UNRESOLVED OUTFLOW -- at galaxy center (ie, AGN unresolved outflow)
# # unresolved_outflow
# # -------------------
#
# vcenter, 0. # Central velocity of the Gaussian in km/s
# fwhm, 1000. # FWHM of the Gaussian in km/s
# amplitude, 1.e12 # Amplitude of the Gaussian, for flux in ~M/L=1 luminosity units
# # with the dimming applied ... roughly ....
# # ********************************************************************************
# # BICONICAL OUTFLOW
# # biconical_outflow
# # -------------------
#
# n, 0.5 # Power law index
# vmax, 500. # Maximum velocity of the outflow in km/s
# rturn, 5. # Turn-over radius in kpc of the velocty profile
# thetain, 30. # Half inner opening angle in degrees. Measured from the bicone axis
# dtheta, 20. # Difference between inner and outer opening angle in degrees
# rend, 10. # Maximum radius of the outflow in kpc
# norm_flux, 8. # Log flux amplitude of the outflow at r = 0.
# # Need to check dimming/flux conventions
# tau_flux, 1. # Exponential decay rate of the flux
# biconical_profile_type, both # Type of velocity profile:
# # 'both', 'increase', 'decrease', 'constant'
# biconical_outflow_dispersion, 80. # Dispersion (stddev of gaussian) of biconical outflow, km/s
# # ********************************************************************************
# # VARIABLE BAR FLOW -- in xhat direction along bar in cartesian coordinates,
# # with bar at an angle relative to galaxy major axis (blue)
# # CAUTION!!!
# # variable_bar_flow
# # -------------------
#
# vbar_func_bar_flow, -90.*np.exp(-R/5.) # Bar flow FUNCTION [km/s]. Positive: Outflow. Negative: Inflow.
# phi, 90. # Azimuthal angle of bar [degrees], counter-clockwise from blue major axis.
# # Default is 90 (eg, along galaxy minor axis)
# bar_width, 2 # Width of the bar perpendicular to bar direction.
# # Bar velocity only is nonzero between -bar_width/2 < ygal < bar_width/2.
# # ********************************************************************************
# # AZIMUTHAL PLANAR RADIAL FLOW -- in Rhat direction in cylindrical coordinates
# # (eg, radial in galaxy midplane), with an added azimuthal term
# # CAUTION!!!
# # azimuthal_planar_radial_flow
# # -------------------
#
# vr_func_azimuthal_planar_flow, -90.*np.exp(-R/5.) # Radial flow [km/s].
# # Positive: Outflow. Negative: Inflow.
# m, 2 # Number of modes in the azimuthal pattern. m=0 gives a purely radial profile.
# phi0, 0. # Angle offset relative to the galaxy angle [deg],
# # so the azimuthal variation goes as cos(m(phi_gal - phi0))
# # ********************************************************************************
# # SPIRAL DENSIY WAVE FLOW -- as in Davies et al. 2009, ApJ, 702, 114
# # Here assuming CONSTANT velocity -- try to match real Vrot...
# # CAUTION!!! NO SPACES IN FUNCTION DEFINITONS!
# # spiral_flow
# # -------------------
#
# Vrot_func_spiral_flow, 150.+0.*R # Unperturbed rotation velocity of the galaxy
# dVrot_dR_func_spiral_flow, 0.*R # Derivative of Vrot(R) -- ideally evaluated analytically, otherwise very slow.
# rho0_func_spiral_flow, 1.e11*np.exp(-R/5.) # Unperturbed midplane density profile of the galaxy
# f_func_spiral_flow, (np.sqrt(m**2-2.)*Vrot(R)/cs)*np.log(R) # Function describing the spiral shape, m*phi = f(R)
# # with k = df/dR
# k_func_spiral_flow, (np.sqrt(m**2-2.)*Vrot(R)/cs)/R # Function for radial wavenumber
#
# m, 2 # Number of photometric/density spiral arms.
# cs, 10. # Sound speed of medium, in km/s.
# epsilon, 1. # Density contrast of perturbation (unitless).
# Om_p, 0. # Angular speed of the driving force, Omega_p
# phi0, 0. # Angle offset of the arm winding, in degrees. Default: 0.
# ********************************************************************************
# ********************************************************************************
# ********************************************************************************
# ********************************************************************************
# ZHEIGHT PROFILE
# ---------------
# Initial Values
sigmaz, 0.9 # Gaussian width of the galaxy in z, in kpc
# Fixed? True if its a fixed parameter, False otherwise
sigmaz_fixed, False
# Parameter bounds. Lower and upper bounds
sigmaz_bounds, 0.1 1.0
# Tie the zheight to the effective radius of the disk?
# If set to True, make sure sigmaz_fixed is False
zheight_tied, True
# GEOMETRY
# --------
# Initial Values
inc, 72. # Inclination of galaxy, 0=face-on, 90=edge-on
pa, 145.
xshift, 0. # pixels
yshift, -0.13 # pixels
vel_shift, 0. # km/s
# Fixed? True if its a fixed parameter, False otherwise
inc_fixed, False
pa_fixed, False
xshift_fixed, False
yshift_fixed, False
vel_shift_fixed, False
# Parameter bounds. Lower and upper bounds
inc_bounds, 42. 82.
pa_bounds, 132. 152.
xshift_bounds, -1.5 1.5 # pixels
yshift_bounds, -1.5 1.5 # pixels
vel_shift_bounds, -100. 100. # km/s
# **************************** Fitting Settings ********************************
fit_method, mpfit # mcmc, nested, or mpfit
do_plotting, True # Produce all output plots?
fitdispersion, True # Simultaneously fit the velocity and dispersion?
fitflux, False # Also fit for the flux?
# MPFIT Settings
#----------------
maxiter, 200 # Maximum number of iterations before mpfit quits
Add some settings for this notebook example:
plot_type = 'png'
Run Dysmalpy fitting: 2D wrapper, with fit method= MPFIT
dysmalpy_fit_single.dysmalpy_fit_single(param_filename=param_filename,
datadir=datadir, outdir=outdir_mpfit,
plot_type='png', overwrite=True)
INFO:DysmalPy:dysmalpy.Galaxy:
********************************************************************
*** INFO ***
instrument.fov[0,1]=(37,37) is being reset to match 2D maps (27, 27)
********************************************************************
INFO:DysmalPy:*************************************
INFO:DysmalPy: Fitting: GS4_43501 using MPFIT
INFO:DysmalPy: obs: OBS
INFO:DysmalPy: velocity file: /Users/jespejo/anaconda3/envs/test_dysmalpy/lib/python3.11/site-packages/dysmalpy/tests/test_data/GS4_43501_Ha_vm.fits
INFO:DysmalPy: dispers. file: /Users/jespejo/anaconda3/envs/test_dysmalpy/lib/python3.11/site-packages/dysmalpy/tests/test_data/GS4_43501_Ha_dm.fits
INFO:DysmalPy: nSubpixels: 1
INFO:DysmalPy: mvirial_tied: False
INFO:DysmalPy:
MPFIT Fitting:
Start: 2024-11-04 15:59:59.073699
INFO:DysmalPy:Iter 1 CHI-SQUARE = 4863.95731 DOF = 404
disk+bulge:total_mass = 11
disk+bulge:r_eff_disk = 5
dispprof_LINE:sigma0 = 39
geom_1:inc = 72
geom_1:pa = 145
geom_1:xshift = 0
geom_1:yshift = -0.13
geom_1:vel_shift = 0
INFO:DysmalPy:Iter 2 CHI-SQUARE = 3391.858655 DOF = 404
disk+bulge:total_mass = 11.02256749
disk+bulge:r_eff_disk = 6.314016605
dispprof_LINE:sigma0 = 33.36135137
geom_1:inc = 80.1969795
geom_1:pa = 144.0786096
geom_1:xshift = -0.3915117315
geom_1:yshift = -0.2053268919
geom_1:vel_shift = 0.470343974
INFO:DysmalPy:Iter 3 CHI-SQUARE = 3070.291053 DOF = 404
disk+bulge:total_mass = 11.17493291
disk+bulge:r_eff_disk = 8.822236059
dispprof_LINE:sigma0 = 28.7835928
geom_1:inc = 79.525108
geom_1:pa = 143.0253885
geom_1:xshift = -0.4590096803
geom_1:yshift = -0.2284259521
geom_1:vel_shift = 0.3610861955
INFO:DysmalPy:Iter 4 CHI-SQUARE = 2124.533719 DOF = 404
disk+bulge:total_mass = 11.23876486
disk+bulge:r_eff_disk = 10.84611558
dispprof_LINE:sigma0 = 26.41974081
geom_1:inc = 78.24159899
geom_1:pa = 144.0666743
geom_1:xshift = -0.5126667736
geom_1:yshift = -0.5584810651
geom_1:vel_shift = 0.8034814807
INFO:DysmalPy:Iter 5 CHI-SQUARE = 1858.95427 DOF = 404
disk+bulge:total_mass = 11.22284368
disk+bulge:r_eff_disk = 10.55852661
dispprof_LINE:sigma0 = 26.95611728
geom_1:inc = 77.02489459
geom_1:pa = 145.2829668
geom_1:xshift = -0.5206432527
geom_1:yshift = -0.8504292774
geom_1:vel_shift = 1.037174844
INFO:DysmalPy:Iter 6 CHI-SQUARE = 1799.964493 DOF = 404
disk+bulge:total_mass = 11.24931363
disk+bulge:r_eff_disk = 11.37911566
dispprof_LINE:sigma0 = 26.63716887
geom_1:inc = 77.25523512
geom_1:pa = 146.549937
geom_1:xshift = -0.3864614764
geom_1:yshift = -0.9338220054
geom_1:vel_shift = 0.9299774842
INFO:DysmalPy:Iter 7 CHI-SQUARE = 1755.790678 DOF = 404
disk+bulge:total_mass = 11.23880745
disk+bulge:r_eff_disk = 11.10273218
dispprof_LINE:sigma0 = 24.51821952
geom_1:inc = 76.93602404
geom_1:pa = 146.1685048
geom_1:xshift = -0.1835254312
geom_1:yshift = -1.073486741
geom_1:vel_shift = 0.7280803198
INFO:DysmalPy:Iter 8 CHI-SQUARE = 1755.033023 DOF = 404
disk+bulge:total_mass = 11.24050004
disk+bulge:r_eff_disk = 11.11478921
dispprof_LINE:sigma0 = 24.69401168
geom_1:inc = 76.63700138
geom_1:pa = 146.7513809
geom_1:xshift = -0.1953901937
geom_1:yshift = -1.070126186
geom_1:vel_shift = 0.7261058258
INFO:DysmalPy:Iter 9 CHI-SQUARE = 1750.217717 DOF = 404
disk+bulge:total_mass = 11.24061494
disk+bulge:r_eff_disk = 11.11607529
dispprof_LINE:sigma0 = 25.01605169
geom_1:inc = 76.67592341
geom_1:pa = 146.7793715
geom_1:xshift = -0.185203532
geom_1:yshift = -1.070160913
geom_1:vel_shift = 0.7289132677
INFO:DysmalPy:Iter 10 CHI-SQUARE = 1750.082449 DOF = 404
disk+bulge:total_mass = 11.24063747
disk+bulge:r_eff_disk = 11.11622614
dispprof_LINE:sigma0 = 25.11113924
geom_1:inc = 76.67072013
geom_1:pa = 146.7774838
geom_1:xshift = -0.1858546538
geom_1:yshift = -1.070151373
geom_1:vel_shift = 0.7287818672
INFO:DysmalPy:Iter 11 CHI-SQUARE = 1750.022426 DOF = 404
disk+bulge:total_mass = 11.24064704
disk+bulge:r_eff_disk = 11.11627646
dispprof_LINE:sigma0 = 25.12610883
geom_1:inc = 76.66746622
geom_1:pa = 146.7783321
geom_1:xshift = -0.1858525364
geom_1:yshift = -1.070192788
geom_1:vel_shift = 0.7289287385
INFO:DysmalPy:Iter 12 CHI-SQUARE = 1749.878728 DOF = 404
disk+bulge:total_mass = 11.24065361
disk+bulge:r_eff_disk = 11.11629826
dispprof_LINE:sigma0 = 25.15805123
geom_1:inc = 76.6668677
geom_1:pa = 146.7786316
geom_1:xshift = -0.1856982592
geom_1:yshift = -1.070196524
geom_1:vel_shift = 0.7289823194
INFO:DysmalPy:Iter 13 CHI-SQUARE = 1749.864321 DOF = 404
disk+bulge:total_mass = 11.24065464
disk+bulge:r_eff_disk = 11.11629646
dispprof_LINE:sigma0 = 25.16055143
geom_1:inc = 76.66647618
geom_1:pa = 146.7787255
geom_1:xshift = -0.1856770347
geom_1:yshift = -1.070198525
geom_1:vel_shift = 0.7289948036
INFO:DysmalPy:Iter 14 CHI-SQUARE = 1749.835605 DOF = 404
disk+bulge:total_mass = 11.24065611
disk+bulge:r_eff_disk = 11.11628035
dispprof_LINE:sigma0 = 25.16788657
geom_1:inc = 76.66603205
geom_1:pa = 146.7790628
geom_1:xshift = -0.185658349
geom_1:yshift = -1.070197881
geom_1:vel_shift = 0.7290073537
INFO:DysmalPy:Iter 15 CHI-SQUARE = 1749.806514 DOF = 404
disk+bulge:total_mass = 11.24065803
disk+bulge:r_eff_disk = 11.11632404
dispprof_LINE:sigma0 = 25.17372288
geom_1:inc = 76.66555138
geom_1:pa = 146.7795419
geom_1:xshift = -0.1856254474
geom_1:yshift = -1.070213453
geom_1:vel_shift = 0.7290327874
INFO:DysmalPy:Iter 16 CHI-SQUARE = 1749.801064 DOF = 404
disk+bulge:total_mass = 11.2406584
disk+bulge:r_eff_disk = 11.1163272
dispprof_LINE:sigma0 = 25.17472982
geom_1:inc = 76.6652504
geom_1:pa = 146.7795455
geom_1:xshift = -0.185617335
geom_1:yshift = -1.070215226
geom_1:vel_shift = 0.7290310108
INFO:DysmalPy:Iter 17 CHI-SQUARE = 1749.793622 DOF = 404
disk+bulge:total_mass = 11.24065907
disk+bulge:r_eff_disk = 11.11632529
dispprof_LINE:sigma0 = 25.17763378
geom_1:inc = 76.66499883
geom_1:pa = 146.7796311
geom_1:xshift = -0.1856239212
geom_1:yshift = -1.070220109
geom_1:vel_shift = 0.7290276926
INFO:DysmalPy:Iter 18 CHI-SQUARE = 1749.782348 DOF = 404
disk+bulge:total_mass = 11.24066023
disk+bulge:r_eff_disk = 11.11631647
dispprof_LINE:sigma0 = 25.18002597
geom_1:inc = 76.66412292
geom_1:pa = 146.779641
geom_1:xshift = -0.1856137391
geom_1:yshift = -1.070225286
geom_1:vel_shift = 0.729034287
INFO:DysmalPy:Iter 19 CHI-SQUARE = 1749.760549 DOF = 404
disk+bulge:total_mass = 11.24066182
disk+bulge:r_eff_disk = 11.11635545
dispprof_LINE:sigma0 = 25.18380674
geom_1:inc = 76.6633357
geom_1:pa = 146.7797417
geom_1:xshift = -0.1855873429
geom_1:yshift = -1.070223568
geom_1:vel_shift = 0.7290825612
INFO:DysmalPy:Iter 20 CHI-SQUARE = 1749.756106 DOF = 404
disk+bulge:total_mass = 11.24066218
disk+bulge:r_eff_disk = 11.11635462
dispprof_LINE:sigma0 = 25.18491728
geom_1:inc = 76.66313876
geom_1:pa = 146.7798613
geom_1:xshift = -0.185584283
geom_1:yshift = -1.070222421
geom_1:vel_shift = 0.7290837185
INFO:DysmalPy:
End: 2024-11-04 16:00:11.117629
******************
Time= 12.04 (sec), 0:12.04 (m:s)
MPFIT Status = 2
MPFIT Error/Warning Message = None
******************
------------------------------------------------------------------
Dysmalpy MPFIT fitting complete for: GS4_43501
output folder: /Users/jespejo/Dropbox/Postdoc/Data/dysmalpy_test_examples/JUPYTER_OUTPUT_2D_FITTING_WRAPPER/MPFIT/
------------------------------------------------------------------
Examine results
Result plots
Read in parameter file
params = utils_io.read_fitting_params(fname=param_filename)
# Override data + output paths:
params['datadir'] = datadir
params['outdir'] = outdir_mpfit
# Add the plot type:
params['plot_type'] = plot_type
f_galmodel = params['outdir'] + '{}_model.pickle'.format(params['galID'])
f_results = params['outdir'] + '{}_{}_results.pickle'.format(params['galID'],
params['fit_method'])
Best-fit plot
# Look at best-fit:
filepath = outdir_mpfit+"{}_{}_bestfit_{}.{}".format(params['galID'],
params['fit_method'],
params['obs_1_name'],
params['plot_type'])
Image(filename=filepath, width=600)
Directly generating result plots
Reload the galaxy, results files:
gal, results = fitting.reload_all_fitting(filename_galmodel=f_galmodel,
filename_results=f_results,
fit_method=params['fit_method'])
Plot the best-fit results:
results.plot_results(gal)
Results reports
We now look at the results reports, which include the best-fit values and uncertainties (as well as other fitting settings and output).
# Print report
print(results.results_report(gal=gal))
###############################
Fitting for GS4_43501
Date: 2024-11-04 16:00:12.513962
obs: OBS
Datafiles:
vel : /Users/jespejo/anaconda3/envs/test_dysmalpy/lib/python3.11/site-packages/dysmalpy/tests/test_data/GS4_43501_Ha_vm.fits
disp: /Users/jespejo/anaconda3/envs/test_dysmalpy/lib/python3.11/site-packages/dysmalpy/tests/test_data/GS4_43501_Ha_dm.fits
fit_velocity: True
fit_dispersion: True
fit_flux: False
moment: False
zcalc_truncate: True
n_wholepix_z_min: 3
oversample: 1
oversize: 1
Fitting method: MPFIT
fit status: 2
pressure_support: True
pressure_support_type: 1
###############################
Fitting results
-----------
disk+bulge
total_mass 11.2407 +/- 0.0135
r_eff_disk 11.1164 +/- 0.4071
n_disk 1.0000 [FIXED]
r_eff_bulge 1.0000 [FIXED]
n_bulge 4.0000 [FIXED]
bt 0.3000 [FIXED]
mass_to_light 1.0000 [FIXED]
noord_flat True
-----------
halo
mvirial 11.0000 [FIXED]
fdm 0.1583 [TIED]
conc 5.0000 [FIXED]
-----------
dispprof_LINE
sigma0 25.1849 +/- 1.1831
-----------
zheightgaus
sigmaz 1.8883 [TIED]
-----------
geom_1
inc 76.6631 +/- 0.5704
pa 146.7799 +/- 0.3628
xshift -0.1856 +/- 0.0104
yshift -1.0702 +/- 0.0212
vel_shift 0.7291 +/- 0.0600
-----------
Adiabatic contraction: False
-----------
Red. chisq: 4.3311
-----------
obs OBS: Rout,max,2D: 11.1201
To directly save the results report to a file, we can use the following:
# Save report to file:
f_report = params['outdir'] + '{}_fit_report.txt'.format(params['galID'])
results.results_report(gal=gal, filename=f_report)
Also note the fitting wrappers automatically save two versions of the report files:
fbase = '{}_{}_bestfit_results'.format(params['galID'], params['fit_method'])
f_report_pretty = params['outdir'] + fbase + '_report.info'
f_report_machine = params['outdir'] + fbase + '.dat'
The “pretty” version, automatically saved as *_best_fit_results_report.info, is formatted to be human-readable, and includes more information on the fit settings at the beginning (for reference).
with open(f_report_pretty, 'r') as f:
lines = [line.rstrip() for line in f]
for line in lines: print(line)
###############################
Fitting for GS4_43501
Date: 2024-11-04 16:00:11.922759
obs: OBS
Datafiles:
vel : /Users/jespejo/anaconda3/envs/test_dysmalpy/lib/python3.11/site-packages/dysmalpy/tests/test_data/GS4_43501_Ha_vm.fits
disp: /Users/jespejo/anaconda3/envs/test_dysmalpy/lib/python3.11/site-packages/dysmalpy/tests/test_data/GS4_43501_Ha_dm.fits
fit_velocity: True
fit_dispersion: True
fit_flux: False
moment: False
zcalc_truncate: True
n_wholepix_z_min: 3
oversample: 1
oversize: 1
Fitting method: MPFIT
fit status: 2
pressure_support: True
pressure_support_type: 1
###############################
Fitting results
-----------
disk+bulge
total_mass 11.2407 +/- 0.0135
r_eff_disk 11.1164 +/- 0.4071
n_disk 1.0000 [FIXED]
r_eff_bulge 1.0000 [FIXED]
n_bulge 4.0000 [FIXED]
bt 0.3000 [FIXED]
mass_to_light 1.0000 [FIXED]
noord_flat True
-----------
halo
mvirial 11.0000 [FIXED]
fdm 0.1583 [TIED]
conc 5.0000 [FIXED]
-----------
dispprof_LINE
sigma0 25.1849 +/- 1.1831
-----------
zheightgaus
sigmaz 1.8883 [TIED]
-----------
geom_1
inc 76.6631 +/- 0.5704
pa 146.7799 +/- 0.3628
xshift -0.1856 +/- 0.0104
yshift -1.0702 +/- 0.0212
vel_shift 0.7291 +/- 0.0600
-----------
Adiabatic contraction: False
-----------
Red. chisq: 4.3311
-----------
obs OBS: Rout,max,2D: 11.1201
The “machine” version, automatically saved as *_best_fit_results.dat, is formatted as a machine-readable space-separated ascii file. It includes key parameter fit information, as well as the best-fit reduced chisq.
with open(f_report_machine, 'r') as f:
lines = [line.rstrip() for line in f]
for line in lines: print(line)
# component param_name fixed best_value l68_err u68_err
disk+bulge total_mass False 11.2407 0.0135 0.0135
disk+bulge r_eff_disk False 11.1164 0.4071 0.4071
disk+bulge n_disk True 1.0000 -99.0000 -99.0000
disk+bulge r_eff_bulge True 1.0000 -99.0000 -99.0000
disk+bulge n_bulge True 4.0000 -99.0000 -99.0000
disk+bulge bt True 0.3000 -99.0000 -99.0000
disk+bulge mass_to_light True 1.0000 -99.0000 -99.0000
halo mvirial True 11.0000 -99.0000 -99.0000
halo fdm TIED 0.1583 -99.0000 -99.0000
halo conc True 5.0000 -99.0000 -99.0000
dispprof_LINE sigma0 False 25.1849 1.1831 1.1831
zheightgaus sigmaz TIED 1.8883 -99.0000 -99.0000
geom_1 inc False 76.6631 0.5704 0.5704
geom_1 pa False 146.7799 0.3628 0.3628
geom_1 xshift False -0.1856 0.0104 0.0104
geom_1 yshift False -1.0702 0.0212 0.0212
geom_1 vel_shift False 0.7291 0.0600 0.0600
mvirial ----- ----- 11.0000 -99.0000 -99.0000
fit_status ----- ----- 2 -99.0000 -99.0000
adiab_contr ----- ----- False -99.0000 -99.0000
redchisq ----- ----- 4.3311 -99.0000 -99.0000
noord_flat ----- ----- True -99.0000 -99.0000
pressure_support ----- ----- True -99.0000 -99.0000
pressure_support_type ----- ----- 1 -99.0000 -99.0000
obs:OBS:moment ----- ----- False -99.0000 -99.0000
obs:OBS:Routmax2D ----- ----- 11.1201 -99.0000 -99.0000