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 2blank) RANGE=5, 2 (IF 2blank 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.