PRO noise_model_ccd, imdata, bpmdata, masterdark, masterflat, ron, gain, coeff2, coeff3, out_type, imout, bpmout, status
; Description: This function uses a standard CCD noise model to calculate the pixel variances
; for a set of calibrated image data "imdata", which should ideally represent the
; image model, taking into account the master dark image "masterdark", the master
; flat field image "masterflat", the CCD readout noise "ron" (ADU), the CCD gain
; "gain" (e-/ADU), and the non-linearity correction coefficients "coeff2" and
; "coeff3". The standard CCD noise model that is used is as follows:
;
; Var(M_ij) = (sig0^2)/(H^2 * F_ij^2) + D_ij/(G * F_ij^2) + M_ij/(G * F_ij)
; H = dg(Y_ij)/dY_ij
; Y_ij = D_ij + (F_ij * M_ij)
;
; where M_ij is the model for the calibrated image data "imdata", D_ij is the
; master dark image (scaled by exposure time) used to calibrate the image data,
; and F_ij is the normalised master flat field image used to calibrate the image
; data. The quantity sig0 represents the CCD readout noise (ADU), G represents
; the CCD gain (e-/ADU) and the function g is the non-linearity relationship for
; the CCD. Finally, Var(M_ij) represents the variance on the model for the
; calibrated image data.
; This noise model can be derived from considering the following model for the
; raw image data X_ij in ADU:
;
; X_ij = B_ij + g(D_ij + (F_ij * M_ij)) (Equation 1)
;
; where B_ij is the master bias image. Note that we are assuming that the
; non-linearity component distorts the detected counts on image readout, which
; is different to a non-linearity effect that changes the effective number of
; counts that are detected. The noise model for X_ij is then given by:
;
; Var(X_ij) = sig0^2 + (H^2 * (D_ij + (F_ij * M_ij)))/G (Equation 2)
;
; The model for the raw image data can be inverted to yield the model for the
; calibrated image data:
;
; M_ij = ((g^(-1)(X_ij - B_ij)) - D_ij)/F_ij (Equation 3)
;
; Noise propagation of Var(X_ij), defined in equation 2, using the relationship
; between X_ij and M_ij described by equation 3, yields the standard CCD noise
; model for Var(M_ij).
; The module returns the CCD noise model for the calibrated image data as an array
; "imout" consisting of either pixel uncertainties sqrt(Var(M_ij)), pixel variances
; Var(M_ij), or inverse pixel variances 1/Var(M_ij), depending on the value of the
; parameter "out_type". The bad pixel mask "bpmdata" is updated to include pixels
; for which the noise model could not be evaluated and this updated bad pixel mask
; is returned via the parameter "bpmout". Note that this module does not modify the
; input parameters in any way.
;
; Input Parameters:
;
; imdata - FLOAT/DOUBLE ARRAY - A two-dimensional image array representing a set of calibrated
; image data. Ideally this image array should represent the MODEL
; for the calibrated image data, since uncertainties correspond
; to the model, not the data. Negative or zero pixel values in
; this image will cause the offending pixels to be added to the
; output bad pixel mask "bpmout", since the Poisson noise is
; undefined or zero for these values.
; bpmdata - BYTE/INTEGER/LONG ARRAY - A bad pixel mask array of the same size as "imdata" which
; flags bad pixels with a value of "0" and good pixels with
; a value of "1". If this array does not have the same
; dimensions as "imdata", then all image pixels will be
; considered good.
; masterdark - FLOAT/DOUBLE ARRAY - The master dark image, scaled by exposure time, used to
; calibrate the image data. The pixel values in this image have
; units of ADU and represent the dark current specific to the
; image data "imdata". Negative pixel values will be ignored in
; the CCD noise model, and will not contribute to the pixel
; variance. If this array does not have the same dimensions as
; "imdata", then the contribution of the dark current to the
; pixel variances will be assumed to be zero.
; masterflat - FLOAT/DOUBLE ARRAY - The normalised master flat field image used to calibrate the
; image data. Negative or zero pixel values in this image will
; cause the offending pixels to be added to the output bad pixel
; mask "bpmout", since the master flat field image is a
; reciprocal component of the CCD noise model. If this array does
; not have the same dimensions as "imdata", then the flat field
; component of the CCD noise model will be ignored.
; ron - FLOAT/DOUBLE - The CCD readout noise (ADU). A zero or negative value of "ron" will lead to
; the readout noise component being ignored in the CCD noise model.
; gain - FLOAT/DOUBLE - The CCD gain (e-/ADU). A zero or negative value of "gain" does not make
; sense, and the gain will be assumed to be "1.0".
; coeff2 - FLOAT/DOUBLE - The quadratic coefficient in the CCD linearisation equation:
;
; g^(-1)(x) = x + (coeff2 * x^2) + (coeff3 * x^3)
;
; where x represents the image counts (ADU) after overscan and bias
; correction, and g^(-1)(x) represents the non-linearity corrected image
; counts (ADU). The non-linearity correction generates a correction to the
; image pixel variances, and any correction to the pixel variances that is
; zero or negative will cause the offending pixels to be added to the
; output bad pixel mask "bpmout", since the correction to the pixel
; variances is a reciprocal component of the CCD noise model.
; coeff3 - FLOAT/DOUBLE - The cubic coefficient in the CCD linearisation equation:
;
; g^(-1)(x) = x + (coeff2 * x^2) + (coeff3 * x^3)
;
; where x represents the image counts (ADU) after overscan and bias
; correction, and g^(-1)(x) represents the non-linearity corrected image
; counts (ADU). The non-linearity correction generates a correction to the
; image pixel variances, and any correction to the pixel variances that is
; zero or negative will cause the offending pixels to be added to the
; output bad pixel mask "bpmout", since the correction to the pixel
; variances is a reciprocal component of the CCD noise model.
; out_type - STRING - The array of pixel variances will be returned in one of four different forms.
; If "out_type" is set to the string "SIGMA", then the output array "imout" is
; returned as an array of pixel uncertainties. If "out_type" is set to the string
; "INV_SIGMA", then the output array "imout" is returned as an array of inverse
; pixel uncertainties. If "out_type" is set to the string "VAR", then the output
; array "imout" is returned as an array of pixel variances. If "out_type" is set
; to the string "INV_VAR", then the output array "imout" is returned as an array
; of inverse pixel variances. If "out_type" is set to any other value, then it is
; reset to the default value "VAR". In all cases, the output array "imout" is set
; to zero where a pixel is flagged as a bad pixel in the output bad pixel mask
; "bpmout".
;
; Output Parameters:
;
; imout - DOUBLE ARRAY - The output array of image pixel uncertainties, inverse pixel uncertainties,
; variances or inverse variances, as requested by the parameter "out_type"
; (see above). This array is of the same dimensions as "imdata", and bad
; pixels in this array flagged by "bpmout" (see below) will have a value of
; "0.0".
; bpmout - INTEGER ARRAY - The output bad pixel mask array of the same size as "imdata" which flags
; bad pixels with a value of "0" and good pixels with a value of "1". This
; bad pixel mask array includes the bad pixels flagged in "bpmdata" along
; with zero or negative pixels from the arrays "imdata" and "masterflat",
; which are also considered as bad pixels. It also includes those pixels
; for which the non-linearity correction generates a zero or negative
; correction to the image pixel variances.
; status - INTEGER - If the module successfully calculated the CCD noise model, then "status" is
; returned with a value of 1, otherwise it is returned with a value of 0.
;
; Author: Dan Bramich (dan.bramich@hotmail.co.uk)
;
; History:
;
; 08/05/2009 - Rewrite of the algorithm to invert the polynomial non-linearity correction to employ
; the newly created modules "solve_realroots_quadratic.pro" and
; "solve_realroots_cubic.pro". This section of the code now runs approximately ten times
; faster, and calculates the exact real roots (dmb).
; 24/04/2009 - Further module optimisation performed (dmb).
; 07/09/2008 - Module created (dmb).
;Check that "imdata" is an image of the correct number type
;If "bpmdata" is not an image of the correct number type and dimensions, then assume that all
;image pixels are good
;If "masterdark" is not an image of the correct number type and dimensions, then assume that
;there is no dark current
;If "masterflat" is not an image of the correct number type and dimensions, then assume that
;all pixels have the same sensitivity
;Check that "ron" and "gain" are numbers of the correct type, and assume standard values if they
;are out of range
;Check that "coeff2" and "coeff3" are numbers of the correct type
;Check that "out_type" is a string with an acceptable value
;Check that there is at least one good pixel
;Determine the effect of the non-linearity correction on the noise model, by inverting the non-linearity
;correction and differentiating. This is only necessary if the readout noise is greater than zero.
;Create a temporary image that is the calibrated image data multiplied by the master flat field image
;plus the dark current, and ensure that the bad pixels are set to sensible values
;Invert the polynomial non-linearity correction
;Calculate the roots of the quadratic equations
;For each quadratic equation, select the positive root closest to the input value. If such a root
;does not exist, then there is a problem with the non-linearity correction that has been used, and
;the module will exit
;Check that there is at least one root
;Determine the positive root closest to the input value
;Calculate the roots of the cubic equations
;For each cubic equation, select the positive root closest to the input value. If such a root
;does not exist, then there is a problem with the non-linearity correction that has been used, and
;the module will exit
;Determine the positive root closest to the input value
;Differentiate the inverted non-linearity correction
;Calculate an array of variance corrections due to the non-linearity correction
;Include zero or negative variance corrections due to the non-linearity correction in the output bad
;pixel mask, and check that there is still at least one good pixel
;Prepare the variance corrections for application in the CCD noise model
;Calculate the variance component from the image data "imdata"
;Add in the variance component from the master dark image
;Add in the variance component from the readout noise
;Prepare the output array of image pixel uncertainties, variances or inverse variances, as requested
;by the parameter "out_type"