Table of Contents

## Name

rotcurshape - fit parameters to a velocity field using shaped rotation curves in a co-planar disk

## Synopsis

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

## Description

rotcurshape derives the rotation curve and kinematic parameters from an observed velocity field of a coplanar galactic disk by fitting a (set of) rotation curve shape(s). Instead of dividing the galactic disk in a set of rings (see rotcur(1NEMO) ) it fits the shape of the rotation curve in an annulus of a planar disk. This is particularly useful for low inclination galaxies and/or for central regions of a galaxy where normally the rotation curve would vary accross the ring/annulus.

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) = -----------------------------------------
r
```
and where the radius r is measured in the plane of the galaxy:

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).

## Parameters

The following parameters are known (those with ** are still here since we cloned this program off rotcur, but are likely to be slashed or (re)implemented):
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.

## Rotcur Functions

A rotcur function needs to provide a routine that returns a rotation curve, as well as all of its partial derivates w.r.t. the parameters. Most rotation curve shapes have two parameters, a velocity and a radial scale parameter, and they are usually the first and second parameter. Note that the user needs to know how many parameters a rotcur function has. Apart from a number of pre-defined ones, the user can write his/her own in the C language and loaded via the load= keyword (see also loadobj(3NEMO) .

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 = r;
return p * r;
}
real rotcur_plummer(real r, int np, real *p, real *d)
{
real x = radius/p;
real y = pow(1+x*x,-0.75);
d = y;
d = -x*p/p*(1-x*x/2)/(1+x*x)/y;
return p * x * y;
}
```

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

```Name:    Parameters:    Formula:                Comments:
-----    -----------    --------                ---------
linear    omega       v=omega*r                -
flat    v0           v=v0                    should return rotcur solution
plummer    v0,a        v=v0*x/(1+x^2)^(3/4)            -
core1    v0,a          v=v0*x/(1+x)                -
core2    v0,a        v=v0*x/(1+x^2)^(1/2)            -
core    v0,a,c      v=v0*x/(1+x^c)^(1/c)            -
arctan    c0,a        v=2*v0/pi*arctan(x)            -
poly    v0,a,p2,..     v=v0*(x+p2*x^2+p3*x^3+.....pN*x^N)     a needs to be fixed !
iso    v0,a        v=v0*sqrt(1-atan(x)/x)            -
exp     v0,a        v=v0*(1-exp(-x))                     -
nfw    v0,a,c      v=v0*sqrt((ln(1+cx)-cx/(1+cx))/x/(ln(1+c)-c/(1+c)))    needs a fixed
’c’
power    v0,a,p      v=v0*x^p                             a needs to be fixed (scale
free)!
```

## Example

First an example of creating a synthetic velocity field with ccdvel(1NEMO) , and analysing it with rotcurshape, by using a simple rotation curve entered by a sneaky construction in shell variables \$r and \$v. The shape function is a core1 with amplitude 200 and core radius 20:
```    % 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:  5658
(this is also the Testfile version)
```
The reason why the error bars are not exactly zero is because ccdvel interpolates from the (r,v) table, and the finite pixel sizes.

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 = r;
return p * 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
```

## Files

```func_rotcur.c    helper 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.

## See Also

rotcurtab(1NEMO) , rotcur(1NEMO) , ccdvel(1NEMO) , tabcomment(1NEMO) , rotcurves(1NEMO) , pvtrace(1NEMO) , ccdmom(1NEMO) , tabnllsqfit(1NEMO) , rotcurshape(5NEMO) , loadobj(3NEMO) , gal(AIPS), nllsqfit(3NEMO)
```GBKFIT: http://supercomputing.swin.edu.au/projects/gbkfit/
```

Peter Teuben

## Update History

```20-jul-02    0.9 cloned off rotcur        PJT
10-sep-02    1.0 implemented an image version for resid=    PJT
19-sep-02    1.0d added exp (vMoorsel&Wells 1985), and nfw    PJT
13-dec-02     1.0h added (corrected) power law rotcur    PJT
11-oct-03    fixed up for adass 2003      PJT
26-may-04    1.2e fixed sqrt(n) bug in sigma, improved nsigma>0 iterations    PJT
13-jun-04    1.3: added fit= option to save a fitted map    PJT
30-jan-08    1.4: with new rotcurtab minor overhaul of code    PJT
28-may-20    added arctan    PJT
```