Program: ELLINT
Purpose: Integrate map in elliptical rings. Write statistics
to log file, Ascii file or GDS table.
Category: ANALYSIS, TABLES
File: ellint.c
Author: M.G.R. Vogelaar
Keywords:
INSET= Give set, subsets:
Maximum number of subsets is 2048.
Dimension of the subsets must be 2.
OPTION= Give OPTION: [1]
OPTION=1 Calculate SUM, MEAN, MEDIAN, AREA etc. in
rings.
OPTION=2 Calculate projected and face-on surface
brightness in mapunits/area for all rings.
OPTION=3 Calculate mass surface density in Msun/pc^2
given the total mass (MASS=) in Msun and
distance (DISTANCE=) in Mpc. The total area
under the density distribution will be
scaled to the total mass.
MASS= Give mass in Msun:
Total mass of your object in solar masses. This value
will be used (together with DISTANCE=) to scale the
results to mass surface densities.
DISTANCE= Give distance in Mpc:
See description at MASS=
** PIXELS= Calculations and plotting in pixels? [N]
Do calculations and plotting in pixels, even when the
conversion to seconds of arc is possible. If PIXELS=Y
subsequent keywords (RADII=, WIDTH= etc) are given in
pixels. There is a status message in the log file about
this conversion.
RADII= Give radii in seconds of arc:
If PIXELS=Y or no conversion to arcsec could be made
then the prompt is changed to:
Give radii in pixels:
The radius is measured in seconds of arc (or pixels)
from the centre (POS=) of the inner ellipse to the
centre of the ring, along the major axis. Maximum
number of radii is 512.
WIDTH= Width of ring (arcsec or pixels):
There can be as many widths as there are radii.
If the number of widths is less, the last value is
copied until the number of widths is equal to the
number of radii. A ring (radius R, Width W) is defined
by the inner ellipse with radius R-(W/2) and an outer
ellipse with radius R+(W/2)
PA= Position angle major axis (deg):
Measured as per astronomical convention from north
to east. For the number of position angles see remarks
at WIDTH=
SEGMENTS= Segments (angles wrt major axis): [0,360]
For options 1 and 2 it is also possible to calculate
statistics in a part of a ring, a so called segment.
Default segment is the entire ring (0,360). Required
to define a valid segment are two angles. The maximum
number of segments is 360. The first pair of segment
angles MUST be 0,360.
See description.
INCL= Inclination (deg):
If i is the inclination in degrees then
minor / major = cos(i).
For the number of inclinations see remarks at WIDTH=
POS= Position of center of ellipses: [0,0]
(Enter position in grids or physical coordinates)
This position applies to ALL ellipses.
RANGE= Range of intensities to include: [all levels]
Input consists of two values. The first value may
be greater than the second. See the description
for the correct use of this keyword.
** OVERLAP= Weight data in overlapping rings: [Y]/N
Data can be weighted in overlapping rings. If at certain
position a pixel with image value X is encountered in
N rings, its image value decreases to X/N if OVERLAP=Y
SUBPIX= Give number of 'sub'pixels in x and y per pixel: [2 2]
Make a better approximation of the true area in a ring
by dividing a pixel into SUBPIX= 'sub'pixels. The
first number is the subdivision in x. The second is the
subdivision in y. If only x is entered, then y is copied
from the value of x.
** ORDER= Table order 1)seg/rad, 2)rad/seg: [1]/2
The generated statistics table can display all
segments for one ring (ORDER=1) or display all radii
for one segment (ORDER=2).
MEDIAN= Do you want to calculate the MEDIAN also! Y/[N]
For OPTION=1 it is possible to calculate the median
of the (sub)pixels in a ring or segment.
OVERLAY= Overlay ellipse(s) in GIDS? [N]/Y
If GIDS displays a set with the same name as INSET=
then you are prompted with this keyword. It is asked
hidden if the names do not match and the keyword is
skipped if GIDS is not in use.
ELLIPSE= Give index of ellipse to plot: [no plot]
Keyword is unhidden if OVERLAY=Y else hidden.
Syntax: ELLIPSE=n m
Indices correspond to the n-th radius and m-th segment.
If m is omitted, m is set to 1 (first segment).
The first number corresponds to a radius index,
the second to a segment index. Indices start at 1.
Plot the selected ellipse and fill the LAST selected
ring & segment with dots on positions where a pixel is
found that is include in the ring. The default ends the
plotting loop.
This plot illustrates the definition of the angles
and the sampling of a ring or segment.
** FILENAME= Write statistics to: [No file]
Results are always written to a GDS table, but
is also possible to write data to an Ascii file. The
name of this file is given in FILENAME=
If the file already exists, the keyword OVERWRITE=
is prompted. The output is formatted according to
FORMAT=
OVERWRITE= File exists, ok to overwrite? [Y]/N
Only asked if FILENAME= is a name of an existing
file on disk.
FORMAT= Enter format for number output: [ffffff.ff]
See description for syntax and examples.
Keywords to initialize plots:
** PGBOX= Corners of box Xl,Yl, Xh,Yh: [calculated]
(In the same units as the radius)
If a plot is made of one or more ellipses, then
the program needs to know what the plot boundaries are.
Default are the corners the same as the box (in pixels)
returned by GIDS, or if GIDS is not used, the default
corners are set by the largest ellipse.
PGBOX= is also prompted (hidden) after keyword PLOTOPT=
to set the limits of the box in which profiles are
plotted.
** PGMOSAIC= View surface subdivisions x,y: [1,1]
It is possible to have more plots on the view surface
(f.i. the plot window).
The number of plots can vary in both directions x and y.
GRDEVICE= Give graphics device to [list of available devices]
plot on:
** PGPAPER= Give width (cm), aspect ratio: [0.0, 1.0]
Change the size of the (output) view surface.
The aspect ratio is defined as height/width.
The default is a calculated size and aspect ratio 1.
** PGWIDTH= Give line width (1-21): [2]
Only required if a hardcopy device is selected.
Note that the plot keywords are asked before plotting
ellipses and are asked again before plotting profiles.
Keywords related to GDS tables:
** GDSTABLE= Store results in table header of INSET= Y/[N]
If GDSTABLE=Y a GDS table is created in the header
of INSET= and most of the results are written to this
table (so that it can extracted by porgram TABLE or
by your private COLA program). However, writing to such
a table can slow down the output. Therefore the
default will NOT create a GDS table.
** TABNAME= Give name of table to store results: [ELLINT]
Columns are created at set level. See 'description'
to see what the names and contents of the columns
is. Tables can be listed/edited/plotted with program
TABLE.
** TABAPPEND= Append to existing columns? Y/[N]
If a table already exists, it is possible to append
to this table if you enter TABAPPEND=Y
The default always creates a new table.
PLOTOPT= OPTION=1 MEDIAN=Y:
1=sum 2=mean 3=rms 4=min 5=max 6=A 7=A-bl 8=med: [quit]
OPTION=1 MEDIAN=N:
1=sum 2=mean 3=rms 4=min 5=max 6=A 7=A-bl: [quit]
OPTION=2:
1=sd 2=sd-t 3=fo-sd 4=fo-sd-t 5=sum 6=A 7=A-bl: [quit]
OPTION=3
1=Msd-t 2=M-t 3=Cu.M 4=Msd 5=M 6=Cu.M 7=S 8=A 9=A-bl:[q]
Plot calculated statistics as function of radius.
At this moment one can alter the default limits
of a plot with PGBOX= (is a hidden keyword).
Example: PLOTOPT=1 PGBOX=0 0 1200 10
Description: GEOMETRY
========
ELLINT can be used to find the radial intensity
distributions in galaxies or, for example, to
find the mean intensity of instrumental rings
in maps. The source of the data is specified by INSET=
The output is a table where the sum, mean, median and
rms, are given in map units (OPTION=1), or as (a sort
of) surface brightness (OPTION=2: map units/arcsec^2
usually for optical work) or, as mass surface density
(OPTION=3: Mo/pc^2). The last option is a scaling and
requires the total mass of your object (galaxy) in
solar masses (MASS=) and the distance in Mpc (DISTANCE=)
(a negative sum will result in a negative mass).
The maximum number of subsets is 2048 and the subsets
must be two dimensional. The program gets the axis
information from the header (can be displayed when
Hermes set to TEST mode) and calculates a factor to
convert from the header axis units to seconds of arc.
However, if you want that all calculations and plotting
is done in pixels, use PIXELS=Y. The program informs you
if it uses the conversion to seconds of arc or not.
At the RADII= prompt, give radii measured in seconds of
arc (or pixels) from the center along the major
axis.
The keyword POS= is used to set the position
of the center of all the ellipses.
For example, if you want the new central position in
input pixel coordinates (-10.680130, -8.700096) and
position (0,0) in your input corresponds to the physical
coordinates 198.468100, 46.002610 deg. you specify:
POS=-10.680130 -8.700096 (pixels) or:
POS= U 198.493667 U 45.983306 (deg, deg) or:
POS= * 13 13 58.48 * 45 58 59.90 (hms, dms)
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. This angle is given in PA=.
If possible, the program automatically corrects for the
rotation angle of the actual axis from the stated header
coordinate. It is assumed that the spatial latitude axis
carries the rotation information. If nothing is found
in the header, 0 degrees is assumed. If you give less
angles than there are radii, the missing angles are copied
from the last one.
A ring is defined after the keyword WIDTH= by an
inner ellipse with radius R-(width/2) and an outer
ellipse with radius R+(width/2). The width is given in
the same units as the radii. For each radius you can
define a width. If you give less widths than there are,
radii, the program copies the last value of
the width until there are as much widths as there are
radii. Note that the width of the defined segments is
not constant in the defined ellipse but is only
constant in its (de)projection to a circle.
Integration in circles is simply accomplished by setting
INCL=0.0 else you integrate in an ellipse with
minor/major = COS(INCL=).
If you want to integrate a complete circle with radius
R, then use:
INCL=0.0 , RADII=R/2 , WIDTH=R
For each radius an inclination can be given, but if
you give less inclinations than there are radii, the
rest of the inclinations are copied from the last one
until there are as much inclinations as there are radii.
You can give an expression instead of the angles.
If, for example, you want to give two ratio's instead of
two angles, try an input like:
INCL=DEG(ACOS([0.3 0.5]))
Between the brackets are the arguments 0.3 and 0.5.
They are converted to the angles 72.5424 and 60.0 deg.
by the input routines.
For all defined rings, the calculations are carried
out for one or more segments (SEGMENTS=).
Two numbers mark one segment. The numbers
[0,360 270,90 90,270] for example, define three
segments. The first one is the complete ring, the second
runs from the minor axis to the opposite minor axis
through the major axis (from -90 deg. incl. to 90 deg.
excl. in the direction of the east) and the third segment
runs from the opposite minor axis to the the minor axis
through the opposite major axis (from 90 deg. incl. to
-90 deg. excl.)
The angles can be negative but then they are converted
to values in the range [0,360> degrees.
If option 3 was selected the first segment must be
0,360, otherwise the scaling would not be correct.
You can visualize the way the program interprets
the axes and the angles with keyword ELLIPSE=
The arguments are the index of a radius and a segment.
ELLIPSE=2,3 plots the third segment of your second
ellipse. The default skips plotting.
The plot-data is always extracted from the first subset.
You can overlay such a plot on an image displayed in GIDS.
The image can be a zoomed image. If it is possible
to apply an overlay then you are prompted with
OVERLAY=
STATISTICS
==========
ELLINT computes in a given elliptical ring with
certain major-, minor axis, width, angle, inclination
and in each defined segment:
OPTION=1 : the sum, mean and rms of these values.
If the values are represented by Xi and there
are Ni valid values, then:
sum = SUM(Xi)
mean = sum / Ni
var = (SUM(Xi^2) - Ni*mean^2) / (Ni-1)
rms = SQRT(var)
A valid value is not a blank and is within
a range set by RANGE= Values outside this range
are treated as blanks. A 'subpixel' (there are
SUBPIX= 'subpixels' in a pixel in x and y) is
treated like a pixel.
Examples of the use of the RANGE= keyword:
RANGE=2, 5
(IF 2blank)
RANGE=5, 2
(IF 2blank ELSE include)
At the RANGE= keyword, the values -INF and INF can
be input also. These values represent the minimum
and maximum values of the current system.
RANGE=5, INF
(IF value>5 THEN include ELSE value==>blank)
If regions do overlap in rings, it is possible
to integrate them separate or weight the common data
over the overlapping regions. The keyword OVERLAP=
sets the action. If OVERLAP=Y a value is weighted by
the number of times the corresponding position is
present in a ring.
OPTION=2: -Mean surface brightness projected on the sky,
averaged over all non blank area.
-Mean surface brightness projected on the sky,
averaged over total area.
-Mean face-on surface brightness averaged over
all non blank area.
-Mean face-on surface brightness averaged over
total area.
The mean surface brightnesses are calculated by
summing and averaging the surface densities per
pixel. This is equivalent to the sum divided by the
area. Here we distinguish two areas. The first
area is equal to the area of all pixels that are
not blank.
The second area is the area that is equal to
to the area of all pixels.
It depends on what you think of what a blank
actually represents whether you use the one
surface brightness or the other.
The face-on areas are increased by a factor
1/COS(inclination)
sb = Xi / [Area of one pixel / COS(inclination)]
The units are map units per second of arc
squared.
OPTION=3: First for each ring the mean face-on surface
density (fo-sd) (derived in the same way as the
mean face-on surface brightness of OPTION=2)
is multiplied by PI*(R(n+1)^2-Rn^2), the area between
two annuli, with index n and n+1, in the segment
0,360 deg. (this is the reason that for OPTION=3
the first segment in a sequence must always be
0,360 degrees).
The values are summed for all rings and the result
is stored ('total sum').
If M is the total mass of the object in solar
masses (MASS+), D its distance (DISTANCE=) in Mpc
and the conversion of an area in seconds of arc to
an area in pc is given by:
x(pc) = 4.848 * D(Mpc) * a(arcsec)
then the mass surface density sigM per ring is
sigM = (face-on sd/total sum)*M/(4.848*D)^2
We distinguish two values for the mass surface
density. This is because the mass surface density
is derived from the surface density (see OPTION=2)
and that value depends on the area that was used
(total area or non-blank area). Program ELLINT
lists both values in a table.
TABLES
======
Examples of output table in log file:
OPTION=1:
radius cum.sum sum mean median rms
arcsec JY/BEAM JY/BEAM JY/BEAM JY/BEAM JY/BEAM
========================================================= ...
0.00 1.49 1.49 0.37 0.37 0.01
30.00 9.19 7.71 0.39 0.39 0.01
subpix unique area area-bl segLO segHI width
# pixels pixels pixels deg. deg. arcsec
====================================================== ...
4 4 4.00 0.00 0.00 360.00 30.00
20 20 20.00 0.00 0.00 360.00 30.00
pa incl
deg. deg.
==============
336.00 30.00
336.00 30.00
cum.sum: The cumulative sum, i.e. the 'running' sum of the
values in the next column 'sum'.
sum: Add the image values represented by each 'subpixel'
in the ring.
The image value of a 'subpixel' is equal to the
image value of the pixel to which it belongs
divided by the number of 'subpixels' in a pixel.
mean: Mean image value in the ring.
median: Only listed if MEDIAN=Y. It the number of
'subpixels' is equal to N, then the median element
of the data is the (N+1)/2 th element of the
sorted data if N is odd. If N is even, its the
arithmetic mean of the N/2 th and N/2+1 th element.
rms: Measure of data's width around its central value.
subpix: Number of 'subpixels' that was found in a ring.
This number depends on the value of SUBPIX=
It will increase if SUBPIX= is increased.
unique: Number of different pixels that contributed to the
ring. If SUBPIX= is increased then this number
will converge to a fixed number. This is because
a larger number for SUBPIX= means a better
approximation of the area of a ring.
area: Area of ring covered by non blank image values, counted
in 'subpixels' and scaled to pixels.
area-bl: Same as area, but now only the blank 'subpixels'
are counted.
segLO: All the statistics described above can also be
applied to segments. If a ring is divided into
segments (SEGMENTS=) then the generated table displays
all statistics for all segments per ring (ORDER=1) or
it displays statistics for all radii per segment
(ORDER=2). This column lists the start angles of the
segments.
segHI: End angles of segments.
width: Width of current ring.
pa: Position angle of major axis of current ring.
incl: Inclination.
OPTION=2:
radius sb sb-t fo-sb fo-sb-t sum
arcsec JY/BEAM/'' JY/BEAM/'' JY/BEAM/'' JY/BEAM/'' JY/BEAM
=============================================================
0.00 0.0017 0.0017 0.0014 0.0014 1.4866
30.00 0.0017 0.0017 0.0015 0.0015 7.7064
sb: Mean surface brightness projected on the sky,
averaged over all non blank area.
sb-t: Mean surface brightness projected on the sky,
averaged over total area (i.e. blank pixel area
included).
fo-sb: Mean face-on surface brightness averaged over
all non blank area.
fo-sb-t: Mean face-on surface brightness averaged over
total area.
OPTION=3:
radius radius Msd-t Mass-t Cum.Mass Msd Mass Cum.Mass
arcsec Kpc Mo/pc^2 10^9 M0 10^9 M0 Mo/pc^2 10^9 M0 10^9 M0
===================================================================
0.00 0.00 5.05 0.00 0.00 4.67 0.00 0.00
30.00 0.99 5.24 0.03 0.04 4.84 0.03 0.03
radius (Kpc): Radius converted from arcsec to kpc:
R(Kpc) = 4.848 * D(Mpc)/1000 * Radius('');
Msd-t: Mass surface density derived from mean face-on
surface density (as in OPTION=2) averaged over
total area.
Mass-t: Scaled MASS= proportional to surface density
times area of a ring. The surface density
is the same as used in 'Msd-t'.
Cum.Mass: The cumulative 'Mass-t'. The value in the
last row of this column, is equal to MASS=
Msd: Mass surface density derived from mean face-on
surface density (as in OPTION=2) averaged over
area in which there are no blank pixels.
Mass: Scaled MASS= proportional to surface density
times area of a ring. The surface density
is the same as used in 'Msd'.
Cum.Mass: The cumulative 'Mass'. The value in the
last row of this column, is equal to MASS=
The result is written to disk in an Ascii file with
name set by FILENAME=
If this file already exists, the user is warned
with OVERWRITE=
Note that the first lines in the file are meant as
comments and are prefixed with an exclamation mark (!).
If you want more precision in the output, use the
FORMAT= keyword.
Results are always stored in a table. The name of the table
is 'ELLINT' by default. You can change this with TABNAME=
The column names are:
CUMSUM, Cumulative sum
SUM, Sum in ring
MEAN, Mean in ring
RMS, Rms in ring
MINVAL, Minimum of data
MAXVAL, Maximum of data
NPTS, Non blank pixels in ring
NBLANKS, Blanks in ring
SEGMLO, Start angle of segment
SEGMHI, End angle of segment
WIDTH, Width of ring
PAMAJ, Pos. angle Major axis
INCL, Inclination of ring
SUBGRIDn Subset grid of n-th non subset axis
FREQ etc. Physical coordinates of non subset axis.
LONGGRID Centre position in grids (longitude)
LATGRID Centre position in grids (latitude)
LONGPHYS Physical coord. of CP (longitude)
LATPHYS Physical coord. of CP (latitude)
Data for all subsets, radii and segments are stored in one
table. If you want to append to an existing table, use
TABAPPEND=Y The application TABLE is used to do the
formatting of numbers and the plotting of data.
Number output (FORMAT=):
FORMAT= is used to print floating point numbers
in an user specified format. Input is a format image
consisting of characters representing the
wanted output format for 'NUMBER'. The formatted number
is returned in 'RESULT'.
The syntax for 'FORMAT'is:
flag(s):
Zero or more flags, in any order, which modify the
meaning of the conversion specification. The flag
characters and their meanings are:
- The result of the conversion is left-
justified within the field.
+ The result of a signed conversion always
begins with a sign, "+" or "-".
string:
Characters, some with special meaning.
If the string (f.i. FFFFF.FF or gggg.gg or wwwww)
contains no dot, the number of characters specify
a minimum field width. For an output field, if the
converted value has fewer characters than the field
width, it is padded on the left (or right, if the
left-adjustment flag, - has been given) to the field
width.
If the string contains a dot, the total number of
characters including the dot is the minimum field width
and the number of characters after the dot is the
precision.
The characters are used to determine the conversion
type. If the string contains an:
'e' or 'E'
The floating-point-number argument is
printed in the style [-]drddde+dd,
where there is one digit before the
radix character, and the number of
digits after it is equal to the
precision. The E conversion character
produces a number with E introducing
the exponent instead of e. The
exponent always contains at least two
digits. However, if the value to be
printed requires an exponent greater
than two digits, additional exponent
digits are printed as necessary.
'g' or 'G'
The floating-point-number argument is
printed in style f or e (or int style E
n the case of a G conversion
character), with the precision
specifying the number of significant
digits. The style used depends on the
value converted; style e is used only
if the exponent resulting from the
conversion is less than -4 or greater
than or equal to the precision.
others
Strings without 'e', 'E', 'g' and 'G'
indicate a floating point conversion.
The floating point number argument is
printed in decimal notation in the
style [-]dddrddd, where the number of
digits after the radix character, r, is
equal to the precision specification.
If the result of a conversion is longer than the
field width, an asterisk is returned. If the
input number is a blank, a 'b' is returned.
Examples: FORMAT= +eeeeee.eeee
Number: 43345.5436
Result: +4.3346e+04
Remark: exponential format
signed conversion
field width: 12
precision: 4
FORMAT= gggg.ggggg
Number: 34.43
Result: 34.430
Remark: Field width is 10
Number of significant digits is 5
FORMAT= +ffff.ff
Number: 23.456
Result: +23.46
Remark: signed conversion
FORMAT= -ffff
Number: 345.87
Result: 346
Remark: left justified
FORMAT= -+ffff.fff
Number: 34.43
Result: +34.430
Remark: left justified
signed conversion
FORMAT= eee.eeee
Number: 234.43
Result: *
Remark: Field width too small
for conversion
FORMAT= ffff.ff
Number: blank
Result: b
Remark: input was a blank
Example:
ellint
ELLINT Version 1.0 (Oct 31 1991)
ellint ellipse=3 1
ELLINT Version 1.0 (Oct 31 1991)
INSET="../M8320" f 30
Set ../M8320 has 3 axes
RA-NCP from -50 to 40
DEC-NCP from -50 to 20
FREQ-OHEL from 1 to 59
Choose on of the following options:
OPTION=1 Calculate SUM, MEAN, #POINTS along ellipse,
in map units
OPTION=2 Calculate surface brightness in map units
OPTION=3 Calculate surface brightness in Msun/pc^2
OPTION=1
Calculations in seconds of arc
RADII=100:200:50
WIDTH=50
PA=60
INCL=deg(acos(0.5))
SEGMENTS=
POS=* 13 14 00 * 46 0 0
GRDEVICE=x11
Calculations in seconds of arc
FILENAME=
====== ELLINT RESULTS from set ../M8320 (31-OCT-1991) =====
etc., etc.
Notes:
Updates: Oct 14, 1991: VOG, Document created.
Feb 17, 1994: JPT, Header lines of file output preceded by "!".
Mar 8, 1995: VOG, Removed bug where pixels at 360 deg are not
counted in segment 0-360 deg.
Jan 24, 2002: VOG, Increased number of segments to 360