;+
; NAME:
;                  AUX_PSFH
;
; PURPOSE:
;                  Plot the SFH resulting from n-SSPs fit   
;
; USAGE:
;                  plot_sfr, <filename>, PS=ps, \OVER, TIMESHIFT=timeshift, MULTCNST=multcnst, \NORMAL, FILL=fill, XTITLE=xtitle, YTITLE=ytitle, YRANGE=yrange, OVER=over, _EXTRA=extra
;
; INPUTS:
;   <filename>     Name of a file containing results from multi-SSP
;   fit with ULYSS, should end on ".res". string               
;
; KEYWORDS:
;   \OVER          Overplot SFH
;
;   TIMESHIFT      [Gyr] shift the age of the components by a constant
;
;   NORMAL         Normilize the total mass to 1. It can be used
;   together with MULTFRAC. In this case the total mass is first
;   normilized to 1 and after multiplied by 
;
;   MULTCNST       To get the total mass in Solar masses, multiply the
;   component wieght by a constant. For example, in the case of
;   standard input ulyss models and SDSS data, MULTCNST =
;   0.3127e7*D^2, where D^2 is the distance in Mpc (see
;   http://ulyss.univ-lyon1.fr/models.html) 
;
;   FILL           color of the area under the SFH curve
;
;   XTIT           x-axis title
;
;   YTIT           y-axis title
; 
;   YRANGE         y-axis range
;
;   _EXTRA         any keywords accepted by IDL PLOT procedure 
;
; OPTIONAL OUTPUTS:
;   PS             name of the output postscript file. Note that
;   '.eps' will be authomatically added 
;
; RESTRICTIONS:    If a component without free parameters is used (for
; example a stellar spectrum) the program will fail because it looks
; for the parameters names to indentify the number of SSPs used for
; the fit
;
;
; PROCEDURE:
;
;
;
; EXAMPLE:
;  aux_sfh, uly_root+'/data/VazMiles_z-0.40t07.94.fits', [uly_ssp(AL=[10,1000], AG=300.),uly_ssp(AL=[1000,13000])], 'test'
;  loadct, 39
;  aux_psfh, 'test.res', color=40, thick=5, fill=220
;
; HISTORY:
;  200X/XX/XX: MK, creation
;  2012/10/29: MK, documentation 
;  
;-
; CATEGORY:
;                  AUX

pro aux_psfh, filename, PS=ps, OVER=over, $
              TIMESHIFT=timeshift, $
              MULTCNST=multcnst, NORMAL=normal, $
              FILL=fill, $
              XTITLE=xtitle, YTITLE=ytitle, YRANGE=yrange, _EXTRA=extra


  if n_elements(filename) eq 0 then begin
     print, 'Usage: AUX_PSFH, filename, ...'
     return
  endif

;; read solution structure, take the last (most recent) solution
  aa = uly_solut_tread(filename)
  aaa = aa[n_elements(aa)-1]

;; shift the age of the components by a constant in Gyr
  if n_elements(timeshift) eq 0 then timeshift = 0

;; find the number of SSP components 
  cmp_n = 0
  ncmp = 0
  for j = 0,n_elements(aaa[0].cmp)-1 do begin
     if ptr_valid((aaa.cmp[j].para)) then begin
        if (*(aaa.cmp[j].para))[0].name eq 'age' then begin
           cmp_n = [cmp_n,j]
           ncmp ++
        endif
     endif else continue 
  endfor
  if ncmp eq 0 then message, 'No valid SSP CMP in the res file'
  cmp_n = cmp_n[1:*] ; remove the first 0, used to define the variable
  aaa = (aaa.cmp[cmp_n]) ; take only the SSP components

;; set the titles of the axes
  if n_elements(xtitle) eq 0 then xtitle = 'Age, Gyr'
  if n_elements(ytitle) eq 0 then ytitle = (keyword_Set(multcnst) or keyword_Set(normal)) ? 'SFR * 1e3 (M!Dsol!N/yr)' : 'SFR'

;; read the age and the fraction
  age = (*(aaa[0]).para)[0].value
  for k=1, ncmp-1 do age = [age,(*(aaa[k]).para)[0].value]
  age = exp(age) / 1000. + timeshift*[replicate(1,ncmp-1),0] ; in Gyr

;; normalize if asked and eventually multiply by a constant
  frac = (n_elements(normal) gt 0) ? aaa.weight / total(aaa.weight, /NaN) : aaa.weight
  if keyword_Set(multcnst) then frac *= multcnst 

  print, filename, ' had produce ', total(frac, /NaN), ' Msol'

; take the middle between the ages
  age_max = [(( alog(age)+alog(age[1:*])) / 2.),alog(12.)]
  age_min = [2*alog(age[0])-age_max[0], age_max[0:ncmp-1]]

  sfr_uly = (keyword_Set(multcnst) or keyword_Set(normal)) ? $
            frac /  ((exp(age_max) - exp(age_min))) / 1.e+6 : $  ; in units [1e3 Msol/yr]
            frac /  ((exp(age_max) - exp(age_min))) / 1.e+9      ; in units [data_unit/cmp_unit/yr]

  x_plot = exp(reform(transpose([[age_min],[age_min]]), 2*n_elements(age_min)))
  y_plot = [0.,reform(transpose([[sfr_uly],[sfr_uly]]),2*ncmp),0.]

  heap_free, aa
  heap_free, aaa

;;--------------- plot the SFR  --------------------

  if keyword_set(ps) then begin

     set_plot, 'ps'

     !p.charsize=2
     !P.charthick=10
     !p.thick=2
     !x.thick=3
     !y.thick=3
     !P.FONT = 0
     !X.MARGIN=[9,1]
     !Y.MARGIN=[4,1]
     !X.TICKLEN = 0.03
     !Y.TICKLEN = 0.03
     !y.minor = 1

     filename = ps & xs = 9 & ys = 8
     device, FILENAME=filename+'.eps', /COLOR, XS=xs, YS=ys, FONT_SIZE=4, /TIMES, /ENC

  endif

  if  not keyword_set(over) then begin
     plot, x_plot, y_plot, TITLE=title, $
           XST=3, XR=[0.01,12.1], XTIT=xtitle, $ ; /XLOG, $
           YST=3, YR=yrange, YTIT=ytitle,$
           /NODATA
  endif
  if keyword_set(fill) then $
     polyfill, [x_plot, reverse(x_plot)], [y_plot, replicate(0.,n_elements(y_plot))], COL=fill;, ORIENT=-45
  oplot, x_plot, y_plot, _EXTRA=extra, COLOR=color ; oplot results ULYSS

  if keyword_set(ps) then begin
     device, /CLOSE
     set_plot, 'x'
  endif

end
;;-- end -------------------------------------------