;
; Computes the PDF of data
; Input(s): data (array)
; Optional inputs: nbins, min_val, max_val, log_binning, bins_in
; Outputs: bins_loc, pdf, hist, cmf, mode_val, mean_val, stddev_val
; written by Christoph Federrath, 2016
;
pro get_pdf, data, NBINS=nbins, MIN_VAL=min_val, MAX_VAL=max_val, LOG_BINNING=log_binning, BINS_IN=bins_in, $
            BINS_LOC=bins_loc, PDF=pdf, HISTOGRAM=histogram, CMF=cmf, DELTA_BINS=delta_bins, $
            MODE_VAL=mode_val, MEAN_VAL=mean_val, STDDEV_VAL=stddev_val, $
            DEBUG=debug, _EXTRA=extra

    if not keyword_set(nbins) then nbins = 100
    if not keyword_set(min_val) then min_val = min(data)
    if not keyword_set(max_val) then max_val = max(data)
    if keyword_set(log_binning) then begin
        ; make log bins
        bins = 10d0^(dindgen(nbins)/(nbins-1)*alog10(max_val/min_val)+alog10(min_val))
    endif else begin
        ; make lin bins
        bins = dindgen(nbins)/(nbins-1)*(max_val-min_val)+min_val
    endelse
    if keyword_set(bins_in) then begin
        bins = bins_in ; user bins
        nbins = n_elements(bins_in) ; overwrite nbins
    endif

    if keyword_set(debug) then DEBUG = 1 else DEBUG = 0

    ; individual bin widths
    delta_bins = shift(bins,-1)-bins
    ; find indices
    located_list = value_locate(bins, data)
    hist = histogram(located_list)

    ; now treat cases of different return values from value_locate()
    ; outside bounds, etc...
    if min(located_list) eq -1 then begin
        if max(located_list) eq nbins-1 then begin
            nbins_valid = n_elements(hist)-2
            is_bins = max(located_list)-nbins_valid
        endif else begin
            nbins_valid = n_elements(hist)-1
            is_bins = max(located_list)-nbins_valid+1
        endelse
        is_hist = 1 & ie_hist = nbins_valid
    endif else begin
        if max(located_list) eq nbins-1 then begin
            nbins_valid = n_elements(hist)-1
            is_bins = max(located_list)-nbins_valid
        endif else begin
            nbins_valid = n_elements(hist)
            is_bins = max(located_list)-nbins_valid+1
        endelse
        is_hist = 0 & ie_hist = nbins_valid-1
    endelse
    ie_bins = is_bins+nbins_valid-1

    ; output is always nbins-1 bins starting at left edge
    pdf = dblarr(nbins-1)
    histogram = dblarr(nbins-1)
    bins = bins[0:nbins-2]
    delta_bins = delta_bins[0:nbins-2]
    bins_loc = bins + delta_bins/2d0 ; get location bins (bin-centered grid)
    cmf = dblarr(nbins-1)

    if (DEBUG) then print, 'nbins_valid = ', nbins_valid

    if nbins_valid gt 0 then begin
        ; PDF normlization (division by bin widths and total number of points in histogram)
        pdf[is_bins:ie_bins] = double(hist[is_hist:ie_hist])
        pdf /= double(total(hist[is_hist:ie_hist]))
        pdf /= delta_bins
        if (DEBUG) then print, 'PDF normalization: ', total(pdf*delta_bins,/double)
        histogram[is_bins:ie_bins] = hist[is_hist:ie_hist]
        cmf = total(pdf*delta_bins,/cumulative,/double)
    endif else begin
        return
    endelse

    ; get the mode, mean, and stddev
    ind = where(pdf eq max(pdf), count)
    mode_val = mean((bins_loc[ind])[0:count-1],/double)
    if (DEBUG) then print, 'mode = ', mode_val
    mean_val = total(pdf*bins_loc*delta_bins,/double)
    if (DEBUG) then print, 'mean = ', mean_val
    ms_val = total(pdf*bins_loc^2*delta_bins,/double)
    stddev_val = sqrt(ms_val-mean_val^2)
    if (DEBUG) then print, 'stddev_val = ', stddev_val

    return

end ; ===============================================