At this point is counterproductive to carry around notation originally found in the Zhong Wang model, which was a wholly different functional form. However, roughly speaking, as in previous reports the new equation to be fit is of the form:

DN(obs) = a t + b t^2 + c t^3

In principle, one could, in the manner of the quadratic function, reduce this to two unknowns and DN(lin), but it doesn't help very much since this just makes it harder to apply the standard root inversion.

I have already written an algorithm that does the fitting. It works in essence much like lincal currently does - iterative sigma rejection and identification of the saturation point based on signal rate. The three parameters a, b, and c look like

Non-linearity Parameters for Ch. 3 CTA non-linearity special test sequence.

Saturation and Fitting Error Map for Ch. 3 CTA non-linearity special test sequence.

The saturation points and errors are also given above. Saturation is typically 48-50,000 DN, and the avg. error of the fit vs. the data is 0.7%.

Root Inversion and Accuracy of Method

There are several ways to invert the cubic equation. The one enumerated by Rick is found in section 5.6 of Numerical Recipes. I examined how well this works by comparing it to an iterative root finder I write myself to solve this specific application. The following piece of IGOR code was used to implement Rick's approach:

| Find roots of function:
| DN(obs) = s = a t + b t^2 + c t^3

function cuberoot(a,b,c,s)
variable a,b,c,s


| recast into Rick's form:
| y = t^3 + a2*t^2 + a1*t + a0

variable xn,yn,delta,h,a2,a1,a0,theta,t

a0=-s/c
a1=a/c
a2=b/c

xn = -a2/3
yn = ((xn+a2)*xn+a1)*xn+a0
delta = sqrt(a2^2-3*a1)/3
h = 2*delta^3
theta = acos(-yn/h)/3
t = xn + 2*delta*cos(4*pi/3+theta)

return(t)

Performing the above inversion 65,536 times on a 300MHz G3 took less than 1 second. In order to check the above root for accuracy, I did the inversion for all possible DN ranges (0.1 - 65536) for a "typical" pixel (pixel 125,125 of the Ch. 4 array). The roots so found were then reapplied to the cubic equation to see if the non-linear DN out were the same as that that went in. The results were almost universally exact (using 32-bit floats). The greatest deviation seen was 1e-5%, which could probably be attributed to small inaccuracies in the trig functions. Hence we can conclude that the above root-finder is more than sufficiently accurate for our usage. The above root also always corresponds to the same root that my own code settled on as the appropriate one.

Alternately, one could use something like Newton-Raphson to find the roots, although I think the above is just as easy.

Error Propogation

Since we have an analytic form for the root, propagating the errors is relatively straightforward, although not particularly pleasant to code. We use the standard error propagation formula, wherein if we have a function

y=f(x1,x2,x3...)

then

sigma(y) = Sqrt( (D[f,x1] sigma(x1))^2 + (D[f,x2] sigma(x2))^2 + ...)

where D[f,x] is the partial derivative of f with respect to x. Analytically, the function f is

                                                               2
                                                           -2 b    a
                                                      -(b (----- + -))
                                                              2    c
                                                           9 c           s
                                                 -27 (---------------- - -)
                                                            3 c          c
                                          ArcCos[--------------------------]
                                                          2
                                                         b    3 a 3/2
                                                      2 (-- - ---)
                      2                                   2    c
                     b    3 a      4 Pi                  c
              2 Sqrt[-- - ---] Cos[---- + ----------------------------------]
                      2    c        3                     3
        -b           c
f =     --- + ---------------------------------------------------------------
        3 c                                  3