Program:       FLAT

Purpose:       Apply a flat field correction to a map by fitting a
               2-D polynomial to the background.

Category:      MODELS, CALCULATION, FITTING, CALCULATION

File:          flat.c

Author:        M. Vogelaar

Keywords:
               
   INSET=      Give set, subsets:
               Maximum number of subsets is 2048.   
               Dimension of the subsets must be 2.               

   BOX=        Area which needs to be flat fielded           [whole map]
   
   OUTSET=     Output Set/Subset(s)  

   ORDER=      Order of polynomial to be fitted to the map           [1]
               highest allowed order is 21


** RANGE=      Give range of levels to include in fit:      [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.



Description:   This program can be used for example to extract a 
               background by fitting a polynomial to it.
               The program uses a set (INSET=) which usually
               is blotted first. The subsets must be two dimensional.
               The area on which FLAT works is given in BOX=
               It writes the results to OUTSET= in the same box.
               Beware if OUTSET= already exists. It only overwrites
               data in the given box. A polynomial plane will be fitted 
               to the data using the least squares method. The order
               of this polynomial is given in ORDER= This number cannot
               exceed 21.
               Values outside the range defined in RANGE= are treated 
               as blanks, i.e. they are not included in the fit and are 
               interpolated after a plane is fitted to the data.
               Examples for the use of the RANGE= keyword:

               RANGE=2, 5
               (IF 2<value<5 THEN include in fit ELSE value==>blank)

               RANGE=5, 2
               (IF 2<value<5 THEN value==>blank ELSE include in fit)
                             
               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)
               
               The fitted background can be subtracted from the 
               original set with the program COMBIN.

               Polynomial:
               In general 2-dim polynomial of order N: 
               f(x,y) = SIG SIG Amn x^m y^n where n+m <= N
               
Example:       Special case, order = 2:
               Replace Amn by ai, where i is the column position of a in 
               the coefficient matrix. Then the polynomial can be written as:
               f(x,y) = z = a0 + a1.x + a2.y + a3.x^2 + a4.y^2 + a5.x.y
               To minimize Chi-squared we calculate the derivatives to the 
               parameters ai of the expression SIG[ (zfit-f(x,y))^2 ] which
               for this special case results in the coefficient matrix:
               
               SIG(z)     = a0.SIG(1) + a1.SIG(x) + a2.SIG(y) + a3.SIG(x^2) + a4.SIG(y^2) + a5.SIG(x.y)
               SIG(z.x)   = a0.SIG(1.x) + a1.SIG(x.x) + a2.SIG(y.x) + a3.SIG(x^2.x) + a4.SIG(y^2.x) + a5.SIG(x.y.x) 
               SIG(z.y)   = a0.SIG(1.y) + a1.SIG(x.y) + a2.SIG(y.y) + a3.SIG(x^2.y) + a4.SIG(y^2.y) + a5.SIG(x.y.y)
               SIG(z.x^2) = a0.SIG(1.x^2) + ....
               SIG(z.y^2) = a0.SIG(1.y^2) + ....
               SIG(z.x.y) = a0.SIG(1.x.y) + ....
               
               If we call the left side vector B and the coefficient matrix A
               then:
               B = A.(a0,a1,a2,a3,a4,a5)T
               
               The unknown vector (a0,a1,a2,a3,a4,a5) is calculated by 
               multiplying the inverse of A with B.
               
               The algorithm is splitted into two parts. First we fill a
               matrix for an arbitrary order N matrix. Secondly we use
               a Gauss-Jordan algorithm to calculate the inverse. In special
               cases one expects better results by applying the 
               Singular Value Decomposition method but until now we are 
               satisfied with the 'Gauss-Jordan' fit results for the data 
               we processed.
               
               The calculations are in double precision. We noticed that the
               conversion of floats in an image with high values resulted 
               in a small increase of minimum chi-squared, but only by a very 
               small amount, too small to be significant.
              
               For the data points listed below the results are:
               From Mathematica:
               2.78061e+12 - 6.73389e+9 x - 1.60955e+10 x^2 - 1.03662e+9 y 
               - 2.83247e+8 xy - 9.30063e9 y^2
               
               flat.c:
               2.78061e+12 - 6.73389e+09 X**1 Y**0 - 1.03662e+09 X**0 Y**1
               - 1.60955e+10 X**2 Y**0 - 9.30063e+09 X**0 Y**2 
               - 2.83247e+08 X**1 Y**1
               chi2=5.90453e+27
              
               Fit from article (without reference, just as an example):
               # Ax^2 + By^2 + Cxy + Dx + Ey + F
               #
               # A = -16.1e9
               # B = -9.3e9
               # C = -0.28e9
               # D = 6.11e5
               # E = 0.32e5
               # F = 0.47e2
               chi2 = 1.07231e+28
X =     5.0,
        397.0,
        1995.0,
        1994.0,
        10.0,
        1177.0,
        600.0,
        1420.0,
        1141.0,
        990.0,
        574.0,
        1740.0,
        1456.0,
        864.0,
        1604.0,
        1332.0,
        554.0,
        1643.0,
        442.0,
        1266.0
             
Y =     3996.0,
   	3736.0,
   	3994.0,
   	7.0,
   	59.0,
   	16.0,
   	1610.0,
   	2070.0,
   	1160.0,
   	3085.0,
   	2389.0,
   	1586.0,
   	3278.0,
   	3846.0,
   	3030.0,
   	3658.0,
   	3198.0,
   	900.0,
   	1937.0,
   	1661.0

Z =    	-1.485E+17,
   	-1.328E+17,
   	-2.147E+17,
   	-6.402E+16,
   	-3.415E+13,
   	-2.231E+16,
   	-3.017E+16,
   	-7.314E+16,
   	-3.385E+16,
   	-1.052E+17,
   	-5.877E+16,
   	-7.291E+16,
   	-1.354E+17,
   	-1.505E+17,
   	-1.282E+17,
   	-1.544E+17,
   	-1.006E+17,
   	-5.141E+16,
   	-3.828E+16,
   	-5.205E+16


Notes:         

Updates:       Oct 14,  1991: VOG, Document created.
               Mar 03,  2004: VOG, Added to documentation, changed calculations
                                   to double precision and tested the 
                                   algorithm for a test series of 20 data 
                                   points and a second order polynomial for 
                                   which fitted parameters where known from the
                                   literature. FLAT produced better results
                                   because it found a smaller value for 
                                   chi-squared compared to the value from
                                   the literature.