Channels 3 and 4 of IRAC are made with a different detector material than Channels 1 and 2. Not unsurprisingly, they have different linearity properties, and in general are much more non-linear than the InSb arrays. As a result, the quadratic in time solution used by LINCAL currently breaks downs (fails to fit) at unacceptably low DN. This has been well-documented in previous reports.
2. The Solution
After examining several alternatives, we have decided that the best fit that can be easily acheived is by extending the current approach to one additional order. That is, a fit should be made of the form
S = A t + B t^2 + C t^3
where S is the observed DN, and t is the exposure time. This is described as method 4 in Nick Collins' document and II.B.1 in Jason's document.
3. Derivation of Roots for FOWLINEARIZE
In order to apply this, the cubic equation must be inverted to solve for t, and the linearized counts S' are given by A t. There are three possible roots, but Rick Arendt has shown that only one is the appropriate solution:
We want a root (solution for R) to the equation
y(R) = D' = BR + aAR^2 + bCR^3
The second derivative of the equation is
y" = 2aA + 6 bCR
Normal detector behavior is that
The only cubic function that fits all these contraints is one with bC < 0, and a local maximum at Rmax > 0. (Additionally, the inflection point defined by y" = 0 will have to be near R=0. This rejects the y" < 0 portions of cubic functions with bC > 0.) If the fitted detector linearity parameters a and b don't result in such a curve, then we are dealing with an odd pixel, and this model is not appropriate for correcting its data. Thus, as long as D' < y(Rmax), there will be 2 positive roots to the equation (and one negative root). If D' > y(Rmax), then are dealing with a saturated datum that can't or shouldn't be corrected. The only root that makes physical sense is the middle root which lies at 0 < R < Rmax. The other positive root at R > Rmax would correspond to a value that is low again because it has dropped off after saturation. This root can be given by:
function cubic, a ; ; y = x^3 + a(2)*x^2 + a(1)*x + a(0) ; ; Finds the middle root for y in the case that ; three real roots exist. ; The solution for three real roots follows the form ; and notation of Nickalls, R. W. D. 1993, Math. Gazette, ; 77, 354. (http://www.m-a.org.uk/eb/mg/) ; WARNING: This code will fail if y does not have three real roots! xn = -a(2)/3.d0 ;yn = a(0)+a(1)*xn+a(2)*xn^2+xn^3 yn = ((xn+a(2))*xn+a(1))*xn+a(0) delta = sqrt(a(2)^2-3.d0*a(1))/3. h = 2.d0*delta^3 theta = acos(-yn/h)/3. z = xn + 2d0*delta*cos(4*!dpi/3.+theta) return,z end
Rick implemented the above IDL root solver and found that it can perform all of the inversions for an entire image in approximately 1 second on an 800MHz PC. I suspect it will be much, much faster when executed in C, and hence it should not represent a significant bottle neck in processing time. The full code (which gives all 3 roots) is:
function cubic, a ; ; 0 = x^3 + a(2)*x^2 + a(1)*x + a(0) ; ; Finds the 3 solutions for x ; The solution for three real roots (within the ; first IF statment) follows the form ; and notation of Nickalls, R. W. D. 1993, Math. Gazette, ; 77, 354. (http://www.m-a.org.uk/eb/mg/) ; The solution for other cases foloows the form and ; notation of Abromowitz, M. & Stegun, I. A. "Handbook ; of Mathematical Functions", Section 3.8.2, p. 17 q = dcomplex(a(1)/3.d0 - a(2)^2/9.d0,0) r = dcomplex((a(1)*a(2)-3d0*a(0))/6.d0 - a(2)^3/27.d0,0) if (float(q^3+r^2) lt 0) then begin ; case of three real roots xn = -a(2)/3.d0 yn = a(0)+a(1)*xn+a(2)*xn^2+xn^3 delta = sqrt(a(2)^2-3.d0*a(1))/3. h = 2.d0*delta^3 theta = acos(-yn/h)/3. z = xn + 2d0*delta*cos(indgen(3)*2*!dpi/3.+theta) return,z endif ; s1 and s2 are a little complicated to get the correct answer when ; they are cube roots of negative real numbers. s1 = r+sqrt(q^3+r^2) s2 = r-sqrt(q^3+r^2) if (abs(s1) ne 0) then begin if ((imaginary(s1) eq 0) and (float(s1) lt 0)) then $ s1 = -1*(-s1)^(1./3) else s1 = s1^(1./3) endif if (abs(s2) ne 0) then begin if ((imaginary(s2) eq 0) and (float(s2) lt 0)) then $ s2 = -1*(-s2)^(1./3) else s2 = s2^(1./3) endif z = dcomplexarr(3) z(0) = s1 + s2 - a(2)/3 z(1) = -1.*(s1 + s2)/2. - a(2)/3 +complex(0,sqrt(3))/2*(s1-s2) z(2) = -1.*(s1 + s2)/2. - a(2)/3 -complex(0,sqrt(3))/2*(s1-s2) return,z end
4. Fowler Number Correction
Finally, there is the impact on the fowler number correction to the linearization. Fortunately, because S can be parameterized as a function of time, the same formalism used for the quadratic function can be used here. This is is given by this document.