;+
; NAME:
;         p3d_wavecal_match_maskwd
;
;         $Id: p3d_wavecal_match_maskwd.pro 79 2010-03-04 14:24:25Z christersandin $
;
; PURPOSE:
;         This routine finds out which lines in the calibration line list
;         (REFWAVE, at positions REFPIXL) are actually found in the spectrum
;         SPEC(WAVELEN).
;
;
; AUTHOR:
;         Christer Sandin
;         Astrophysikalisches Institut Potsdam (AIP)
;         An der Sternwarte 16
;         D-14482 Potsdam, GERMANY
;
; COPYRIGHT:
;         p3d: a general data-reduction tool for fiber-fed IFSs
;
;         Copyright 2009,2010 Astrophysikalisches Institut Potsdam (AIP)
;
;         This program is free software; you can redistribute it and/or modify
;         it under the terms of the GNU General Public License as published by
;         the Free Software Foundation; either version 3 of the License, or
;         (at your option) any later version.
;
;         This program is distributed in the hope that it will be useful, but
;         WITHOUT ANY WARRANTY; without even the implied warranty of
;         MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
;         General Public License for more details.
;
;         You should have received a copy of the GNU General Public License
;         along with this program; if not, see <http://www.gnu.org/licenses>.
;
;         Additional permission under GNU GPL version 3 section 7
;
;         If you modify this Program, or any covered work, by linking or
;         combining it with IDL (or a modified version of that library),
;         containing parts covered by the terms of the IDL license, the
;         licensors of this Program grant you additional permission to convey
;         the resulting work.
;
; CATEGORY:
;         p3d :: wavelength calibration
;
; CALLING SEQUENCE:
;         p3d_wavecal_match_maskwd,spec,wavelen,refpixl,refwave,fwhm, $
;             identline,identpixl,identwave,noise=,poldeg=,niterate=, $
;             findwidth=,topwid=,logunit=,verbose=,error=,/debug,/help
;
; INPUTS:
;         spec            - A one-dimensional array containing the spectrum
;                           against which the matching is done. It has to be
;                           of floating point type.
;         wavelen         - A one-dimensional wavelength array corresponding
;                           to SPEC.
;         refpixl         - A one-dimensional array that specifies the
;                           location of every line given in REFWAVE in the
;                           array WAVELEN.
;         refwave         - A one-dimensional array that specifies the
;                           wavelengths of known emission lines; used together
;                           with REFPIXL.
;         fwhm            - A scalar decimal value (integer is also ok)
;                           specifying the full width at half maximum value
;                           that is used when summing up the contribution at
;                           the location of a potential emission line; given
;                           in pixels.
;
; KEYWORD PARAMETERS:
;         noise           - A scalar value specifying the noise. If this
;                           variable is not given then it is determined from
;                           regions in SPEC where there are no calibration
;                           lines specified.
;         poldeg [1]      - The degree of the polynomial, that is used to fit
;                           REFWAVE(REFPIXL).
;         niterate [1]    - A scalar integer specifying how many times the
;                           solution is iterated.
;         findwidth [fwhm] - The width of the region, that is used to search
;                           for the emission lines given in REFWAVE(REFPIXL);
;                           given in pixels.
;         topwid          - If set, then error messages are displayed using
;                           DIALOG_MESSAGE, using this widget id as
;                           DIALOG_PARENT, instead of MESSAGE.
;         logunit         - Messages are saved to the file pointed to by this
;                           logical file unit, if it is defined.
;         verbose         - Show more information on what is being done.
;         error           - Returns an error code if set.
;         debug           - The error handler is not setup if debug is set.
;         help            - Show this routine documentation, and exit.
;
; OUTPUTS:
;         identline       - Returns an array with the indices of those lines in
;                           REFWAVE, which have been found in SPEC, based on
;                           the input array WAVELEN. The array is the section
;                           of WAVELEN which is completely bounded by REFWAVE.
;         identpixl       - Returns an array with the pixel number of
;                           the arrays returned in IDENTLINE.
;         identwave       - An array with the wavelength of the positively
;                           matched lines.
;
; COMMON BLOCKS:
;         none
;
; SIDE EFFECTS:
;         none
;
; RESTRICTIONS:
;         IDL version 6.2 or higher is required.
;
; MODIFICATION HISTORY:
;         15.10.2008 - Converted from original routine findlines of
;                      Thomas Becker. /CS
;-
PRO p3d_wavecal_match_maskwd,spec,wavelen,refpixl_,refwave_,fwhm_, $
        identline,identpixl,identwave,noise=noise,poldeg=poldeg, $
        niterate=niter,findwidth=fwidth,topwid=topwid,logunit=logunit, $
        verbose=verbose,error=error,debug=debug,help=help
  compile_opt hidden,IDL2

  if !version.release lt 6.2 then message,'IDL Version <6.2. Cannot continue.'
  error=0 & rname='p3d_wavecal_match_maskwd: '
  if ~n_elements(verbose) then verbose=0
  debug=keyword_set(debug)

  if keyword_set(help) then begin
    doc_library,'p3d_wavecal_match_maskwd'
    return
  endif

  ;;========================================------------------------------
  ;; Setting up an error handler:

  if ~debug then begin
    catch,error_status
    if error_status ne 0L then begin
      p3d_misc_errors,error_status,rname=rname,topwid=topwid
      catch,/cancel
      error=-1
      return
    endif
  endif ;; ~debug

  ;;========================================------------------------------
  ;; Checking input arguments:

  s=size(spec)
  if s[0L] ne 1L or $
    (s[s[0L]+1L] ge 6L and s[s[0L]+1L] le 11L) then begin
    errmsg='SPEC [1] must be set; to a one-dimensional array of floating ' + $
           'point type.'
    goto,error_handler
  endif
  idx=where(spec ne 0.0,count)
  if ~count then begin
    errmsg='SPEC [1] must be set; to a one-dimensional array of floating ' + $
           'point type - with non-zero elements!'
  endif ;; ~count

  sb=size(wavelen)
  if sb[0L] ne 1L or $
    (sb[sb[0L]+1L] ge 6L and sb[sb[0L]+1L] le 11L) then begin
    errmsg='WAVELEN [2] must be set; to a one-dimensional array of deci' + $
           'mal type.'
    goto,error_handler
  endif
  n=sb[sb[0L]+2L]

  sb=size(refpixl_)
  if sb[0L] ne 1L or $
    (sb[sb[0L]+1L] ge 6L and sb[sb[0L]+1L] le 11L) then begin
    errmsg='REFPIXL [3] must be set; to a one-dimensional array of decimal' + $
           ' type.'
    goto,error_handler
  endif
  refpixl=refpixl_

  sw=size(refwave_)
  if sw[0L] ne 1L or $
    (sw[sw[0L]+1L] ge 6L and sw[sw[0L]+1L] le 11L) then begin
    errmsg='REFWAVE [4] must be set; to a one-dimensional array of decimal' + $
           ' type.'
    goto,error_handler
  endif
  refwave=refwave_

  if sb[sb[0L]+2L] ne sw[sw[0L]+2L] then begin
    errmsg='REFPIXL ['+strtrim(sb[sb[0L]+2L],2L)+'] and REFWAVE ['+ $
           strtrim(sw[sw[0L]+2L],2L)+'] must have the same number of elem' + $
           'ents.'
    goto,error_handler
  endif

  s=size(fwhm_)
  if s[s[0L]+2L] ne 1L or $
    (s[s[0L]+1L] ge 6L and s[s[0L]+1L] le 11L) then begin
    errmsg='FWHM [5] must be set; to a scalar integer.'
    goto,error_handler
  endif
  fwhm=long(fwhm_)

  if ~n_elements(fwidth) then fwidth=fwhm
  s=size(fwidth)
  if s[s[0L]+2L] ne 1L or $
    (s[s[0L]+1L] ge 6L and s[s[0L]+1L] le 11L) then begin
    errmsg='FINDWIDTH must be a scalar integer.'
    goto,error_handler
  endif
  fwidth=long(fwidth)

  if ~n_elements(poldeg) then poldeg=1L
  s=size(poldeg)
  if s[s[0L]+2L] ne 1L or $
    (s[s[0L]+1L] ge 4L and s[s[0L]+1L] le 11L) then begin
    errmsg='POLDEG must be a scalar integer.'
    goto,error_handler
  endif
  poldeg=poldeg<(n_elements(refwave)-1L)

  if ~n_elements(niter) then niter=1L
  s=size(niter)
  if s[s[0L]+2L] ne 1L or $
    (s[s[0L]+1L] ge 4L and s[s[0L]+1L] le 11L) then begin
    errmsg='NITERATE must be a scalar integer.'
    goto,error_handler
  endif

  if ~n_elements(noise) then begin
    noisestr=' (determined from the data, see below)'
  endif else begin
    sb=size(noise)
    if sb[sb[0L]+2L] ne 1L or $
      (sb[sb[0L]+1L] ge 6L and sb[sb[0L]+1L] le 11L) then begin
      errmsg='NOISE must be a decimal scalar.'
      goto,error_handler
    endif
    noisestr=string(noise,format='(f9.3)')
  endelse

  s=size(spec)

  ;;========================================------------------------------
  ;; Logging the performed operations:

  msg=['Matching lines in the mask with the spectrum image.', $
       '  Using the following parameters:', $
       '               fwhm='+string(  fwhm,format='(i9)'), $ 
       '          findwidth='+string(fwidth,format='(i9)'), $
       '              noise='+noisestr, $
       '  polynomial degree='+string(poldeg,format='(i9)'), $
       '           niterate='+string(niter,format='(i9)')]
  error=p3d_misc_logger(msg,logunit,rname=rname,verbose=verbose ge 1)

  ;;========================================------------------------------
  ;; At first the calibration line reference wavelength array is fitted with a
  ;; polynomial, and thereafter WAVELEN is trimmed to use only this range:

  s=s[1L]
  x=dindgen(s)
  identline_in=lindgen(n)

  for jj=1L,niter do begin

    ;; Printing information about the iteration:
    msg='Iteration '+strtrim(jj,2L)+'/'+strtrim(niter,2L)+'.'
    error=p3d_misc_logger(msg,logunit,loglevel=2,rname=rname, $
        verbose=verbose ge 2)

    params=poly_fit(refpixl,refwave,poldeg,/double,status=status)
    if status ne 0L then begin
      errmsg='[iteration '+strtrim(jj,2L)+'/'+strtrim(niter,2L)+'] POLY_FIT: '
      case status of
        1: errmsg+='REFPIXL is singular.'
        2: errmsg+='A small pivot element was used, accuracy is likely lost.'
        3: errmsg+='Undefined error estimate (NaN) was encountered.'
      endcase
      goto,error_handler
    endif ;; status ne 0L

    pixwave=dblarr(s) & for k=0L,poldeg do pixwave+=x^k*params[k]

    cutpos=where(wavelen ge min(pixwave) and wavelen lt max(pixwave),nswave)
    if ~nswave then begin
      errmsg='None of the elements of WAVELEN are within the boundaries of' + $
             ' REFWAVE ['+strtrim(min(pixwave),2L)+'-'+ $
                          strtrim(max(pixwave),2L)+'].'
      goto,error_handler
    endif ;; nswave
    swave=wavelen[cutpos]

    wpix=dblarr(nswave)
    for k=0L,nswave-1L do begin
      tmp=where(pixwave ge swave[k])
      wpix[k]=tmp[0L]
    endfor ;; k=0L,nswave-1L

    ;;========================================------------------------------
    ;; Logging the performed operations:

    if ~(jj-1L) then begin
      case poldeg of
        1: tmp='1st'
        2: tmp='2nd'
        3: tmp='3rd'
        else: tmp=strtrim(poldeg,2L)+'th'
      endcase
      msg=['Fitted a '+tmp+'-order polynomial using the ' + $
           'known wavelengths of the line mask as ordinates:']
      for i=0L,n_elements(refpixl)-1L do msg=[msg,'  '+string(i,refwave[i], $
          refpixl[i],format='(i4," ",f9.3," [pixel ",f9.3,"]")')]
      error=p3d_misc_logger(msg,logunit,rname=rname,verbose=verbose ge 2)
    endif ;; ~(jj-1L)

    ;;========================================------------------------------
    ;; Estimating a scalar value on the noise using regions in the input
    ;; spectrum where there are no lines:

    noline=lonarr(s)+1L
    for k=0L,nswave-1L do $
       noline[(wpix[k]-fwidth)>0L:(wpix[k]+fwidth)<(s-1L)]=0L
    nolpos=where(noline)

    if ~n_elements(noise) then begin
      pos=where(spec[nolpos] ne 0.0)
      noise=median(abs(spec[nolpos[pos]]-median(spec[nolpos[pos]])),/double)

      msg=['The noise is calculated to be: '+string(noise,format='(f9.3)')]
      error=p3d_misc_logger(msg,logunit,rname=rname,verbose=verbose ge 1)
    endif ;; ~n_elements(noise)

    ;;========================================------------------------------
    ;; Looping through all elements of the trimmed wavelength array in order
    ;; to examine if there is a line at any position:

    lines=dblarr(nswave)
    offst=dblarr(nswave)
    for k=0L,nswave-1L do begin

      ;;==============================--------------------
      ;; Finding the position of the maximum value of the spectrum in the
      ;; wavelength range of the current line:

      lindex=floor(wpix[k]-2L*fwidth)>   0L
      uindex= ceil(wpix[k]+2L*fwidth)<(s-1L)
      curve=spec[lindex:uindex]
      dummy=max(curve,pos)
      off=-2L*fwidth+pos

      ;; Recalculating the maximum if the offset is non-zero:
      if off ne 0L then begin
        lindex=floor(wpix[k]+off-2L*fwidth)>   0L
        uindex= ceil(wpix[k]+off+2L*fwidth)<(s-1L)
        curve=spec[lindex:uindex]
        dummy=max(curve,pos)
      endif ;; off ne 0L

      ;;==============================--------------------
      ;; Checking if the maximum is above 3 times the noise level => real peak:

      t1a=(wpix[k]+off-fwidth*3L/2L)>0L
      t1b=(wpix[k]+off-fwidth      )>0L
      t2a=(wpix[k]+off+fwidth      )<(s-1L)
      t2b=(wpix[k]+off+fwidth*3L/2L)<(s-1L)
      msp=~t1b?spec[t2a:t2b]:(t2a eq s-1L?spec[t1a:t1b]: $
                              [spec[t1a:t1b],spec[t2a:t2b]])
      tmp=median(msp,/double)
      if curve[pos]-tmp gt 3d0*noise and $
         lindex+pos gt 0L and lindex+pos lt s-1L then begin
        if (spec[lindex+pos] gt spec[lindex+pos-1L]) and $
           (spec[lindex+pos] gt spec[lindex+pos+1L]) then begin
          low=(pos-fwhm)>0L
          upp=(pos+fwhm)<(uindex-lindex-1L)
             sum=total(curve[low:upp])
          weight=total(curve[low:upp]*dindgen(upp-low+1L))
          lines[k]=weight/sum+lindex+low
          offst[k]=jj eq 1L?abs(lines[k]-wpix[k]):abs(lines[k]-plines[k])
        endif
      endif ;; curve[pos]-tmp gt 3d0*noise
    endfor ;; k=0L,nswave-1L

    ;;==============================--------------------
    ;; Removing entries of different wavelength at the same pixel,
    ;; keeping the entry with the smallest offset:

    pidx=where(lines ge 1d0,count,complement=nidentline,ncomplement=ncount)
    pmsk=lonarr(nswave)+1L & if ncount ne 0L then pmsk[nidentline]=0L
    for k=0L,count-1L do begin
      if ~pmsk[pidx[k]] then continue
      idx=where(lines[pidx] eq lines[pidx[k]],count)
      if count ge 2L then begin
        tmp1=min(offst[pidx[idx]],pos)
        tmp2=where(idx ne idx[pos])
        lines[pidx[idx[tmp2]]]=0d0
        pmsk[pidx[idx]]=0L
      endif ;; count ge 2L
    endfor ;; k=0L,nswave-1L

    ;; Iterating the solution:
    idx=where(lines ge 1d0,count,complement=nidentline,ncomplement=ncount)
    mask=lonarr(n)
    if count ne 0L then begin
      mask[cutpos[idx]]=1L
      identline=cutpos[idx]
      refpixl=lines[idx]
      refwave=swave[idx]
    endif ;; count ne 0L

    plines=lines
  endfor ;; jj=1L,niter

  ;;========================================------------------------------
  ;; An emission line is present where LINES>=1.0:

  if ~count then begin
    msg='Could not match any entries in the line mask with the input spectrum.'
  endif else begin
    identpixl=refpixl
    identwave=refwave

    ;; Logging the performed operations:
    msg=['Matched the following '+strtrim(count,2L)+' entries of the line ' + $
         'mask using the input spectrum [wavelength, pixel position]:']
    for i=0L,count-1L do msg=[msg,'  '+string(i+1L,identwave[i], $
       identpixl[i],format='(i4," ",f9.3," [",f9.3,"]")')]

    idx=where(~mask,ncount)
    if ncount ne 0L then begin
      nidentwave=refwave_[idx]
      msg=[msg,'Discarded the following '+strtrim(ncount,2L)+' entries of' + $
           ' the lines mask using the input spectrum [wavelength, pixel p' + $
           'osition]:']
      for i=0L,ncount-1L do msg=[msg,'  '+string(count+i+1L,nidentwave[i], $
         format='(i4," ",f9.3)')]
    endif ;; ncount ne 0L

  endelse ;; ~count 
  error=p3d_misc_logger(msg,logunit,rname=rname,verbose=verbose ge 1)

  return

error_handler:
  error=p3d_misc_logger(errmsg,logunit,rname=rname,topwid=topwid, $
      verbose=verbose,/error)
  return
END ;;; procedure: p3d_wavecal_match_maskwd