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.