Program:       MOMENTS2D

Purpose:       Calculate FWHM and position angle of an ellipse, that represents
               a 2-dim gaussian, using zero-th, first and second moments
               of a distribution. The program can find asymmetrical features 
               in a data set.

Category:      ANALYSIS, CALCULATION

File:          moments2d.c

Author:        M.G.R. Vogelaar

Keywords:

   INSET=      Give set, subset:
               
               Maximum number of subsets is 1.


   BOX=        Give box in .....                        [entire subset]


   CLIPLO=     Enter one or more lower clip values:
               
               To select a range of data which you want to include in
               the moments analysis, use keywords CLIPLO= and CLIPHI=
               The entries in both keywords correspond so e.g. 
               CLIPLO=10 15 CLIPHI=20 40 means that the first range 
               includes data with levels between 10 and 20 (in units of
               the input data) and the second range includes data with 
               levels between 15 and 40.
               
              
   CLIPHI=     Enter exactly ... values for the upper clip:
               
               See description at CLIPLO=
               You have to specify exactly the same amount of levels
               as in CLIPLO=
               
              
   ZEROLEVEL=  Enter a value for the zero level of the data:      [calc]
               
               For the Full Width Half Maximum values that are 
               calculated using moments analysis, you need to define
               the intensity offset of your data. The default is that
               the program estimates this level by averaging the piixels
               on the border of the box in BOX=. If it cannot find any
               non-blank pixels on the border, the zero level is set 
               to 0.0.
               
              
   MARKRANGE=  Enter the clip range to mark in the plot:          [none]
               
               If you enter two values here then all data displayed
               in gids with a value within this range, will be marked 
               so that one can see which data is include in a moments
               analysis. If you enter just one value, then all pixels 
               with values greater than this one will be marked.
               

   FILENAME=   Enter a name or press return for log file only: [no file]
               
               If a name is specified, an ASCII file is created to
               store data. If you press carriage return, there will 
               be no output to an ASCII file only the GIPSY screen and
               log file.





Description:   Using moments analysis, one can return an estimated center,
               the Full Width Half Maximum in x and y and the position 
               angle of 2dim gaussian represented by data in a box that 
               is set by BOX=
               Only data is used with levels in ranges defined by CLIPLO=
               and CLIPHI=
                                                                          
               The analysis method uses the zero-th, first and second
               moment of the data to estimate the parameters of the       
               gaussian. The function can be represented by:
                                                                          
               z(x,y) = z0.exp[-{a(x-xc)^2+b(y-yc)^2+2c(x-xc)(y-yc)}]+o 
                  
               where 'o' is an offset in the data entered in ZEROLEVEL=
               xc and yc are the estimated centers and a,b,c are parameters
               that are converted to FWHM and position angle.
                                                                       
               Definition of the moments:
                                                                          
               First moment M0: M0 = SjSi Tji   (S == Sigma character)    
               For the peak position in grid coordinates we have:          
               x_av = 1/M0.SjSi[xi.Tji]   and  y_av = 1/M0.SjSi[yi.Tji]   
                                                                          
               Second moments Mxx = 1/M0.SjSi[xi^2.Mji]  - x_av^2         
                              Myy = 1/M0.SjSi[yi^2.Mji]  - y_av^2         
                              Mxy = 1/M0.SjSi[xi.yi.Mji] - y_av.x_av      

               Note that calculations are done both in grids and arcsec. 
               the grid parameters are used to plot an ellipse. The
               other parameters are only written to log file or file on disk. 

               The conversions to FWHM and angle are given by the formulas:
               
                           theta = 0.5 * atan(2.0*c/(a-b)); 
               
               A dummy variable 'p' is defined as:
               
                           p = sqrt( (a-b)*(a-b)+4.0*c*c );
                           
               then:
               
                    FWHM in x = 2.0 * sqrt( 2.ln(2)/(a+b+p) );
                    FWHM in y = 2.0 * sqrt( 2.ln(2)/(a+b-p) );
                    
               or v.v. if b > a.

               The moments analysis can be used for instance to find 
               asymmetric distributions in galaxies. The method is 
               described in more detail in the documentation of GAUFIT2D.
               
               
               Overlays:
               ========
               
               Ellipses are plotted in GIDS. If set. subset loaded in GIDS
               are the same as set and subset for this program, then 
               ellipses are plotted over the map. Else, program VIEW
               is started to display the right set and subset.


Notes:         

Example:       

Updates:       Oct 09, 2000: VOG, Document created.