Program: RECTIFY
Purpose: Review the geometrical component of a model galaxy
by fitting surface densities in (user supplied) tilted
rings.
Category: ANALYSIS, MODELS
File: rectify.c
Author: M. Vogelaar
Keywords:
INSET= Give (sub)set with data to fit:
Maximum number of subsets is 1. The dimension of the
input structure is 2.
This set contains to data used to fit the rings.
It must be possible to convert the units of its subset
axes to seconds of arc.
BOX= Give BOX in ..., .... [entire subset]
Limit the size of the input data from INSET=
The dots correspond to the axis names of the input
(sub)set. A box bigger than the frame of the subset
is not allowed.
OUTSET= Give output set (,subset):
With the fitted model parameters, a set is created
with the same structure as the input set. Only the
part defined with BOX= is written to the output set.
Data outside any rings will be blank in the output
(unless the output set already existed and contained
non blank values at these positions).
PROSET= Give projected output set (,subset): [no set]
The fitted model can be projected onto the plane
of the central ring. The name of this projection
is given with PROSET= The grid spacings of this set
are equal in both directions. The position angle
of the central ring is aligned with the y-axis.
KEYWORDS RELATED TO THE MODEL:
POSITION= Central position (of galaxy):
All positions in the model will be calculated wrt.
this position so it will be the center of all rings.
The position is entered as grids or physical
coordinates.
A two dimensional subset with axes RA and DEC:
POSITION=10 5 Grid (RA,DEC)=(10,5)
POSITION=U 45.8 U 20.02 Phys. coord. at RA=45.8,
DEC=20.02 in
units of the RA, DEC axes.
POSITION=PC RA, DEC of projection centre
POSITION=AC RA, DEC of axis centre
POSITION=* 10 12 8 * -67 8 9.6
(RA=10h12m8s, DEC=-67d8m9.6s
in epoch of set)
RADII= Give (outer) radius of rings (arcsec):
Maximum number of rings is 2048, but this maximum is
also limited by the available memory.
SEGMENTS= Give number of segments per ring for n rings: [calc.]
If you assume axial symmetry for your model, there is
no need to define segments in a ring: enter SEGMENTS=1.
Else, give for each defined ring the number of segments
that you want to use. If you give less than 'n' values
the remaining number of segments will be copied from
the last one entered. The default for ring k with
k=1..n is: segments = int((k-1)/4)*4
POSANG= Give n position angles (deg. N-->E):
(See description)
n is the number of rings defined with RADII=
If you give less than 'n' angles, the remaining angles
all get the value of the last specified angle.
Values for the position angles and inclinations could
(for example) be obtained with program ROTCUR.
INCLINATION= Give n inclinations (degree):
(See description) If less than n items are entered, then
the missing items are copied from the last one entered.
Z0= Give n scale heights (arcsec):
(See description) If less than n items are entered, then
the missing items are copied from the last one entered.
PROFILE= Give dens.prof. perpend. to plane of the rings: [list]
Complete tilted ring model with a vertical structure
for the HI layer.
The options are:
1 gaussian layer ( dispersion = 1*Z0(ring) )
2 sech2 layer.
3 exponential layer.
4 Lorentzian layer.
5 box layer.
ISEED= Seed (should be negative): [-1]
Different seeds generate different random number
sequences.
MAXCLOUDS= Max. number of Monte Carlo clouds: [100000]
Use MAXCLOUDS= random positions in a ring for
the Monte Carlo simulation.
See also description at EQDENS=
EQDENS= Equal density of M.C. points in all rings? Y/[N]
The number in MAXCLOUDS= is used to set the number of
Monte Carlo points that cover the largest ring area.
If EQDENS=Y and if the area of a ring is smaller
than the area of the largest ring, then scale the number
of generated M.C. points so that equal areas in different
rings will approximately contain the same number of
clouds.
KEYWORDS RELATED TO PLOTTING
============================
GRDEVICE= Plot device: [List of devices]
Destination of plot, Screen or Hardcopy.
** MOSAIC= View surface subdivisions x,y: [1,1]
View surface can contain a number of plot pages in
in X and Y direction (mosaic). Default is 1 plot in
both X- and Y direction. With this keyword, you can
force to display more than one perspective view of
your model on one physical device page.
** LINEWIDTH= Give line width (1-21): [1]
The width of lines etc. that are plotted, can be
controlled by this keyword. For a hardcopy, the default
is 2.
RADPLOT= Plot radius vs. pos.ang/inclination? Y/[N]
Plot input radius vs. input position angle and
radius vs. inclination in separate plots.
PLOTMODEL= Draw rings in perspective view? Y/[N[
Display model in a perspective view, using
parameters defined in DISTANCE= RHO=, PHI= & THETA=
MCPLOT= Monte Carlo plot? Y/[N]
Create a 2d-plot with positions generated by the
random number generator.
TRPLOT= Tilted ring overview? Y/[N]
For an overview of the tilted ring, 1) plot the
components of the spin vectors in the xy plane of the
central disk (first ring) and 2) plot the intersection
of all rings with the central disk. Just before
plotting, keyword GRDEVICE= is asked, so you are able
to change the output destination for this plot.
DRPLOT= Plot surface density vs radius? Y/[N]
Plot the fitted parameters agains radius. Draw a
horizontal line from inner radius to outer radius and
mark outer radius. Just before plotting, keyword
GRDEVICE= is asked, so you are able to change the
output destination for this plot.
COMPVALS= Values to include in plot: [None]
If you plotted the fitted ring parameter vs. radius
(DRPLOT=Y) it can be handy to include some data
for comparison. The data is given as a list of
reals (entered manually or via a default file or a
recall file). The maximum number of values that can
be specified is equal to the maximum number of rings.
VIEW RESULTS OBTAINED IN PREVIOUS RUN
=====================================
To view previously obtained results stored in a GDS
table in a GIPSY descriptor file, use the following
keywords:
** TABSET= View table results from set: [skip plot]
Name (no subset level) of set containing
table data. This must be a set previously made
by RECTIFY.
TABNAME= Give name of GDS table: [stop]
The name of the table (not a column name) can be
obtained by running the task TABLE. In most cases the
table name will be RECTIFY, but if you used a default
file for RECTIFY, it will be the name of that default
file. After TABSET= and TABNAME= a plot will appear
with the fitted values (densities) as function
of radius and the keyword PLOTSEG= will be prompted.
PLOTSEG= Ring number, max. segments: [stop]
This keyword is asked in a loop that can be aborted
by pressing carriage return. The first number is the
ring number. The second is the maximum number of
segments that you want to display in your plot.
The default for the second number is the total number
of segments that was defined for this ring.
KEYWORDS RELATED CREATING GDS TABLES
====================================
TABNAME= Give name of GDS table to store results: [Name of task]
A table will be stored in the output set OUTSET=
The maximum length of the table name is less than 9
characters. Columns in a table are created on set level.
The tables can be viewed and manipulated with
program TABLE. These tables are not the tables that
are written to your screen and GIPSY log file but have
similar contents.
OVERTAB= Overwrite table? [Y]/N
If a table exists, ask user permission to overwrite.
If OVERTAB=N, a new name is asked.
KEYWORDS RELATED TO THE VIEWING TRANSFORMATION
==============================================
It is possible to view your rings in a 3-plot. The
viewing transformation is controlled by the foll
DISTANCE= Distance of the eye from the screen: [calculated]
Control amount of perspective with this
unitless number.
** RHO= Distance between viewpoint and origin. [calculated]
THETA= Angle view vector wrt. pos. x-axis (deg): [30]
Angle between view vector and positive x-axis wrt.
the positive x-axis (degrees).
x-axis is equivalent to the first subset axis.
PHI= Angle view vector wrt. pos. z-axis (deg): [60]
Angle between view vector and positive z-axis wrt.
the positive z-axis (degrees).
PHI= ranges from 0 to 180 deg.
z-axis is equivalent to the image value axis.
REPEAT= Repeat perspective view? [Y}/N
Replot same subset with new values for DISTANCE=
RHO=, PHI=, THETA=. The defaults are the previous values
of these keywords.
Description: RECTIFY
Introduction
============
RECTIFY is used as part of a series of programs to
create a model galaxy which is a 'best fit' to an
existing HI data cube. You can review the geometrical
component of your model galaxy after RECTIFY fitted
the surface densities in the tilted (and perhaps
segmented) rings that are part of that model.
So RECTIFY determines fitted surface densities in a
tilted ring model as function of radius and azimuth
(SEGMENTS=). The program uses data from an existing
HI map (INSET=) to fit the parameters of the model,
which is a set of rings with radius, position angle,
inclination, scale height (Z0=) and a definition for the
assumed vertical structure for the HI layer (PROFILE=).
RECTIFY fits the surface densities per ring or per segment
and transforms the model into a set (OUTSET=) so that it
will be possible to compare the model with the original
set. The fitted segments can also be projected onto the
plane of the central disk, the name of this set is given
in PROSET=
Model parameters
================
A (simplified) tilted ring model consists of a number
of concentric rings with arbitrary spatial orientations.
A ring is characterized by an inner and outer radius.
The outer radii are given by RADII=. The inner radii
are calculated. The first disk has inner radius equal
to 0 and outer radius equal to the first number in RADII=.
The orientation of a ring is determined by two angles.
1) POSANG=. The position angle of the major
axis of a galaxy is defined as the angle taken in
anti-clockwise direction between the north direction
in the sky and the major axis of the receding half
of that galaxy (Rots 1975) astron, astrophys 45, 43.
2) INCLINATION= The inclination of the ring with respect
to the plane of the sky.
The tilted ring model is completed with a choice
of the width of the HI layer (Z0=) for each ring and
a vertical structure (PROFILE=) for all rings. The
options for this vertical structure are:
1 gaussian layer ( dispersion = 1*Z0(ring) )
2 sech2 layer.
3 exponential layer.
4 Lorentzian layer.
5 box layer.
Each of these distributions can be plotted and tested
with the GIPSY program RANDOM.
If all ring parameters are given, a table is written in
the Log file.
(Non) axial symmetry
====================
If you assume axial symmetry, there is no need to define
segments in a ring. You must tell to program to ignore
the existence of segments with SEGMENTS=0
However, for non axial symmetry, segments can be defined
for each ring. The last value that you enter with
SEGMENTS= is copied for all remaining radii that did not
get a number of segments with SEGMENTS= There is a
suitable default for this keyword. Assuming that the
width of the rings is one bundle and the index number
of a ring is k (k=1..nrings) then the default number of
segments for that ring is: s = int((k-1)/4)*4 (if k = 1
then s = 1)
resulting in a series s = 1 1 1 1 4 4 4 4 8 8 8...32
with a maximum of 32.
Examples: (You entered: RADII= 1 2 3 4 5 6 7 8)
SEGMENTS= 2 2 4 6 8 ==> 2 2 4 6 8 8 8 8 .... 8
SEGMENTS= 1 ==> 1 1 1 1 1 1 1 1 .... 1
SEGMENTS= ==> 1 1 1 1 4 4 4 4 .... 32
In practice the total number of segments will be limited
by the available memory. This amount will be different on
different machines.
Method
======
The surface densities that we want to fit are calculated
using a Monte Carlo simulation. A certain number of
'clouds' given in MAXCLOUDS= are randomly distributed
over a ring. Each position corresponds to a position
of a pixel in the input set. If we count the number of
hits in each pixel and divide that number by the number
of generated clouds in the ring, it is possible to use
this value as a matrix element containing weights in a
system with linear equations:
A(1,1)m(1) + A(1,2)m(2) + ... + A(1,NR)m(NR) = u(1)
A(2,1)m(1) + A(2,2)m(2) + ... + A(2,NR)m(NR) = u(2)
......................................................
A(i,1)m(1) + A(i,2)m(2) + ... + A(i,NR)m(NR) = u(i)
......................................................
A(NP,1)m(1) + A(NP,2)m(2) + ... + A(NP,NR)m(NR) = u(NP)
NP = Number of pixels in the input box (BOX=)
NR = Number of rings
A = Matrix element containing weights
m = Wanted surface densities
u(i) = Pixel value of pixel from input with index i
A cloud has a random position determined by a random
radius between the inner radius (included) and the
outer radius (excluded) and a random azimuth between 0
and PI. But if the radius increases, also the area in-
creases, so we have to correct the uniform random radius:
To fill a circle uniformly, the correction would be:
(rnd is an uniform random number between 0 and 1)
Radius = Rhi*sqrt( rnd )
For a ring with inner radius Rlo and outer radius Rhi
we find the correction formula:
Radius = sqrt( (Rhi^2-Rlo^2)*rnd + Rlo^2 )
A second correction can be applied for the decreasing
density of points for increasing areas. The correction
is initiated with EQDENS=Y. The number in MAXCLOUDS=
is used to cover the largest area. For other areas
the number of random clouds is smaller. If EQDENS=N,
the number in MAXCLOUDS= will be used for all rings.
Radius and angle correspond to a position in the
xy plane (in the xyz system of a ring). The value of
z is randomly selected from the distribution given
in PROFILE=. The characteristic width of these
distributions is equal to the value of Z0 of the ring
that is examined.
Because NP >> NR, the system is 'overdetermined' and
we use a least squares method to solve m(i). The method
is described in Numerical Recipes in C (Press cs., 2nd ed.)
$15.4, page 674. The Matrix solution is obtained with
a Gauss-Jordan elimination, Num. Rec. $2.1 pag 39.
Segments are treated as rings. The number NR has to be
interpreted as the number of segments in your model and
not as the number of rings anymore.
Output
======
The output set OUTSET= is a set made with the fitted
parameters, i.e. the matrix elements m(j) are known
after the least squares fit, the normalized matrix A(i,j)
is left unchanged and therefore the values of u[i] can
be calculated. This vector u[i] is in fact the output
set in OUTSET=
A set which contains an image of the input set as seen
from a position perpendicular to the plane of the central
disk (ring with first input radius, position angle and
inclination) is created with PROSET= Random positions
in the rings are converted to positions in the plane
of the central ring and so these rings are sampled in
that plane. The matrix A(i,j) stores the hits at all
possible positions and after normalization of 'A' a new
vector 'u' is calculated using the fitted parameters 'm'.
The number of generated random positions is equal related
to the number in MAXCLOUDS=
Log file
========
Results are written in the GIPSY log file. A table shows
for each segment the fitted parameters (surface densities)
and their estimated errors and the sum off all pixels in
a segment (and the errors).
The minimized Chi square, the reduced Chi square and an
estimated measurement error per pixel are also listed.
Tables
======
The results of the fitting are saved in a GDS table
with the name of the current task as default name.
This name is changed with TABNAME= An existing table can
be overwritten but you are always prompted with OVERTAB=
before you can do so. The table has the following column
names:
LORAD -Inner radius in arcsec
HIRAD -Outer radius in arcsec
WIDTH -Width in arcsec
POSANG -Position angle of ring in degrees
INCLIN -Inclination of ring in degrees
AREA -Area of ring in arcsec^2
SURFDENS -Fitted surface density in map units of INSET=
SURFDERR -Errors in covariance matrix corresponding to
-fitted densities
Z0 -Scale height of ring in arcsec
Contents of the columns can be viewed with program TABLE.
But RECTIFY also has a possibility to view previously
stored data. To enable this feature start RECTIFY with
TABSET= and TABNAME= You get a plot of all stored
surface densities. In a loop it is possible to plot
the values of the segments in a ring. The keyword
PLOTSEG= accepts a ring index number and a second number
that specifies the maximum number of segments that you
want in your plot. The loop (and the program) is aborted
by pressing carriage return.
Plots
=====
1) Position angles and inclinations
Keyword RADPLOT=Y shows two plots on one page. The
first is a plot of radius versus position angle and
the second plot shows radius versus inclinations.
The data used in this plot is written in the log file.
2) The model
A view of the defined model (set of circles) from an
arbitrary angle can be obtained with PLOTMODEL=Y.
First a viewpoint has to be specified (RHO=, THETA=,
PHI= & DISTANCE=, angles in degrees, distance in grids)
The 'viewpoint' transformation
converts positions in world coordinates (== ring coordinates)
into 'eye' coordinates expressed in a coordinate system
centered at the viewpoint. A perspective transfor-
mation produces the actual 2-dim. 'screen' coordinates.
It has a single vanishing point and the screen axes
are parallel to the 'eye' coordinates.
If you want to view the projection of the model on the
sky, use THETA=90 and PHI=180.
In the plot, a coordinate system that is aligned with the
sky (XY) and with the line of sight (Z) is displayed.
The plot also shows the xyz coordinate system
of the central disk and labels it with p0, q0, and s0.
The vector s is the spin vector of a ring. It is aligned
with the rotation axis of the rings. The vector p
is a unit vector aligned with the major axis of the
(projected) ring. The unit vector q is perpendicular
to p and s. For all rings the vectors p and s are plotted.
3) Overview tilted ring
With TRPLOT=Y, a 2d-plot of the tilted ring model is
created. The first plot shows the end points of the
spin vectors of each ring with respect to the central
spin vector (positioned at 0,0 ). A second plot shows
the intersections of the rings with the central disk.
4) Monte Carlo plot
If you want to view the random process of selecting
positions in a ring, use MCPLOT=Y. You will see how
the rings are build up. The plotting of positions is
slow so take a small number for MAXCLOUDS= if you
want to use this feature.
5) Surface density vs. radius
DRPLOT=Y will plot the fitted densities vs. radius.
A line is drawn between the inner radius to the
outer radius along the X-axis at height equal to the
fitted density.
If you want to compare the plotted data with an
array of numbers of another source (maybe some input
densities for another program), use COMPVALS=
to read this data. Then also these values will be plotted.
Make sure that the densities are in the same units.
Example: An example of a COLA file where a set of rings was fitted to
data in NGC3198 p 7 (subset param 7 == total HI map)
! The COLA file can also be used to start nhermes (or ngipsy
! if plots are created) as a batch.
! The syntax is: nhermes -lmylog batch
! The file 'mylog' is the log file for the RECTIFY output.
!
!===========================================================
"RECTIFY
INSET=NGC3198 7
OUTSET=rect_fit
PROSET=rect_depro
BOX=-101 -87 102 106
POSITION=-7.8359895314883 8.8999911586812
SEGMENTS=1::4 4::4 8::4 12::4 16::4 20::4 24
PROFILE=1
ISEED=-10
MAXCLOUDS=
EQDENS=
GRDEVICE=
MOSAIC=
LINEWIDTH=
RADPLOT=
PLOTMODEL=
MCPLOT=
TRPLOT=
DRPLOT=
COMPVALS=
TABSET=
TABNAME=
PLOTSEG=
TABNAME=
OVERTAB=
DISTANCE=
RHO=
THETA=
PHI=
REPEAT=
POSANG= 2.092830e+02 2.119370e+02 2.136150e+02 2.149110e+02
2.162340e+02 2.169150e+02 2.172880e+02 2.175830e+02 2.176790e+02
2.176830e+02 2.176370e+02 2.173810e+02 2.170590e+02 2.166300e+02
2.162810e+02 2.159940e+02 2.158050e+02 2.157010e+02 2.156520e+02
2.155100e+02 2.153090e+02 2.151090e+02 2.149150e+02 2.146560e+02
2.143530e+02 2.140620e+02 2.138230e+02 2.136580e+02 2.136440e+02
2.136890e+02 2.137610e+02 2.138460e+02 2.140110e+02 2.142470e+02
2.146760e+02
RADII= 1.800000e+01 3.600000e+01 5.400000e+01 7.200000e+01
9.000000e+01 1.080000e+02 1.260000e+02 1.440000e+02 1.620000e+02
1.800000e+02 1.980000e+02 2.160000e+02 2.340000e+02 2.520000e+02
2.700000e+02 2.880000e+02 3.060000e+02 3.240000e+02 3.420000e+02
3.600000e+02 3.780000e+02 3.960000e+02 4.140000e+02 4.320000e+02
4.500000e+02 4.680000e+02 4.860000e+02 5.040000e+02 5.220000e+02
5.400000e+02 5.580000e+02 5.760000e+02 5.940000e+02 6.120000e+02
6.300000e+02
INCLINATION= 7.645950e+01 7.349020e+01 7.204080e+01 7.132620e+01
7.107400e+01 7.086810e+01 7.073370e+01 7.054340e+01 7.044120e+01
7.031490e+01 7.027940e+01 7.033710e+01 7.035760e+01 7.025710e+01
7.005100e+01 6.993850e+01 6.990560e+01 6.992830e+01 7.004700e+01
7.015080e+01 7.015080e+01 7.042400e+01 7.064860e+01 7.103810e+01
7.149680e+01 7.163160e+01 7.202490e+01 7.224920e+01 7.247960e+01
7.268960e+01 7.286270e+01 7.295460e+01 7.305270e+01 7.327160e+01
7.393910e+01
Z0= 1.000000e+01 1.000000e+01 1.000000e+01 1.000000e+01
1.000000e+01 1.000000e+01 1.000000e+01 1.000000e+01 1.000000e+01
1.000000e+01 1.000000e+01 1.000000e+01 1.000000e+01 1.000000e+01
1.000000e+01 1.000000e+01 1.000000e+01 1.000000e+01 1.000000e+01
1.000000e+01 1.000000e+01 1.100000e+01 1.200000e+01 1.400000e+01
1.600000e+01 1.800000e+01 2.000000e+01 2.200000e+01 2.400000e+01
2.600000e+01 2.800000e+01 3.000000e+01 3.200000e+01 3.400000e+01
3.600000E+01
COMPVALS=
3.443390e+20/(1.7114975189653e19*cos(RAD(7.645950e+01)))
4.024280e+20/(1.7114975189653e19*cos(RAD(7.349020e+01)))
4.690980e+20/(1.7114975189653e19*cos(RAD(7.204080e+01)))
5.445430e+20/(1.7114975189653e19*cos(RAD(7.132620e+01)))
6.381640e+20/(1.7114975189653e19*cos(RAD(7.107400e+01)))
6.785450e+20/(1.7114975189653e19*cos(RAD(7.086810e+01)))
6.926840e+20/(1.7114975189653e19*cos(RAD(7.073370e+01)))
6.770860e+20/(1.7114975189653e19*cos(RAD(7.054340e+01)))
6.459870e+20/(1.7114975189653e19*cos(RAD(7.044120e+01)))
5.997140e+20/(1.7114975189653e19*cos(RAD(7.031490e+01)))
5.454890e+20/(1.7114975189653e19*cos(RAD(7.027940e+01)))
4.962500e+20/(1.7114975189653e19*cos(RAD(7.033710e+01)))
4.592390e+20/(1.7114975189653e19*cos(RAD(7.035760e+01)))
4.217650e+20/(1.7114975189653e19*cos(RAD(7.025710e+01)))
4.217650e+20/(1.7114975189653e19*cos(RAD(7.005100e+01)))
4.217650e+20/(1.7114975189653e19*cos(RAD(6.993850e+01)))
4.217650e+20/(1.7114975189653e19*cos(RAD(6.990560e+01)))
3.233590e+20/(1.7114975189653e19*cos(RAD(6.992830e+01)))
3.045750e+20/(1.7114975189653e19*cos(RAD(7.004700e+01)))
2.760480e+20/(1.7114975189653e19*cos(RAD(7.015080e+01)))
2.448930e+20/(1.7114975189653e19*cos(RAD(7.015080e+01)))
2.171780e+20/(1.7114975189653e19*cos(RAD(7.042400e+01)))
1.903990e+20/(1.7114975189653e19*cos(RAD(7.064860e+01)))
1.610130e+20/(1.7114975189653e19*cos(RAD(7.103810e+01)))
1.344120e+20/(1.7114975189653e19*cos(RAD(7.149680e+01)))
1.148650e+20/(1.7114975189653e19*cos(RAD(7.163160e+01)))
1.044350e+20/(1.7114975189653e19*cos(RAD(7.202490e+01)))
9.934400e+19/(1.7114975189653e19*cos(RAD(7.224920e+01)))
9.245870e+19/(1.7114975189653e19*cos(RAD(7.247960e+01)))
8.053860e+19/(1.7114975189653e19*cos(RAD(7.268960e+01)))
6.675810e+19/(1.7114975189653e19*cos(RAD(7.286270e+01)))
5.385430e+19/(1.7114975189653e19*cos(RAD(7.295460e+01)))
4.159810e+19/(1.7114975189653e19*cos(RAD(7.305270e+01)))
3.017970e+19/(1.7114975189653e19*cos(RAD(7.327160e+01)))
2.336700e+19/(1.7114975189653e19*cos(RAD(7.393910e+01)))
"
! NOTE LAST QUOTE
Updates: Jul 14, 1994: MV, Document created.