; docformat = 'rst'
;
; NAME: 
;    satForces
; PURPOSE:
;    Compute and display order of magnitude accelerations acting on a satellite
;
;+
; :Description:
;    Compute and display order of magnitude accelerations acting on a satellite
;
; :Categories:
;    Ephemeris
;
; :Params:
;    ID: in, required, type=integer/string
;        Target identifyer ('*', Name, Provisional designation, Number) 
;
; :Returns: A structure containing the ephemeris
;
; :Keywords:
;
; :Examples:
;
; :Uses:
;    
;
; :Author:
;    B.Carry (IMCCE)
;
; :History:
;     Change History::
;        Original Version written in December 2014, B. Carry (IMCCE)
;-
function satForces, id, a;, perturbers=perturbers, moons=moons, planets=planets, sun=sun, GR=GR, dump=dump

;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---    I -- Initialization And Input Verification              -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
  COMPILE_OPT hidden

 ;--I.1-- Input Target ID  --------------------------------------------------------
  if not keyword_set(id) then begin
    message, /ioError, 'Syntax - val = satForces( id [, a, perturbers=, moons=, planets=, sun=, /GR])'
    return, -1
  endif

;tbd: ID can be a:xxx, c:..., p:..., 


 ;--I.2-- Pertubing Moons  --------------------------------------------------------

 ;--I.3-- Pertubing Planets  ------------------------------------------------------
  pINFO = replicate({name:'', a:0.d, e:0., M:0.}, 8)
    pINFO(0).name ='Mercury'
    pINFO(0).M    = 3.301e23
    pINFO(0).a    = 0.7
           
    pINFO(3).name ='Mars'
    pINFO(3).M    = 6.4e23
    pINFO(3).a    = 1.52
    pINFO(3).e    = 0.0934

    pINFO(4).name ='Jupiter'
    pINFO(4).M    = 1.899e27
    pINFO(4).a    = 5.204
    pINFO(4).e    = 0.048775

  if not keyword_set(planets) then begin
    planets=intarr(8)+1
  endif

 ;--I.4-- Pertubing Sun  ----------------------------------------------------------
  M_sun = 1.98855e30


  G = 6.67384e-11 ;SI: m3.kg-1.s-2
  au=149597871e3  ;SI: m

;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---   II -- Computatino of  Accelerations                      -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;


 ;-log
  r = [1,1e8];findgen(10000e3)+1

 ;-linear
  r = [1, $
       1+findgen(100), $
       100+2*(findgen(100)+1), $
       300+4*(findgen(100)+1), $
       700+8*(findgen(100)+1), $
       1500+16*(findgen(100)+1), $
       3100+32*(findgen(100)+1), $
       6300+64*(findgen(100)+1), $
       12700+128*(findgen(100)+1)  ]*1e3


 ;--huenna
  m = 5.598+17  ;-kg
  a=3559.
  e=0.281
  J2 = 0.1
  J4 = 0

 ;-alauda
;  m = 6.450943e+18
;  a=1252.100
;  e=0.04741
;  J2 = 0.187
;  J4 = 0.151
;
; ;-elektra
;  m = 6.6e+18
;  a=1325
;  e=0.1489
;  J2 = 0.13
;  J4 = 0
;  orbit={a:3.12,e:0.208}


  gC = G*m/(R*R)



  gJ2 = J2*gC
  gJ4 = J4*gC

  print, min(gJ2), max(gJ2)

  orbit = astElements(id)

  help, orbit




  gSun = 2*G*M_sun*r / (orbit.a*au)^3.
  gSunMax = 2*G*M_sun*r / ((orbit.a*(1-orbit.e))*au)^3.
  gSunMin = 2*G*M_sun*r / ((orbit.a*(1+orbit.e))*au)^3.


  aph  = (orbit.a*(1+orbit.e))*au
  peri = (orbit.a*(1-orbit.e))*au

  qMar = (pINFO(3).a*(1+pINFO(3).e))*au
  pMar = (pINFO(3).a*(1-pINFO(3).e))*au
  qJup = (pINFO(4).a*(1+pINFO(4).e))*au
  pJup = (pINFO(4).a*(1-pINFO(4).e))*au


  gMar    = 2*G*pINFO(3).m*r / (abs(orbit.a-pINFO(3).a)*au)^3.
  gMarMax = 2*G*pINFO(3).m*r / (abs(pMar-aph))^3.
  gMarMin = 2*G*pINFO(3).m*r / (abs(qMar-peri))^3.

  gJup    = 2*G*pINFO(4).m*r / (abs(orbit.a-pINFO(4).a)*au)^3.
  gJupMax = 2*G*pINFO(4).m*r / (abs(pJup-aph))^3.
  gJupMin = 2*G*pINFO(4).m*r / (abs(qJup-peri))^3.



  rkm = r/1e3

  cgPlot, rkm, gc, $
;          /xlog, $
          /yLog, $
          yr=[1e-13,1e-7], $
          xr=[10,10000], $
          xtitle='Distance from barycenter (km)', $
          ytitle='Acceleration (m.s!U-2!N)', $
          xminor=5

  cgPlot, rkm, gSun,    /over, color='Orange'
  cgPlot, rkm, gSunMin, /over, color='Orange'
  cgPlot, rkm, gSunMax, /over, color='Orange'

  cgPlot, rkm, gJ2,     /over, color='Grey'
  cgPlot, rkm, gJ4,     /over, color='Grey'

  cgPlot, rkm, gMar,    /over, color='Red'
  cgPlot, rkm, gMarMax, /over, color='Red'
  cgPlot, rkm, gMarMin, /over, color='Red'

  cgPlot, rkm, gJup,    /over, color='Green'
  cgPlot, rkm, gJupMax, /over, color='Green'
  cgPlot, rkm, gJupMin, /over, color='Green'
  

;print, max(gJupMax), max(gJupMin)



  cgPlot, /over, [a,a], [1e-20,1e10]     , color='Cornflower Blue'
  cgPlot, /over, [a,a]*(1-e),[1e-20,1e10], color='Cornflower Blue', linestyle=2
  cgPlot, /over, [a,a]*(1+e),[1e-20,1e10], color='Cornflower Blue', linestyle=2






  return, 0

end