PRO mk_xrt_composite, l_idx, l_da, s_idx, s_da, c_idx, c_da, $
                      despike=despike, lsat_map=lsat_map, verbose=verbose

; =========================================================================
;+
; PROJECT:
;       Solar-B / XRT
;
; NAME:
;
;       MK_XRT_COMPOSITE
;
; CATEGORY:
;
;       Data graphics.
;
; PURPOSE:
;
;       Combine a SINGLE long/short exposure pair of images to increase the
;       dynamic range and to replace saturation with unsaturated data.
;       See Note #1.
;
; CALLING SEQUENCE:
;
;       MK_XRT_COMPOSITE, l_idx, l_da, s_idx, s_da, c_idx, c_da
;                         [,despike=despike] [,/verbose]
;
; INPUTS:
;
;       L_IDX    - [Mandatory] (structure scalar)
;                  This is the "index" for the long exposure image.
;                  Can be Level-0 or Level-1, as long as it matches L_DA.
;       L_DA     - [Mandatory] (2D number array, [Nx,Ny])
;                  This is the 2D data array for the long exposure image.
;                  Can be Level-0 or Level-1, as long as it matches L_IDX.
;       S_IDX    - [Mandatory] (structure scalar)
;                  This is the "index" for the short exposure image.
;                  Can be Level-0 or Level-1, as long as it matches S_DA.
;       S_DA     - [Mandatory] (2D number array, [Nx,Ny])
;                  This is the 2D data array for the short exposure image.
;                  Can be Level-0 or Level-1, as long as it matches S_IDX.
;
; KEYWORDS:
;
;       DESPIKE  - [Optional] (Boolean or integer) Set if 1 pass of 
;                  despike wanted to remove radiation-belt/cosmic-ray 
;                  spikes, or set to 1-3 passes of despike. (Two passes 
;                  works well even for heavy spiking.) Default = no 
;                  despiking if the keyword is not set. 
;                  (See <xrt_prep.pro> for more information.)
;       LSAT_MAP - [Optional] (Byte array) Set to a boolian map of the 
;                  location of saturated pixels to be replaced. Use this 
;                  if the data have already been xrt_prepped and the map
;                  was generated or if xrt_prep is not to be used. Array 
;                  must be same size as images to be composited. 1 represents
;                  the pixel to be replaced.
;       /VERBOSE - [Optional] (Boolean) If set, print out extra messages.
;
; OUTPUTS:
;
;       C_IDX    - [Mandatory] (structure scalar)
;                  This is the "index" for the composite image. It will
;                  have the same content as L_IDX (after it has been 
;                  through <xrt_prep.pro> and normalized), except that
;                  there will be an update to the HISTORY field.
;       C_DA     - [Mandatory] (2D number array, [Nx,Ny])
;                  This is the 2D data array for the composite image.
;
; EXAMPLES:
;
;       Basic usage:
;       IDL> mk_xrt_composite, l_idx, l_da, s_idx, s_da, c_idx, c_da
;       IDL> tvscl, xrtdisp(c_da, /log)
;
;       Recommended to do despiking on Level-0 data:
;       IDL> mk_xrt_composite, l_idx, l_da, s_idx, s_da, c_idx, c_da, $
;       IDL>                   despike=2
;       IDL> tvscl, xrtdisp(c_da, /log)
;
; COMMON BLOCKS:
;
;       none
;
; NOTES:
;
;      1) A simple proceedure for combining long and short XRT images. 
;         This procedure simply generates a composite image for XRT
;         data.  Only makes 1 composite at a time. It checks to see
;         if the images have been xrt_prepped, and preps them if
;         they have not been. The index structure for the long image
;         will be used for the combined index structure with the
;         history updated with the short image information.
;      2) If xrt_prep is not used to process the images both the long
;         and short images must be normalized for exposure before they 
;         are inputted into this routine. Exposure normalization ensures
;         there is no intensity variation for the combined image.
;
; CONTACT:
;
;       Comments, feedback, and bug reports regarding this routine may be
;       directed to this email address:
;                xrt_manager ~at~ head.cfa.harvard.edu
;
; MODIFICATION HISTORY:
; 
progver = 'v2006-May-09' ;--- (P.Jibben) Written.
progver = 'v2007-Jun-13' ;--- (P.Jibben, S.Saar) This will handle the 
;                             bloom pixels in a special way without 
;                             altering the other data. Also changed
;                             the way the history is saved, as well
;                             as changed the data_lev keyword to 2.
progver = 'v2007-Jun-28' ;--- (P.Jibben) This will allow users to pass
;                             their own saturated pixel map so they are
;                             not forced to use xrt_prep.
progver = 'v2007-Jul-13' ;--- (M.Weber) Do not replace saturated
;                             pixels with values below the
;                             minimum value of the long exposure.
;                             I.e., preserve minimum level of long
;                             exposure.
;
;-
; =========================================================================

  despike=keyword_set(despike)
  q_lsat_map=keyword_set(lsat_map)
  
  loud=keyword_set(verbose)
  if (not keyword_set(verbose)) then quiet=1

  ; Added this bit so that inputs are not overwritten.
  l_idx1 = l_idx
  l_da1  = l_da
  s_idx1 = s_idx
  s_da1  = s_da
  if (q_lsat_map eq 1) then l_map=lsat_map

  ; determine if prepping is necessary.
  l_prep=get_history(l_idx1,'XRT_PREP',found=l_found)
  s_prep=get_history(s_idx1,'XRT_PREP',found=s_found)

  ; If you prep one you must prep both.
  if (l_found ne s_found) then begin
     print,'Long/Short images must have same type of processing... Returning.'
     return
  endif

  ; If you prep the long you must provide the saturated pixel map.
  if (l_found eq 1 and q_lsat_map eq 0) then begin
     print,'Must use the lsat_map keyword for xrt_prepped data... Returning.'
     return
  endif

  ; Make sure that l_da & lsat_map dimensions match. 
  if (l_found eq 1 or q_lsat_map eq 1) then begin
        sz_map=size(l_map,/dimensions)
        sz_lda=size(lda1,/dimensions)
        ndim1=n_elements(l_da1)
        ndim2=n_elements(l_map)
        if (ndim1 ne ndim2) then begin
           print,'Error long image and sat_map array must have same dimensions... Returning.'
           return
        endif
     endif
  
  ;Need to make sure that images are normalized for exposure if xrt_prep is used.
  ; If xrt_prep not used, assume user normalized images themselves.
  l_norm=get_history(l_idx1,'XRT_RENORMALIZE',found=nl_found)
  s_norm=get_history(s_idx1,'XRT_RENORMALIZE',found=ns_found)

  ; Perhaps they normalized them themselves.
  if (l_found eq 1 and nl_found eq 0) then begin
   l_da1=temporary(l_da1)/get_xrt_expdur(l_idx1, /sec)
   l_idx1.e_etim=1000000 ; [microsec] for consistancy with XRT_PREP & returned index
   if (loud eq 1) then print,'Renormalizing the exposures.'
  endif

  if (s_found eq 1 and nl_found eq 0) then begin
   s_da1=temporary(s_da1)/get_xrt_expdur(s_idx1, /sec)
  endif

  ; Use xrt_prep on images... If needed.
  if (s_found eq 0) then begin
   if (loud eq 1) then print,'Prepping short for composite'
   s_map=1
   ;ns_pix=1
  ; xrt_prep, temporary(s_idx1), temporary(s_da1), s_idx1, s_da1, despike=despike, $
  ;           /float, /normalize, sat_map=s_map, n_sat_pixels=ns_pix, $
  ;           verbose=verbose, quiet=quiet
   xrt_prep, temporary(s_idx1), temporary(s_da1), s_idx1, s_da1, despike=despike, $
             /float, /normalize, n_sat_pixels=ns_pix, $
             verbose=verbose, quiet=quiet
endif

  if (l_found eq 0) then begin
   if (loud eq 1) then print,'Prepping long for composite'
   l_map=1
   nl_pix=1
   xrt_prep, temporary(l_idx1), temporary(l_da1), l_idx1, l_da1, despike=despike, $
             /float, /normalize, sat_map=l_map, n_sat_pixels=nl_pix, $
             verbose=verbose, quiet=quiet
  endif

  if (q_lsat_map eq 1) then nl_pix = 1
  ; Use sat_map to replace saturated pixels of long exposure.
  if (nl_pix gt 0) then begin
   ss=where(l_map eq 1)
   l_da1[ss]=s_da1[ss] > min(l_da1)
   c_da=l_da1
  endif else begin
   c_da=l_da1
   message,'No saturated pixels to replace returning long image',/continue
  endelse

  ;; ----- Added 11-June-2007 -----;;
	; handeling the edge pixels
	; using Saar's method of filtering whole image and then replacing, not just filtering bloom pixels.
	; P.R. Jibben -He also calcuated shift better, only in vertical direction.
	; Only doing this for long image
	siz=size(c_da)
	ny=siz(2)
 	donut=[[1,1,1],[1,0,1],[1,1,1]] / 8.
	Result=CONVOL(c_da,donut)
	sl_map1=shift(l_map,[0,1])
	sl_map1[*,0]=0
	sl_map2=shift(fix(l_map),[0,-1])
	sl_map2[*,ny-1]=0
	bloomers=temporary(sl_map1)+temporary(sl_map2)
	bloomers=(bloomers eq 1)	
	replace_pix=where(bloomers eq 1,cnt)
	bloom_tag=''
	if (cnt gt 0) then begin
	 nbloom=strcompress(cnt,/rem)
	 c_da[replace_pix]=Result[replace_pix]
	 bloom_tag='Replaced ' + nbloom + ' bloom pixels in Composite.'
	endif

  ;generating the c_idx

  c_idx=l_idx1
  ;; ----- Added 11-June-2007 ----- ;;
  ;Updating the index data level
  ;level2 data
  c_idx.data_lev=2

  tagval0=''
  change_tag_value,c_idx,tagval0,'HISTORY'
  tagval0=[tagval0,'XRT_MK_COMPOSITE ' + progver + ': ']
  long_hist=get_history(l_idx1)
  short_hist=get_history(s_idx1)
  tagval0=[tagval0,bloom_tag,'Long image history --',long_hist,'Short image history --',short_hist]
  update_history,c_idx,tagval0,/noroutine
  ;update_history, c_idx, long_hist, /noroutine
  ;update_history, c_idx, 'XRT_MK_COMPOSITE ' + progver + $
  ;                ': Short image history --', $
  ;               /noroutine
  ;update_history, c_idx, short_hist, /noroutine

  RETURN

END ;======================================================================