ROTCUR(1NEMO) manual page

This HTML automatically generated with rman for NEMO
Table of Contents

Name

rotcur - fit parameters to velocity field using tilted rings

Synopsis

rotcur in=velfie [parameters=values ...]

Description

rotcur derives the kinematical parameters from the observed velocity field by fitting a tilted-rings model to the velocity field. It does a non-linear least-squares-fitting to the function:

         v(x,y) = VSYS + VROT * cos(theta) * sin(INC)
                      - (x-XPOS) * sin(PA) -/+ (y-YPOS) * cos(PA) ?
where:   cos(theta) = -----------------------------------------
                                       r
                      - (x-XPOS) * cos(PA) - (y-YPOS) * sin(PA)
and:     sin(theta) = -----------------------------------------
                                   r * cos(INC)

where the radius r is measured in the plane of the galaxy:

         r = sqrt( (x-XPOS)**2 + (y-YPOS)**2/cos(INC)**2 )

In the above formula v(x,y) denotes the radial velocity at rectangular sky coordinates x and y, VSYS the systemic velocity, VROT the rotational velocity, INC the inclination angle and theta the azimuthal distance from the major axis in the plane of the galaxy. In itself, theta is a function of the inclination (INC) and the position angle (PA) of the major axis. XPOS and YPOS denote the position of the rotation center in pixels w.r.t. 0,0 being the lower left corner of the map. rotcur can fit for each ring the 6 parameters VSYS, VROT, INC, PA, XPOS and YPOS, though any combination of them can be fixed (see fixed=). The position angle PA of the major axis is defined as the angle, taken in anti-clockwise direction between the north direction on the ‘‘sky’’ and the major axis of the receding half (positive radial velocity) of the galaxy.

Values equal the undefined value (currently 0.0) are ignored in the fit. See ccdmath(1NEMO) or fitsccd(1NEMO) on how to create a velocity field with such undefined values.

By setting the keyword fitmode= appropriately, rotcur can also be used to fit a model of pure expansional velocities:

         v(x,y) = VSYS + VROT * sin(theta) * sin(INC)

use fitmode=sin,1 in this case. currently XPOS and YPOS must be held fixed in this case.

For a detailed discussion of this method see K. Begeman, Astr. & Astrophys. 223, 47. (1989) and references therein.

See also velfit(1NEMO) on the classic WWB73 (Warner, Wright and Baldwin, 1973, MNRAS, 163,163) method.

Parameters

The following parameters are recognized:

in=image_vel: Input velocity field map, either in image(5NEMO) format, or table(5NEMO) format. See also imagemode= below. No default.
radii=r0,r1,r2,...rN: Inner and outer radii for N touching rings (in arcsec). Remember that the number of radii is hence one more than the number of rings. The units= keyword can be used to scale your physical units in the header to ‘‘human readable’’ units (arcsec). Radii must be either all increasing or decreasing. Note that the corresponding keyword in GIPSY has a different meaning, where they denote the ring centers, and a different keyword WIDTHS= controls the width of each ring. No default.
vrot=v1,v2,...vN: Array of rotation (or expansion) velocities, one for each ring. If less than the number of rings is given, the trailing ones are filled with the last supplied value.
pa=p1,p2,...pN: Array of position angles (degrees), one for each ring. As is the convention, this is the PA of the receding side of the galaxy. See also vrot=.
inc=i1,i2,...iN: Array of inclinations (degrees), one for each ring. See also vrot=.
vsys=v0: Systemic velocity. Only one number applies, though VSYS can vary accross rings.
center=x0,y0: Rotation center (grid units w.r.t. lower left corner being 0,0). Two numbers are required. Default: center of map.
frang=: free angle around minor axis (degrees), in the plane of the galaxy, from which data is excluded from the fit (thus the total cone size of ignored data around the minor axis is 2*frang) [Default: 20.0].
side=: Choose the side of the galaxy to fit the velocity field to. Valid options are receding, approaching or both side(s). [Default: both].
weight=: Choice of geometric weighting function with which points are weighed into the least squares solution as a function of galactic angle away from the major axis. Valid options are: uniform, cosine, and cos-squared. [Default: cosine].
fixed=: List of parameters, separated by commas, to be kept fixed during the fit. Choose any of the following: vsys, vrot, pa, inc, xpos, ypos, although at least one parameters should be kept free. [Default: none, i.e. all parameters free.
ellips=: The names of two parameters for which to calculate an error ellips. (see fixed=). For the two parameters it shows the major and minor axis, position angle of the one sigma deviation ellipse. [Default: not used]
beam=: The beam size (FWHM, in arcsec) for beam correction. One or two numbers required. Currently these are only used to correct error bars for the number of independant points per beam. If not given, each point is assumed independant. [no correction].
dens=image_den|true|false: Image containing the density. From this local derivatives (dN/dx)/N and (dN/dy)/N are computed numerically, and used for an estimate of beam smearing corrections. See also testbsc(1NEMO) . Note that for bizarre reasons this keyword can be used to trigger reading a weight column in tabular input mode (Bimagemode=). [Default: required if a beam is supplied].
wtmap=wt_image: If given, the weights (like 1/RMS^2) of the pixels.
tab=: If specified, this output table is used in append mode! This table can be piped through tabcomment(1NEMO) and fed to ccdvel(1NEMO) to create model velocity fields. [Default: not used].
resid=: If specified, this output will either contain a residual image (OBS-FIT) if the input file was an image, or a table with X, Y, Velocity and Residual Velocities in each ring. Each ring will be separated from the previous one with a simple commented line specifying which ring it refers to. [Default: not used].
tol=: Tolerance for convergence of nllsqfit [Default: 0.001].
lab=: Mixing parameter for nllsqfit [Default: 0.001]
itmax=: Maximum number of allowed nllsqfit iterations [Default: 50]
units=: Units of input axes for radius and velocity. Valid options are deg, arcmin, arcsec, rad for radius. A numeric value can also be given, in which case your image pixel separation from the image header is multiplied by this number to get to the ‘‘arcsec’’ that will be quoted in the tables. The units for velocity can only be numeric, and will be the factor by which the velocities in the map are multiplied. [Default: deg]
blank=: Value of the blank pixel that needs to be ignored. [Default: 0.0].
inherit=t|f: Logical denoting if the initial conditions for subsequent fitted rings should be inherited from the previous successfully fitted ring. The fixed parameters keep of course their fixed value. [Default: t]
reuse=t|f: Reuse pixels between rings. If neighboring rings have a different geometries, it can occur that pixels will be reused. This flag will prevent that. Note, if you don’t reuse pixels, rotcur will more likely produce a different rotation curve if you start at outer rings,e.g. radii=100:0:-10. [Default: t]
fitmode=cos|sin,1
nsigma=: Reject outlier points will fall outside nsigma times the dispersion away from the mean velocity in a ring. By default, it will not reject any outliers.
imagemode=t|f: Image input file mode? By default the input file is an image, alternatively a simple ascii table with X and Y positions in columns 1 and 2, and radial velocities in column 3, and optional errors in the radial velocity in column 4 (activated by setting dens=t). [Default: t]
wwb73=t|f: Use a simpler WWB73 (Warner, Wright, Baldwin 1973) linear method of fitting? [false]

AWK

The standard output is normally not very useful; it displays, for each iteration, the run of parameters plus the number of points and mean error in the ring. The following awk(1) scripts may be useful to extract information per ring: iter. number, vsys, vrot, pa, icn, xpos, ypos, npoints, sigma_vel.

BEGIN{count=0;line="";}
{
  if ($1 == "radius"){
    if (count != 0){
      printf("%s  %s0,rad,line);
      rad=$4;
    }else{
      count=1; 
      rad=$4;
    }
  }else{
    line=$0
  }
}
END{printf("%s  %s0,rad,line);}

Example

Here is an example of creating a synthetic velocity field with ccdvel, and analysing it with rotcur:

    % set r=‘nemoinp 0:100:5‘
    % set v=‘nemoinp 0:100:5 | tabmath - - "100*%1/(20+%1)" all‘
    % ccdvel out=map1.vel rad="$r" vrot="$v" pa=30 inc=60
    % rotcur in=map1.vel radii=0:100:5 vrot=0:100:5 pa=30 inc=60 vsys=0
tab=map1.rotcur units=arcsec,1
    % head map1.rotcur
  radius   systemic   error  rotation   error position    error   inclination
 error   x-position   error   y-position   error
          velocity           velocity           angle               angle
              of center            of center
 (arcsec)   (km/s)   (km/s)   (km/s)    (km/s)(degrees)  (degrees)  (degrees)
  (degrees) (grids w.r.t. (0,0))  (grids w.r.t. (0,0))
     7.50      0.00     0.44     27.08    0.50    30.45       1.09      59.72 
     2.47      63.50      0.20      63.50      0.30 96
    12.50     -0.00     0.24     38.54    0.28    30.66       0.41      59.27  
    0.94      63.50      0.13      63.50      0.20 150
    ...
    % tail map1.rotcur
 average inclination      :    59.93  (   0.043)  degrees
 average position angle   :    30.05  (   0.044)  degrees
 average systemic velocity:    -0.00  (   0.001)  km/s
 average x-position       :    63.50  (   0.000)  grids
 average y-position       :    63.50  (   0.001)  grids

Bugs

Failures in nllsqfit are not handled gracefully, and may error out the program. Usage of the error= system keyword can be used to bypass such bad rings, use with caution though and study the output.

Errorbars quoted in the table are only an estimate since the beam size is not known. Multiply these numbers by the square root of the number of pixels per beam to get a more realistic estimate. If the beamsize, beam=, is given, the formulae of Sicking (1997) was previously used to correct the errors for:

    factor = sqrt(4.PI.B_x.B_y/(D_x.D_y))

but this is clearly too large, instead

    factor = sqrt(PI/(4ln2).B_x.B_y/(D_x.D_y))

is now used.

Sign of the pixelsize Dx,Dy in the CCD header is ignored, and an astronomical image is assumed. See also the reuse= keyword.

For calculations of residuals in overlapping rings (e.g. warps) only the last ring velocity will be used.

Author

K. Begeman (original GIPSY Sheltran version, now also available in C), P. J. Teuben (NEMO C version)

Update History

19/jul/83    original program                         KGB
9/mar/85    revision of program                     KGB
23/may/86    migrated to VAX-VMS                      KGB
27/nov/88    UNIX version                               pjt
8-feb-91    flushed buffers ROTCUR.DAT each write    pjt
30-apr-91    moved to NEMO into C                       pjt
10-sep-91    documentation improved               pjt
17-oct-91    added Npts to table output          pjt
21-may-92    added Bob Gruendl’s rotcur awk scripts    PJT
12-jun-92    added inherit=t as default          PJT
13-aug-92    implemented fitmode= without XPOS,YPOS    PJT
15-oct-99    compute residuals and add resid=    PJT
14-mar-01    V2.5: clarifications, added nsigma=        PJT
9-may-01    V2.6a: corrected error correction factor    PJT
10-aug-01    clarified some differences between NEMO and GIPSY versions    PJT
26-jan-02    Added unit (scale factor) for velocity too    PJT
26-jun-02    V2.8: added tabular input for irregular spaced data, fixed example    PJT
11-sep-02    V2.9: implemented map residual velocity field    PJT
30-jan-03    V2.10: allow tables to use error in velocity    PJT
2-jun-04    V2.12: finally implemented the reuse= option    PJT
6-jun-20    V2.13: added wtmap=    PJT
18-jan-21    V2.14: beam error factor back to standard, not Sicking    PJT

Table of Contents

Name
Synopsis
Description
Parameters
AWK
Example
Bugs
See Also
Author
Update History