
PRO FLAT_KUHN_LIN_I, con, gain, bitmap, cnt, p, disp
;+
; NAME:
;    FLAT_KUHN_LIN_I
;
; PURPOSE:
;    Iterate Kuhn-Lin style gain table (see FLAT_KUHN_LIN.PRO).
;
; CATEGORY:
;    BBSO Archiving System.
;    Image Processing.
;
; CALLING SEQUENCE:
;    FLAT_KUHN_LIN_I, con, gain, bitmap, cnt, p, disp
;
; INPUTS:
;    CON:    Constant term of the gain table.
;    BITMAP: Bitmap containing good pixel for each input flat field frame.
;    CNT:    Valid pixel pair count.
;    DISP:   Array of displacements with respect to the centered image, also
;            reduced in size.
;    P:      STRUCTURE image specs.
;
; OUTPUTS:
;    GAIN:   Reduced size gain table.
;
; KEYWORDS:
;    None.
;
; MODIFICATION HISTORY:
;    1997-10-27: Anders Johanneson, Big Bear Solar Observatory
;        Original verion of DO_ITER.PRO.
;    1999-12-15: Carsten Denker, Big Bear Solar Observatory
;        Intergration in the new BBSO archiving system, variables are passed
;        as structures.
;
;-

    ;definitions
loopct = 0
gain_n = con

    ;----- recompute gain
FOR iq = 1, p.n - 1 DO BEGIN
    mskiq = bitmap / ( 2l^iq OR bitmap )
    FOR ir = 0, iq - 1 DO BEGIN
        mskir = bitmap / ( 2l^iq OR bitmap )
        dx = disp( 0, iq ) - disp( 0, ir )
        dy = disp( 1, iq ) - disp( 1, ir )
        jxl = MAX( [ 0, -dx ] )
        jxh = MIN( [ 0, -dx ] ) + p.nxx - 1
        jyl = MAX( [ 0, -dy ] )
        jyh = MIN( [ 0, -dy ] ) + p.nyy - 1
        ixl = MAX( [ 0,  dx ] )
        ixh = MIN( [ 0,  dx ] ) + p.nxx - 1
        iyl = MAX( [ 0,  dy ] )
        iyh = MIN( [ 0,  dy ] ) + p.nyy - 1
        msk = ( mskiq( jxl : jxh, jyl : jyh ) AND mskir( ixl : ixh, iyl : iyh ) )
        gain_n( jxl : jxh, jyl : jyh ) = gain_n( jxl : jxh, jyl : jyh ) $
                                         + gain( ixl : ixh, iyl : iyh ) * msk
        gain_n( ixl : ixh, iyl : iyh ) = gain_n( ixl : ixh, iyl : iyh ) $
                                         + gain( jxl : jxh, jyl : jyh ) * msk

        loopct = loopct + 1
    ENDFOR
ENDFOR
gain_n = gain_n / ( cnt > 1 )
index = WHERE( cnt EQ 0, ct )
IF ct GT 0 THEN gain_n( index ) = 0.
index = WHERE( cnt GT 0, n )

    ;----- calculate average of gain table for normalization (Caution: must use double
    ;----- floats here!)
sum2 = TOTAL( DOUBLE( gain_n( index ) ) )
sum3 = TOTAL( DOUBLE( gain_n( index )^2 ) )

    ;----- ignore pixels more than 5*sigma away from mean
ave2 = sum2 / n
five_sigma = 5 * SQRT( sum3 / n - ave2 * ave2 )

index = WHERE( ABS( gain_n - ave2 ) GT five_sigma, ct )
IF ct GT 0 THEN BEGIN
    sum2 = sum2 - TOTAL( DOUBLE( gain_n ( index ) ) )
    n = n - ct
ENDIF

    ;----- normalize this iteration's gain table
ave2 = sum2 / n

IF ct GT 0 THEN BEGIN
   gain_n = gain_n - ave2
   sum3 = TOTAL( DOUBLE( ( gain_n - gain )^ 2 ) )
ENDIF
gain = gain_n



RETURN

END



