FUNCTION remove_vis_xtalk, im_in, $
                           xtalk_im, $
                           mode = mode, $
                           instrument = instrument

; Use the signal in im_in to determine the crosstalk signal in
; the other quadrants
; Assumes that the image is ALREADY bias subtracted

; assume a 5th order crosstalk gain for possible expansion.
; Will start with a linear gain

gain1 = FLTARR( 6, 3 ) ; by quadrant
gain2 = FLTARR( 6, 3 )
gain3 = FLTARR( 6, 3 )
gain4 = FLTARR( 6, 3 )

xshift1 = INTARR( 3 )
xshift2 = INTARR( 3 )
xshift3 = INTARR( 3 )
xshift4 = INTARR( 3 )

CASE instrument of
    1 : BEGIN                   ; HRI

        gain1(1,0) = 4.0E-4     ;Done
        gain1(1,1) = 3.5E-4     ;Done         ;5.0E-4 (TV4)
        gain1(1,2) = 3.5e-4     ;Done         ;4.0E-4 (TV4)
        gain2(1,0) = 3.5E-4     ;Done
        gain2(1,1) = 7.0E-4     ;Done         ;8.0E-4 (TV4)
        gain2(1,2) = 4.0e-4     ;Done         ;5.0E-4 (TV4)
        gain3(1,0) = 3.0E-4     ;Done
        gain3(1,1) = 8.0E-4     ;Done
        gain3(1,2) = 5.0e-4     ;Done         ;4.5E-4 (TV4)
        gain4(1,0) = 3.3E-4     ;Done
        gain4(1,1) = 4.0E-4     ;Done
        gain4(1,2) = 6.0E-4     ;Done         ;6.0E-4 (TV4) 
        
        xshift1 = [-1,1,-1]     ;Done
        xshift2 = [1,1,-1]      ;done
        xshift3 = [1,-1,-1]     ;Done
        xshift4 = [1,-1,1]      ;Done


    END
    2 : BEGIN                   ; MRI 


        gain1(1,0) = 6.0E-4    ;Done
        gain1(1,1) = 5.0E-4    ;Done
        gain1(1,2) = 8.0e-4    ;Done         ;2.0E-4 (TV4) 
        gain2(1,0) = 5.0E-4    ;Done 
        gain2(1,1) = 3.0E-4    ;Done
        gain2(1,2) = 3.0E-4    ;Done         ;4.0E-4 (TV4) 
        gain3(1,0) = 4.0E-4    ;Done 
        gain3(1,1) = 4.0E-4    ;Done
        gain3(1,2) = 3.5e-4    ;Done         ;2.0E-4 (TV4)
        gain4(1,0) = 9.0E-4    ;Done
        gain4(1,1) = 3.5E-4    ;Done 
        gain4(1,2) = 3.0E-4    ;Done

        xshift1 = [-1,1,-1]    ;Done
        xshift2 = [1,1,-1]     ;Done
        xshift3 = [1,-1,-1]    ;Done
        xshift4 = [1,-1,1]     ;Done


    END
    3 : BEGIN                  ;ITS

        gain1(1,0) = 6.5E-4    ;Done
        gain1(1,1) = 3.0E-4    ;Done
        gain1(1,2) = 14.0E-4   ;Done
        gain2(1,0) = 6.0E-4    ;Done 
        gain2(1,1) = 5.0E-4    ;Done
        gain2(1,2) = 5.0E-4    ;Done
        gain3(1,0) = 3.0E-4    ;Done
        gain3(1,1) = 6.0E-4    ;Done
        gain3(1,2) = 10.0E-4   ;Done
        gain4(1,0) = 15.0E-4   ;Done
        gain4(1,1) = 5.5E-4    ;Done
        gain4(1,2) = 8.0E-4    ;Done

        xshift1 = [-1,1,-1]    ;Done
        xshift2 = [1,1,-1]     ;Done
        xshift3 = [1,-1,-1]    ;Done
        xshift4 = [1,-1,1]     ;Done

    END
    ELSE : BEGIN
        PRINT, '^^^ remove_vis_xtalk: Instrument ID unknown'
        im_out = -1
        GOTO, Finish
    END
ENDCASE

imsize = SIZE( im_in )
IF imsize( 0 ) NE 2 THEN BEGIN
    PRINT, '^^^ remove_vis_xtalk: Image is not 2-D...'
    im_out = -1
    GOTO, Finish
ENDIF

xsize = imsize( 1 )
ysize = imsize( 2 )

CASE mode OF
    1 : BEGIN
        xtalk_im = FLTARR( 1024, 1024 )
        q1 = im_in(   9: 511, 512:1014 )
        q2 = im_in( 512:1014, 512:1014 )
        q3 = im_in(   9: 511,   9: 511 )
        q4 = im_in( 512:1014,   9: 511 )
    END
    2 : BEGIN
        xtalk_im = FLTARR( 512, 512 )
        q1 = im_in(    5:255, 256:506 )
        q2 = im_in(  256:506, 256:506 )
        q3 = im_in(    5:255,   5:255 )
        q4 = im_in(  256:506,   5:255 )
    END
    ELSE : BEGIN
        PRINT, 'Modes other than 1 and 2 are not supported yet'
        im_out = -1
        GOTO, Finish
    END
ENDCASE

q1x_2 = POLY( q1, gain2( *, 0 ) )
q1x_3 = POLY( q1, gain3( *, 0 ) )
q1x_4 = POLY( q1, gain4( *, 0 ) )

q2x_1 = POLY( q2, gain1( *, 0 ) )
q2x_3 = POLY( q2, gain3( *, 1 ) )
q2x_4 = POLY( q2, gain4( *, 1 ) )

q3x_1 = POLY( q3, gain1( *, 1 ) )
q3x_2 = POLY( q3, gain2( *, 1 ) )
q3x_4 = POLY( q3, gain4( *, 2 ) )

q4x_1 = POLY( q4, gain1( *, 2 ) )
q4x_2 = POLY( q4, gain2( *, 2 ) )
q4x_3 = POLY( q4, gain3( *, 2 ) )

; Produce the rotated images, and produce INTEGER images

q1x_2 = SHIFT( ROTATE( q1x_2, 5), xshift1( 0 ) )
q1x_3 = SHIFT( ROTATE( q1x_3, 7), xshift1( 1 ) )
q1x_4 = SHIFT( ROTATE( q1x_4, 2), xshift1( 2 ) )

q2x_1 = SHIFT( ROTATE( q2x_1, 5), xshift2( 0 ) )
q2x_3 = SHIFT( ROTATE( q2x_3, 2), xshift2( 1 ) )
q2x_4 = SHIFT( ROTATE( q2x_4, 7), xshift2( 2 ) )

q3x_1 = SHIFT( ROTATE( q3x_1, 7), xshift3( 0 ) )
q3x_2 = SHIFT( ROTATE( q3x_2, 2), xshift3( 1 ) )
q3x_4 = SHIFT( ROTATE( q3x_4, 5), xshift3( 2 ) )

q4x_1 = SHIFT( ROTATE( q4x_1, 2), xshift4( 0 ) )
q4x_2 = SHIFT( ROTATE( q4x_2, 7), xshift4( 1 ) )
q4x_3 = SHIFT( ROTATE( q4x_3, 5), xshift4( 2 ) )

; Produce the crosstalk image

CASE mode OF
    1 : BEGIN
        xtalk_im(   9: 511, 512:1014 ) = q2x_1 + q3x_1 + q4x_1
        xtalk_im( 512:1014, 512:1014 ) = q1x_2 + q3x_2 + q4x_2
        xtalk_im(   9: 511,   9: 511 ) = q1x_3 + q2x_3 + q4x_3
        xtalk_im( 512:1014,   9: 511 ) = q1x_4 + q2x_4 + q3x_4
    END
    2 : BEGIN
        xtalk_im(    5:255, 256:506 ) = q2x_1 + q3x_1 + q4x_1
        xtalk_im(  256:506, 256:506 ) = q1x_2 + q3x_2 + q4x_2
        xtalk_im(    5:255,   5:255 ) = q1x_3 + q2x_3 + q4x_3
        xtalk_im(  256:506,   5:255 ) = q1x_4 + q2x_4 + q3x_4
    END
    ELSE : BEGIN
        PRINT,' How did you get here???'
    END
ENDCASE

im_out = im_in - xtalk_im

;print,mean(im_out(8:511   ,512:1014)),median(im_out(8:511   ,512:1014))
;print,mean(im_out(512:1014,512:1014)),median(im_out(512:1014,512:1014))
;print,mean(im_out(8:511   ,9:511)),median(im_out(8:511   ,9:511))
;print,mean(im_out(512:1014,9:511)),median(im_out(512:1014,9:511))

Finish:

RETURN, im_out
END