PRO eis_add_wd_to_spec, wd, mask, wpix, spec, quiet=quiet, iexp=iexp, $
                        sum=sum

;+
; NAME
;
;     EIS_ADD_WD_TO_SPEC
;
; PROJECT
;
;     Hinode/EIS
;
; EXPLANATION
;
;     This routine is intended to work with the EIS_PIXEL_MASK and
;     EIS_MASK_SPECTRUM routines. Spatial pixels that are identified
;     in the pixel mask (MASK) are averaged to yield a 1D spectrum
;     that goes into the output SPEC structure. The format of the SPEC
;     structure is defined in EIS_MASK_SPECTRUM; EIS_ADD_WD_TO_SPEC
;     fills in the tags of this structure.
;
;     A minority of EIS data-sets take multiple exposures at each
;     raster position and there are two options for how the spectrum
;     is created in this case:
;
;     (i) (Default) The routine identifies the exposures with the
;         largest exposure times, and then averages only these. The
;         largest exposure times are defined to be those > 90% of the
;         maximum exposure time. E.g., if the exposure times are 5s,
;         10s, 95s and 100s, then the 95 and 100s exposures will be
;         averaged, and the shorter exposures ignored.
;     (ii) A specific exposure time can be specified using the IEXP=
;         keyword. If there are five different exposure times, then
;         IEXP can take values between 0 and 4.
;
;     In some circumstances it may be preferable to sum the specta at
;     the different spatial locations instead of averaging them. In
;     this case, the /SUM keyword should be set. Applications include
;     the case where a small bright point is not resolved by the
;     instrument. The true intensity of the bright point is then the
;     sum of spectra that have been dispersed by the instrument.
;
; INPUTS
;
;     WD     A structure in the format returned by EIS_GETWINDATA()
;
;     MASK   A pixel mask in the format returned by EIS_PIXEL_MASK.
;
;     WPIX   A vector of pixel numbers that correspond to where the
;            EIS window data appear on the CCD.
;
;     SPEC   A spectrum structure.
;
; OPTIONAL INPUTS
;
;     IEXP   If the data-set has multiple exposures per raster position,
;            then setting IEXP allows a particular exposure to be
;            selected. 
;
; KEYWORDS
;
;     QUIET  If set then routine will not generate text warnings.
;
;     SUM    If set, then the spectra will be summed instead of averaged
;
; OUTPUTS
;
;     The spectrum structure, SPEC, is modified to include the
;     averaged spectra derived from WD and the pixel mask.
;
; CALLS
;
;     EIS_ADJUST_MASK
;
; HISTORY
;
;     Ver.1, 27-Oct-2008, Peter Young
;     Ver.2, 15-Dec-2009, Peter Young
;       added /QUIET keyword; updated header; now prints warning if
;       pixel mask has been truncated.
;     Ver.3, 10-May-2010, Peter Young
;       I now pass the observation time to eis_adjust_mask.
;     Ver.4, 25-Aug-2010, Peter Young
;       Routine now deals with data-sets with multiple exposures per
;       raster position (nexp_prp>1).
;     Ver.5, 2-Apr-2012, Peter Young
;       Added IEXP= keyword to provide an alternative option for
;       data-sets with nexp_prp>1.
;     Ver.6, 30-May-2012, Peter Young
;       Added keyword /SUM.
;     Ver.7, 12-Aug-2025, Peter Young
;       For /sum keyword, I now scale the intensity and error by the
;       data quality to prevent dips in the spectrum due to missing
;       data.
;-

nexp_prp=wd.hdr.nexp_prp

;
; For data-sets with nexp_prp>1, the following code checks the
; exposure times of each exposure. The idea is that if there is a
; short exposure and a long exposure (5s and 50s, say), then we only
; want to use the long exposure. However, if there are two exposures
; with the same exposure time then we want to use both exposures.
;
; Another option has been added whereby a specific exposure can be
; selected (IEXP).
;
IF nexp_prp GT 1 THEN BEGIN
  IF n_elements(iexp) NE 0 THEN BEGIN
    k_nexp_prp=iexp
    n_nexp_prp=n_elements(iexp)
    exposure_times=wd.exposure_time[*,iexp]
    exp_time=median(exposure_times[where(exposure_times NE 0.)])
  ENDIF ELSE BEGIN
    exp_times=fltarr(nexp_prp)
    FOR i=0,nexp_prp-1 DO BEGIN
      exposure_times=wd.exposure_time[*,i]
      exp_times[i]=median(exposure_times[where(exposure_times NE 0.)])
    ENDFOR 
    exp_time_max=max(exp_times)
    k=where(exp_times/exp_time_max GT 0.90,nk)
   ;
    k_nexp_prp=k    ; position(s) of largest exposure times
    n_nexp_prp=nk   ; number of exposures with large exposure times
    IF keyword_set(sum) THEN BEGIN
      exp_time=total(exp_times)
    ENDIF ELSE BEGIN
      exp_time=average(exp_times)
    ENDELSE 
  ENDELSE 
ENDIF ELSE BEGIN
  k_nexp_prp=0
  n_nexp_prp=1
  exposure_times=wd.exposure_time
  exp_time=median(exposure_times[where(exposure_times NE 0.)])
ENDELSE

IF spec.exposure_time EQ -1. THEN spec.exposure_time=exp_time

int=wd.int
err=wd.err
;
siz=size(int)

qual_max=n_elements(where(mask.image EQ 1))
spec.qual_max=qual_max*n_nexp_prp

;
; For each wavelength I average the image over the selected pixels.
;
info_swtch=0
FOR j=0,siz[1]-1 DO BEGIN
 ;
 ; Create the pixel mask for this wavelength
 ;
  wvl_check=spec.wvl[wpix[j]]
  newmask=eis_adjust_mask(mask,wvl_check,info=info,time=wd.hdr.date_obs)
  IF info EQ 1 THEN info_swtch=1
  swtch=newmask.image
 ;
 ; Case 1
 ; ------
 ; Only one exposure is used to derive the average spectum.
 ;
  IF n_nexp_prp EQ 1 THEN BEGIN
    IF nexp_prp EQ 1 THEN BEGIN
      int_img=reform(int[j,*,*]) ; 2D image for wavelength pixel j
      err_img=reform(err[j,*,*]) ; 2D error image for wavelength pixel j
    ENDIF ELSE BEGIN
      int_img=reform(int[j,*,*,k_nexp_prp]) ; 2D image for wavelength pixel j
      err_img=reform(err[j,*,*,k_nexp_prp]) ; 2D error image for wavelength pixel j
    ENDELSE
  ;
  ; Remove any missing pixels from the mask
  ;
   k=where(err_img EQ wd.missing,nk)
   IF nk NE 0 THEN swtch[k]=0b
  ;
   k=where(swtch EQ 1,nk)
   spec.qual[wpix[j]]=nk
   ;
   IF nk NE 0 THEN BEGIN
     IF keyword_set(sum) THEN BEGIN 
       spec.int[wpix[j]]=total(int_img[k]) * $
                         float(spec.qual_max) / float(spec.qual[wpix[j]])
       spec.err[wpix[j]]=sqrt( total( err_img[k]^2 ) )  * $
                         float(spec.qual_max) / float(spec.qual[wpix[j]])
     ENDIF ELSE BEGIN
       spec.int[wpix[j]]=average(int_img[k])
       spec.err[wpix[j]]=sqrt( total( err_img[k]^2 ) ) / float(nk)
     ENDELSE 
    ;
   ENDIF ELSE BEGIN
     spec.int[wpix[j]]=wd.missing
     spec.err[wpix[j]]=wd.missing
   ENDELSE
 ;
 ;
 ; Case 2
 ; ------
 ; The following code deals with the case when more than one exposure
 ; is used to generate the averaged spectrum.
 ;
  ENDIF ELSE BEGIN
    int_all=-1
    err_all=-1
    FOR m=0,n_nexp_prp-1 DO BEGIN
      int_img=reform(int[j,*,*,k_nexp_prp[m]]) ; 2D image for wavelength pixel j
      err_img=reform(err[j,*,*,k_nexp_prp[m]]) ; 2D error image for wavelength pixel j
     ;
      swtch_m=swtch
     ;
      k=where(err_img EQ wd.missing,nk)
      IF nk NE 0 THEN swtch_m[k]=0b
     ;
      k=where(swtch_m EQ 1,nk)
      IF nk NE 0 THEN BEGIN
        IF int_all[0] EQ -1 THEN BEGIN
          int_all=int_img[k]
          err_all=err_img[k]
        ENDIF ELSE BEGIN
          int_all=[int_all,int_img[k]]
          err_all=[err_all,err_img[k]]
        ENDELSE 
        spec.qual[wpix[j]]=spec.qual[wpix[j]]+nk
      ENDIF
    ENDFOR
    ;
    IF spec.qual[wpix[j]] NE 0 THEN BEGIN
      n_all=n_elements(int_all)
      IF keyword_set(sum) THEN BEGIN 
        spec.int[wpix[j]]=total(int_all) * $
                          float(spec.qual_max) / float(spec.qual[wpix[j]])
        spec.err[wpix[j]]=sqrt( total( err_all^2 ) ) * $
                          float(spec.qual_max) / float(spec.qual[wpix[j]])
      ENDIF ELSE BEGIN
        spec.int[wpix[j]]=average(int_all)
        spec.err[wpix[j]]=sqrt( total( err_all^2 ) ) / float(n_all)
      ENDELSE 
    ENDIF ELSE BEGIN
    ENDELSE 
  ENDELSE 
ENDFOR

IF info_swtch EQ 1 AND NOT keyword_set(quiet) THEN BEGIN
  message,/info,/cont,'WARNING: mask truncated for some pixels in this window.'
ENDIF 

END