;+
; NAME:
;         p3d_wavecal_fit_maskwd
;
;         $Id: p3d_wavecal_fit_maskwd.pro 140 2010-03-28 17:51:09Z christersandin $
;
; PURPOSE:
;         This is the routine where the actual wavelength calibration
;         takes place. The wavelengths in the mask are fitted with a
;         polynomial of order DISPDEG against the corrected positions of the
;         lines in the spectrum image, that is assumed to contain
;         these lines.
;
;         Additional parameters allow the procedure to be modified:
;         If REFDIST is set to any integer above 1 then the fitting is made
;         with REFDIST as the stride. I.e., with REFDIST==2 every second
;         spectrum is fitted. If NROWS>0 then the 2*NROWS+1 most adjacent
;         spectra are used in the fitting of every spectrum (boundaries are
;         different). If CROSSDEG>=0 then the fitted parameters are fitted once
;         again, across all spectra for every order in the fitted parameters
;         separately (use is not recommended).
;
;         In order to check the quality of the fit the maximum residual (or all
;         residuals if LOGLEVEL==2) is output to a logfile or the screen (if
;         VERBOSE is set). Additionally the Chi-square value and the status of
;         POLY_FIT are also shown.
;
;         The polynomial parameters are provided for every spectrum of the
;         image in an array upon completion.
;
; 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.
;
; AUTHOR:
;         Christer Sandin
;         Astrophysikalisches Institut Potsdam (AIP)
;         An der Sternwarte 16
;         D-14482 Potsdam, GERMANY
;
; CATEGORY:
;         p3d :: wavelength calibration
;
; CALLING SEQUENCE:
;         p3d_wavecal_fit_maskwd,image,nlines,linepos,wavelength,refdist, $
;             dispdeg,crossdeg,out,dispmid,residualcut=,chisq=,maxspecnum=, $
;             refrow=,deadfibers=,sdeadfibers=,topwid=,logunit=, $
;             verbose=,error=,/debug,/help
;
; INPUTS:
;         image           - A two-dimensional image with emission lines of an
;                           arc (calibration) lamp.
;         nlines          - A scalar integer specifying the number of
;                           calibration lines.
;         linepos         - A two-dimensional array holding the positions of
;                           every matched calibration line for all spectra;
;                           given in pixels.
;         wavelength      - A one-dimensional array with the corresponding
;                           wavelength of every emission line specified in
;                           LINEPOS.
;         refdist         - A scalar integer specifying if any spectrum is to
;                           be skipped. If REFDIST==1 1 then every spectrum is
;                           used, if REFDIST==2 then every second, etc.
;         dispdeg         - A scalar integer specifying the order of the
;                           polynomial that is used to fit every spectrum on
;                           the dispersion axis.
;         crossdeg [-1]   - A scalar integer specifying the order of the
;                           polynomial that is used to fit the separate
;                           parameters of the dispersion-direction fit.
;
; KEYWORD PARAMETERS:
;         residualcut [-1]- Residuals are calculated between the polynomial fit
;                           of the wavelength function and WAVELENGTH. If
;                           RESIDUALCUT is set to a positive value then the
;                           weights are set to zeros where the residual is
;                           larger than RESIDUALCUT.
;         chisq           - A one-dimensional array that returns the chi²
;                           values of the fits for every spectrum.
;         maxspecnum      - Returns the index of the spectrum which has the
;                           largest residual.
;         refrow [s[2]/2] - The starting spectrum.
;         deadfibers      - A one-dimensional array of integers specifying
;                           which fiber positions are to be interpolated
;                           instead of fitted.
;         sdeadfibers     - A string array with as many elements as DEADFIBERS
;                           specifying the type of fiber in DEADFIBER (dead,
;                           low transmission, unused,...). This array is used
;                           in order to only skip [D]ead or [U]nused fibers.
;         stawid          - If set to a valid ID then a log message is written
;                           using this ID for relevant actions.
;         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:
;         out             - A two-dimensional array that contains the
;                           parameters of the polynomial fit [2nd dim.] for
;                           every spectrum in IMAGE.
;         dispmid         - .
;
; COMMON BLOCKS:
;         none
;
; SIDE EFFECTS:
;         none
;
; RESTRICTIONS:
;         IDL version 6.2 or higher is required.
;
; MODIFICATION HISTORY:
;         15.10.2008 - Converted from original routine wavepos of
;                      Thomas Becker. /CS
;         29.06.2009 - Added full documentation. /CS
;-
PRO p3d_wavecal_fit_maskwd,image,nlines,linepos,wavelength,refdist,dispdeg, $
        crossdeg,out,dispmid,residualcut=residualcut,chisq=chisq, $
        maxspecnum=maxspecnum,refrow=refrow,deadfibers=deadfibers_, $
        sdeadfibers=sdeadfibers,stawid=stawid,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_fit_maskwd: '
  if ~n_elements(verbose) then verbose=0
  debug=keyword_set(debug)
  usestawid=~n_elements(stawid)?0L:widget_info(stawid,/valid_id)

  if keyword_set(help) then begin
    doc_library,'p3d_wavecal_fit_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 the input arguments:

  s=size(image)
  if s[0L] ne 2L or $
    (s[s[0L]+1L] ge 6L and s[s[0L]+1L] le 11L) then begin
    errmsg='IMAGE [1] must be set; to a two-dimensional array of decimal type.'
    goto,error_handler
  endif

  sr=size(nlines)
  if sr[sr[0L]+2L] ne 1L or $
    (sr[sr[0L]+1L] ge 4L and sr[sr[0L]+1L] le 11L) then begin
    errmsg='NLINES must be set to a scalar integer.'
    goto,error_handler
  endif
  n=nlines

  sl=size(linepos)
  if sl[0L] ne 2L or $
    (sl[sl[0L]+1L] ge 6L and sl[sl[0L]+1L] le 11L) then begin
    errmsg='LINEPOS [2] must be set; to a two-dimensional array of decimal' + $
           ' type.'
    goto,error_handler
  endif

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

  if sl[1L] ne sw[sw[0L]+2L] then begin
    errmsg='LINEPOS [# wavelength bins='+strtrim(sl[1L],2L)+'] and WAVELEN' + $
           'GTH ['+strtrim(sw[sw[0L]+2L],2L)+'] must have the same number ' + $
           'of elements.'
    goto,error_handler
  endif

  sr=size(refdist)
  if sr[sr[0L]+2L] ne 1L or $
    (sr[sr[0L]+1L] ge 4L and sr[sr[0L]+1L] le 11L) then begin
    errmsg='REFDIST [4] must be set to a scalar integer; >=0.'
    goto,error_handler
  endif
  if refdist lt 0L then begin
    errmsg='REFDIST [4] must be set to a scalar integer; >=0.'
    goto,error_handler
  endif

  sr=size(dispdeg)
  if sr[sr[0L]+2L] ne 1L or $
    (sr[sr[0L]+1L] ge 4L and sr[sr[0L]+1L] le 11L) then begin
    errmsg='DISPDEG [6] must be set; to a scalar integer.'
    goto,error_handler
  endif
  ldispdeg=dispdeg<(n_elements(lwavelength)-1L)

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

  sr=size(crossdeg)
  if sr[sr[0L]+2L] ne 1L or $
    (sr[sr[0L]+1L] ge 4L and sr[sr[0L]+1L] le 11L) then begin
    errmsg='CROSSDEG [7] must be set; to a scalar integer.'
    goto,error_handler
  endif

  if ~n_elements(residualcut) then residualcut=-1L
  sr=size(residualcut)
  if sr[sr[0L]+2L] ne 1L or $
    (sr[sr[0L]+1L] ge 6L and sr[sr[0L]+1L] le 11L) then begin
    errmsg='RESIDUALCUT, if set, must be a decimal scalar; RESIDUALCUT>0.0.'
    goto,error_handler
  endif
  if residualcut ne -1L and residualcut lt 0d0 then begin
    errmsg='RESIDUALCUT, if set, must be a decimal scalar; RESIDUALCUT>0.0.'
    goto,error_handler
  endif

  farr=lindgen(s[2L])
  usedeadfibers=0L
  spc=size(deadfibers_)
  dfmask=lonarr(s[2L])+1L
  if spc[spc[0L]+2L] ge 1L then begin
    if spc[spc[0L]+1L] ge 4L and spc[spc[0L]+1L] le 11L then begin
      errmsg='DEADFIBERS must, if set, be an array of integer type.'
      goto,error_handler
    endif
    if min(deadfibers_) lt 0L or max(deadfibers_) gt s[2L] then begin
      errmsg='DEADFIBERS must, if set, be an array of integer type.'
      goto,error_handler
    endif
    usedeadfibers=1L
    deadfibers=deadfibers_-1L

    spc=size(sdeadfibers)
    if spc[spc[0L]+2L] ne n_elements(deadfibers) or $
       spc[spc[0L]+1L] ne 7L then begin
      errmsg='SDEADFIBERS must be an array of string type with as' + $
             'many elements as DEADFIBERS.'
       goto,error_handler
    endif

    count=0L
    for i=0L,n_elements(deadfibers)-1L do begin
      tmp=where(farr eq deadfibers[i],count_)
      if count_ eq 1L and $
         (strlowcase(strmid(sdeadfibers[i],0L,1L)) eq 'd' or $
          strlowcase(strmid(sdeadfibers[i],0L,1L)) eq 'u') then begin
        dfmask[tmp[0L]]=0L
        count++
      endif
    endfor ;; i=0L,n_elements(deadfibers)
    if count eq s[2L] then begin
      errmsg='All spectra are masked as dead, cannot continue.'
      goto,error_handler
    endif
;   if count gt 0L then farr=farr[where(dfmask)]
  endif ;;spc[spc[0L]+2L] ge 1L

  ;;========================================------------------------------
  ;; Logging the input parameters:

  msg='Creating a dispersion mask; fitting the input spectra with polynomials.'
  msg=[msg, $
       'Using the following parameters:', $
       '    refdist='+string(refdist,format='(i9)')+', specifying the stri' + $
       'de of fitted spectra.', $
       '    dispdeg='+string(dispdeg,format='(i9)')+', dispdeg is the orde' + $
       'r of the polynomial that is used on the dispersion axis.', $
       '   crossdeg='+string(crossdeg,format='(i9)')+', crossdeg is the or' + $
       'der of the polynomial that is used on the cross-dispersion axis.', $
       'residualcut='+string(residualcut,format='(f9.3)')+', residuals lar' + $
       'ger than residualcut are ignored in the second fitting attempt.']
  error=p3d_misc_logger(msg,logunit,rname=rname,verbose=verbose ge 1)

;  ;;========================================------------------------------
;  ;; Creating a reference image, with the dispersion direction in dimension 1,
;  ;; which contains all spectra (REFDIST==1,NROWS=0), or every REFDIST
;  ;; spectrum, or an average f(NROWS):
;
  s=size(image)
  lowsize=refrow/refdist
  uppsize=(s[2L]-1L-refrow)/refdist
  refsize=lowsize+uppsize+1L

  if usestawid then begin
    msg='Fitting a polynomial Wavelength(arc line position) for every spec' + $
        'trum.'
    widget_control,stawid,set_value=msg
  endif

  ;;========================================------------------------------
  ;; Looping through all spectra in REFIMAGE. Fitting a polynomial
  ;; of order DISPDEG to the function WAVELENGTH(IMLINEPOS):

  weight  =lonarr(n,refsize)+1L ;; The weight for every line is 1 per default
  params  =dblarr(refsize,ldispdeg+1L)
  chisq   =dblarr(refsize)
  status  =lonarr(refsize)
  residual=dblarr(n,refsize)
  mresidual=dblarr(refsize)

  for k=0L,refsize-1L do begin
    params[k,*]=poly_fit(linepos[*,k],lwavelength,ldispdeg,chisq=chisq_, $
        status=status_,/double)
    chisq[k]=chisq_ & status[k]=status_

    ;; Assembling the polynomial function:
    curve=dblarr(n)
    for m=0L,ldispdeg do curve+=double(linepos[*,k])^m*params[k,m]
    ;; Calculating the residual of the wavelength polynomial from
    ;; the input data:
    residual[*,k]=lwavelength-curve
    mresidual[k]=max(abs(residual[*,k]))

    ;; Assigning a 0 weight to those wavelength bins where the residual
    ;; is larger than RESIDUALCUT:
    if residualcut ne -1L then weight[*,k]=abs(residual[*,k]) lt residualcut

  endfor ;; k=0L,refsize-1L
  tmp=max(mresidual*dfmask,maxspecnum)

  ;;========================================------------------------------
  ;; Logging the input parameters:

  stride=10L & ttstr=string(replicate(32b,40L))
  xrows=(n-1L)/stride & xrmod=n mod stride & if xrmod eq 0L then $
     xrmod=stride-1L
  xrfst=(n<stride)-1L

  tmp=logunit[1L] lt 2L?'':' :: residuals of individual lines'
  msg=['Results of the fitting procedure (status==0 [poly_fit] is O.K.):', $
       '   Spec.  max(res)     chisq  status'+tmp]
  if logunit[1L] ge 2L then begin
    msg=[msg,'                   (Wavelengths [Å])    '+ $
         strjoin(string(lwavelength[0L:xrfst],format='(f8.2)'),' ')]
    for i=1L,xrows do begin
      x1=i*stride
      x2=x1+(((i eq xrows)?xrmod:stride)-1L)
      msg=[msg,'                                        '+ $
           strjoin(string(lwavelength[x1:x2],format='(f8.2)'),' ')]
    endfor
  endif ;; logunit[1L] ge 2L

  for k=0L,refsize-1L do begin
    usetmpx=0L
    if logunit[1L] lt 2L then begin
      tmp=''
    endif else begin
      tmp=' :: '+strjoin(string(residual[0L:xrfst,k],format='(e8.1)'),' ')
      if xrows ge 1L then begin
        tmpx=strarr(xrows) & usetmpx=1L
        for i=1L,xrows do begin
          x1=i*stride
          x2=x1+(((i eq xrows)?xrmod:stride)-1L)
          tmpx[i-1L]=ttstr+strjoin(string(residual[x1:x2,k], $
                                          format='(e8.1)'),' ')
        endfor ;; i=1L,xrows
      endif     ;; xrows ge 1L
    endelse      ;; logunit[1L] lt 2L
    maxrc=max(abs(residual[*,k])) lt residualcut?' ':'*'
    msg=[msg,'  '+string(k+1L,format='(i6)')+'  '+ $
         string(max(abs(residual[*,k])),maxrc,format='(e8.2,a1)')+'  '+ $
         string(chisq[k],status[k],format='(e7.1,2x,i6)')+tmp]
    if usetmpx then for i=0L,xrows-1L do msg=[msg,tmpx[i]]
  endfor
  error=p3d_misc_logger(msg,logunit,rname=rname,verbose=verbose ge 1)

  ;;========================================------------------------------
  ;; There is one weight value per emission line. It is set to 1 if the
  ;; emission line was identified through at least half of the spectra:

  weight=total(weight,2L) ge refsize/2L
  ok=where(weight,count)

  ;;========================================------------------------------ 
  ;; Re-doing the fitting, taking the now known weights into account:

  if count lt n then begin
    ldispdeg=ldispdeg<(count-1L)
    params  =dblarr(refsize,ldispdeg+1L)
    residual=dblarr(count,refsize)

    ;;========================================------------------------------
    ;; Looping through all spectra in REFIMAGE. Fitting a polynomial
    ;; of order DISPDEG to the function WAVELENGTH(LINEPOS):

    for k=0L,refsize-1L do begin
      params[k,*]=poly_fit(linepos[ok,k], $
          lwavelength[ok],ldispdeg,chisq=chisq_,status=status_,/double)
      chisq[k]=chisq_ & status[k]=status_

      ;; Assembling the polynomial function:
      curve=dblarr(count)
      for m=0L,ldispdeg do curve+=double(linepos[*,k])^m*params[k,m]

      ;; Calculating the residual of the wavelength polynomial from
      ;; the input data:
      residual[*,k]=lwavelength-curve
      mresidual[k]=max(abs(residual[*,k]))
    endfor ;; k=0L,refsize-1L
    tmp=max(mresidual*dfmask,maxspecnum)

    ;;========================================------------------------------
    ;; Logging the input parameters:

    xrows=(count-1L)/stride & xrmod=count mod stride
    if xrmod eq 0L then xrmod=stride-1L
    xrfst=(count<stride)-1L

    tmp=logunit[1L] lt 2L?'':' :: residuals of individual lines'
    msg=['All lines could not be fitted with small residuals (<residualc' + $
         'ut). Procedure redone.', $
         '   Spec.  max(res)     chisq  status'+tmp]
    if logunit[1L] ge 2L then begin
      msg=[msg,'                   (Wavelengths [Å])    '+ $
           strjoin(string(lwavelength[0L:xrfst],format='(f8.2)'),' ')]
      for i=1L,xrows do begin
        x1=i*stride
        x2=x1+(((i eq xrows)?xrmod:stride)-1L)
        msg=[msg,'                                        '+ $
             strjoin(string(lwavelength[x1:x2],format='(f8.2)'),' ')]
      endfor
    endif ;; logunit[1L] ge 2L

    for k=0L,refsize-1L do begin
      usetmpx=0L
      if logunit[1L] lt 2L then begin
        tmp=''
      endif else begin
        tmp=' :: '+strjoin(string(residual[0L:xrfst,k],format='(e8.1)'),' ')
        if xrows ge 1L then begin
          tmpx=strarr(xrows) & usetmpx=1L
          for i=1L,xrows do begin
            x1=i*stride
            x2=x1+(((i eq xrows)?xrmod:stride)-1L)
            tmpx[i-1L]=ttstr+strjoin(string(residual[x1:x2,k], $
                                            format='(e8.1)'),' ')
          endfor
        endif
      endelse
      maxrc=max(abs(residual[*,k])) lt residualcut?' ':'*'
      msg=[msg,'  '+string(k+1L,format='(i6)')+'  '+ $
           string(max(abs(residual[*,k])),maxrc,format='(e8.2,a1)')+'  '+ $
           string(chisq[k],status[k],format='(e7.1,2x,i6)')+tmp]
      if usetmpx then for i=0L,xrows-1L do msg=[msg,tmpx[i]]
    endfor
    error=p3d_misc_logger(msg,logunit,rname=rname,verbose=verbose ge 1)

  endif ;; count lt n

  out=dblarr(refsize,dispdeg+1L)
  out[*,0L:ldispdeg]=params

  dispmid=mean(params[*,1L],/double)

  ;;========================================------------------------------
  ;; Making an additional polynomial fit across separate spectra,
  ;; if crossdeg>=0: This fitting is made separately for every polynomial
  ;; parameter of the fitting on the dispersion axis:

  if crossdeg ge 0L then begin
    fit_params=dblarr(crossdeg+1L,dispdeg+1L)
    for k=0L,dispdeg do $
       fit_params[*,k]=poly_fit(refpos,params[*,k],crossdeg,/double)

    par=dblarr(s[2L],ldispdeg+1L)
    for k=0L,dispdeg do for L=0L,crossdeg do $
       par[*,k]+=fit_params[L,k]*dindgen(s[2L])^L
    out=par
  endif ;; crossdeg ge 0L

  return

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