; docformat = 'rst'
;
; NAME: 
;    ssoModelSFD
; PURPOSE:
;    Create simple power-law H-Distributions
;+
; :Description:
;    Create simple power-law H-Distributions
;
; :Categories:
;    Database, SSO
;
; :Params:
;    h: in, required, type=fltarr
;      An array of absolute magnitudes
;    bright: in, required, type=float
;      The power-law index for a single-slope distribution, or of its
;      bright end only in double/divot distributions.
;    pivot: in, optional, type=float
;      The absolute magnitude of the knee in double/divot distributions.
;    faint: in, optional, type=float
;      The power-law index for the faint end ofdouble/divot distributions.
;    contrast: in, optional, type=float
;      The ratio between the bright and faint end in divot differential distributions
;    reference: in, optional, type=struct
;      Give a reference value to scale the differential power-laws.
;      
;
; :Returns: A structure with the absolute magnitude array, the
;    differential and cumulative distribution, and a description of
;    the model. The fields are::
;
; :Keywords:
;    single: in, optional, type=boolean, default=1
;      Model the differential H-distribution by a single slope
;    knee: in, optional, type=boolean, default=0
;      Model the differential H-distribution by two slopes
;    divot: in, optional, type=boolean, default=0
;      Model the differential H-distribution by two slopes with a drop at the knee
;
; :Examples:
;
; :Uses:
;
; :Author:
;    B.Carry (OCA)
;
; :History:
;    Change History::
;      Original Version written in February 2017, B. Carry (OCA)
;-
function ssoModelSFD, h, bright, pivot, faint, contrast, reference=reference, $
                      single=single, knee=knee, divot=divot
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---    I -- Initialization And Input Verification              -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
  COMPILE_OPT hidden, idl2

 ;--I.1-- Check Inputs  -------------------------------------------------------------------------;
  if n_params() lt 2 then begin
    message, /ioError, 'Syntax: sfd = ssoModelSFD( h, bright [pivot, faint, contrast, reference=...])'
    return, -1
 endif


 ;--I.2-- Determine Which Function to Use  ------------------------------------------------------;
  if keyword_set(pivot) and keyword_set(faint) then knee=1
  if keyword_set(contrast) then begin
     knee=0
     divot=1
  endif

  modelType= [keyword_set(single), keyword_set(knee), keyword_set(divot) ]
  if total(modelType) eq 0 then modelType=[1,0,0] 

  if total(modelType) gt 1 then begin
    if not keyword_set(constrast) then divot=0
    if not keyword_set(pivot) or not keyword_set(faint) then knee=0
    modelType= [keyword_set(single), keyword_set(knee), keyword_set(divot) ]
    if total(modelType) eq 0 then modelType=[1,0,0] 
  endif

  
 ;--I.3-- H-Distribution Model Structure  -------------------------------------------------------;
  emptyModel={type:'',bright:0.,faint:0.,pivot:0.,contrast:0.}



;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---   II -- Single Slope H-distribution                        -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
  if modelType[0] eq 1 then begin

   ;--II.1-- Store Model Information  ----------------------------------------------------------
    model=emptyModel
    model.type  = 'single'
    model.bright= bright

   ;--II.2-- Compute Power-Law H-Distribution  -------------------------------------------------
    histo = 10d^(h * bright)

   ;--II.3-- Scale H-Distribution to Input Reference  ------------------------------------------
    if keyword_set(reference) then begin
      C = reference.N * 10d^( -bright * reference.H )
      histo *= C
    endif
    
   ;--II.4-- Cumulative H-Distribution  --------------------------------------------------------
    cumul = total(histo,/Cumulative)

   ;--II.5-- Return  ---------------------------------------------------------------------------
    return, {h:h, histo:{obs:histo, model:histo}, cumul:{obs:cumul, model:cumul}, model:model}

  endif


;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---  III -- Double Slope (Knee) H-distribution                 -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
  if modelType[1] eq 1 then begin

   ;--III.1-- Store Model Information  ---------------------------------------------------------
    model=emptyModel
    model.type  = 'knee'
    model.bright= bright
    model.faint = faint
    model.pivot = pivot
    
   ;--III.2-- Compute Power-Law H-Distribution  ------------------------------------------------
   ;--III.2.1-- Bright End
    histo = 10d^(h * bright)

   ;--III.2.2-- Faint End
    fRange = where( h ge pivot, nbF, comp=bRange, ncomp=nbB )
    histo[fRange] = 10d^(h[fRange] * faint) 

   ;--III.2.3-- Scale Faint End to the Maximum of the Bright End
    contrast = 1.*histo[bRange[nbB-1]] / histo[fRange[0]]
    histo[fRange] *= contrast
    
   ;--III.3-- Scale H-Distribution to Input Reference  -----------------------------------------
    if keyword_set(reference) then begin
      C = reference.N * 10d^( -bright * reference.H )
      histo *= C
    endif
    
   ;--III.4-- Cumulative H-Distribution  -------------------------------------------------------
    cumul = total(histo,/Cumulative)

   ;--III.5-- Return  --------------------------------------------------------------------------
    return, {h:h, histo:{obs:histo, model:histo}, cumul:{obs:cumul, model:cumul}, model:model}

  endif



  
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---   IV -- Double Slope with Drop (Divot) H-distribution      -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
  if modelType[2] eq 1 then begin

   ;--IV.1-- Store Model Information  ----------------------------------------------------------
    model=emptyModel
    model.type     = 'divot'
    model.bright   = bright
    model.faint    = faint
    model.pivot    = pivot
    model.contrast = contrast
    
   ;--IV.2-- Compute Power-Law H-Distribution  -------------------------------------------------
   ;--IV.2.1-- Bright End
    histo = 10d^(h * bright)

   ;--IV.2.2-- Faint End
    fRange = where( h ge pivot, nbF, comp=bRange, ncomp=nbB )
    histo[fRange] = 10d^(h[fRange] * faint) 

   ;--IV.2.3-- Scale Faint End to the Maximum of the Bright End
    fac = (1.*histo[bRange[nbB-1]] / histo[fRange[0]])/contrast
    histo[fRange] *= fac
    
   ;--IV.3-- Scale H-Distribution to Input Reference  ------------------------------------------
    if keyword_set(reference) then begin
      C = reference.N * 10d^( -bright * reference.H )
      histo *= C
    endif
    
   ;--IV.4-- Cumulative H-Distribution  --------------------------------------------------------
    cumul = total(histo,/Cumulative)

;    cgplot, h, cumul, /ylog, xr=[6,12]
;    cgplot, h, histo, /over, lineStyle=2
;
;    h6 = where( h eq 6)
;    h8=wherE( h eq 8)
;    print, cumul[[h6,h8]]

;    stop


   ;--IV.5-- Return  ---------------------------------------------------------------------------
    return, {h:h, histo:{obs:histo, model:histo}, cumul:{obs:cumul, model:cumul}, model:model}

  endif
  
end