Thin disk

Back to the main page

Thindisk_fit. Version 2021

Thindisk_fit is a program to model the intensity of a geometrically thin disk and to fit the intensity of an spectral line tracing its kinematics. The data to be fitted are read from a FITS file with 3 axes: the position along the major axis of the projection of the disk onto the plane of the sky, the position along the minor axis, and the velocity. The inner and outer radius of the disk are r1 and r2, and the angle between the disk axis and the line of sight is the inclination i (i=0 for a face-on disk; i=90 deg for a edge-on disk). Two different motions are considered in the disk:

• rotation around its axis, with a rotation velocity given by a power law vr(r/r0)^qr

• infall toward the disk center, with an infall velocity also given by a power law  vi(r/r0)^qi

where r0 is an arbitrary reference radius, and vr and vi are the rotation and infall velocities at the reference radius.

The kinematics of the disk is computed as follows. We consider a cube of velocity channels and, for each channel, a grid of points onto the plane of the disk. For each point of the grid, the projection of its rotation and infall velocities onto the line of sight, vz, is calculated. A Gaussian line profile of width dv, centered on vz, and with an intensity proportional to (r/r0)^qd, is added to the channels associated with the grid point. After the cube onto the plane of the disk is computed, its projection onto a cube with a grid onto the plane of the sky is calculated. Finally, each channel map of the plane-of-the-sky cube is convolved with a Gaussian beam.

The final cube contains all the information on the kinematics of the model disk. However, the intensity scale is arbitrary. The scaling factor, the same for all channel maps, is obtained by minimizing the sum of the squared differences between the data cube and the model cube.

The fitting procedure consists in sampling the m-dimensional parameter space (m is the number of parameters fitted) using a pseudo-random sequence. For a description of the method, see Estalella (2017). The range explored for each parameter can be chosen. Any parameter can be kept fixed by setting its range to zero. The procedure also finds the uncertainties of the values found for the fitted parameters.

The input data file has to be a FITS data cube with three axes. The first axis has to be along the projected-disk major axis, the second axis along the minor axis, and the third axis has to be the radial velocity. The best fit parameters obtained from the fit are used to generate a model FITS file with the same geometry as the data file.

Model parameters 

The calculation uses 3 parameters defining the iterative fitting procedure, 4 fixed parameters and a maximum of 12 fit parameters.

• Iteration parameters:

  • Niter, number of iterations per loop;

  • Nloop, number of loops;

  • LoopFct, factor of decrease of the search ranges for each loop;

• Fixed parameters:

  • D, source distance (pc), only used in the central mass estimation;

  • bmaj, bmin, bpa, major and minor axes width of the beam (arcsec), and position angle of the beam major axis (degrees), from north to east;

• Fit parameters:

  • dv, linewidth (km/s);

  • x0, y0, position of the disk center, in offsets from the map center (arcsec);

  • v0, disk central velocity (km/s);

  • i, disk inclination angle (i=0 face-on, i=90 edge-on) (deg);

  • r1, r2, disk inner and outer radii (arcsec);

  • qd, power-law index of the disk intensity as a function of radius;

  • vr sin i, qr, disk rotation velocity at the reference radius r0, projected onto the line of sight (km/s), and power-law index of the dependence on radius, vr(r/r0)^qr;

  • vi sin i, qi, disk infall velocity at the reference radius r0, projected onto the line of sight (km/s), and power-law index of the dependence on radius, vi(r/r0)^qi.

Some of the parameters can be known beforehand, such as dv, x0, y0, v0, or i. Some other can be guessed from physical grounds (i.e., qd=-1; qr=qi=-0.5). The inspection of the data usually suggests whether the dominant motion is rotation or infall. These considerations leave only a few free parameters, which can be estimated through model fitting to the data.

Once the final values of the inclination i and the projected rotation and infall velocities vr sin(i) and vi sin(i) are obtained, the deprojected velocities vr and vi are calculated. Finally, the central mass necessary for producing the rotation (assumed to be Keplerian) or infall (assumed to be free-fall) are estimated as

M*(rot)= vr^2 r0 / G,

M*(inf)= vi^2 r0 / 2G.

In practical units,

[M*(rot) / Msol]= 1.127E-3 [vr / km/s]^2 [r0 / arcsec] [D / pc],

[M*(inf) / Msol]= 5.636E-4 [vi / km/s]^2 [r0 / arcsec] [D / pc].

The program generates plots of the zeroth order (velocity-integrated intensity), first order (intensity-weighted velocity) and second order (intensity-weighted velocity dispersion) moments, position-velocity cuts along the major and minor axes, and channel maps for the data, the model, and the residual data minus model (see the Figures below).

The program code (Fortran) can be downloaded here. The publication you can cite if you use this program is:

Zapata, L. A.; Garay, G.; Palau, A.; Rodríguez, L. F.; Fernández-López, M.; Estalella, R.; Guzmán, A. 2019, ApJ, 872, 176

If you want to use the program and need support, you can contact me at robert.estalella(at)ub.edu

Requisites

• Linux or Mac OS X system.

• Fortran compiler, e.g. gfortran.

• PGplot Graphics Subroutine Library (www.astro.caltech.edu/~tjp/pgplot/) compiled with gfortran. For the installation of PGplot with gfortran, follow the installation procedure described in pgplot_gfortran. Mac users can find a guide for installing PGplot on Mac OS X in mingus.as.arizona.edu/~bjw/software/pgplot_fix.html. You can find a local copy of the same page in pgplot_macosx.

Installation

Once downloaded, untar the tar file and compile the program with the PGplot and X11 libraries

gfortran -o thindisk_fit -l X11 -l pgplot thindisk_fit.f90

You may need to indicate the location of the PGplot or X11 libraries, adding the options, for instance,

-L /usr/local/pgplot -L /usr/X11/lib

Figures 1 - 2. Rotating ring

Fig. 1. Position-velocity cuts along the major (top) and minor (bottom) axes for different values of the inclination of a ring of radius r1=r2=30" with Keplerian rotation (qr=-0.5). The projected rotation velocity is vr sin(i)=10 km/s at r0=1", the beam is circular, bmaj=bmin=5", the linewidth is dv=3 km/s, and the intensity power-law index is qd=-1.

Fig. 2. Position-velocity cuts along the major (top) and minor (bottom) axes for different values of the radius of a ring with an inclination i=80 deg, with Keplerian rotation (qr=-0.5). The projected rotation velocity is vr sin i=10 km/s at r0=1", the beam is circular, bmaj=bmin=5", the linewidth is dv=3 km/s, and the intensity power-law index is qd=-1.

Figures 3 - 6. Rotating disk

Fig. 3 Position-velocity cuts along the major (top) and minor (bottom) axes for different values of the inclination of a disk of radii r1=0, r2=30" with Keplerian rotation (qr=-0.5). The projected rotation velocity is vr sin i=10 km/s at r0=1", the beam is circular, bmaj=bmin=5", the linewidth is dv=3 km/s, and the intensity power-law index is qd=-1.

Fig. 4. Position-velocity cuts along the major (top) and minor (bottom) axes for different values of the radius of a disk with inclination i=80 deg and Keplerian rotation (qr=-0.5). The projected rotation velocity is vr sin i=10 km/s at r0=1", the beam is circular, bmaj=bmin=5", the linewidth is dv=3 km/s, and the intensity power-law index is qd=-1.

Fig. 5. Position-velocity cuts along the major (top) and minor (bottom) axes for different values of the intensity power-law index of a disk with inclination i=80 deg and Keplerian rotation (qr=-0.5). The projected rotation velocity is vr sin i=10 km/s at r0=1", the beam is circular, bmaj=bmin=5", and the linewidth is dv=3 km/s.

Figures 7 - 8. Infalling ring

Fig. 7. Position-velocity cuts along the major (top) and minor (bottom) axes for differebt values of the inclination of a ring of radius r1=r2=30" with free-fall infall (qi=-0.5). The projected infall velocity is vi sin i=10 km/s at r0=1", the beam is circular, bmaj=bmin=5", the linewidth is dv=3 km/s, and the intensity power-law index is qd=-1.

Fig. 8. Position-velocity cuts along the major (top) and minor (bottom) axes for different values of the radius of a ring with an inclination i=80 deg, with free-fall infall (qi=-0.5). The projected infall velocity is vi sin i=10 km/s at r0=1", the beam is circular, bmaj=bmin=5", the linewidth is dv=3 km/s, and the intensity power-law index is qd=-1.

Figures 9 - 12. Infalling disk

Fig. 9. Position-velocity cuts along the major (top) and minor (bottom) axes for different values of the inclination of a disk of radii r1=0, r2=30" with free-fall infall (qi=-0.5). The projected infall velocity is vi sin i=10 km/s at r0=1", the beam is circular, bmaj=bmin=5", the linewidth is dv=3 km/s, and the intensity power-law index is qd=-1.

Fig. 10. Position-velocity cuts along the major (top) and minor (bottom) axes for different values of the radius of a disk with inclination i=80 deg and free-fall infall (qi=-0.5). The projected infall velocity is vi sin i=10 km/s at r0=1", the beam is circular, bmaj=bmin=5", the linewidth is dv=3 km/s, and the intensity power-law index is qd=-1.

Fig. 11. Position-velocity cuts along the major (top) and minor (bottom) axes for different values of the intensity power-law index of a disk with inclination i=80 deg and free-fall infall (qi=-0.5). The projected infall velocity is vi sin i=10 km/s at r0=1", the beam is circular, bmaj=bmin=5", the linewidth is dv=3 km/s.