;+
; Project : IRIS
;
; Name    : fit_iris_xput
;
; Purpose : Calculates coefficients of best-fit time function for throughput
;           Two modes: 1. perform fit: supply xput and single element coef
;                      2. apply fit: supply full coef
;
; Use     : IDL> fit   = fit_iris_xput(tvec,tcoef,coef)
;           IDL> resid = fit_iris_xput(tvec,tcoef,coef,xput_in=xput) - xput
;
; Input  :
;    tvec   = vector of times of datapoints
;    tcoef  = dblarr(2,n). start and end times of intervals of constant coef[i]
;
; Input/Output :
;    coef = fltarr(3*n).  Coefficients of best-fit function.
;           fit=c[1,i]*exp(c[2,i]*t)+c[0,i]
;           with c[2,i] negative, in units of 1/yr.
;           output when performing fit, input when applying fit
;
; Keywords input:
;    xput_in = fltarr(n_elements(tvec)).  Data values to be fitted
;              Presence of proper number of values indicates fitting mode
;              Do not supply if just applying fit.
;
; Output:
;    fit = values of fit for times tvec
;
; Common  :
;
; Restrictions:
;
; Side effects:
;
; Category   : Instrument response
;
; Prev. Hist.:
; Written    : JPW, 27-Jan-2015
;              JPW, 28-Jan-2015  changed coefficient set to multiple offsets
;              JPW, 30-Mar-2015  multipley decay constants, all moved into coef
;                                and trying to roughly optimize decay constants
;-

function fit_iris_xput,tt0,tcc0,ccc,xput_in=xx

tex = 1.5         ; exponent for transition between exp.decay intervals
m = size(tcc0)
if m[1] ne 2 then begin
   print,'incorrect number of elements in tcoef'
   return,0
endif else m=m[m[0]+2]/2

; times in years since tcc[0,0]
yr2sec = 3.6525d2*2.4d1*3.6d3
tcc = anytim(tcc0) / yr2sec
tt  = anytim(tt0) / yr2sec
ntt = n_elements(tt)
if n_elements(xx) eq ntt then dofit=1 else dofit=0

if dofit then begin  ; -> perform least squares fit
 ; fit each interval separately, optimize decay exponent for each
 ccc = fltarr(3,m)
 for j=0,m-1 do begin
   ww=where(tt ge tcc[0,j] and tt le tcc[1,j],nww)
   nee = 15
   ee = -(10.0^(findgen(nee)/10.0))
   resid = fltarr(nee)
   ; repeat for a series of decay exponents
   for k=0,nee-1 do begin
     a = fltarr(2,nww)
     a[0,*] = 1.0
     a[1,*] = exp(ee[k]*(tt[ww]-tcc[0,j]))
     svdc,a,w,u,v
     cc = svsol(u,w,v,xx[ww])
     resid[k] = total((reform(a##cc)-xx[ww])^2.0)
   endfor
   dmy = min(resid,wmi)
   ; repeat for  a series of decay exponents on a finer grid
   nee = 11
   ee = -(10.0^(wmi/10.0+findgen(nee)/100.0-0.05))
   resid = fltarr(nee)
   for k=0,nee-1 do begin
     a = fltarr(2,nww)
     a[0,*] = 1.0
     a[1,*] = exp(ee[k]*(tt[ww]-tcc[0,j]))
     svdc,a,w,u,v
     cc = svsol(u,w,v,xx[ww])
     resid[k] = total((reform(a##cc)-xx[ww])^2.0)
   endfor
   dmy = min(resid,wmi)
   ccc[2,j] = ee[wmi]
 endfor
endif else begin
 if n_elements(ccc) ne (3*m) then begin
   print,'incorrect number of elements in coef'
   return,0
 endif
endelse

; calculation of output fit
ee = transpose(ccc[2,*])
a  = fltarr(2*m,ntt)                   ; base vector
if ntt eq 1 then a=reform(a,2*m,ntt)   ; needs to be 2D!!
for j=0,m-1 do begin
  ; base vector for interval of constant multiplier
  if dofit or j gt 0 then ww=where(tt ge tcc[0,j] and tt le tcc[1,j],nww) $
                     else ww=where(tt le tcc[1,j],nww)
  if nww gt 0 then begin
    a[2*j,ww]   = 1.0
    a[2*j+1,ww] = exp(ee[j]*(tt[ww]-tcc[0,j]))
  endif
  ; base vector for interval when multiplier is linear function of time
  if j gt 0 then begin
    ww = where(tt gt tcc[1,j-1] and tt lt tcc[0,j],nww)
    if nww gt 0 then begin
      dtt = (tt[ww]-tcc[1,j-1]) / (tcc[0,j]-tcc[1,j-1])
      a[2*j,ww]   = dtt^tex
      a[2*j+1,ww] = dtt^tex * exp(ee[j]*(tt[ww]-tcc[0,j]))
    endif
  endif
  if j lt (m-1)  then begin
    ww = where(tt gt tcc[1,j] and tt lt tcc[0,j+1],nww)
    if nww gt 0 then begin
      dtt = (tt[ww]-tcc[1,j]) / (tcc[0,j+1]-tcc[1,j])
      a[2*j,ww]   = (1.0-dtt^tex)
      a[2*j+1,ww] = (1.0-dtt^tex) * exp(ee[j]*(tt[ww]-tcc[0,j]))
    endif
  endif
endfor

if dofit then begin  ; -> perform least squares fit
  svdc,a,w,u,v
  cc = svsol(u,w,v,xx)
  ;;; safety feature: ccc[1,m-1] shall not be negative:
  if cc[2*m-1] lt 0.0 then begin
     cc[2*m-2] = cc[2*m-2]+cc[2*m-1]
     cc[2*m-1] = 0.0
  endif
  ;;;
  ccc[0:1,*] = reform(cc,2,m)
endif else cc=reform(ccc[0:1,*],2*m)

f = reform(a##cc)                 ; fit results over tres range

return,f

end