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