JamPy: Jeans Anisotropic Models for Galactic Dynamics

## Project description

## The JamPy Package

**Jeans Anisotropic Modelling for Galactic Dynamics**

JamPy is a Python implementation of the Jeans Anisotropic Modelling (JAM) formalism for the dynamical modelling of galaxies.

This software can be used e.g. to measure the mass of supermassive black holes in galaxies, to infer their dark-matter content or to measure galaxy masses and density profiles.

The method calculates all the first and second velocity moments, for both the intrinsic and the projected kinematics, in spherical and axisymmetric geometry.

The JAM solution assuming a cylindrically-oriented velocity ellipsoid was introduced in Cappellari (2008), while the solution assuming a spherically-oriented velocity ellipsoid was introduced in Cappellari (2020)

### Attribution

If you use this software for your research, please cite Cappellari (2008) for the cylindrically-aligned JAM solution and Cappellari (2020) for the spherically-aligned JAM solution.

The BibTeX entry for the two main JAM papers are respectively:

@ARTICLE{Cappellari2008, author = {{Cappellari}, Michele}, title = "{Measuring the inclination and mass-to-light ratio of axisymmetric galaxies via anisotropic Jeans models of stellar kinematics}", journal = {MNRAS}, eprint = {0806.0042}, year = 2008, volume = 390, pages = {71-86}, doi = {10.1111/j.1365-2966.2008.13754.x} } @ARTICLE{Cappellari2020, author = {{Cappellari}, Michele}, title = "{Efficient solution of the anisotropic spherically-aligned axisymmetric Jeans equations of stellar hydrodynamics for galactic dynamics}", journal = {MNRAS}, eprint = {1907.09894}, year = 2020, volume = 494, pages = {4819-4837}, doi = {10.1093/mnras/staa959} }

### Installation

install with:

pip install jampy

Without writing access to the global site-packages directory, use:

pip install --user jampy

To upgrade JamPy to the latest version use:

pip install --upgrade jampy

### Documentation

Full documentation is contained in the individual files docstrings.

Usage examples are contained in the directory jampy/examples which is copied by pip within the global folder site-packages.

What follows is the documentation of the two main procedures of the JamPy package, extracted from their Python docstrings. The other procedures are documented in their respective docstrings.

The older JamPy routines from version < v6.0 are generally redundant, but can be imported from jampy.legacy and used like before.

## jam_axi_proj

### Purpose

This procedure calculates a prediction for all the projected first or second velocity moments for an anisotropic (three-integral) axisymmetric galaxy model.

Any of the three components of the first velocity moment or any of the six components of the symmetric velocity dispersion tensor are supported. These include the line-of-sight velocities and the components of the proper motion.

Two assumptions for the orientation of the velocity ellipsoid are supported:

The cylindrically-aligned (R, z, phi) solution was presented in Cappellari (2008)

The spherically-aligned (r, th, phi) solution was presented in Cappellari (2020)

### Calling Sequence

```
jam = jam_axi_proj(
surf_lum, sigma_lum, qobs_lum, surf_pot, sigma_pot, qobs_pot,
inc, mbh, distance, xbin, ybin, align='cyl', analytic_los=True,
beta=None, data=None, epsrel=1e-2, errors=None, flux_obs=None,
gamma=None, goodbins=None, interp=True, kappa=None, ml=None,
moment='zz', nang=10, nlos=150, nodots=False, normpsf=1.,
nrad=20, pixang=0., pixsize=0., plot=True, quiet=False,
rbh=0.01, sigmapsf=0., step=0., vmax=None, vmin=None)
vrms = jam.model # with moment='zz' the output is the LOS Vrms
jam.plot() # Generate data/model comparison when data is given
```

See more examples in the jampy/examples folder inside site-packages.

### Input Parameters

- surf_lum: array_like with shape (n,)
peak surface values of the Multi-Gaussian Expansion (MGE) Gaussians describing the surface brightness of the tracer population for which the kinematics is derived.

The units are arbitrary as they cancel out in the final results.

EXAMPLE: when one obtains the kinematics from optical spectroscopy, surf_lum contains the galaxy optical surface brightness, which has typical units of Lsun/pc^2 (solar luminosities per parsec^2).

- sigma_lum: array_like with shape (n,)
dispersion (sigma) in arcseconds of the MGE Gaussians describing the distribution of the kinematic-tracer population.

- qobs_lum: array_like with shape (n,)
observed axial ratio (q’) of the MGE Gaussians describing the distribution of the kinematic-tracer population.

- surf_pot: array_like with shape (m,)
peak value of the MGE Gaussians describing the galaxy total-mass surface density in units of Msun/pc^2 (solar masses per parsec^2). This is the MGE model from which the model gravitational potential is computed.

EXAMPLE: with a self-consistent model, one has the same Gaussians for both the kinematic-tracer and the gravitational potential. This implies surf_pot = surf_lum, sigma_pot = sigma_lum and qobs_pot = qobs_lum. The global M/L of the model is fitted by the routine when passing the data and errors keywords with the observed kinematics.

- sigma_pot: array_like with shape (m,)
dispersion in arcseconds of the MGE Gaussians describing the galaxy total-mass surface density.

- qobs_pot: array_like with shape (m,)
observed axial ratio of the MGE Gaussians describing the galaxy total-mass surface density.

- inc: float
inclination in degrees between the line-of-sight and the galaxy symmetry axis (0 being face-on and 90 edge-on).

- mbh: float
Mass of a nuclear supermassive black hole in solar masses.

IMPORTANT: The model predictions are computed assuming surf_pot gives the total mass. In the self-consistent case, one has surf_pot = surf_lum and if requested (keyword ml) the program can scale the output model to best fit the data. The scaling is equivalent to multiplying

*both*surf_pot and mbh by a factor M/L. To avoid mistakes, the actual mbh used by the output model is printed on the screen.- distance: float
the distance of the galaxy in Mpc. When the distance is derived from redshift one should use the angular diameter distance D_A here.

- xbin: array_like with shape (p,)
X coordinates in arcseconds of the bins (or pixels) at which one wants to compute the model predictions. The X-axis is assumed to coincide with the galaxy projected major axis. The galaxy center is at (0,0).

In general the coordinates (xbin, ybin) have to be rotated to bring the galaxy major axis on the X-axis, before calling jam_axi_proj.

When no PSF/pixel convolution is performed (sigmapsf=0 or pixsize=0) there is a singularity at (0,0) which should be avoided by the user in the input coordinates.

- ybin: array_like with shape (p,)
Y coordinates in arcseconds of the bins (or pixels) at which one wants to compute the model predictions. The Y-axis is assumed to coincide with the projected galaxy symmetry axis.

### Optional Keywords

- align: {‘cyl’, ‘sph’}, optional.
Assumed alignment for the velocity ellipsoid during the solution of the Jeans equations.

align='cyl' assumes a cylindrically-aligned velocity ellipsoid using the solution of Cappellari (2008)

align='sph' assumes a spherically-aligned velocity ellipsoid using the solution of Cappellari (2020)

- analytic_los: bool, optional
This is True (default) if the line-of-sight integral is performed analytically and False if it is done via numerical quadrature.

An analytic integral is only possible with align='cyl' and only for the second velocity moments. For this reason, when comparing the two second-moment solutions with align='cyl' and align='sph', it may be preferable to set analytic_los=False to ensure that numerical interpolation error is exactly the same in both cases.

When align='sph', or when the user requests a first velocity moment, this keyword is automatically set to False.

- beta: array_like with shape (n,)
radial anisotropy of the individual kinematic-tracer MGE Gaussians (Default: beta=np.zeros(n)):

beta = 1 - (sigma_th/sigma_r)^2 # with align=`sph` beta = 1 - (sigma_z/sigma_R)^2 # with align=`cyl`

- data: array_like with shape (p,), optional
observed first or second velocity moment used to fit the model.

EXAMPLE: In the common case where one has only line-of-sight velocities the second moment is given by:

Vrms = np.sqrt(velBin**2 + sigBin**2)

at the coordinates positions given by the vectors xbin and ybin.

If data is set and ml is negative or None, then the model is fitted to the data, otherwise, the adopted ml is used and just the chi**2 is returned.

- epsrel: float, optional
Relative error requested for the numerical computation of the intrinsic moments (before line-of-sight quadrature). (Default: epsrel=1e-2)

- errors: array_like with shape (p,), optional
1sigma uncertainty associated with the data measurements.

EXAMPLE: In the case where the data are given by the Vrms = np.sqrt(velBin**2 + sigBin**2), from the error propagation:

errors = np.sqrt((dVel*velBin)**2 + (dSig*sigBin)**2)/Vrms,

where velBin and sigBin are the velocity and dispersion in each bin and dVel and dSig are the corresponding 1sigma uncertainties. (Default: constant errors = 0.05*np.median(data))

- flux_obs: array_like with shape (p,), optional
Optional mean surface brightness of each bin for plotting.

- gamma: array_like with shape (n,)
tangential anisotropy of the individual kinematic-tracer MGE Gaussians (Default: gamma=np.zeros(n)):

gamma = 1 - (sigma_phi/sigma_r)^2 # with align=`sph` gamma = 1 - (sigma_phi/sigma_R)^2 # with align=`cyl`

- goodbins: array_like with shape (p,)
Boolean vector with values True for the bins/spaxels which have to be included in the fit (if requested) and in the chi**2 calculation. (Default: fit all bins).

- interp: bool, optional
This keyword is for advanced use only! Set interp=False to force no-interpolation on the sky plane. In this way jam.vel and jam.vel2 contain all the first and second velocity moments at the input coordinates (xbin, ybin), without PSF convolution. By default interp=True and one should generally not change this.

IMPORTANT: If sigmapsf=0 or pixsize=0 or interp=False then PSF convolution is not performed.

This keyword is mainly useful for testing against analytic results or to compute all moments, including proper motions, simultaneously.

- kappa: float, optional
When kappa=None (default) the first velocity moments are scaled in such a way that the projected angular momentum of the data and model is the same [equation 52 of Cappellari (2008)]. When kappa=1 the model first velocity moments are output without any scaling.

- ml: float, optional
Mass-to-light ratio (M/L) to multiply the values given by surf_pot. Setting this keyword is completely equivalent to multiplying the output model by np.sqrt(M/L) after the fit. This implies that the BH mass is also scaled and becomes mbh*ml.

If ml=None (default) the M/L is fitted from the data and the best-fitting M/L is returned in output. The BH mass of the model is also scaled and becomes mbh*ml.

- moment: {‘x’, ‘y’, ‘z’, ‘xx’, ‘yy’, ‘zz’, ‘xy’, ‘xz’, ‘yz’}, optional
String specifying the component of the velocity first or second moments requested by the user in output. All values ar in km/s.

moment='x' gives the first moment <V_x'> of the proper motion in the direction orthogonal to the projected symmetry axis.

moment='y' gives the first moment <V_y'> of the proper motion in the direction parallel to the projected symmetry axis.

moment='z' gives the first moment Vlos = <V_z'> of the line-of-sight velocity.

moment='xx' gives sqrt<V_x'^2> of the component of the proper motion dispersion tensor in the direction orthogonal to the projected symmetry axis.

moment='yy' gives sqrt<V_y'^2> of the component of the proper motion dispersion tensor in the direction parallel to the projected symmetry axis.

moment='zz' (default) gives the usual line-of-sight Vrms = sqrt<V_z'^2>.

moment='xy' gives the mixed component <V_x'V_y'> of the proper motion dispersion tensor.

moment='xz' gives the mixed component <V_x'V_z'> of the proper motion dispersion tensor.

moment='yz' gives the mixed component <V_y'V_z'> of the proper motion dispersion tensor.

- nang: int, optional
The number of linearly-spaced intervals in the eccentric anomaly at which the model is evaluated before interpolation and PSF convolution. (default: nang=10)

- nlos: int (optional)
Number of values used for the numerical line-of-sight quadrature. (default nlos=150)

- nodots: bool, optional
Set to True to hide the dots indicating the centers of the bins in the linearly-interpolated two-dimensional map (default False).

- normpsf: array_like with shape (q,)
fraction of the total PSF flux contained in the circular Gaussians describing the PSF of the kinematic observations. The PSF will be used for seeing convolution of the model kinematics. It has to be np.sum(normpsf) = 1.

- nrad: int, optional
The number of logarithmically spaced radial positions at which the model is evaluated before interpolation and PSF convolution. One may want to increase this value if the model has to be evaluated over many orders of magnitude in radius (default: nrad=20).

- pixang: float, optional
Angle between the observed spaxels and the galaxy major axis X. This angle only rotates the spaxels around their centers,

*not*the whole coordinate system (xbin, ybin), which must be rotated independently by the user before calling jam_axi_proj. Using the keyword is generally unnecessary.- pixsize: float, optional
Size in arcseconds of the (square) spatial elements at which the kinematics is obtained. This may correspond to the side of the spaxel or lenslets of an integral-field spectrograph. This size is used to compute the kernel for the seeing and aperture convolution.

IMPORTANT: If sigmapsf=0 or pixsize=0 or interp=False then PSF convolution is not performed.

- plot: bool
When data is not None setting this keyword produces a plot with the data/model comparison at the end of the calculation.

- quiet: bool
Set this keyword to avoid printing values on the console.

- rbh: float, optional
This scalar gives the sigma in arcsec of the Gaussian representing the central black hole of mass MBH [See Section 3.1.2 of Cappellari (2008)] The gravitational potential is indistinguishable from a point source for radii > 2*rbh, so the default rbh=0.01 arcsec is appropriate in most current situations.

When using different units as input, e.g. pc instead of arcsec, one should check that rbh is not too many order of magnitude smaller than the spatial resolution of the data.

- sigmapsf: array_like with shape (q,)
dispersion in arcseconds of the circular Gaussians describing the PSF of the kinematic observations.

IMPORTANT: If sigmapsf=0 or pixsize=0 or interp=False then PSF convolution is not performed.

IMPORTANT: PSF convolution is done by creating a 2D image, with pixels size given by step=max(sigmapsf, pixsize/2)/4, and convolving it with the PSF + aperture. If the input radii are very large compared to step, the 2D image may require a too large amount of memory. If this is the case one may compute the model predictions at small radii with a first call to jam_axi_proj with PSF convolution, and the model predictions at large radii with a second call to jam_axi_proj without PSF convolution.

- step: float, optional
Spatial step for the model calculation and PSF convolution in arcsec. This value is automatically computed by default as step=max(sigmapsf,pixsize/2)/4. It is assumed that when pixsize or sigmapsf are large, high-resolution calculations are not needed. In some cases, however, e.g. to accurately estimate the central Vrms in a very cuspy galaxy inside a large aperture, one may want to override the default value to force smaller spatial pixels using this keyword.

- vmax: float, optional
Maximum value of the data to plot.

- vmin: float, optional
Minimum value of the data to plot.

### Output Parameters

Stored as attributes of the jam_axi_proj class.

- .chi2: float
Reduced chi**2, namely per degree of freedom, describing the quality of the fit:

d, m = (data/errors)[goodbins], (model/errors)[goodbins] chi2 = ((d - m)**2).sum()/goodbins.sum()

When no data are given in input, this is returned as np.nan.

- .flux: array_like with shape (p,)
PSF-convolved MGE surface brightness of each bin in Lsun/pc^2, used to plot the isophotes of the kinematic-tracer on the model results.

- .kappa: float
Ratio by which the model was scaled to fit the observed velocity [defined by equation 52 of Cappellari (2008)]

- .ml: float
Best fitting M/L by which the mass was scaled to fit the observed moments.

- .model: array_like with shape (p,)
Model predictions for the selected velocity moments for each input bin (xbin, ybin). This attribute is the main output from the program.

Any of the six components of the symmetric proper motion dispersion tensor {'xx', 'yy', 'zz', 'xy', 'xz', 'yz'}, or any of the three first velocity moments {'x', 'y', 'z'}` can be returned in output. The desired model output is selected using the moment keyword. See the moment documentation for details.

- .vel: array_like with shape (3, p)
This attribute generally contains an intermediate result of the calculation and should not be used. Instead, the output kinematic model predictions are contained in the .model attribute.

However, for advanced use only, when setting interp=False and analytic_los=False, this attribute contains the first velocity moments for all the x, y and z components,

*not*PSF convolved, at the sky coordinates (xbin, ybin).- .vel2: array_like with shape (3, 3, p)
This attribute generally contains an intermediate result of the calculation and should not be used. Instead, the output kinematic model predictions are contained in the .model attribute.

However, for advanced use only, when setting interp=False and analytic_los=False, this attribute contains the full 3x3 second velocity moment tensor,

*not*PSF convolved, at the sky coordinates (xbin, ybin).

## jam_axi_intr

### Purpose

This procedure calculates all the intrinsic first and second velocity moments for an anisotropic axisymmetric galaxy model.

This program is useful e.g. to model the kinematics of galaxies like our Milky Way, for which the intrinsic moments can be observed directly, or to compute starting conditions for N-body numerical simulations of galaxies.

Two assumptions for the orientation of the velocity ellipsoid are supported:

The cylindrically-aligned (R, z, phi) solution was presented in Cappellari (2008)

The spherically-aligned (r, th, phi) solution was presented in Cappellari (2020)

### Calling Sequence

```
jam = jam_axi_intr(
dens_lum, sigma_lum, qintr_lum, dens_pot, sigma_pot, qintr_pot,
mbh, Rbin, zbin, align='cyl', beta=None, data=None,
epsrel=1e-2, errors=None, gamma=None, goodbins=None,
interp=True, ml=None, nang=10, nodots=False, nrad=20,
plot=True, proj_cyl=False, quiet=False, rbh=1)
# The meaning of the output is different depending on `align`
sig2R, sig2z, sig2phi, v2phi = jam.model # with align='cyl'
sig2r, sig2th, sig2phi, v2phi = jam.model # with align='sph'
jam.plot() # Generate data/model comparison
```

### Input Parameters

- dens_lum: array_like with shape (n,)
vector containing the peak value of the MGE Gaussians describing the intrinsic density of the tracer population for which the kinematics is derived. The units are arbitarary as they cancel out in the final results. Typical units are e.g. Lsun/pc^3 (solar luminosities per parsec^3)

- sigma_lum: array_like with shape (n,)
vector containing the dispersion (sigma) in pc of the MGE Gaussians describing the galaxy kinematic-tracer population.

- qintr_lum: array_like with shape (n,)
vector containing the intrinsic axial ratio (q) of the MGE Gaussians describing the galaxy kinematic-tracer population.

- surf_pot: array_like with shape (m,)
vector containing the peak value of the MGE Gaussians describing the galaxy total-mass density in units of Msun/pc^3 (solar masses per parsec^3). This is the MGE model from which the model gravitational potential is computed.

- sigma_pot: array_like with shape (m,)
vector containing the dispersion in pc of the MGE Gaussians describing the galaxy total-mass density.

- qintr_pot: array_like with shape (m,)
vector containing the intrinsic axial ratio of the MGE Gaussians describing the galaxy total-mass density.

- mbh: float
Mass of a nuclear supermassive black hole in solar masses.

- Rbin: array_like with shape (p,)
Vector with the R coordinates in pc of the bins (or pixels) at which one wants to compute the model predictions. This is the first cylindrical coordinate (R, z) with the galaxy center at (0,0).

There is a singularity at (0, 0) which should be avoided by the user in the input coordinates.

- zbin: array_like with shape (p,)
Vector with the z coordinates in pc of the bins (or pixels) at which one wants to compute the model predictions. This is the second cylindrical coordinate (R, z), with the z-axis coincident with the galaxy symmetry axis.

### Optional Keywords

- align: {‘cyl’, ‘sph’} optional
If align='cyl' the program computes the solution of the Jeans equations with cylindrically-aligned velocity ellipsoid, presented in Cappellari (2008). If align='sph' the spherically-aligned solution of Cappellari (2020) is returned.

- beta: array_like with shape (n,)
Vector with the axial anisotropy of the individual kinematic-tracer MGE Gaussians (Default: beta=np.zeros(n)):

beta = 1 - (sigma_th/sigma_r)^2 # with align=`sph` beta = 1 - (sigma_z/sigma_R)^2 # with align=`cyl`

- data: array_like of shape (4, p), optional
Four input vectors with the observed values of:

[sigR, sigz, sigphi, vrms_phi] in km/s, when align='cyl' (or align='sph' and proj_cyl=True).

vrms_phi is the square root of the velocity second moment in the tangential direction. If the velocities vphi_j are measured from individual stars then vrms_phi = sqrt(mean(vphi_j^2)). One can also use the relation vrms_phi = sqrt(vphi^2 + sigphi^2), where vphi = mean(vphi_j) and sigphi = std(vphi_j)

[sigr, sigth, sigphi, vrms_phi] in km/s, when align='sph', where vrms_phi is defined above.

- epsrel: float, optional
Relative error requested for the numerical quadratures, before interpolation (Default: epsrel=1e-2).

- errors: array_like of shape (4, p), optional
1sigma uncertainties on data, in the same format (default 5 km/s).

- gamma: array_like with shape (n,)
Vector with the tangential anisotropy of the individual kinematic-tracer MGE Gaussians (Default: gamma=np.zeros(n)):

gamma = 1 - (sigma_phi/sigma_r)^2 # with align=`sph` gamma = 1 - (sigma_phi/sigma_R)^2 # with align=`cyl`

- goodbins: array_like with shape (4, p), optional
Boolean vector of the same shape as data with values True for the bins which have to be included in the fit (if requested) and chi^2 calculation (Default: fit all bins).

- interp: bool, optional
If interp=False no interpolation is performed and the model is computed at every set of input (R, z) coordinates. If interp=True (default), the model is interpolated if the number of requested input (R, z) coordinates is larger than nang*nrad.

- ml: float, optional
Mass-to-light ratio M/L. If ml=None (default) the M/L is fitted to the data and the best-fitting value is returned in output. The mbh is also scaled and becomes mbh*ml. If ml=1 no scaling is applied to the model.

- nang: int, optional
The number of linearly-spaced intervals in the eccentric anomaly at which the model is evaluated before interpolation (default: nang=10).

- nodots: bool, optional
Set to True to hide the dots indicating the centers of the bins in the two-dimensional map (default False).

- nrad: int, optional
The number of logarithmically spaced radial positions at which the model is evaluated before interpolation. One may want to increase this value if the model has to be evaluated over many orders of magnitude in radius (default: nrad=20).

- plot: bool, optional
If plot=True (default) and data is not None, produce a plot of the data-model comparison at the end of the calculation.

- proj_cyl: bool, optional
If align='sph' and proj_cyl=True, the function projects the spherically-aligned moments to cylindrical coordinates and returns the [sig2R, sig2z, sig2phi, v2phi] components as in the case align='cyl'. This is useful for a direct comparison of results with either the spherical or cylindrical alignment, as it allows one to fit the same data with both modelling assumptions.

- quiet: bool, optional
If quiet=False (default), print the best-fitting M/L and chi2 at the end for the calculation.

- rbh: float, optional
This scalar gives the sigma in pc of the Gaussian representing the central black hole of mass mbh [See Section 3.1.2 of Cappellari (2008)]. The gravitational potential is indistinguishable from a point source for radii > 2*rbh, so the default rbh=1 pc is appropriate for observations taken with current telescopes.

### Output Parameters

Returned as attributes of the jam_axi_intr class.

- .chi2: float
Reduced chi^2 (chi^2/DOF) describing the quality of the fit:

d = (data/errors)[goodbins] m = (model/errors)[goodbins] chi2 = ((d - m)**2).sum()/goodbins.sum()

- .flux: array_like with shape (p,)
Vector with the MGE luminosity density at each (R, z) location in Lsun/pc^3, used to plot the isophotes on the model results.

- .ml: float
Best fitting M/L. This value is fitted while ignoring sigphi and it is strictly independent of the adopted tangential anisotropy gamma.

- .model: array_like with shape (4, p)
Contains [sig2R, sig2z, sig2phi, v2phi] with align='cyl'

Contains [sig2r, sig2th, sig2phi, v2phi] with align='sph'

where the above quantities are defined as follows:

- sig2R (sig2r): array_like with shape (p,)
squared intrinsic dispersion in (km/s)^2 along the R (r) direction at each (R, z) location.

- sig2z (sig2th): array_like with shape (p,)
squared intrinsic dispersion in (km/s)^2 along the z (th) direction at each (R, z) location.

- sig2phi: array_like with shape (p,)
squared intrinsic dispersion in (km/s)^2 along the tangential phi direction at each (R, z) location.

- v2phi: array_like with shape (p,)
the second velocity moment in (km/s)^2 along the tangential phi direction at each (R, z) location.

The mean velocity along the tangential direction can be computed as vphi = np.sqrt(v2phi - sig2phi)

NOTE: I return squared velocities instead of taking the square root, to allow for negative values (unphysical solutions).

## License

Other/Proprietary License

Copyright (c) 2003-2022 Michele Cappellari

This software is provided as is without any warranty whatsoever. Permission to use, for non-commercial purposes is granted. Permission to modify for personal or internal use is granted, provided this copyright and disclaimer are included in all copies of the software. All other rights are reserved. In particular, redistribution of the code is not allowed.

## Changelog

- V6.4.0: MC, Oxford, 3 October 2022
jam_sph_proj: Created this new function by renaming the procedure legacy.jam_sph_rms and changing its interface to be consistent with the axisymmetric version.

jam_sph_proj: Included special isotropic formula for testing.

jam_sph_proj: Included Osipkov-Merritt anisotropy for testing.

jam_sph_proj: Made quadrature limits insensitive to scaling.

jam_sph_proj: Simplified integrand with formulas of Cappellari (2020) and using recurrence relations to reduce calls to special functions.

jam_sph_proj: More efficient TANH transformation of the integration variable following Cappellari (2020).

jam_sph_intr: New function to compute the intrinsic moments in spherical symmetry.

jam_axi_proj: Removed fixed minimum radius limit in pc for the interpolation without PSF convolution. This avoids the risk of artificial truncation when using small arbitrary spatial coordinates for testing.

New procedure examples.jam_dark_halo_bayes_example.py.

Renamed quadva as quad1d with modified interface and new singular keyword to skip transforming the integration variable.

- V6.3.3: MC, Oxford, 7 July 2021
jam_axi_proj: Clarified meaning of interp keyword in dosctring. Thanks to Kai Zhu (nao.cas.cn) for the feedback.

jam_axi_proj: print “No PSF/pixel convolution” when interp == False.

- V6.3.2: MC, Oxford, 28 April 2021
Use the new jam_axi_proj instead of legacy software in the examples.

Removed redundant legacy examples.

- V6.3.1: MC, Oxford, 11 November 2020
jam_axi_proj: New keyword analytic_los to chose between numeric or analytic line-of-sight integral for the second velocity moment, when align='cyl'.

jam_axi_proj: Increased default value of nlos keyword.

jam_axi_proj: Raise an error if rbh is too small.

jam_axi_proj and jam_axi_intr: Removed **kwargs argument and included new nodots keyword passed to plot_velfield.

- V6.2.1: MC, Oxford, 15 September 2020
jam_axi_proj: Fixed program stop when data == ml == None. Thank to Bitao Wang (pku.edu.cn) for reporting.

- V6.2.0: MC, Oxford, 17 August 2020
jam_axi_proj: Avoid possible division by zero after convolution, when the tracer MGE is much smaller than the field of view.

jam_axi_proj: Fully broadcasted vmom_proj.

jam_axi_proj: Removed minimum-radius clipping in vmom_proj.

jam_axi_proj: New interp keyword to force no-interpolation when using the full first and second velocity moments simultaneously.

Made jam.plot() callable after jam_axi_proj or jam_axi_intr.

New axisymmetric analytic vs MGE test in mge_vcirc_example.

mge_vcirc: Upgraded formalism.

Fixed Numpy 1.9 VisibleDeprecationWarning.

Updated documentation.

- V6.1.5: MC, Oxford, 23 July 2020
Fixed program stop in first velocity moment without input data, introduced in V6.1.2. Thanks to Bitao Wang (pku.edu.cn) for reporting.

Implemented the kappa input keyword as scalar.

- V6.1.4: MC, Oxford, 16 July 2020
Added kappa to the returned parameters of jam_axi_proj.

Compute both velocity and Vrms in jam_axi_proj_example.

- V6.1.3: MC, Oxford, 13 July 2020
Fixed program stop in legacy.jam_axi_vel due to a variable name typo introduced in V6.1.2.

- V6.1.2: MC, Oxford, 20 June 2020
jam_axi_proj: Fixed input ml being ignored. Thanks to Sabine Thater (univie.ac.at) and Takafumi Tsukui (grad.nao.ac.jp) for reporting.

jam_axi_rms: I reduced the interpolation error before the PSF convolution for all the rotines in the legacy sub-folder, as already implemented in the new jam_axi_proj. Thanks to Takafumi Tsukui (grad.nao.ac.jp) for reporting differences.

jam_axi_intr: Request input data = [sigR, sigz, sigphi, vrms_phi] instead of data = [sigR, sigz, sigphi, vphi].

jam_axi_intr: exclude sigphi from ml fitting. These two changes make the fitted ml strictly independent of the adopted tangential anisotropy gamma.

- V6.0.1: MC, Oxford, 23 April 2020
Fixed model output when fitting ml. Thanks to Selina Nitschai (mpia-hd.mpg.de) for reporting.

- V6.0.0: MC, Oxford, 22 April 2020
Major changes to the whole jampy package: from this version I include the new spherically-aligned solution of the Jeans equations from Cappellari (2020, MNRAS).

Two new functions jam_axi_intr and jam_axi_proj now provide either the intrinsic or the projected moments, respectively, for both the spherically-aligned and cylindrically-aligned JAM solutions.

I moved the previous procedures jam_axi_rms, jam_axi_vel and jam_sph_rms to the jampy.legacy folder.

- V5.0.23: MC, Oxford, 31 October 2019
Use analytic mge_surf in convolution.

- V5.0.22: MC, Oxford, 21 March 2019
Reformatted documentation of all procedures.

- V5.0.21: MC, Oxford, 14 February 2019
Significant speedup of mge_vcirc.

Formatted documentation.

Created package-wide CHANGELOG: before this version, the CHANGELOG file only refers to the procedure jam_axi_rms.

- V5.0.16: MC, Oxford, 27 September 2018
Fixed clock DeprecationWarning in Python 3.7.

- V5.0.15: MC, Oxford, 12 May 2018
Dropped Python 2.7 support.

- V5.0.14: MC, Oxford, 17 April 2018
Fixed MatplotlibDeprecationWarning in Matplotlib 2.2.

Changed imports for jam as a package.

Removed example.

- V5.0.13: MC, Oxford, 7 March 2018
Check that PSF is normalized.

- V5.0.12: MC, Oxford, 22 January 2018
Print a message when no PSF convolution was performed.

Broadcast kernel and MGE convolution loops.

Fixed missing tensor in assertion test.

- V5.0.11: MC, Oxford, 10 September 2017
Make default step depend on sigmapsf regardless of pixsize.

- V5.0.10: MC, Oxford, 10 August 2017
Raise an error if goodbins is all False.

- V5.0.9: MC, Oxford, 17 March 2017
Included flux_obs keyword. Updated documentation.

Fixed DeprecationWarning in Numpy 1.12.

- V5.0.8: MC, Oxford, 17 February 2017
Use odd kernel size for convolution.

Fixed corner case with coordinates falling outside the interpolation region, due to finite machine precision.

- V5.0.7: MC, Oxford, 23 February 2016
Scale rmsModel by the input M/L also when rms is not given. Thanks to Alex Grainger (Oxford) for pointing out the inconsistency.

Pass **kwargs for plotting.

- V5.0.6: MC, Oxford, 18 September 2015
Plot bad bins on the data.

- V5.0.5: MC, Oxford, 23 May 2015
Changed meaning of goodbins to be a boolean vector.

- V5.0.4: MC, Sydney, 5 February 2015
Introduced further checks on matching input sizes.

- V5.0.3: MC, Oxford, 31 October 2014
Modified final plot layout.

- V5.0.2: MC, Oxford, 25 May 2014
Support both Python 2.7 and Python 3.

- V5.0.1: MC, Oxford, 24 February 2014
Plot bi-symmetrized V_rms as in IDL version.

- V5.0.0: MC, Paranal, 11 November 2013
Translated from IDL into Python.

- V4.1.5: MC, Paranal, 8 November 2013
Use renamed CAP_* routines to avoid potential naming conflicts.

- V4.1.4: MC, Oxford, 12 February 2013
Include _EXTRA and RANGE keywords for plotting.

- V4.1.3: MC, Oxford, 1 February 2013
Output FLUX in Lsun/pc^2.

- V4.1.2: MC, Oxford, 28 May 2012
Updated documentation.

- V4.1.1: MC, Oxford, 8 December 2011
Only calculates FLUX if required.

- V4.1.0: MC, Oxford 19 October 2010
Included TENSOR keyword to calculate any of the six components of the symmetric proper motion dispersion tensor (as in note 5 of the paper).

- V4.0.9: MC, Oxford, 15 September 2010
Plot and output with FLUX keyword the PSF-convolved MGE surface brightness.

- V4.0.8: MC, Oxford, 09 August 2010
Use linear instead of smooth interpolation. After feedback from Eric Emsellem.

- V4.0.7: MC, Oxford, 01 March 2010
Forces q_lum && q_pot < 1.

- V4.0.6: MC, Oxford, 08 February 2010
The routine TEST_JAM_AXISYMMETRIC_RMS with the usage example now adopts a more realistic input kinematics.

Updated documentation.

- V4.0.5: MC, Oxford, 6 July 2009
Skip unnecessary interpolation when computing a few points without PSF convolution. After feedback from Eric Emsellem.

- V4.0.4: MC, Oxford, 29 May 2009
Compute FLUX even when not plotting.

- V4.0.3: MC, Oxford 4 April 2009
Added keyword RBH.

- V4.0.2: MC, Oxford, 21 November 2008
Added keywords NRAD and NANG. Thanks to Michael Williams for reporting possible problems with too coarse interpolation.

- V4.0.1: MC, Windhoek, 29 September 2008
Bug fix: when ERMS was not given, the default was not properly set. Included keyword STEP. The keyword FLUX is now only used for output: the surface brightness for plotting is computed from the MGE model.

- V4.0.0: MC, Oxford, 11 September 2008
Implemented PSF convolution using interpolation on a polar grid. Dramatic speed-up of calculation. Further documentation.

- V3.2.0: MC, Oxford, 14 August 2008
Updated documentation.

- V3.1.3: MC, Oxford, 12 August 2008
First released version.

- V2.0.0: MC, Oxford, 20 September 2007
Introduced new solution of the MGE Jeans equations with constant anisotropy sig_R = b*sig_z.

- V1.0.0: Michele Cappellari, Vicenza, 19 November 2003
Written and tested

## Project details

## Release history Release notifications | RSS feed

## Download files

Download the file for your platform. If you're not sure which to choose, learn more about installing packages.