AIPS HELP file for CUBIT in 31DEC24
As of Tue Apr 23 10:29:35 2024
CUBIT: Task to model a galaxy's density & velocity dist'ns
INPUTS
INNAME Input image name (name)
INCLASS Input image name (class)
INSEQ 0.0 9999.0 Input image name (seq. #)
INDISK 0.0 9.0 Input image disk unit #
OUTNAME Output image name (name)
OUTCLASS Output image name (class)
OUTSEQ 0.0 9999.0 Output image name (seq. #)
OUTDISK 0.0 9.0 Output image disk unit #.
BLC 0.0 4096.0 Bottom left corner of input
TRC 0.0 4096.0 Top right corner of input
FACTOR 0.0 1.0 Convergence criterion
(0 => .001)
VCODE Velocity code. Allowed forms:
'CV','SB','BR','RC'. Any-
thing else => 'CV'
DCODE Dens. code in plane.Allowed:
'CD','EX',GS','RG'. Anything
else => 'CD'
FPARM 0.0 Initial guesses of parameters
1,2. X,Y pixel position of
center in RA, DEC plane res-
pectively (specify)
3. position angle of receding
major axis (degrees)
4. inclination (1 - 89 deg)
5. systemic velocity (km/s)
(should be positive)
6. Vmax - max rotation vel
(km/s)(OR-scaling for RC
0=>1)
7. Rmax -radius at which Vmax
occurs (") 0 => DPARM(2)
8. m - Brandt index 0=>1
9. max in-plane density,
n0 (cm**-3)
10. radial density scale lngth
(outer) (")
11. perp EX dens scale length
(") 0 => DPARM(7) if
cparm(4)=2
12. radial inner scale length
(")RG only 0=>FPARM(10)
13. perp GS dens scale length
(") 0 => DPARM(7) if
cparm(4)=3
14. nf1 - fractional EX perp
dens (0=>1 if cparm(4)=1,2)
(0=>0.5 if cparm(4)=4)
15. nf2 - fractional GS perp
dens (0=>1 if cparm(4)=3)
(0=>0.5 if cparm(4)=4)
CPARM 0.0 65534.0 1. 0 => write residuals
1 => write model
2. Distance to galaxy (Mpc)
(specify)
3. Parameters to be held
fixed (see HELP CUBIT)
0 => all parms free
4. perp density dist'n #
(1,2,3,4) see explain
anything else =>1
5. R0 -- center of ring dens
dist'n (RG only)
6. FWHM (km/s) of gaussian
for velocity smoothing
(see EXPLAIN CUBIT)
7. For RG only: 2=exp
3=gaus, 0==>3
DPARM -1.0 8000.0 1. min lim to radius (")
2. max limit to radius
in plane of galaxy
3,4: min/max limits to
cos(psi)
5,6: min/max limits to
sin(psi)
psi measured from receding
major axis in RA,DEC plane
7. upper limit to 1/2 disk
thickness. (specify)
VPARM 0.0 1000.0 User specified rotation curve
(km/s) for VCODE='RC'
RPARM 0.0 5000.0 User specified radii (")
corresponding to VPARMS
ICUT Use only data > in abs value
PIXSTD 0.0 100.0 Estimated rms uncertainty in
input map pixels (Jy/ba)
0 => .001
HELP SECTION
CUBIT
Type: Task
Use: CUBIT models the HI density and velocity distributions of a
galaxy simultaneously. The model assumes that the galaxy is a
disk of finite thickness in circular rotation. Several density
distributions and rotation curves are allowed (see EXPLAIN
CUBIT) and spatial smoothing is performed. A non-linear least
squares fit solves for up to 15 parameters and either the final
model or residuals can be written to the output cube. CUBIT
expects as input, a cube of data with VELOCITY in M/S as the
FIRST AXIS, RA and DEC as the 2nd and 3rd axes respectively,
and pixel values in JY/BEAM. Blanking is supported. The
output cube size and resolution are specified by BLC, TRC and
the resolution of the input cube, respectively.
If you use this task to analyze your data cube and make
reference to the results in a paper, you should reference the
original paper describing this task:
Irwin, Judith A. 1994, "Arcs and bridges in the interacting
Galaxies NGC 5775/NGC 5774", ApJ, 429, 618-633.
Adverbs:
INNAME.....Input image name (name). Standard defaults.
INCLASS....Input image name (class). Standard defaults.
INSEQ......Input image name (seq. #). 0 => highest.
INDISK.....Disk drive # of input image. 0 => any.
OUTNAME....Output image name (name). Standard defaults.
OUTCLASS...Output image name (class). Standard defaults.
OUTSEQ.....Output image name (seq. #). 0 => highest unique.
OUTDISK....Disk drive # of output image. 0 => highest disk
number with sufficient space.
BLC........Bottom left corner in input image of desired
subimage. 0,0,0 => 1,1,1
TRC........Top right corner in input image of desired
subimage. Default is entire image.
VCODE......Code specifying form of rotation curve. 'CV' (constant
velocity), 'SB' (solid body), 'BR' (Brandt), and 'RC'
(user-specified curve) are recognized. Anything else =>
'CV'.
DCODE......Code specifying form of density distribution in plane.
'CD' (constant), 'EX' (exponential), 'GS' (Gaussian), and
'RG' (offset asymmetric ring) distributions are
recognized. Anything else => 'CD'
FPARM......First guesses to parameters to be fitted.
1. X pixel position of center (RA). (input map)
2. Y pixel position of center (DEC).
3. Position angle of receding major axis (deg)
4. Inclination of galaxy plane (1 - 89 degrees)
5. Systemic velocity (km/s)
6. Vmax - Maximum rotation velocity (km/s).
If VCODE=RC, FPARM(6) = scaling factor of VPARMS
where 0=>1.
7. Rmax - Radius at which Vmax occurs (")
0 => DPARM(2) (ignored for CV,RC)
8. m - Brandt rotation curve index - 0=>1 (BR only)
9. n0 - density from in-plane distn(cm**-3)
10. Density scale length in plane of galaxy (").
Ignored for CD.
11. Density scale length perpendicular to plane for
exponential ("). (ignored if CPARM(4)=1)
12. Radial density scale length for RCPARM(4)=1)
14. nf1 - fraction of n0 belonging to h distn 1
0 => 1
15. nf2 - fraction of n0 belonging to h distn 2
CPARM......1. 0 => write residuals. ie. model - data
1 => write model
2. Distance to galaxy (Mpc) - (Specify)
3. Parameters to be held fixed. Sum of 2**i, where
i is the parameter number. eg. parameters 1,2
and 5 fixed => CPARM(3)=2**1 + 2**2 + 2**5=38.
N.B.***Different from GAL***
4. Vertical density distribution code - 1,2,3,4 for
'CD','EX','GS','EG'; anything else => 1
5. R0 - radial location of gaussian center (RG only)
arcsec
6. FWHM of gaussian for velocity smoothing km/s (see
explain CUBIT)
7. For RG only: Specifies form of density distribution
in disk, 2=Exp, 3=Gaus, 0==> 3
DPARM......1. Lower limit to radius in galaxy plane (").
2. Upper limit to radius ...
3. Min limit to cos(psi). Angle in plane of sky
measured from receding major axis.
4. Max limit to cos(psi)
5. Min limit to sin(psi)
6. Max limit to sin(psi)
DPARM(3 to 6)=0 => all quadrants modelled.
7. Max limit to 1/2 disk thickness measured from
center of plane ("). Specify.
All points outside of limits DPARM(1->7) are zeroed
VPARM......Array of up to 30 velocities (km/s) specifying a rotation
curve (RC only)
RPARM......Radii in plane of galaxy from center (") corresponding to
each VPARM (need not be uniformly spaced).
ICUT.......A flux cuttof (Jy/beam). If abs(data) < abs(icut) then
the residuals are set to zero in the least-squares fit
PIXSTD.....Estimated rms uncertainty associated with input image
pixels (Jy/beam) 0=>.001
EXPLAIN SECTION
CUBIT: Task which models the HI density and velocity distributionns
of a galaxy disk.
DOCUMENTOR: J. A. Irwin, U. of Toronto for CUBIT
RELATED PROGRAMS: GAL, AIPS tasks which apply to cubes of data.
DESCRIPTION
CUBIT uses input parameters to calculate intensities (Jy/beam)
for each non-blanked pixel position (V, RA, DEC) from a cube of data.
The model assumes a smooth density distribution in a circular disk of
finite thickness, and circular velocities. Presently, the allowed
density distributions are constant, exponential, gaussian, and offset
asymmetric ring (exp or gauss) distributions. The allowed velocity
distributions are constant, solid body and Brandt curves. The user
also has the option of specifying a rotation curve numerically, which
is allowed to vary in magnitude (but not shape). Segments of the sky
plane or rings in the galaxy plane can also be isolated for separate
modelling (through the DPARMS). The model intensities are calculated
for each pixel and then smoothed spatially to the resolution of the
input map beam. To solve for 'best fit' parameters, the input data
are subtracted from the model and, given the partial derivatives of
the intensity with respect to each free parameter, a non-linear least
squares fit is performed using the Levenberg-Marquardt algorithm
(modeled on subroutines LMSTR, LMSTR1). Any number of parameters may
be held fixed. The solution (with internal errors) is then printed
and an output cube containing either the final model or residuals is
written. Velocity smoothing is supported; see CPARM(6) for details.
The vertical column density (integrated over the distribution in z
specified by the model and fit parameters) is displayed at the end.
(Skip to the bottom for Quick Start Examples)
CALCULATION OF THE MODEL
For each input pixel, its position in the sky plane (x,y) is
calculated with respect to the current x(0), y(0) position of the
galaxy centre, using the receding major axis as the positive X axis.
The radial velocity at the pixel relative to the current systemic
velocity (Vr) is also determined. The fundamental equation:
I(x,y,Vr) = K n(R,h) dl/dVr
is then used to calculate the intensity, I (Jy/ba) where:
dl/dVr (arcsec/(km/s))
is the change in LINE OF SIGHT distance, l, (measured +ve or
-ve from the center of the galaxy's disk) for a given channel
separation, dVr (from header). dl is calculated from l3 - l1
where li is the line of sight distance corresponding to Vi, Vi
being the velocity at the channel center (i=2) or at either
'edge' (1,3). Given Vi, the position, x,y, and the assumed
rotation curve, li can be calculated. The intensity is zeroed
if all li (i = 1 => 3) lie outside the galaxy disk as specified
by Dparm 1, 2 and 7. The intensity is calculated according to
the above equation if any or all li occur within the disk, with
li adjusted to the boundary value if it exceeds the boundary
value. Note that two values each of li (i=1,2,3) may result
for a given x,y,Vr if the disk is sufficiently edge-on and/or
thick. In such cases the computed intensity is the sum of the
intensities contributed from each l position (optically thin
conditions are assumed).
For the 'solid body' rotation curve, l cannot be determined
from the given Vr and rotation curve since l is independent of
V. Instead, all gas along the line of sight occurs at a single
velocity (ie. within a single channel). In this case, the
model velocity is calculated from the pixel position of
interest, Vmax and Rmax. Then, if the model velocity falls
within the velocity channel of the pixel, dl is set to the the
entire line of sight disk width, with the intensity calculated
as before. If the model velocity falls outside of the channel,
the intensity is set to zero.
Along the minor axis, l cannot be determined from Vr, either.
Here all the gas should be within one channel at Vsys. If the
y pixel position lies within 1/20 of a HPBW from the current
minor axis, then the calculated intensity is non-zero only if
the pixel velocity is within a half channel width from Vsys.
The intensity is calculated as before with dl set to the total
line of sight disk width.
If VCODE=RC, l cannot be analytically determined from x, y, and
Vr, so the radius of the point in the galaxy plane, R, is not
known. Therefore R must be determined iteratively with l
initially assumed to be zero until a value of R and l are found
which agrees with the specified rotation curve, Vr and the
position x,y. The iterations stop when (the change in V
(circular) is less than half the channel separation AND the
change in R is less than half the average beam size (for less
than 50 iterations)) OR (the change in V (circular) is less
than the channel separation AND the change in R is less than
the average beam size (for more than 50 iterations)). If
greater than 200 iterations are required to find R (and
therefore l) the last two values of R are simply averaged and
that value is used. The routine prints a message indicating
how many points converged and how many were approximated (some
of these points are later zeroed for other reasons -- e.g.
because they physically fall outside of the galaxy disk -- so
the total number of points listed will be more than the number
of non-zero points printed in a later message). Once R has
been found, l is calculated, then dl, and the intensity is
computed as before.
n(R,h) (cm**-3)
is the density at the position R (in arcsec, the distance from
the center of the galaxy along the plane) and h (the distance
perpendicular to the center of the plane) with R and h found
from x, y and l. Within a single velocity channel, however,
there may be a range of densities, since the line of sight
distance, l, may vary considerably from one edge of the channel
to the other. Consequently, the calculated model density for a
given channel is the (integral of the densities from l1 to l3)
divided by (l1 - l3). The integration is computed numerically
in adjustable steps (in l) with a calculated value being
accepted if it meets one of two criteria: either a) the change
in the function is less than or equal to 5 percent or b) the change in
the function is less than 0.01 X PIXSTD (see below). In cases
where the numerical computation could be compared with an exact
analytical solution, the rms difference between results was
better than 1 percent. Note that by including PIXSTD as one
acceptance criterion, the user can speed up the algorithm by
increasing PIXSTD and accepting cruder density calculations; of
course the final listed errors will then be incorrect.
K
is a constant which converts from cm**-3 arcsec/(km/s) to
Jy/beam. Folded into 'K' is a conversion from arcsec to cm
(through CPARM(2)), a conversion from cm**-2 to degrees K
(through the well-known equation relating the integral of
brightness temperature over velocity to column density), and a
conversion from degrees K to Jy/beam (through the beam size and
velocity (=>frequency) in the header). It is assumed that the
beam size specified in the header corresponds to the reference
velocity in the header. Calculations are done in km/s for
header velocities in m/s.
DESCRIPTION OF ADVERBS
BLC,TRC
These are the usual bottom left and top right corner pixel
positions. The 31DEC07 version of CUBIT can handle any size
of image if the computer has enough memory. Of course,
larger images are slower.
FACTOR
Determines when the iterations should stop. A smaller value
of FACTOR results in increasingly finer adjustments to the
parameters but more computing time. The quality of the fit
is indicated by the 'INFORMATION PARAMETER' which is printed
in the message file. The codes are:
INFORMATION PARAMETER =
0 Improper input parameters
1 Algorithm estimates that the relative error in
the sum of the squares is at most, FACTOR (ie.
change in sum of squares of residuals from the
last iteration to the current one, divided by
sum of squares of residuals from last iteration
is less than or equal to FACTOR)
2 Algorithm estimates that the relative error in
the approximate solution is at most, FACTOR
3 Conditions for INFO=1 and INFO=2 both hold
4 The vector containing the residuals is orthogonal
to the columns of the Jacobian, to machine precision.
5 Number of iterations has reached 100*(N+1), where N
is the number of free parameters.
6 FACTOR is too small. No further reduction in
the sum of squares is possible.
7 FACTOR is too small. No further improvement in
the approximate solution is possible.
VCODE
Form of the rotation curve. It is assumed that V(R,h) = V(R,0)
where R is the radial coordinate measured from the center of
the galaxy outwards along the galaxy's plane
h is the coordinate measured from the center of the
plane in a direction perpendicular to the plane.
'CV' - constant velocity
V(R) = V(max) where V(max) = FPARM(6)
'SB' - solid body
V(R) = (V(max)/R(max)) R
where R(max) = FPARM(7)
'BR' - Brandt curve
V(max) R / R(max)
V(R) = -----------------------
m 3/(2m)
( 1/3 +2/3 (R/(R(max)) )
where m = FPARM(8)
'RC' - User specified rotation curve
V(Ri) = f (Vi)
where Vi are specified through
VPARM's
Ri are specified through
RPARM's
f = FPARM(6)
DCODE
Form of the density distribution in the plane of the galaxy
n(R,h) where h=0.
'CD' - constant density
n(R,0) = n0
where n0 = FPARM(9)
'EX' - exponential density
n(R,0) = n0 exp(-R/a)
n0 = FPARM(9)
a = FPARM(10)
'GS' - gaussian density
2 2
n(R,0) = n0 exp(-R /a )
n0 = FPARM(9)
a = FPARM(10)
'RG' - offset asymmetric ring
If cparm(7)=2, then
n(R,0) = n0 exp(-(R-R0)/a)
If cparm(7)=3, then
2 2
n(R,0) = n0 exp(-(R-R0) /a )
n0 = FPARM(9)
a = FPARM(10) (R >= R0)
a = FPARM(12) (R < R0)
R0 = Cparm(5)
FPARM
1,2 x,y pixel positions of galaxy's center in RA, DEC plane
respectively of input map. If outside of BLC, TRC window,
an error message is printed and the program terminates.
3 Position angle of receding major axis measured from N
through E from 0 to 360 degrees.
4 Inclination of galaxy to line of sight from 1 to 89 degrees.
If outside of this range, an error message is printed and
the program terminates. The optically thin model is
independent of which 'edge' of the galaxy is the near
side.
5 Systemic velocity - km/s.
6 Rmax (as explained under VCODE) arc seconds
7 Vmax (as explained under VCODE) km/s
8 m (as explained under VCODE) Brandt exponent
9 n0 (as explained under DCODE) 1/cm^3
10 a (as explained under DCODE) arc seconds
11 Perpendicular EX density distribution scale height, b1 (see
CPARM(4) below) arc seconds
12 a for R < R0 (RG density code only) arc seconds
13 Perpendicular GS density distribution scale height, b2 (see
CPARM(4) below) arc seconds
14 nf1 -- fraction of n0 needed for perpendicular EX distn
(see CPARM(4) below)
15 nf2 -- fraction of n0 needed for perpendicular GS distn
(see CPARM(4) below)
CPARM
1 1 => write the final model in the output cube
anything else => write the residuals in the output cube
2 Distance to the galaxy (Mpc) (needed to calculate the con-
version constant, K). If not specified, an error message
is printed and the program terminates. Note that an increase
in D gives *higher* values of model intensity, I. This is
because, for a specified density distribution, n(R,h), at a
position, R, h, (in arcsec), the line of sight coordinate
dl contains more cm/arcsec.
3 Code for fixing parameters.
Sum of 2**i, where i is the parameter number. eg.
if the 1st, 6th, and 9th parameters are to be fixed, then
CPARM(3) = 2**1 + 2**6 + 2**9 = 2 + 64 + 512 = 578. If
CPARM(3)=65534, then all parameters are fixed and no
least squares solution is found - either the residuals
or model specified by the input parameters are written.
4 Density distribution number to specify perpendicular density
distribution, where 1 => 'CD' (constant density), 2 => 'EX'
(exponential distribution), 3 => 'GS' (Gaussian
distribution), and 4 => 'EG' (sum of an exponential and
Gaussian). 0 => 1
The density distribution at galactocentric radius, R, and
perpendicular height, h, is specified by:
n(R,h) = n(R,0) x n(h)
where n(h) = 1 if CPARM(4)=1 and if CPARM(4) > 1
n(h) = fn1 exp(-h/b1)+ fn2 exp(-h^2/b2^2)
where b1 = FPARM(11)
b2 = FPARM(13)
fn1 = FPARM(14)
fn2 = FPARM(15)
If CPARM(4)=2, then fn1, b1 are required and 2nd part of
equation is ignored (not added)
if fn1, b1 are zero, they are set to 1 and
DPARM(7), respectively
If fn1 is greater than 1.0, it is set to 1.0
If CPARM(4)=3, then fn2, b2 are required and 1st part of
equation is ignored (not added)
if fn2, b2 are zero, they are set to 1 and
DPARM(7), respectively
If fn1 is greater than 1.0, it is set to 1.0
If CPARM(4)=4, then both fn1 and fn2 are required, fn1+fn2 must
equal 1. If these 2 values are zero, then
fn1 is set to 0.5 and fn2 is set to 0.5
also b1 and b2 are required. If b1,b2 are
zero, then they are set to DPARM(7)
If fn1 or fn2 > 1.0, they are set to 1.0
A warning is issued if fn1+fn2 .ne. 1.0
If only one of these parameters is fixed,
the routine will automatically fix the other
5 Center of offset gaussian when Dcode = RG, arcsec
6 FWHM of gaussian smoothing fcn in velocity space (km/s).
This applies velocity smoothing to the model in the
line-of-sight only. It is left to the observer to translate
this into a dispersion in or perpendicular to the plane of
the galaxy. If CPARM(6)=0, no velocity smoothing is done.
Note that velocity smoothing should, in principle, not
affect the determination of the density, n. Since Integral
of I dVr = Integral of n dl, when dVr increases, I decreases,
and n dl is not affected, i.e. the disk is "heating up" but
not "expanding". If you try to fit a smoothed distribution
with an unsmoothed model, however, n will be in error by
an amount which is not a simple scale of velocity widths
because CUBIT is attempting to minimize residuals.
**NB> The current version of CUBIT does not apply velocity
smoothing in the least squares subroutine (whereas spatial
smoothing *is* included there). This means that a non-zero
CPARM(6) will velocity smooth the model to give a correct
output model or residual cube, *but* will not search
correctly through parameter space to find the best-fit
solution.** If velocity smoothing is required, it is
necessary to run several trials of CUBIT with no velocity
smoothing first to get a general idea what the parameters
are. After that, turn CPARM(6) on, fixing all parameters,
and start systematically working through parameter space
to see if the addition of velocity smoothing can improve
the model (i.e. check the residuals for the smallest rms).
The velocity smoothing subroutine is taken from the aips
task XSMTH, and has been checked against XSMTH for correct
output. The smoothing here corresponds to DPARM(1)=1,
DPARM(3)=6, and DPARM(4)=0.01 in XSMTH, i.e. a gaussian
smoothing function, with a support radius of 3 times the
HWHM, with smoothing over blanks allowed.
7 Switch to determine whether the ring density distribution
will fall off as an exponential or as a gaussian from its
center at R0=cparm(5). Must be exponential both inside
and outside of R0 or gaussian both inside and outside of
R0 (i.e. can't mix a gaussian with an exponential).
DPARM Sets limits to the size of the disk or sections of it for
the least squares fit. All unblanked points outside of
these limits are zeroed.
1,2 Angular extent of the inner and outer disk radius in the plane
of the galaxy.
3,4,5,6 Permits isolation of sections of the galaxy in the plane of
the sky.
Eg. DPARM(3) ~ -1, 0, -1, 0 implies that fitting is done for
the third quadrant only (receding major axis is the positive
X axis). All zero => keep all quadrants
7 1/2 thickness of the disk as measured from the center of the
galaxy plane 'above' or 'below' it.
VPARM Array of up to 30 user specified velocities defining a rotation
curve. All values should be positive, and a Vparm should be
specified for each non-zero rparm. CUBIT does a LINEAR inter-
polation between specified VPARMS to calculate intermediate
velocities. If the radius of the pixel considered is less than
the first Rparm specified,then the the rotation curve velocity
is set to the interpolated value between 0 and Vparm(1). If
the radius is greater than the last Rparm specified, then the
velocity is set to Vparm(last non-zero value).
RPARM Array of up to 30 user specified radii corresponding to the
velocities specified by VPARM. All values should be positive
and increasing in magnitude with array subscript. The inter-
vals need not be evenly spaced.
ICUT If abs(data) < abs(icut) (Jy/beam), then the residuals are
set to zero during the least squares search through parameter
space, i.e. these points are not included in the search for
the best fit. However, during the final pass through,
(final computation of residual or model cube), all these data are
retained (provided they are in bounds as specified by the dparms)
so that the final model appears continuous or so that the user
can retain all low level emission in the residual cube for study.
Use of ICUT does not speed up the algorithm.
PIXSTD The estimated rms error in the input map pixels (Jy/ba). The
errors listed for the solution depend upon the error associated
with each pixel. Note that the error calculation assumes that
every pixel constitutes an independent data point. Usually,
this will not be the case, so the real errors will be somewhat
higher than those calculated. CUBIT should be run several
times with initial input parameters both higher and lower than
the 'true' values; variation in the final iterated value may
be a better representation of the real errors than the calculated
values. Otherwise, it is left to the user to determine the
relationship between the error per pixel (computed by CUBIT)
and the error per independent data point.
GETTING STARTED
The input cube requires beam dimensions and position angle in the
header. If the cube has not been cleaned, these quantities must be
inserted with PUTHEAD.
The 31DEC07 version of CUBIT has no limits on the V, X, or Y axes'
number of pixels. The speed of the program will be more or less
proportional to the total number of non-blanked pixels. You could
consider using SUBIM or LGEOM to reduce the number of pixels during
the early exploratory executions of CUBIT. When good estimates of the
parameters are found, then apply them to the full data cube.
**CUBIT does not recognize a "rotation" in the header** If you line
up your galaxy with the major axis along the 'x' axis to reduce the
total number of pixels, for example, the position angle of the
receding major axis (FPARM3) will be either 90 or 180 degrees. For
faster execution, reduce the cube size and BLANK all regions with no
emission (STRONGLY RECOMMENDED but leave a few pixels around real
emission, both in RA, DEC, and VEL). Also, several (ie. 3-4) pixels
per beam are required for the smoothing algorithm to function
accurately.
HINTS
Try blanking regions of no emission to make the routine run faster.
Do a first trial, putting in realistic values for the input
parameters from other knowledge, e.g. center and vsys from the
1st moment, etc. Once you are satisfied with a model for the
whole galaxy, try each side (advancing and receding) separately,
but keep X_0, Y_0 and V_sys fixed when operating on each side
separately. The best 1st rotation curve to use is the Brandt curve,
which can be made to look similar to a solid body curve by specifying
a large enough R_max, V_max and truncating the distribution if
necessary, before that maximum is reached. If it looks like you
need a ring/torus, rather than a disk (i.e. RG), you will have to
systematically set the ring center in trial after trial until
you're satisfied that the residuals are minimized. The best thing
to do is set up a batch job.
If a user specified rotation curve is desired, this can only be
run on each side separately, so it's best to start with a Brandt
curve anyway so that at least you can find X_0, Y_0, and V_sys
in order to keep them fixed when this is attempted.
Dparm(7) (the thickness of the disk) and Dparm(2) (the maximum
radius) should be large enough that all emission will be contained
within these boundaries. Specifying them too large is still okay
when using the exponential or gaussian distributions because the
scale lengths of these distributions will automatically reduce the
emission to zero at the appropriate places within these limits.
The last thing you should do when satisfied with all parameters is
to include a velocity dispersion (if the line widths are still too
narrow). A velocity dispersion is again set by hand and must be
varied systematically through different trials until the residuals
are minimized. Ensure that the other parameters are allowed to
vary freely when doing this and don't just fix the parameters to
values modelled previously when a velocity dispersion wasn't used
because there will be new best fit parameters.
APPLICABILITY TO CO DATA
There is no reason why this routine cannot be used for modelling
CO data, provided you believe the N_H_2 to CO conversion. The
input cube should be in the same format, i.e. Vel, RA, DEC, with
pixels in Jy/beam. Then the only output which will be in error
will be the calculated value of n_max. To correct this value, do
the following:
n_H_2 = X nu(HI)^2
------------- -------- n_HI
1.823 x 10^18 nu(CO)^2
where n_H_2 is the molecular hydrogen density desired from the CO
intensity, X is the CO to H_2 conversion factor (e.g. 3 x 10^20),
nu(HI) is the frequency of emission at the band
center if it were emitting in HI (e.g. 1416.5 x 10^6 Hz for a
recessional velocity of 820 km/s), nu(CO)^2 is the frequency that
the emission of emission at the band center for CO (e.g. 114.96
x 10^9 Hz for V = 820 km/s) and n_HI is the value of n returned
by the model (i.e. n_max). Using the rest frequency for HI and
CO in the above is probably okay for most circumstances.
The other requirement for the code to work for CO data is that the
molecular clouds do not shadow each other in position velocity space.
RUN TIME AND SPACE REQUIREMENTS
The run time depends upon the size of the input cube and the
number of unblanked pixels it contains, how close the input trial
parameters are to the solution, and the number of free parameters.
A model (all parameters fixed) can be generated in a few minutes for
a 60,000 data point cube on a SUN 3/180. The same conditions,
but with all parameters free, may require 5 or 6 hours of cpu time.
The routine is both "cpu hungry" and "memory hungry" -- not ideal,
since data cubes can be very large. However, modern computers are
very much larger and faster than the aforementioned SUN 3/180.
QUICK START EXAMPLES
1) I want to create a model galaxy that is inclined to the line of
sight and then determine its 0th moment. The galaxy should
have a density distribution that declines as a gaussian in the plane
of the galaxy and declines as an exponential perpendicular to the
plane. The rotation curve should be smooth -- a Brandt curve
is sufficient. Let the galaxy be centered in the cube, and have
parameters as indicated by the sample parameters, below.
DO THE FOLLOWING:
a) TASK 'CUBIT'; GETN n (Specify the task and the cube that you
already have as the input so that the model cube can get and
set the header info. *NB* input cube must be in VRD format --
if not, use TRANS and, if needed, ALTSW to switch from frequency
to velocity)
b) OUTNA 'outputcub'; OUTCL 'VRD' (give a name for the output cube)
c) CPARM(3)=65534 (fix all parameters -- no iterating)
d) CPARM(1)=1 (model output, rather than residuals)
e) CPARM(2)=12 (specify the distance of the galaxy in Mpc)
f) CPARM(4)=2 (exponential vertical z distribution)
g) VPARM 'BR' (use the Brandt rotation curve)
h) DCODE 'GS' (Gaussian density distribution in plane of galaxy)
i) FPARM(1)=57; FPARM(2)=39 (set the RA, DEC, pixel positions of
center of galaxy in the cube)
j) FPARM(3)=45; FPARM(4)=30 (set the position angle and inclination
-- note, need to avoid exactly edge-on or face-on orientations but
can get within a few degrees)
k) FPARM(5)=695; FPARM(6)=220 (set systemic velocity -- see your
input cube -- want to ensure that the systemic velocity is actually
in the cube, and the maximum of the rotation curve, in km/s)
l) FPARM(7)=200; FPARM(8)=0.1 (set the radius (arcsec) at which the
rotation curve will peak and how 'curvy' the rotation
curve is, e.g. m=0 for a flattish rotation curve, larger values
for one that tend to turn over)
m) FPARM(9)=0.1; FPARM(10)=150 (set the maximum density (cm^-3)
corresponding, in this case, to the peak of the gaussian at the
galaxy center as well as the corresponding radial scale length
(arcsec))
n) FPARM(11)=15; DPARM(7)=1000 (set the vertical (z) scale height
and set an upper limit to the region over which the vertical height
values are calculated -- dparm(7) should be any value much greater
than fparm(11), except for the case of a constant vertical density
distribution which isn't very physical.)
o) GO
A model cube is then created. You can then treat this cube like any
other cube of data in aips, finding the 0th or 1st moments (MOMNT),
etc.
2) I want to see how well the model galaxy that I just created
compares to the galaxy in the input cube.
DO THE FOLLOWING:
Repeat #1 but now set the switch, CPARM(1)=0 which will compute
the residuals, rather than the model.
3) I want to find the best fit parameters that suit the galaxy I
have in the input cube. What's the best way to proceed?
DO THE FOLLOWING:
a) (Optional) Start by rotating the galaxy to align it with the major
axis (LGEOM) and trimming regions (BLANK) that clearly have no
emission. The program will run faster. Leave several cells around
the emission (in v-r-d space) so that smoothing algorithms don't
have trouble.
b) Assign trial values of FPARMs as described in #1. We recommend
fixing FPARM(1) and FPARM(2) (the x, y position of the center), as
well as FPARM(5) (Vsys) at first since these are fairly easy to
determine without running CUBIT and help the algorithm along
i.e. set CPARM(3) = 2**1+2**2+2**5 = 38.
Don't forget that the FPARM(3) refers to the position angle of
the *receding* major axis. If this is wrong, the program won't
iterate. If IFLAG = 1 warnings are generated, it means that
the search through parameter space is approaching some limits, such
as the center of the galaxy, and the program is adjusting those
limits to keep parameters in the right range. This is not fatal
but the result will be best if the *last* iteration does not have
any warnings.
c) Once a good solution is found, repeat using CPARM(3)=0 so that
fine adjustments in central position and Vsys can be found.
Set CPARM(1)=0 to create a residual cube.
d) Measure the rms of the residual cube (IMEAN).
e) Repeat steps b through d for different density distributions
and/or including a velocity smoothing. The Brandt curve is the
most well-tested and we recommend it over the other velocity
distributions. We have also found it useful to try a variety
of initial inputs just in case the program is iterating to a
local minimum in parameter space, rather than the best
global minimum. The best fit result is considered to be
the one that gives the lowest rms in the residual cube.
f) Once the best parameters and distribution become more certain,
do a run of CUBIT for which PIXSTD is specified, so that the formal
error bars on the final parameters have their proper meanings.
g) (Optional) We have found that, after a few trials with CUBIT,
the above process can be batched. Overnight, a number of
residual cubes are generated and the rms of each can be measured in
the morning to see what the best solution is.
h) (Optional) Adjusting FACTOR to a smaller value will fine-tune
the results but take longer to iterate to a solution. In
practice, we have found that the result is not significantly
better than leaving FACTOR at the default, but there could be
some improvement depending on the galaxy being modeled.
i) (Optional) You might want to experiment with each side of
the galaxy, separately or with rings that isolate certain parts
of the galaxy, if this is scientifically desirable.
***********THE 31DEC07 VERSION OF CUBIT IS COMPATIBLE WITH***********
****************AIPS VERSIONS 31DEC04 AND LATER**************