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.