; docformat = 'rst'
;
; NAME: 
;    ssoDistProb
; PURPOSE:
;    Compute the distribution of heliocentric distance from (a,e)
;+
; :Description:
;    Compute the distribution of heliocentric distance from (a,e)
;
; :Categories:
;    Database, SSO
;
; :Params:
;    a: in, required, type=fltarr
;      Semin-major axis
;    e: in, required, type=float
;      Orbital eccentricity
;      
; :Returns: An array of the heliocentric distance
;
; :Keywords:
;    min: in, optional, type=float, default=0
;      The minimum heliocentric distance in the output probability array
;    max: in, optional, type=float, default=500
;      The maximum heliocentric distance in the output probability array
;    step: in, optional, type=float, default=0.01
;      The bin size of heliocentric distance in the output probability array
;    delta: out, optional, type=float
;      The array of heliocentric distances associated with the probability
;
; :Examples:
;    Compute the distribution of heliocentric distance for a randomn Trojan::
;      IDL> proba = ssoDistProb( 5.2, 0.01, delta=dSun )
;      IDL> cgPlot, dSun, proba, xRange=[4,6]
;      IDL> cgPlot, dSun, proba, /OverPlot, psym='OpenCircle'
;
; :Uses:
;
; :Author:
;    B.Carry (OCA)
;
; :History:
;    Change History::
;      Original Version written in February 2017, B. Carry (OCA)
;-
function ssoDistProb, a, e, min=min, max=max, step=step, delta=delta
  COMPILE_OPT hidden, idl2

 ;--I-- Initialization and Input Validation  ---------------------------------------------------;
 ;--I.1-- Syntax
  if n_params() lt 2 then begin
    message, /ioError, 'Syntax: distance = ssoDistProb(a, e [N=])'
 endif

 ;--I.2-- Default Value
  if not keyword_set(min)  then min = 0
  if not keyword_set(max)  then max = 500
  if not keyword_set(step) then step= 0.05

 ;--I.3-- Heliocentric Array
  nbDelta = (max-min)/step + 1
  delta = min + (findgen(nbDelta)/(nbDelta-1))*(max-min)

 ;--I.4-- Singular Null Eccentricity
  if e eq 0 then begin 
    proba=delta*0.
    trash = min( abs(delta-a), pA )
    proba[pA]=1
    return, proba
  endif

 ;--II-- Compute SSO Distance, and Residence Time  ---------------------------------------------;
  N = 360.
  theta = 2*!PI * indgen( N )/ N
  r = a*( 1 - e*e ) / (1 + e*cos(theta) )
  time = 0.5*r*r
  
 ;--III-- Interpolate to Regular Distances  ----------------------------------------------------;
  inside = where( delta ge min(r) and delta le max(r), comp=outside )
  tempo = interpol( time, r, delta[inside] )
  proba=delta*0.
  proba[inside] = tempo
  
 ;--IV-- Normalized the Probability  -----------------------------------------------------------;
  proba /= total(proba)
  return, proba

end