HTML automatically generated with rman
Table of Contents
rotcurshape - fit parameters to a velocity field using shaped rotation
curves in a co-planar disk
rotcurshape in=velfie [parameters=values
...]
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).
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.
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[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:
----- ----------- -------- ---------
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)!
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[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.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.
rotcurtab(1NEMO)
,
rotcur(1NEMO)
, ccdvel(1NEMO)
, tabcomment(1NEMO)
, rotcurves(1NEMO)
, pvtrace(1NEMO)
,
ccdmom(1NEMO)
, tabnllsqfit(1NEMO)
, rotcurshape(5NEMO)
, loadobj(3NEMO)
,
gal(AIPS), nllsqfit(3NEMO)
, rotcur(5NEMO)
GBKFIT: http://supercomputing.swin.edu.au/projects/gbkfit/
Peter Teuben
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
Table of Contents