(*
  simple ephemeris from
  http://iau-comm4.jpl.nasa.gov/XSChap8.pdf
  Orbital Ephemerides of the Sun, Moon, and Planets 
  presented by E. Myles Standish and James G. Williams
  angles for ecliptic J2000 in degree
  elements values with time changing:
               a(AU)      e         i(degree)       L(degree)      \Om+\om       \Omega
Юпитер     5.20248019  0.04853590  1.29861416     34.33479152  14.27495244 100.29282654
          -0.00002864  0.00018026 -0.00322699   3034.90371757   0.18199196   0.13024619
Сатурн     9.54149883  0.05550825  2.49424102     50.07571329  92.86136063 113.63998702
          -0.00003065 -0.00032044  0.00451969   1222.11494724   0.54179478  -0.25015002
Уран      19.18797948  0.04685740  0.77298127    314.20276625 172.43404441  73.96250215
          -0.00020455 -0.00001550 -0.00180155    428.49512595   0.09266985   0.05739699
Нептун    30.06952752  0.00895439  1.77005520    304.22289287  46.68158724 131.78635853
           0.00006447  0.00000818  0.00022400    218.46515314   0.01009938  -0.00606302
            b            c             s              f
Юпитер    -0.00012452   0.06064060   -0.35635438    38.35125000
Сатурн     0.00025899  -0.13434469    0.87320147    38.35125000
Уран       0.00058331  -0.97731848    0.17689245     7.67025000
Нептун    -0.00041348   0.68346318   -0.10162547     7.67025000
   b   -   градусы в столетие за столетие,
   c   -   градусы,
   s   -   градусы,
   f   -   градусы за столетие.
   \omega=(\Omega+\omega)-\Omega
   M=L-(\Omega+\omega)+b*T*T+c*cos(f*T)+s*sin(f*T)
   T in century from J2000.0
*)
unit ungephem ; { to obtain simple elements for giant planets }
interface
uses
  unsagcot ; { type tvect3 some constants }
{ simple heliocentric ecliptic
  Kepler elements for giant planets }
procedure forgianel ( n  : byte ;      { number of planets }
                      t  : extended ;  { current moment in julian day }
                      var  p  : tvect3 ; { positional elements a e i }
                      var  u  : tvect3 ) ; { angle elements \O \o M }
{ simple heliocentric equatorial position for giant planets }
procedure forgiapos ( n  : byte ;      { number of planets }
                      t  : extended ;  { current moment in julian day }
                      var  r  : tvect3 ) ; { position x y z equatorial }
implementation
uses
  unforkep , { procedure fromkeplerpos }
  unforequ ; { procedure fromeclipequa }
{ for the Jupiter simple middle elements }
procedure elmupi ( d  : extended ;  { current moment in julian day }
                   var  p  : tvect3 ; { positional elements a e i }
                   var  u  : tvect3 ) ; { angle elements \O \o M }
var
  q,v,w  : extended ; { temporary angles }
begin
  p[1]:= 5.20248019-0.00002864*d; { semimajor axis in AU }
  p[2]:= 0.04853590+0.00018026*d; { eccentricity }
  p[3]:= 1.29861416-0.00322699*d; { inclination in degree }
  v:= 14.27495244+0.18199196*d;   { longitude of perihelium degree }
  u[1]:= 100.29282654+0.13024619*d; { the ascending node in degree }
  u[2]:= v-u[1]; { the argument of perihelium in degree }
  w:= 34.33479152+3034.90371757*d-v; { the mean anomaly in degree }
  q:=grarad*38.35125000*d; { for auxiliary periodic term in radian }
  u[3]:=w-0.00012452*d*d   { the mean anomaly with secular and }
         +0.06064060*cos(q)-0.35635438*sin(q); { periodic terms }
end;
{ for the Saturn simple middle elements }
procedure elmatu ( d  : extended ;  { current moment in julian day }
                   var  p  : tvect3 ; { positional elements a e i }
                   var  u  : tvect3 ) ; { angle elements \O \o M }
var
  q,v,w  : extended ; { temporary angles }
begin
  p[1]:= 9.54149883-0.00003065*d; { semimajor axis in AU }
  p[2]:= 0.05550825-0.00032044*d; { eccentricity }
  p[3]:= 2.49424102+0.00451969*d; { inclination in degree }
  v:= 92.86136063+0.54179478*d;   { longitude of perihelium degree }
  u[1]:= 113.63998702-0.25015002*d; { the ascending node in degree }
  u[2]:= v-u[1]; { the argument of perihelium in degree }
  w:= 50.07571329+1222.11494724*d-v; { the mean anomaly in degree }
  q:=grarad*38.35125000*d; { for auxiliary periodic term in radian }
  u[3]:=w+0.00025899*d*d   { the mean anomaly with secular and }
         -0.13434469*cos(q)+0.87320147*sin(q); { periodic terms }
end;
{ for the Uranus simple middle elements }
procedure elmran ( d  : extended ;  { current moment in julian day }
                   var  p  : tvect3 ; { positional elements a e i }
                   var  u  : tvect3 ) ; { angle elements \O \o M }
var
  q,v,w  : extended ; { temporary angles }
begin
  p[1]:=19.18797948-0.00020455*d; { semimajor axis in AU }
  p[2]:= 0.04685740-0.00001550*d; { eccentricity }
  p[3]:= 0.77298127-0.00180155*d; { inclination in degree }
  v:=172.43404441+0.09266985*d;   { longitude of perihelium degree }
  u[1]:= 73.96250215+0.05739699*d; { the ascending node in degree }
  u[2]:= v-u[1]; { the argument of perihelium in degree }
  w:= 314.20276625+428.49512595*d-v; { the mean anomaly in degree }
  q:=grarad*7.67025000*d; { for auxiliary periodic term in radian }
  u[3]:=w+0.00058331*d*d   { the mean anomaly with secular and }
         -0.97731848*cos(q)+0.17689245*sin(q); { periodic terms }
end;
{ for the Neptune simple middle elements }
procedure elmnep ( d  : extended ;  { current moment in julian day }
                   var  p  : tvect3 ; { positional elements a e i }
                   var  u  : tvect3 ) ; { angle elements \O \o M }
var
  q,v,w  : extended ; { temporary angles }
begin
  p[1]:=30.06952752+0.00006447*d; { semimajor axis in AU }
  p[2]:= 0.00895439+0.00000818*d; { eccentricity }
  p[3]:= 1.77005520+0.00022400*d; { inclination in degree }
  v:= 46.68158724+0.01009938*d;   { longitude of perihelium degree }
  u[1]:=131.78635853-0.00606302*d; { the ascending node in degree }
  u[2]:= v-u[1]; { the argument of perihelium in degree }
  w:= 304.22289287+218.46515314*d-v; { the mean anomaly in degree }
  q:=grarad*7.67025000*d; { for auxiliary periodic term in radian }
  u[3]:=w-0.00041348*d*d   { the mean anomaly with secular and }
         +0.68346318*cos(q)-0.10162547*sin(q); { periodic terms }
end;
{ simple heliocentric ecliptic
  Kepler elements for giant planets }
procedure forgianel ( n  : byte ;      { number of planets }
                      t  : extended ;  { current moment in julian day }
                      var  p  : tvect3 ; { positional elements a e i }
                      var  u  : tvect3 ) ; { angle elements \O \o M }
var
  d  : extended ; { moment from J2000 in century }
begin
  d:=(t-jd2000)/julianc; { from J2000 in century const unsagcot }
  case  n  of
    5 : elmupi(d,p,u);  { the Jupiter }
    6 : elmatu(d,p,u);  { the Saturn }
    7 : elmran(d,p,u);  { the Uranus }
    8 : elmnep(d,p,u);  { the Neptune }
  end; { case for giant planet }
end;
{ simple heliocentric equatorial position for giant planets }
procedure forgiapos ( n  : byte ;      { number of planets }
                      t  : extended ;  { current moment in julian day }
                      var  r  : tvect3 ) ; { position x y z equatorial }
var
  p  : tvect3 ; { positional elements a e i }
  u  : tvect3 ; { angle elements \O \o M }
  q  : tvect3 ; { heliocentric ecliptic position }
begin
  forgianel(n,t,p,u);   { ecliptic heliocentric Kepler elements }
  fromkeplerpos(p,u,q); { unforkep heliocentric ecliptic position }
  fromeclipequa(q,r);   { unforequ to heliocentric equator position }
end;
end.