; @(#)poly_diffim.pro	1.2 11/12/98 - NRL IDL LIBRARY
;
; Modified:
; 981104	James Tappin 	- changed hdr-> structure
;
pro Poly_diffim, dindex, dname=dname, in_place=in_place, $
                 divide=divide, cr_mask=cr_mask, fix_sum=fix_sum

;+
; POLY_DIFFIM
;	Multiple image differencer.
;
; Usage:
;	poly_diffim, {dindex | dname=dname}
;
; Argument:
;	dindex	int	input	The number of the subtrahend image.
;
; Keyword:
;	dname	string	input	The name of the subtrahend image
;				(which must be loaded).
;	in_place ??	input	If set, then do the subtraction in
;				place and overwrite the old image in
;				memory.
;	divide	int	input	Whether to subtract images (divide=0
;				or unset), divide images (1) or divide
;				and subtract 1 (2).
;	cr_mask	??	input	If set, then attempt to use
;				SIGMA_FILTER to remove cosmic rays.
;	fix_sum	??	input	If set, then divide the input image by
;				the onboard summing.
;
; Restrictions:
;	The DINDEX argument and the DNAME keyword are exclusive. Only
;	the currently selected images are processed (the subtrahend
;	isn't processed whether or not it is selected).
;
; Effects:
;	A set of difference images is generated and selected (the
;	previously selected images are deselected.
;
; History:
;	Original: 21/3/96; SJT
;	Modify to cope with assorted summings: Nov 96; SJT
;-

@chandle.com
@wload.com

if (n_params() eq 1 and keyword_set(dname)) then  $
  message, "May only use one means of defining the subtrahend."

if (n_params() eq 0 and not keyword_set(dname)) then  $
  message, "Must specify the subtrahend."

if (n_params() eq 0) then begin ; If specified by name find the image
                                ; to be subtracted.
    dindex = (where(dname eq string_name, imatch))(0)
    if (imatch eq 0) then $
      message, "Specified image not found in loaded set." $
    else if (imatch gt 1) then  $
      message, /inform, "Multiple matches, using first."
endif

Get_subt:

ghandle, dindex                 ; Get subtrahend

sr = sxpar(head, 'R1ROW') - 1
sc = (sxpar(head, 'R1COL')/2)*2 - 20

pxf = sxpar(head, 'R2ROW')-sxpar(head, 'R1ROW')+1
pxa = sz(2)
scale = float(pxf)/float(pxa)

subt_raw = fltarr(1024, 1024)
subt_raw(sc, sr) = rebin(image, sz(1)*scale, sz(2)*scale, /samp)

dname = name
shdr = head
crpixs = [sxpar(shdr, 'CENTER_X'),  $
          sxpar(shdr, 'CENTER_Y')]
print, crpixs
if (crpixs(0) ne 1) then begin
    scentre = 1b
    rads = sxpar(shdr, 'RADIUS')
endif else scentre = 0b

smask_raw = finite(subt_raw) and subt_raw ne 0

sel_img(dindex) = 0             ; Deselect the subtrahend.
sel_array = where(sel_img)

nimages = n_elements(sel_array)

for j = 0, nimages-1 do begin
    ghandle, sel_array(j)
    index = sel_array(j)
    ir = sxpar(head, 'R1ROW') - 1
    ic = (sxpar(head, 'R1COL')/2)*2 - 20
    n_psum = (sxpar(head, 'SUMROW') > sxpar(head, 'LEBYSUM'))* $
      (sxpar(head, 'SUMCOL') > sxpar(head, 'LEBXSUM')) > 1
    
    imscl = float(sz(2)) / $
      float(sxpar(head, 'R2ROW') - sxpar(head, 'R1ROW') + 1)
    
    if (imscl ne 1.0) then begin
        subt = rebin(/samp, subt_raw, 1024*imscl, 1024*imscl)
        smask = rebin(/samp, smask_raw, 1024*imscl, 1024*imscl)
        ir = ir*imscl
        ic = ic*imscl
    endif else begin
        subt = subt_raw
        smask = smask_raw
    endelse
    
    if (keyword_set(cr_mask)) then  $
      image = sigma_filter(image, 15, n_sigma = 4., /all) 
    
    if (keyword_set(fix_sum)) then image = image/float(n_psum)
    
    if (keyword_set(divide)) then  begin
        name = name + '-:-' + dname
        divok = where(smask(ic:ic+sz(1)-1, ir:ir+sz(2)-1), ndivok)
        if (ndivok ne 0) then begin
            image(divok) = image(divok) / $
              (subt(ic:ic+sz(1)-1, Ir:ir+sz(2)-1))(divok)
            image = image*smask(ic:ic+sz(1)-1, ir:ir+sz(2)-1)
            if (divide eq 2) then image = image-1.0
        endif
    endif else begin
        name = name + '-' + dname
        image = (image-subt(ic:ic+sz(1)-1, ir:ir+sz(2)-1))* $
          smask(ic:ic+sz(1)-1, ir:ir+sz(2)-1)
    endelse        
    
    print, j, sel_array(j), name, format = "(2I4,2x,a)"
    
    sxaddpar, head, 'FILENAME', name
    sxaddpar, head, 'HISTORY', 'Differences by POLY_DIFFIM'
    
    crpixi = [sxpar(head, 'CENTER_X'),  $
              sxpar(head, 'CENTER_Y')]
    if ((!Err eq -1 or crpixi(0) eq 1) and scentre) then begin
        sxaddpar, head, 'CENTER_X', crpixs(0)-ic+sc, 'From Subtrahend'
        sxaddpar, head, 'CENTER_Y', crpixs(1)-ir+sr, 'From Subtrahend'
        sxaddpar, head, 'RADIUS', rads, 'From Subtrahend'
    endif
    
    
    if (keyword_set(in_place)) then begin
        handle_value, h_array(index), image, /set
        handle_value, h_name(index), name, /set
        handle_value, h_head(index), head, /set
        h = lasco_fitshdr2struct(head)
        handle_value, h_hs(index), h, /set
        a3 = rebin(image, sz(1)/8, sz(2)/8) ; Borut's redfined A3
;        a3(*) = 0
;        a3r = (h.r1row-1)/8
;        a3c = (h.r1col-21)/8
;        a3(a3c, a3r) = tmp
        handle_value, h_small(index), a3, /set
    endif else begin
        sel_img(sel_array(j)) = 0
        chandle
        sel_img = [sel_img, 1]
    endelse
endfor

;	Make sure the names array is up to date.

s_handle = N_ELEMENTS(h_name)
string_name = strarr(s_handle)
 
for i = 0, s_handle-1 do begin
    HANDLE_VALUE, h_name(i), string_temp
    string_name(i) = string_temp
endfor

sel_array = where(sel_img)

end