; CUBIT ;--------------------------------------------------------------- ;! Model a galaxy's density and velocity distribution from full cube ;# TASK MODELING ;----------------------------------------------------------------------- ;; Copyright (C) 2007 ;; Associated Universities, Inc. Washington DC, USA. ;; Developed by Judith Irwin, University of Toronto 1989-2000 ;; ;; This program is free software; you can redistribute it and/or ;; modify it under the terms of the GNU General Public License as ;; published by the Free Software Foundation; either version 2 of ;; the License, or (at your option) any later version. ;; ;; This program is distributed in the hope that it will be useful, ;; but WITHOUT ANY WARRANTY; without even the implied warranty of ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ;; GNU General Public License for more details. ;; ;; You should have received a copy of the GNU General Public ;; License along with this program; if not, write to the Free ;; Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, ;; MA 02139, USA. ;; ;; Correspondence concerning AIPS should be addressed as follows: ;; Internet email: aipsmail@nrao.edu. ;; Postal address: AIPS Project Office ;; National Radio Astronomy Observatory ;; 520 Edgemont Road ;; Charlottesville, VA 22903-2475 USA ;----------------------------------------------------------------------- CUBIT LLLLLLLLLLLLUUUUUUUUUUUU CCCCCCCCCCCCCCCCCCCCCCCCCCCCC CUBIT Task to model a galaxy's density & velocity dist'ns 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 ---------------------------------------------------------------- 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 R 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 ---------------------------------------------------------------- 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% 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%. 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**************