;+
; NAME:
;         p3d_tracing_calculate_lprofs_fits
;
;         $Id: p3d_tracing_calculate_lprofs_fits.pro 151 2010-04-16 08:02:14Z christersandin $
;
; PURPOSE:
;         Fit a cross-dispersion cut with a sum of Gaussians using MPFIT.
;
;         If the keyword FULL is specified then the profile center positions
;         are fitted as well.
;
; AUTHOR:
;         Christer Sandin and Joris Gerssen
;         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 :: tracing of spectra on the CCD
;
; CALLING SEQUENCE:    
;         p3d_tracing_calculate_lprofs_fits,x,y,dy,traces,par,out,dpar, $
;             mpmaxiter=,proffun=,/full,foffset=,topwid=,logunit=,error=, $
;             verbose=,/debug,/help
;
; INPUTS:   
;         x             - A one-dimensional array with the abscissae of Y.
;         y             - A one-dimensional array with the data to be fitted.
;         dy            - A one-dimensional array with the variance of Y.
;         traces        - A one-dimensional array with the traces along Y.
;         par           - A one-dimensional array with fitting initial values.
;                         Overwritten with the fitted parameter values upon
;                         return.
;
; KEYWORDS: 
;         mpmaxiter     - A scalar integer specifying the maximum number of
;                         allowed iterations for MPFIT.
;         proffun       - A scalar string with the name of the function to
;                         use when (re-)calculating the line profile.
;         full          - If set, then the profile center positions are fitted,
;                         instead of using the fixed trace mask positions. The
;                         allowed offset in this case is specified by FOFFSET.
;         foffset [0.5] - If FULL is set then this value specifies the maximum
;                         offset of the fitted line center position.
;         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.
;         error         - Returns an error code if set
;         verbose       - Show more information on what is being done
;         debug         - The error handler is not setup if debug is set.
;         help          - Show this routine documentation, and exit
;
; OUTPUTS:
;         out           - A one-dimensional array with the resulting fitted
;                         function of Y.
;         dpar          - The error of PAR (if the fit was good, chisq==1).
;
; COMMON BLOCKS:
;         none
;
; SIDE EFFECTS:
;         none
;
; RESTRICTIONS:
;         IDL version 6.2 or higher is required.
;
;-
PRO p3d_tracing_calculate_lprofs_fits_ngaussmodel,x,a,fx, $
        x0=x0,fun=proffun,full=full
  compile_opt hidden,IDL2

  nx=n_elements(x) & ng=n_elements(x0)

  ;; The function is given in a compact form without using loops:
  case proffun of
    'lorentzian': begin
      ;; The Lorentzian function is defined as:
      ;;     fx=norm/!DPI*0.5*FWHM/((x-x0)^2+(0.5*FWHM)^2)
      ;;
      ;;     a[0]  :: zero-point
      ;;     a[1]  :: slope
      ;;     a[2]  :: FWHM
      ;;     a[3   :  NG+2] :: normalization
      ;;     a[3+NG:2*NG+2] :: line center positions (/FULL)

      if keyword_set(full) then begin
        tmp1=((rebin(x,nx,ng)-rebin(transpose(a[(3L+ng):(2L*ng+2L)]),nx,ng))/ $
              (0.5d0*a[2L]))^2
      endif else begin
        tmp1=((rebin(x,nx,ng)-rebin(transpose(x0),nx,ng))/(0.5d0*a[2L]))^2
      endelse ;; keyword_set(full)

      tmp2=rebin(transpose(a[3L:2L+ng]),nx,ng)/(0.5d0*a[2L]*!DPI)
      fx=total(tmp2/(tmp1+1),2L)+a[0L]+a[1L]*x
    end ;; 'lorentzian'
    'gauss/lorentz': begin
      ;; The Gaussian/Lorentzian function is defined as: 
      ;;   fx=norm*exp((m-1)*fac)/(1+m/ln2*fac), where fac=0.5*((x-x0)/sigma)^2
      ;;
      ;;     a[0L]  :: zero-point
      ;;     a[1L]  :: slope
      ;;     a[2L]  :: sigma
      ;;     a[3L]  :: m
      ;;     a[4   :  NG+3] :: normalization
      ;;     a[4+NG:2*NG+3] :: line center positions (/FULL)

      if keyword_set(full) then begin
        tmp1=0.5d0*((rebin(x,nx,ng)- $
                     rebin(transpose(a[(4L+ng):(2L*ng+3L)]),nx,ng))/a[2L])^2
      endif else begin
        tmp1=0.5d0*((rebin(x,nx,ng)-rebin(transpose(x0),nx,ng))/a[2L])^2
      endelse ;; keyword_set(full)

      tmp2=rebin(transpose(a[4L:3L+ng]),nx,ng)
      fx=total(tmp2*exp(-tmp1)/(1d0+a[3L]*tmp1*1.44269504088896d0),2L)+ $
             a[0L]+a[1L]*x
    end ;; 'gauss/lorentz'
    'doublegauss': begin
      ;; The Double Gaussian function is defined as:
      ;;   fx=norm/(sqrt(2*!DPI))*
      ;;        (exp(-0.5*fac1^2)/sigma1+A*exp(-0.5*fac2^2)/sigma2),
      ;;     where fac1=(x-x0)/sigma1, fac2=(x-x0)/sigma2
      ;;
      ;;     a[0]  :: zero-point
      ;;     a[1]  :: slope
      ;;     a[2]  :: sigma1
      ;;     a[3]  :: sigma2
      ;;     a[4]  :: normalization-constant A
      ;;     a[5   :  NG+4] :: normalization
      ;;     a[5+NG:2*NG+4] :: line center positions (/FULL)

      if keyword_set(full) then begin
        tmp=rebin(x,nx,ng)-rebin(transpose(a[(5L+ng):(2L*ng+4L)]),nx,ng)
      endif else begin
        tmp=rebin(x,nx,ng)-rebin(transpose(x0),nx,ng)
      endelse ;; keyword_set(full)

      tmp1a=-0.5d0*(tmp/a[2L])^2
      tmp1b=-0.5d0*(temporary(tmp)/a[3L])^2
      tmp2=rebin(transpose(a[5L:4L+ng]),nx,ng)/sqrt(2d0*!DPI)
      fx=total(tmp2*(exp(tmp1a)/a[2L]+a[4L]*exp(tmp1b)/a[3L]),2L)+a[0L]+a[1L]*x
    end ;; 'doublegauss'
    else: begin
      ;; The Gaussian function is defined as:
      ;;   fx=norm/(sigma*sqrt(2*!DPI))*exp(-0.5*fac^2), where fac=(x-x0)/sigma
      ;;
      ;;     a[0]  :: zero-point
      ;;     a[1]  :: slope
      ;;     a[2]  :: sigma 
      ;;     a[3   :  NG+2] :: normalization
      ;;     a[3+NG:2*NG+2] :: normalization  

      if keyword_set(full) then begin
        tmp1=-0.5d0*((rebin(x,nx,ng)- $
                      rebin(transpose(a[(3L+ng):(2L*ng+2L)]),nx,ng))/a[2L])^2
      endif else begin 
        tmp1=-0.5d0*((rebin(x,nx,ng)-rebin(transpose(x0),nx,ng))/a[2L])^2
      endelse ;; keyword_set(full)

      tmp2=rebin(transpose(a[3L:2L+ng]),nx,ng)/(a[2L]*sqrt(2d0*!DPI))
      fx=total(tmp2*exp(tmp1),2L)+a[0L]+a[1L]*x
    end ;; else
  endcase ;; proffun

  return
end ;; procedure: p3d_tracing_calculate_lprofs_ngaussmodel


PRO p3d_tracing_calculate_lprofs_fits,x,y,dy,traces,par,out,dout,fixw=fixw, $
        show=show,mpmaxiter=mpmaxiter,proffun=proffun,full=full, $
        foffset=foffset,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_tracing_calculate_lprofs_fits: '
  if ~n_elements(verbose) then verbose=0
  debug=keyword_set(debug)

  if keyword_set(help) then begin
    doc_library,'p3d_tracing_calculate_lprofs_fits'
    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

  ;;========================================------------------------------
  ;; Number of parameters in the model:

  fixw=keyword_set(fixw)
  full=keyword_set(full)
  if full then if ~n_elements(foffset) then foffset=0.5d0

  ntr=n_elements(traces)
  case proffun of
    'lorentzian': begin
      npar=3L+ntr*(full?2L:1L)

      ;; Constrain parameters (no parameter can be <0):
      parinfo=replicate({fixed:0L,limited:[1L,0L],limits:[0d0,1d0]},npar)

      ;; 0      - zero-level
      ;; 1      - slope
      ;; 2      - FWHM
      ;; 3:ng+2      - intensities
      ;; 3+ng:2*ng+2 - line center (/FULL)

      ;; The slope is not limited:
      parinfo[1L].limited=[0L,0L]

      ;; Limit the range of FWHM (in pixels):
      if fixw then begin
        parinfo[2L].fixed=1L
      endif else begin
        parinfo[2L].limited=[1L,1L]
        parinfo[2L].limits=[0.8*par[2L],1.2*par[2L]]
      endelse

      ;; Limit the range of the line center positions (in pixels; /FULL):
      if full then begin
        parinfo[(3L+ntr):(2L*ntr+2L)].limited=[1L,1L]
        tmp=transpose(par[(3L+ntr):(2L*ntr+2L)])
        parinfo[(3L+ntr):(2L*ntr+2L)].limits=[tmp-foffset,tmp+foffset]
      endif ;; full

    end
    'gauss/lorentz': begin
      npar=4L+ntr*(full?2L:1L)

      ;; Constrain parameters (no parameter can be <0):
      parinfo=replicate({fixed:0L,limited:[1L,0L],limits:[0d0,1d0]},npar)

      ;; 0  - zero-level
      ;; 1  - slope
      ;; 2  - sigma
      ;; 3  - m
      ;; 4:ng+3      - intensities
      ;; 4+ng:2*ng+3 - line center (/FULL)

      ;; The slope is not limited:
      parinfo[1L].limited=[0L,0L]

      if fixw then begin
        parinfo[2L].fixed=1L
        parinfo[3L].fixed=1L
      endif else begin
        ;; Limit the range of sigma (in pixels):
        parinfo[2L].limited=[1L,1L]
        parinfo[2L].limits=[0.8*par[2L],1.2*par[2L]]

        ;; Limit the range of m:
        parinfo[3L].limited=[1L,1L]
        parinfo[3L].limits=[0d0,1d0]
      endelse

      ;; Limit the range of the line center positions (in pixels; /FULL):
      if full then begin
        parinfo[(4L+ntr):(2L*ntr+3L)].limited=[1L,1L]
        tmp=transpose(par[(4L+ntr):(2L*ntr+3L)])
        parinfo[(4L+ntr):(2L*ntr+3L)].limits=[tmp-foffset,tmp+foffset]
      endif ;; full

    end
    'doublegauss': begin
      npar=5L+ntr*(full?2L:1L)

      ;; Constrain parameters (no parameter can be <0):
      parinfo=replicate({fixed:0L,limited:[1L,0L],limits:[0d0,1d0]},npar)

      ;; 0  - zero-level
      ;; 1  - slope
      ;; 2  - sigma1
      ;; 3  - sigma2
      ;; 4  - A
      ;; 5:ng+4      - intensities
      ;; 5+ng:2*ng+4 - line center (/FULL)

      ;; The slope is not limited:
      parinfo[1L].limited=[0L,0L]

      ;; Limit the range of sigma1 (in pixels):
      if fixw then begin
        parinfo[2L].fixed=1L
        parinfo[3L].fixed=1L
        parinfo[4L].fixed=1L
      endif else begin
        parinfo[2L].limited=[1L,1L]
        parinfo[2L].limits=[0.7*par[2L],1.3*par[2L]]

        ;; Limit the range of sigma2 (in pixels):
        parinfo[3L].limited=[1L,1L]
        parinfo[3L].limits=parinfo[2L].limits*[1.5d0,5d0]

        ;; Limit the range of A:
        parinfo[4L].limited=[1L,1L]
        parinfo[4L].limits=[0d0,1d0]
      endelse ;; fixw

      ;; Limit the range of the line center positions (in pixels; /FULL):
      if full then begin
        parinfo[(5L+ntr):(2L*ntr+4L)].limited=[1L,1L]
        tmp=transpose(par[(5L+ntr):(2L*ntr+4L)])
        parinfo[(4L+ntr):(2L*ntr+4L)].limits=[tmp-foffset,tmp+foffset]
      endif ;; full

    end
    else: begin ;; 'gaussian'
      npar=3L+ntr*(full?2L:1L)

      ;; Constrain parameters (no parameter can be <0):
      parinfo=replicate({fixed:0L,limited:[1L,0L],limits:[0d0,1d0]},npar)

      ;; 0  - zero-level
      ;; 1  - slope
      ;; 2  - Gauss sigma
      ;; 3:ng+2      - intensities
      ;; 3+ng:2*ng+2 - line center (/FULL)

      ;; The slope is not limited:
      parinfo[1L].limited=[0L,0L]

      ;; Limit the range of sigma (in pixels):
      if fixw then begin
        parinfo[2L].fixed=1L
      endif else begin
        parinfo[2L].limited=[1L,1L]
;       parinfo[2L].limits=[0.8*par[2L],1.2*par[2L]] ;; works better with VIRUS
        parinfo[2L].limits=[0.7*par[2L],1.3*par[2L]]
      endelse ;; fixw

      ;; Limit the range of the line center positions (in pixels; /FULL):
      if full then begin
        parinfo[(3L+ntr):(2L*ntr+2L)].limited=[1L,1L]
        tmp=transpose(par[(3L+ntr):(2L*ntr+2L)])
        parinfo[(3L+ntr):(2L*ntr+2L)].limits=[tmp-foffset,tmp+foffset]
      endif ;; full
    end
  endcase ;; proffun

  ;; Calling the fitting routine MPCURVEFIT:
  out=mpcurvefit(x,y,1/dy,par,sigpar,/autoderivative,chisq=chisq, $
          dof=dof,errmsg=errmsg,functargs={x0:traces,fun:proffun,full:full}, $
          function_name='p3d_tracing_calculate_lprofs_fits_ngaussmodel', $
          iter=niter,itmax=mpmaxiter,nfev=nfev,parinfo=parinfo, $
          quiet=verbose lt 3,status=status)

  if status le 0L then begin
    msg='Failure when executing MPCURVEFIT.'
    errmsg=[msg,errmsg]
    goto,error_handler
  endif

  dpar=sigpar*sqrt(chisq/(n_elements(par)))

  if ~n_elements(nfev) then nfev=-1L
  if ~n_elements(dof) then dof=-1L

  ;; Printing output information:
  msg=     '    MPFIT status:       '+string(status,format='(i5)')
  msg=[msg,'    Number of iterations and function evaluations: '+ $
       string(niter,nfev,format='(i5,", ",i5)')]
  error=p3d_misc_logger(msg,logunit,loglevel=3,rname=rname,topwid=topwid, $
            verbose=verbose ge 3)
  msg='    Chi-square value and degrees-of-freedom:       '+ $
      strtrim(chisq,2L)+', '+strtrim(dof,2L)
  error=p3d_misc_logger(msg,logunit,loglevel=2,rname=rname,topwid=topwid, $
            verbose=verbose ge 2)

  if keyword_set(show) then begin ;; use only for testing/developing
    device,decomposed=0
     plot,x,y,/xstyle,psym=10L
    oplot,x,y,psym=7,thick=2.0
    oplot,x,out,color=2
    device,decomposed=1
  endif ;; keyword_set(show)

  return

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