Galaxy models¶

Code author: Wilfried Mercier - IRAP <wilfried.mercier@irap.omp.eu>

Useful functions for galaxy modelling and other related computation.

Table of contents

1D profiles
2D profiles

1D profiles ¶

A collection of 1D surface brightness profiles.

galaxy.models.bulge(r, re, b4=None, Ie=None, mag=None, offset=None)[source]¶

Code author: Wilfried Mercier - IRAP <wilfried.mercier@irap.omp.eu>

Computes the value of the intensity of a de Vaucouleur bulge at position r defined as

\[\Sigma(r) = I_{\rm{e}} e^{\left [ \left (r/R_{\rm{e}} \right )^{1/4} - 1 \right ]},\]

with \(R_{\rm{e}}\) the effective radius and \(I_{\rm{e}}\) the surface brightness at the effective radius.

Parameters

r (float) – position at which the profile is computed
re (float) – half-light radius
b4 (float) – (Optional) b4 factor appearing in the Sersic profile. If None, its value will be computed.
Ie (float) – (Optional) surface brightness at half-light radius
mag (float) – (Optional) total integrated magnitude used to compute Ie if not given
offset (float) – (Optional) magnitude offset in the magnitude system used

Returns

surface brightness

Return type

float

galaxy.models.exponential_disk(r, re, b1=None, Ie=None, mag=None, offset=None)[source]¶

Code author: Wilfried Mercier - IRAP <wilfried.mercier@irap.omp.eu>

Computes the value of the intensity of an exponential disk at position r defined as

\[\Sigma(r) = I_{\rm{e}} e^{\left [ r/R_{\rm{e}} - 1 \right ]},\]

with \(R_{\rm{e}}\) the effective radius and \(I_{\rm{e}}\) the surface brightness at the effective radius.

Parameters

r (float) – position at which the profile is computed
re (float) – half-light radius
b1 (float) – (Optional) b1 factor appearing in the Sersic profile. If None, its value will be computed.
Ie (float) – (Optional) surface brightness at half-light radius
mag (float) – (Optional) total integrated magnitude used to compute Ie if not given
offset (float) – (Optional) magnitude offset in the magnitude system used

Returns

surface brightness

Return type

float

galaxy.models.hernquist(r, a, M)[source]¶

Code author: Wilfried Mercier - IRAP <wilfried.mercier@irap.omp.eu>

Hernquist profile defined as

\[\rho(r) = \frac{M_{\rm{b}}}{2\pi} \frac{a}{r} (r + a)^{-3},\]

with \(M_{\rm{b}}\) the total mass and \(a\) the scale radius.

Parameters

a (int or float) – scale radius
M (int or float) – total mass
r (int or float or ndarray[int] or ndarray[float]) – radial distance(s) where to compute the Hernquist profile. Unit must be the same as a.

Returns

Hernquist profile evaluated at the given distance(s). Unit is that of M/a^3.

Return type

float or ndarray[float]

Raises

TypeError – if r, M and a are neither int, nor float
ValueError – if np.any(r) < 0, if a <= 0 or if M < 0

galaxy.models.nfw(r, c, rs)[source]¶

Code author: Wilfried Mercier - IRAP <wilfried.mercier@irap.omp.eu>

NFW profile defined as

\[\rho(r) = \delta_{\rm{c}} \rho_{\rm{crit}} (r/r_{\rm{s}})^{-1} (1 + r/r_{\rm{s}})^{-2},\]

with \(r_{\rm{s}} = r_{200} / c\) the halo scale radius, with \(r_{200}\) the virial radius where the mean overdensity is equal to 200, \(c\) the halo concentration, \(\rho_{\rm{crit}} = 3 H_0^2 / (8\pi G)\) the Universe closure density, and \(\delta_{\rm{c}}\) the halo overdensity.

Parameters

c (int or float) – halo concentration
r (int or float or ndarray[int] or ndarray[float]) – radial distance(s) where to compute the profile. Unit must be the same as rs.
rs (int or float) – scale radius

Returns

NFW profile evaluated at the given distance. Unit is that of a 3D mass density in SI (i.e. kg/m^3).

Return type

float or ndarray[float]

Raises

TypeError – if r, c and rs are neither int, nor float
ValueError – if np.any(r)<0, if c<=0, or if rs<=0

galaxy.models.sersic_profile(r, n, re, Ie=None, bn=None, mag=None, offset=None)[source]¶

Code author: Wilfried Mercier - IRAP <wilfried.mercier@irap.omp.eu>

General Sersic profile defined as

\[\Sigma(r) = I_{\rm{e}} e^{\left [ \left (r/R_{\rm{e}} \right )^{1/n} - 1 \right ]},\]

with \(R_{\rm{e}}\) the effective radius, \(I_{\rm{e}}\) the surface brightness at the effective radius and \(n\) the Sersic index.

Note

Compute it with:

n, re and Ie

n, re, mag and offset

Parameters

r (float) – position at which the profile is computed
re (float) – half-light radius
bn (float) – (Optional) bn factor appearing in the Sersic profile. If None, its value will be computed.
Ie (float) – (Optional) surface brightness at half-light radius
mag (float) – (Optional) total integrated magnitude used to compute Ie if not given
offset (float) – (Optional) magnitude offset in the magnitude system used

Returns

surface brightness

Return type

float

2D profiles ¶

A collection of 2D models of surface brightness profiles with or without inclination effects and PSF convolution.

galaxy.models.bulgeDiskOnSky(nx, ny, Rd, Rb, x0=None, y0=None, Id=None, Ib=None, magD=None, magB=None, offsetD=None, offsetB=None, inclination=0, PA=0, combine=True, PSF={'FWHMX': 0.8, 'FWHMY': 0.8, 'name': 'Gaussian2D', 'sigmaX': None, 'sigmaY': None, 'unit': 'arcsec'}, noPSF=False, arcsecToGrid=0.03, fineSampling=1, samplingZone={'dx': 2, 'dy': 2, 'where': 'centre'}, skipCheck=False, verbose=True)[source]¶

Code author: Wilfried Mercier - IRAP <wilfried.mercier@irap.omp.eu>

Generate a bulge + (sky projected) disk 2D model (with PSF convolution).

How to use

Apart from the mandatory inputs, it is necessary to provide

an intensity at Re for each profile
a total magnitude value for each profile and a corresponding magnitude offset per profile (to convert from magnitudes to intensities)

Infos about sampling

fineSampling parameter can be used to rebin the data. The shape of the final image will depend on the samplingZone used

if the sampling is performed everywhere (‘where’ keyword in samplingZone equal to ‘all’), the final image will have dimensions (nx*fineSampling, ny*fineSampling)
if the sampling is performed around the centre (‘where’ equal to ‘centre’), the central part is over-sampled, but needs to be binned in the end so that pixels have the same size in the central part and around. Thus, the final image will have the dimension (nx, ny).

Warning

Rd and Rb should be given in pixel units.

If you provide them in arcsec, you must update the arcsecToGrid value to 1 (since 1 pixel will be equal to 1 arcsec).

Parameters

nx (int) – size of the model for the x-axis
ny (int) – size of the model for the y-axis
Rb (float) – bulge half-light radius. Best practice is to provide it in pixels.
Rd (float) – disk half-light radius. Best practice is to provide it in pixels.
arcsecToGrid (float) – (Optional) pixel size conversion in arcsec/pixel, used to convert the FWHM/sigma from arcsec to pixel
Ib (float) – (Optional) bulge intensity at (bulge) half-light radius. If not provided, magnitude and magnitude offset must be given instead.
Id (float) – (Optional) disk intensity at (disk) half-light radius. If not provided, magnitude and magnitude offset must be given instead.
inclination ((Optional) int or float) – (Optional) disk inclination on sky. Generally given between -90° and +90°. Value must be given in degrees.
magB (float) – (Optional) bulge total magnitude
magD (float) – (Optional) disk total magnitude
offsetB (float) – (Optional) bulge magnitude offset
offsetD (float) – (Optional) disk magnitude offset
noPSF (bool) – (Optional) whether to not perform PSF convolution or not
PA (int or float) – disk position angle (in degrees)
fineSampling (int(>0)) – fine sampling for the pixel grid used to make high resolution models. For instance, a value of 2 means that a pixel will be split into two subpixels.
PSF (dict) – (Optional) Dictionnary of the PSF (and its parameters) to use for the convolution. For now, only 2D Gaussians are accepted as PSF.
samplingZone (dict) –
where to perform the over sampling. Dictionnaries should have the following keys:
- ’where’ (type str) -> either ‘all’ to perform everywhere or ‘centre’ to perform around the centre
- ’dx’ (type int) -> x-axis maximum distance from the centre coordinate. A sub-array with x-axis values within [xpos-dx, xpos+dx] will be selected. If the sampling is performed everywhere, ‘dx’ does not need to be provided.
- ’dy’ (type int) -> y-axis maximum distance from the centre coordinate. A sub-array with y-axis values within [ypos-dy, ypos+dy] will be selected. If the sampling is performed everywhere, ‘dy’ does not need to be provided.
skipCheck (bool) – whether to skip the checking part or not
x0 (int or float) – x-axis centre position. Default is None so that nx//2 will be used.
y0 (int or float) – y-axis centre position. Default is None so that ny//2 will be used.
verbose (bool) – whether to print text on stdout or not

Returns

X, Y grids and the total (sky projected + PSF convolved) model of the bulge + disk decomposition

Return type

2D ndarray, 2D ndarray, 2D ndarray

galaxy.models.Sersic2D(nx, ny, listn, listRe, x0=None, y0=None, listIe=None, listMag=None, listOffset=None, listInclination=None, listPA=None, combine=True, fineSampling=1, samplingZone={'dx': 5, 'dy': 5, 'where': 'centre'}, skipCheck=False, verbose=True)[source]¶

Code author: Wilfried Mercier - IRAP <wilfried.mercier@irap.omp.eu>

Generate a (sky projected) 2D model (image) of a sum of Sersic profiles. Neither PSF smoothing, nor projections onto the sky whatsoever are applied here.

How to use

Apart from the mandatory inputs, it is necessary to provide

an intensity at Re for each profile
a total magnitude value for each profile and a corresponding magnitude offset per profile (to convert from magnitudes to intensities)

fineSampling parameter can be used to rebin the data. The shape of the final image will depend on the samplingZone used

if the sampling is performed everywhere (‘where’ keyword in samplingZone equal to ‘all’), the final image will have dimensions (nx*fineSampling, ny*fineSampling)
if the sampling is performed around the centre (‘where’ equal to ‘centre’), the central part is over-sampled, but needs to be binned in the end so that pixels have the same size in the central part and around. Thus, the final image will have the dimension (nx, ny).

Parameters

listn (list[int] or list[float]) – list of Sersic index for each profile
listRe (list[float]) – list of half-light radii for each profile
nx (int) – size of the model for the x-axis
ny (int) – size of the model for the y-axis
combine (bool) – (Optional) whether to combine (sum) all the components and return a single intensity map, or to return each component separately in lists
listIe (list[float]) – (Optional) list of intensities at re for each profile
listInclination (list[int] or list[float]) – (Optional) list of inclination of each Sersic component on the sky in degrees
listPA (list[int] or list[float]) – (Optional) list of position angle of each Sersic component on the sky in degrees. Generally, these values are given between -90° and +90°.
fineSampling (int(>0)) – (Optional) fine sampling for the pixel grid used to make high resolution models. For instance, a value of 2 means that a pixel will be split into two subpixels.
samplingZone (dict) –
(Optional) where to perform the sampling. Default is everywhere. Dictionnaries should have the following keys:
- ’where’ (type str) -> either ‘all’ to perform everywhere or ‘centre’ to perform around the centre
- ’dx’ (type int) -> x-axis maximum distance from the centre coordinate. An sub-array with x-axis values within [xpos-dx, xpos+dx] will be selected. If the sampling is performed everywhere, ‘dx’ does not need to be provided.
- ’dy’ (type int) -> y-axis maximum distance from the centre coordinate. An sub-array with y-axis values within [ypos-dy, ypos+dy] will be selected. If the sampling is performed everywhere, ‘dy’ does not need to be provided.
skipCheck (bool) – (Optional) whether to skip the checking part or not
verbose (bool) – (Optional) whether to print info on stdout or not
x0 (int or float) – (Optional) x-axis centre position. Default is None so that nx//2 will be used.
y0 (int or float) – (Optional) y-axis centre position. Default is None so that ny//2 will be used.

Returns

X, Y grids and the intensity map if combine is True
X, Y grids and a listof intensity maps for each component if combine is False

Raises

TypeError –
- if ‘dx’ and ‘dy’ keys are not in samplingZone
- if fineSampling, nx and ny are neither int, nor np.integer
ValueError –
- if ‘where’ key value in samplingZone is neither ‘all’, nor ‘centre’
- if nx, ny or arcsecToGrid are < 0
- if at least one n or one Re is < 0
- if fineSampling < 1
- if at least one PA is not in the range [-90, 90] deg
- if Ie and mag and offset are None

galaxy.models.bulge2D(nx, ny, Rb, x0=None, y0=None, Ib=None, mag=None, offset=None, inclination=0, PA=0, PSF={'FWHMX': 0.8, 'FWHMY': 0.8, 'name': 'Gaussian2D', 'sigmaX': None, 'sigmaY': None, 'unit': 'arcsec'}, noPSF=False, arcsecToGrid=0.03, fineSampling=1, samplingZone={'dx': 5, 'dy': 5, 'where': 'centre'}, skipCheck=False, verbose=True)[source]¶

Code author: Wilfried Mercier - IRAP <wilfried.mercier@irap.omp.eu>

Generate a 2D model for a de Vaucouleur bulge. This model can be sky projected and PSF convolved.

How to use

You must provide the size of the image with nx and ny as well as a bulge effective radius (usually in pixel unit) with Rb. Additionally, one must provide one of the following

surface brightness at Rb with Ib
total magnitude and magnitude offset to go from magnitude to Ie using the parameters mag and offset

The PSF convolution only accepts 2D Gaussians for now. You can either provide

a FWHM in the X and Y directions
a dispersion in the X and Y directions

You must also provide a unit for the FWHM or sigma values. This unit must be recognised by astropy. Since images have dimensions in pixel unit, if the FWHM or sigma values are not given in pixel unit, you must also provide arcsecToGrid to convert from physical unit to pixel unit.

If you do not want to convolve with the PSF, provide noPSF=True.

fineSampling parameter can be used to rebin the data. The shape of the final image will depend on the samplingZone used

if the sampling is performed everywhere (‘where’ keyword in samplingZone equal to ‘all’), the final image will have dimensions (nx*fineSampling, ny*fineSampling)
if the sampling is performed around the centre (‘where’ equal to ‘centre’), the central part is over-sampled, but needs to be binned in the end so that pixels have the same size in the central part and around. Thus, the final image will have the dimension (nx, ny).

Parameters

nx (int) – size of the model for the x-axis
ny (int) – size of the model for the y-axis
Rb (float) – bulge half-light radius. Best practice is to provide it in pixels.
x0 (int or float) – (Optional) x-axis centre position. Default is None so that nx//2 will be used.
y0 (int or float) – (Optional) y-axis centre position. Default is None so that ny//2 will be used.
Ib (float) – (Optional) bulge intensity at (bulge) half-light radius. If not provided, magnitude and magnitude offset must be given instead.
magB (float) – (Optional) bulge total magnitude
offsetB (float) – (Optional) bulge magnitude offset
inclination ((Optional) int or float) – (Optional) inclination on sky. Generally given between -90° and +90°. Value must be given in degrees.
PA (int or float) – position angle (in degrees)
PSF (dict) – (Optional) Dictionnary of the PSF (and its parameters) to use for the convolution. For now, only 2D Gaussians are accepted as PSF.
noPSF (bool) – (Optional) whether to not perform PSF convolution or not
arcsecToGrid (float) – (Optional) pixel size conversion in arcsec/pixel, used to convert the FWHM/sigma from arcsec to pixel
fineSampling (int(>0)) – fine sampling for the pixel grid used to make high resolution models. For instance, a value of 2 means that a pixel will be split into two subpixels.
samplingZone (dict) –
where to perform the over sampling. Dictionnaries should have the following keys:
- ’where’ (type str) -> either ‘all’ to perform everywhere or ‘centre’ to perform around the centre
- ’dx’ (type int) -> x-axis maximum distance from the centre coordinate. A sub-array with x-axis values within [xpos-dx, xpos+dx] will be selected. If the sampling is performed everywhere, ‘dx’ does not need to be provided.
- ’dy’ (type int) -> y-axis maximum distance from the centre coordinate. A sub-array with y-axis values within [ypos-dy, ypos+dy] will be selected. If the sampling is performed everywhere, ‘dy’ does not need to be provided.
skipCheck (bool) – whether to skip the checking part or not
verbose (bool) – (Optional) whether to print info on stdout or not

Returns

X coordinate array, Y coordinate array and the 2D bulge model

Return type

2D ndarray[float], 2D ndarray[float], 2D ndarray[float]

Example

Comparing different models:

The first does not use fineSampling, thus it shows an intensity profile, rather than a flux profile
The second one uses fineSampling=81 so that it matches more the ouput of Galfit
The last one also uses fineSampling but is convolved by a Gaussian with \({\rm{FWHM}} = 0.087~\rm{arcsec}\)

from   matplotlib.colors   import LogNorm
from   matplotlib          import rc
from   matplotlib.gridspec import GridSpec
from   wilfried.galaxy     import models as mod
import matplotlib.pyplot   as     plt
import matplotlib          as     mpl

# Bulge model using fine sampling and without PSF
X, Y, bulge1 = mod.bulge2D(100, 100, 35, mag=20, offset=30, noPSF=True)

# Bulge mode with fine sampling but without PSF
X, Y, bulge2 = mod.bulge2D(100, 100, 35, mag=20, offset=30, noPSF=True, fineSampling=81)

# Bulge model with fine sampling and with PSF convolution (FWHM=0.8 arcsec)
X, Y, bulge3 = mod.bulge2D(100, 100, 35, mag=20, offset=30, noPSF=False, fineSampling=81,
                          PSF={'name':'Gaussian2D', 'FWHMX':0.8, 'FWHMY':0.8, 'unit':'arcsec'}, arcsecToGrid=0.03)

###############################
#          Plot part          #
###############################

# Setup figure and axes
rc('font', **{'family': 'serif', 'serif': ['Times']})
rc('text', usetex=True)
mpl.rcParams['text.latex.preamble'] = r'\usepackage{newtxmath}'

f            = plt.figure(figsize=(9, 3))
gs           = GridSpec(1, 3, figure=f, wspace=0, hspace=0, left=0.2, right=0.8, top=0.9, bottom=0.3)

ax1          = f.add_subplot(gs[0])
ax2          = f.add_subplot(gs[1])
ax3          = f.add_subplot(gs[2])

ax1.set_title(r'No fine sampling',      size=15)
ax2.set_title(r'Fine sampling = $9^2$', size=15)
ax3.set_title(r'Fine sampling \& PSF',  size=15)

for a in [ax1, ax2, ax3]:
   a.set_xticklabels([])
   a.set_yticklabels([])
   a.tick_params(axis='x', which='both', direction='in')
   a.tick_params(axis='y', which='both', direction='in')
   a.yaxis.set_ticks_position('both')
   a.xaxis.set_ticks_position('both')

# Show bulges
ret1  = ax1.imshow(bulge1, origin='lower', norm=LogNorm(), cmap='plasma')
ret2  = ax2.imshow(bulge2, origin='lower', norm=LogNorm(), cmap='plasma')
ret3  = ax3.imshow(bulge3, origin='lower', norm=LogNorm(), cmap='plasma')

# Add colorbar
cb_ax = f.add_axes([0.2, 0.2, 0.6, 0.025])
cbar  = f.colorbar(ret1, cax=cb_ax, orientation='horizontal')
cbar.set_label(r'Surface brightness [arbitrary unit]', size=15)
cbar.ax.tick_params(labelsize=15)

plt.show()

(Source code, png, hires.png, pdf)

Comparison with Galfit models

We can also compare the ouput from this function with a similar model made with Galfit:

from   matplotlib.colors   import LogNorm, Normalize
from   matplotlib          import rc
from   matplotlib.gridspec import GridSpec
from   astropy.io          import fits
from   wilfried.galaxy     import models as mod
import matplotlib.pyplot   as     plt
import matplotlib          as     mpl
import numpy               as     np

##############################
#     Generate 2D models     #
##############################

X, Y, bulge1 = mod.bulge2D(100, 100, 15, mag=21, offset=30, noPSF=True,  fineSampling=81, samplingZone={'where':'centre', 'dx':5, 'dy':5})
X, Y, bulge2 = mod.bulge2D(100, 100, 15, mag=21, offset=30, noPSF=False, fineSampling=81, samplingZone={'where':'centre', 'dx':5, 'dy':5},
                          PSF={'name':'Gaussian2D', 'FWHMX':0.087, 'FWHMY':0.087, 'unit':'arcsec'}, arcsecToGrid=0.03)

#####################################
#     Load GALFIT output images     #
#####################################

with fits.open('/home/wilfried/wilfried_libs/galaxy/.test/bulge_withoutPSF.fits') as hdul:
   bulge_nopsf = hdul[0].data

with fits.open('/home/wilfried/wilfried_libs/galaxy/.test/bulge_withPSF.fits') as hdul:
   bulge_psf   = hdul[0].data

##########################
#       Plot parts       #
##########################

rc('font', **{'family': 'serif', 'serif': ['Times']})
rc('text', usetex=True)
mpl.rcParams['text.latex.preamble'] = r'\usepackage{newtxmath}'

f            = plt.figure(figsize=(9, 6))
gs           = GridSpec(2, 3, figure=f, wspace=0, hspace=0, left=0.1, right=0.9, top=0.95, bottom=0.15)

axs          = []
for i in gs:
   axs.append(f.add_subplot(i))

axs[0].set_ylabel(r'without PSF', size=15)
axs[0].set_title(r'bulge2D',      size=15)
axs[1].set_title(r'Galfit bulge', size=15)
axs[2].set_title(r'Residuals',    size=15)
axs[3].set_ylabel(r'with PSF',    size=15)


for a in axs:
   a.set_xticklabels([])
   a.set_yticklabels([])

   a.tick_params(axis='x', which='both', direction='in')
   a.tick_params(axis='y', which='both', direction='in')
   a.yaxis.set_ticks_position('both')
   a.xaxis.set_ticks_position('both')

# First line
mmax  = np.nanmax([bulge1, bulge2, bulge_nopsf, bulge_psf])
mmin  = np.nanmin([bulge1, bulge2, bulge_nopsf, bulge_psf])

diff1 = 100 - 100*bulge_nopsf/bulge1
diff2 = 100 - 100*bulge_psf/bulge2
dmax  = np.nanmax([diff1, diff2])
dmin  = np.nanmin([diff1, diff2])

ret11 = axs[0].imshow(bulge1,      origin='lower', cmap='plasma',  norm=LogNorm(  vmin=mmin, vmax=mmax))
ret12 = axs[1].imshow(bulge_nopsf, origin='lower', cmap='plasma',  norm=LogNorm(  vmin=mmin, vmax=mmax))
ret13 = axs[2].imshow(diff1,       origin='lower', cmap='viridis', norm=Normalize(vmin=dmin, vmax=dmax))

# Second line
ret21 = axs[3].imshow(bulge2,      origin='lower', cmap='plasma',  norm=LogNorm(  vmin=mmin, vmax=mmax))
ret22 = axs[4].imshow(bulge_psf,   origin='lower', cmap='plasma',  norm=LogNorm(  vmin=mmin, vmax=mmax))
ret23 = axs[5].imshow(diff2,       origin='lower', cmap='viridis', norm=Normalize(vmin=dmin, vmax=dmax))

# Add colorbar 1
cb_ax1 = f.add_axes([0.1, 0.1, 0.53, 0.025])
cbar1  = f.colorbar(ret11, cax=cb_ax1, orientation='horizontal')
cbar1.set_label(r'Surface brightness [arbitrary unit]', size=15)
cbar1.ax.tick_params(labelsize=15)

# Add colorbar 2
cb_ax2 = f.add_axes([0.64, 0.1, 0.26, 0.025])
cbar2  = f.colorbar(ret13, cax=cb_ax2, orientation='horizontal')
cbar2.set_label(r'Relative difference (\%)', size=15)
cbar2.ax.tick_params(labelsize=15)

plt.show()

(Source code, png, hires.png, pdf)

Galaxy models¶

1D profiles¶

2D profiles¶

1D profiles ¶

2D profiles ¶