Table of Contents

Much like *rotcur*,
this program does a non-linear least-squares-fit to the function:

v(x,y) = VSYS + VROT(r) * cos(theta) * sin(INC) where again - (x-XPOS) * sin(PA) + (y-YPOS) * cos(PA) cos(theta) = ----------------------------------------- rand where the radius

The parameterized rotation curve VROT(r;p1,p2,...pN) can be selected from a number of predefined shapes or dynamically loaded from a C-coded routine (see ROTCUR FUNCTIONS below). Up to 5 functions can be added in quadrature (but be aware of degeneracies)

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 (note that MIRIAD and GIPSY
all have different conventions where 0,0 is). *rotcurshape*

can fit the 5 geometric parameters *VSYS*, *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 (see also
EXAMPLES below).

**in=***v_image|xyv_tab|rv_tab*- Input velocity
field map, normally in
*image(5NEMO)*format. Input data can also be a table, see**imagemode=**and**rotcurmode=**below how to control the input table. No default. **radii=***r0,r1*- Inner and outer radius of the disk to fit (in arcsec).
For images, the
**units=**keyword can be used to scale your physical units in the header to ‘‘human readable’’ units (arcsec). By default the whole disk will be used. **pa=***pa*- Kinematic position angle (in degrees) of the disk (i.e. receding side). Angles can be 0..360 or -180..180.
**inc=***inc*- Inclination (degrees) of the disk. Angles can be 0..90.
**vsys=***v0*- Systemic velocity.
**center=***x0,y0*- Rotation center (grid units w.r.t. lower left corner being 0,0 for NEMO maps). Two numbers are required. For tabular input the units must be the same as those used in the table. 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:**0.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**. If a density is given (see**dens=**), the weight function is multiplied by this geometric weight factor. [Default:**uniform**]. **fixed=**- List of disk
geometry parameters to be kept fixed during the fit. Choose any of the following:
**vsys,pa,inc,xpos,ypos**, although**all**can also be use to fix all geometric parameters. [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) 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]. ** See CAVEATS below **
**dens=***image_den*- Image containing the density. Currently only used as an additional weight factor. If a map of velocity dispersions would be available, a density map would need to be contstructed from its inverse square, i.e. dens=1/sigma^2. Default: not used.
**tab=**- If specified, this
output table is used in
*append*mode! **deprecated** [Default: not used]. **resid=**- Residuals (OBS-FIT) if specified. If the input file was a table (see
also
**imagemode=**), this will contain a listing of X, Y, Velocity and Residual Velocities in the disk, else it will be an*image(5NEMO)*of the residuals. Note for**nsigma>0**the residual map will be two maps, and you will need*csf(1NEMO)*to get access to the 2nd (supposedly improved) residual map. [Default: not used]. **fit=t|f**- Controls if the residual output map (
**resid=**should contain the residuals [the default], or the fitted map. If the fit is selected, the mask values are ignored and the map is filled with the fitted values. [Default:**f**] **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. For tabular input the spatial scaling factor is ignored, since all spatial coordinates need to be in the same coordinate system. [Default:**deg**] **blank=**- Value of the blank pixel that needs to be ignored.
[Default:
**0.0**]. **nsigma=**- Reject outlier points will fall outside nsigma times the dispersion away from the mean velocity in the disk. By default, it will not reject any outliers. Use with care, only useful if close enough to a solution and just a few outliers need to be removed.
**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 can be used by the
*xyv_tab*(see**in=**). The units of the spatial coordinates now need to be the same as**center=**, and the**units=**factor is ignored in this case. Future expansion will likely allow weight factors to be added in other columns. [Default: t] **rotcurmode=t|f**- Input table
is a rotation curve (R,V) in columns 1 and 2. Radius is allowed to be negative,
as this option is implemted as a special version of XYV where we fix XPOS=0,INC=30,PA=0.
The center can therefore be fitted via YPOS. Note that
*tabnllsqfit(1NEMO)*is also quite efficient to use, except the rotation curve functions would need to be re-written in their function interface. See also func_rotcur.c for a useful helper routine. **load=***so_file*- Name of a shared object file containing
rotation curve(s). The function names must be
**rotcur_***name*, where*name*is the identifier name of the rotation curve used in the subsquent**rotcur#=**keywords. **rotcur1=***name1,p1,p2,...pN,m1,m2,..mN*- Name of first rotation curve, followed by the initial estimates of its parameters (you do need to know how many there are), followed by an equal number of 1s (free) and 0s (fixed) to denote which parameters are free or fixed during the fitting process.
**rotcur2=***name2...*- see rotcur1
**rotcur3=***name3...*- see rotcur1
**rotcur4=***name4...*- see rotcur1
**rotcur5=***name5...*- see rotcur1. The final composite rotation curve will
be the sum (in quadrature) of up to these 5 components.

Here are two examples,
a simple linear rotation curve with one parameter, and a slightly more
involved Plummer disk/sphere rotation curve with two parameters:

#include <nemo.h> real rotcur_linear(real r, int n, real *p, real *d) { d[0] = r; return p[0] * r; } real rotcur_plummer(real r, int np, real *p, real *d) { real x = radius/p[1]; real y = pow(1+x*x,-0.75); d[0] = y; d[1] = -x*p[0]/p[1]*(1-x*x/2)/(1+x*x)/y; return p[0] * x * y; }

Here is a list of the builtin rotation curves, where x=r/a is the dimensionless
radius:

Name:Parameters:Formula:Comments: --------------------------------- linearomegav=omega*r- flatv0v=v0should return rotcur solution plummerv0,av=v0*x/(1+x^2)^(3/4)- core1v0,av=v0*x/(1+x)- core2v0,av=v0*x/(1+x^2)^(1/2)- corev0,a,cv=v0*x/(1+x^c)^(1/c)- arctanc0,av=2*v0/pi*arctan(x)- polyv0,a,p2,..v=v0*(x+p2*x^2+p3*x^3+.....pN*x^N)a needs to be fixed ! isov0,av=v0*sqrt(1-atan(x)/x)- expv0,av=v0*(1-exp(-x))- nfwv0,a,cv=v0*sqrt((ln(1+cx)-cx/(1+cx))/x/(ln(1+c)-c/(1+c)))needs a fixed ’c’powerv0,a,pv=v0*x^pa needs to be fixed (scale free)!

% set r=‘nemoinp 0:60‘ % set v=‘nemoinp 0:60 | tabmath - - "100*%1/(20+%1)" all‘ % ccdvel out=map1.vel rad="$r" vrot="$v" pa=30 inc=60 % rotcurshape in=map1.vel radii=0,60 pa=30 inc=60 vsys=0 units=arcsec,1 \ rotcur1=core1,100,20,1,1 tab=- VSYS: 2.36846e-18 0.00110072 XPOS: 63.5 0.000759475 YPOS: 63.5 0.00100543 PA: 30 0.0010416 INC: 60.0001 0.00229122 P1: 100.392 0.00757645 P2: 20.2883 0.0045192 NPT: 5658The reason why the error bars are not exactly zero is because ccdvel interpolates from the (r,v) table, and the finite pixel sizes.(this is also the Testfile version)

Here is an example to
write your own C code with a rotation curve, and load it in during runtime:
(examples are in $NEMO/src/image/rotcur/shape):

% cat mylinear.c #include <nemo.h> real rotcur_linear(real r, int n, real *p, real *d) {d[0] = r;return p[0] * r; } % bake mylinear.so % rotcurshape in=map1.vel radii=0,10 load=mylinear.so rotcur1=linear,10,1

Here is a contrived example of creating a velocity model field with rotcurshape
by supplying a zero map, fixing all parameters, and computing the residual
velocity field (OBS-FIT). Of course you will get -1 times the velocity field,
but still. It is an alternative to *ccdvel(1NEMO)*

% ccdmath out=map0.vel fie=0 size=128,128 % rotcurshape map0.vel 0,40 30 45 0 blank=-999 resid=map.vel \rotcur1=plummer,200,10,0,0 fixed=all units=arcsec,1

When **nsigma>0** is used to find an improved solution, the residual map now
contains 2 maps, and it is the 2nd map which contains the supposed improved
map. Here is an example how to extract and display that 2nd map using *csf(1NEMO)*

% rotcurshape map0.vel 0,40 30 45 0 blank=-999 resid=map.vel \rotcur1=plummer,200,10,0,0 fixed=all units=arcsec,1 nsigma=2 % csf map.vel map2.vel item=Image select=2 % nds9 map2.vel

func_rotcur.chelper routine for tabnllsqfit to use rotcur functions $NEMO/src/image/rotcur/shape/directory with other example shape functions

CaveatParameters and fix/free masks to rotation curve parameters should be all set, in order for the respective routines to figure out the correct number of parameters. For example, the poly rotation curve can only determine the correct order of the polynomial by counting the number of arguments given in that option, e.g. rotcur1=poly,100,10,2,2,1,1,1,1 would use a 3th order polynomial. rotcurshape does surprisingly bad on exact data, and often complains about taking the sqrt of a negative number. Adding a little noise will speed up convergence! rotcurshape sometimes needs inititial conditions pathetically close to a minimum to converge, and more than often complains with the message ### Warning [rotcurshape]: nllsqfit=-4: must find better solution (n=225) ### Warning [rotcurshape]: ROTCUR: problems with matrix inversion

if **beam=** is used, the map is also used to estimate beam smearing corrections.
This is still a totally untested feature of the code.

GBKFIT: http://supercomputing.swin.edu.au/projects/gbkfit/

20-jul-020.9 cloned off rotcurPJT 10-sep-021.0 implemented an image version for resid=PJT 19-sep-021.0d added exp (vMoorsel&Wells 1985), and nfwPJT 13-dec-021.0h added (corrected) power law rotcurPJT 11-oct-03fixed up for adass 2003PJT 26-may-041.2e fixed sqrt(n) bug in sigma, improved nsigma>0 iterationsPJT 13-jun-041.3: added fit= option to save a fitted mapPJT 30-jan-081.4: with new rotcurtab minor overhaul of codePJT 28-may-20added arctanPJT