tpsystools_4.04/source/StAstroP.pas

// Upgraded to Delphi 2009: Sebastian Zierer

(* ***** BEGIN LICENSE BLOCK *****
 * Version: MPL 1.1
 *
 * The contents of this file are subject to the Mozilla Public License Version
 * 1.1 (the "License"); you may not use this file except in compliance with
 * the License. You may obtain a copy of the License at
 * http://www.mozilla.org/MPL/
 *
 * Software distributed under the License is distributed on an "AS IS" basis,
 * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
 * for the specific language governing rights and limitations under the
 * License.
 *
 * The Original Code is TurboPower SysTools
 *
 * The Initial Developer of the Original Code is
 * TurboPower Software
 *
 * Portions created by the Initial Developer are Copyright (C) 1996-2002
 * the Initial Developer. All Rights Reserved.
 *
 * Contributor(s):
 *
 * ***** END LICENSE BLOCK ***** *)

{*********************************************************}
{* SysTools: StAstroP.pas 4.04                           *}
{*********************************************************}
{* SysTools: Astronomical Routines (general Planetary)   *}
{*********************************************************}

{$I StDefine.inc}

{ ************************************************************** }
{ Sources:                                                       }
{   1. Astronomical Algorithms, Jean Meeus, Willmann-Bell, 1991. }
{                                                                }
{   2. Planetary and Lunar Coordinates (1984-2000), U.S. Govt,   }
{         1983.                                                  }
{                                                                }
{   3. Supplement to the American Ephemeris and Nautical Almanac,}
{        U.S. Govt, 1964.                                        }
{                                                                }
{   4. MPO96-MPO98 source files, Brian D. Warner, 1995-2000.     }
{                                                                }
{ ************************************************************** }

unit StAstroP;

interface

const
  StdDate = 2451545.0;            {Ast. Julian Date for J2000 Epoch}
  OB2000  = 0.409092804;          {J2000 obliquity of the ecliptic (radians)}

type
  TStEclipticalCord = packed record
    L0,
    B0,
    R0   : Double;
  end;

  TStRectangularCord = packed record
    X,
    Y,
    Z   : Double;
  end;

  TStPlanetsRec = packed record
    RA,
    DC,
    Elong  : Double;
  end;
  TStPlanetsArray = array[1..8] of TStPlanetsRec;


procedure PlanetsPos(JD : Double; var PA : TStPlanetsArray);


implementation

uses
  Windows,
  StDate, StMerc, StVenus, StMars, StJup, StSaturn, StUranus, StNeptun,
  StPluto, StMath;

var
  PlanEC   : TStEclipticalCord;
  PlanRC,
  SunRC    : TStRectangularCord;
  SunEQ    : TStPlanetsRec;


{--------------------------------------------------------------------------}

function RealAngle(Value2, Value1, Start : Double) : Double;
begin
  Result := Start;

  if (Value1 = 0) then begin
    if Value2 > 0 then
      Result := Pi / 2.0
    else
      Result := 3.0 * Pi / 2.0;
  end else begin
    if (Value2 > 0.0) then begin
      if (Value1 < 0.0) then
        Result := Start + Pi
      else
        Result := Start;
    end else begin
      if (Value2 = 0) then begin
        if Value1 > 0 then
          Result := 0
        else
          Result := Pi;
      end else begin
        if (Value2 < 0) then begin
          if (Value1 < 0) then
            Result := Start + Pi
          else
            Result := Start + (2.0 * Pi)
        end;
      end;
    end;
  end;
end;

{--------------------------------------------------------------------------}

function SunOfDate(JD : Double) : TStRectangularCord;
  {-compute J2000 XYZ coordinates of the Sun}
var
  T0,
  A,
  L,
  B,
  RV,
  TX,
  TY,
  TZ  : Double;

begin
  T0 := (JD - StdDate) / 365250;

{solar longitude}
  L := 175347046
     +   3341656 * cos(4.6692568 +  6283.07585*T0)
     +     34894 * cos(4.6261000 + 12566.1517*T0)
     +      3497 * cos(2.7441000 +  5753.3849*T0)
     +      3418 * cos(2.8289000 +     3.5231*T0)
     +      3136 * cos(3.6277000 + 77713.7715*T0)
     +      2676 * cos(4.4181000 +  7860.4194*T0)
     +      2343 * cos(6.1352000 +  3930.2097*T0)
     +      1324 * cos(0.7425000 + 11506.7698*T0)
     +      1273 * cos(2.0371000 +   529.6910*T0)
     +      1199 * cos(1.1096000 +  1577.3435*T0)
     +       990 * cos(5.2330000 +  5884.9270*T0)
     +       902 * cos(2.0450000 +    26.1490*T0)
     +       857 * cos(3.5080000 +   398.149*T0)
     +       780 * cos(1.1790000 +  5223.694*T0)
     +       753 * cos(2.5330000 +  5507.553*T0)
     +       505 * cos(4.5830000 + 18849.228*T0)
     +       492 * cos(4.2050000 +   775.523*T0)
     +       357 * cos(2.9200000 +     0.067*T0)
     +       317 * cos(5.8490000 + 11790.626*T0)
     +       284 * cos(1.8990000 +   796.298*T0)
     +       271 * cos(0.3150000 + 10977.079*T0)
     +       243 * cos(0.3450000 +  5486.778*T0)
     +       206 * cos(4.8060000 +  2544.314*T0)
     +       205 * cos(1.8690000 +  5573.143*T0)
     +       202 * cos(2.4580000 +  6069.777*T0)
     +       156 * cos(0.8330000 +   213.299*T0)
     +       132 * cos(3.4110000 +  2942.463*T0)
     +       126 * cos(1.0830000 +    20.775*T0)
     +       115 * cos(0.6450000 +     0.980*T0)
     +       103 * cos(0.6360000 +  4694.003*T0)
     +       102 * cos(0.9760000 + 15720.839*T0)
     +       102 * cos(4.2670000 +     7.114*T0)
     +        99 * cos(6.2100000 +  2146.170*T0)
     +        98 * cos(0.6800000 +   155.420*T0)
     +        86 * cos(5.9800000 +161000.690*T0)
     +        85 * cos(1.3000000 +  6275.960*T0)
     +        85 * cos(3.6700000 + 71430.700*T0)
     +        80 * cos(1.8100000 + 17260.150*T0);

  A := 628307584999.0
     + 206059 * cos(2.678235 +  6283.07585*T0)
     +   4303 * cos(2.635100 + 12566.1517*T0)
     +    425 * cos(1.590000 +     3.523*T0)
     +    119 * cos(5.796000 +    26.298*T0)
     +    109 * cos(2.966000 +  1577.344*T0)
     +     93 * cos(2.590000 + 18849.23*T0)
     +     72 * cos(1.140000 +   529.69*T0)
     +     68 * cos(1.870000 +   398.15*T0)
     +     67 * cos(4.410000 +  5507.55*T0)
     +     59 * cos(2.890000 +  5223.69*T0)
     +     56 * cos(2.170000 +   155.42*T0)
     +     45 * cos(0.400000 +   796.30*T0)
     +     36 * cos(0.470000 +   775.52*T0)
     +     29 * cos(2.650000 +     7.11*T0)
     +     21 * cos(5.340000 +     0.98*T0)
     +     19 * cos(1.850000 +  5486.78*T0)
     +     19 * cos(4.970000 +   213.30*T0)
     +     17 * cos(2.990000 +  6275.96*T0)
     +     16 * cos(0.030000 +  2544.31*T0);
  L := L + (A * T0);

  A := 8722 * cos(1.0725 +  6283.0758*T0)
     +  991 * cos(3.1416)
     +  295 * cos(0.437  + 12566.1520*T0)
     +   27 * cos(0.050  +     3.52*T0)
     +   16 * cos(5.190  +    26.30*T0)
     +   16 * cos(3.69   +   155.42*T0)
     +    9 * cos(0.30   + 18849.23*T0)
     +    9 * cos(2.06   + 77713.77*T0);
  L := L + (A * sqr(T0));

  A := 289 * cos(5.842 +  6283.076*T0)
     +  21 * cos(6.05  + 12566.15*T0)
     +   3 * cos(5.20  +   155.42*T0)
     +   3 * cos(3.14);
  L := L + (A * sqr(T0) * T0);
  L := L / 1.0E+8;


{solar latitude}
  B := 280 * cos(3.199 + 84334.662*T0)
     + 102 * cos(5.422 +  5507.553*T0)
     +  80 * cos(3.88  +  5223.69*T0)
     +  44 * cos(3.70  +  2352.87*T0)
     +  32 * cos(4.00  +  1577.34*T0);
  B := B / 1.0E+8;

  A := 227778 * cos(3.413766 + 6283.07585*T0)
     +   3806 * cos(3.3706 + 12566.1517*T0)
     +   3620
     +     72 * cos(3.33 + 18849.23*T0)
     +      8 * cos(3.89 +  5507.55*T0)
     +      8 * cos(1.79 +  5223.69*T0)
     +      6 * cos(5.20 +  2352.87*T0);
  B := B + (A * T0 / 1.0E+8);

  A := 9721 * cos(5.1519 + 6283.07585*T0)
     +  233 * cos(3.1416)
     +  134 * cos(0.644  + 12566.152*T0)
     +    7 * cos(1.07   + 18849.23*T0);
  B := B + (A * sqr(T0) / 1.0E+8);

  A := 276 * cos(0.595 + 6283.076*T0)
     +  17 * cos(3.14)
     +   4 * cos(0.12  + 12566.15*T0);
  B := B + (A * sqr(T0) * T0 / 1.0E+8);


{solar radius vector (astronomical units)}
  RV := 100013989
     +   1670700 * cos(3.0984635 +  6283.07585*T0)
     +     13956 * cos(3.05525   + 12566.15170*T0)
     +      3084 * cos(5.1985    + 77713.7715*T0)
     +      1628 * cos(1.1739    +  5753.3849*T0)
     +      1576 * cos(2.8649    +  7860.4194*T0)
     +       925 * cos(5.453     + 11506.770*T0)
     +       542 * cos(4.564     +  3930.210*T0)
     +       472 * cos(3.661     +  5884.927*T0)
     +       346 * cos(0.964     +  5507.553*T0)
     +       329 * cos(5.900     +  5223.694*T0)
     +       307 * cos(0.299     +  5573.143*T0)
     +       243 * cos(4.273     + 11790.629*T0)
     +       212 * cos(5.847     +  1577.344*T0)
     +       186 * cos(5.022     + 10977.079*T0)
     +       175 * cos(3.012     + 18849.228*T0)
     +       110 * cos(5.055     +  5486.778*T0)
     +        98 * cos(0.89      +  6069.78*T0)
     +        86 * cos(5.69      + 15720.84*T0)
     +        86 * cos(1.27      +161000.69*T0)
     +        65 * cos(0.27      + 17260.15*T0)
     +        63 * cos(0.92      +   529.69*T0)
     +        57 * cos(2.01      + 83996.85*T0)
     +        56 * cos(5.24      + 71430.70*T0)
     +        49 * cos(3.25      +  2544.31*T0)
     +        47 * cos(2.58      +   775.52*T0)
     +        45 * cos(5.54      +  9437.76*T0)
     +        43 * cos(6.01      +  6275.96*T0)
     +        39 * cos(5.36      +  4694.00*T0)
     +        38 * cos(2.39      +  8827.39*T0)
     +        37 * cos(0.83      + 19651.05*T0)
     +        37 * cos(4.90      + 12139.55*T0)
     +        36 * cos(1.67      + 12036.46*T0)
     +        35 * cos(1.84      +  2942.46*T0)
     +        33 * cos(0.24      +  7084.90*T0)
     +        32 * cos(0.18      +  5088.63*T0)
     +        32 * cos(1.78      +   398.15*T0)
     +        28 * cos(1.21      +  6286.60*T0)
     +        28 * cos(1.90      +  6279.55*T0)
     +        26 * cos(4.59      + 10447.39*T0);
  RV := RV / 1.0E+8;

  A := 103019 * cos(1.107490 +  6283.075850*T0)
     +   1721 * cos(1.0644   + 12566.1517*T0)
     +    702 * cos(3.142)
     +     32 * cos(1.02     + 18849.23*T0)
     +     31 * cos(2.84     +  5507.55*T0)
     +     25 * cos(1.32     +  5223.69*T0)
     +     18 * cos(1.42     +  1577.34*T0)
     +     10 * cos(5.91     + 10977.08*T0)
     +      9 * cos(1.42     +  6275.96*T0)
     +      9 * cos(0.27     +  5486.78*T0);
  RV := RV + (A * T0 / 1.0E+8);

  A := 4359 * cos(5.7846 +  6283.0758*T0)
     +  124 * cos(5.579  + 12566.152*T0)
     +   12 * cos(3.14)
     +    9 * cos(3.63   + 77713.77*T0)
     +    6 * cos(1.87   +  5573.14*T0)
     +    3 * cos(5.47   + 18849.23*T0);
  RV := RV + (A * sqr(T0) / 1.0E+8);

  L := (L + PI);
  L := Frac(L / 2.0 / PI) * 2.0 * Pi;
  if L < 0 then
    L := L + (2.0*PI);
  B := -B;

  TX := RV * cos(B) * cos(L);
  TY := RV * cos(B) * sin(L);
  TZ := RV * sin(B);

  Result.X :=               TX +     4.40360E-7 * TY - 1.90919E-7 * TZ;
  Result.Y := -4.79966E-7 * TX + 0.917482137087 * TY - 0.397776982902 * TZ;
  Result.Z :=                    0.397776982902 * TY + 0.917482137087 * TZ;
end;

{--------------------------------------------------------------------------}

function EclipticToRectangular(Longitude, Latitude,
                               RadiusVector : Double) : TStRectangularCord;
var
  var1,
  var2,
  var3 : Double;
begin
  var1 := RadiusVector * cos(Longitude) * cos(Latitude);
  var2 := RadiusVector * sin(Longitude) * cos(Latitude);
  var3 := RadiusVector * sin(Latitude);

  Result.X := var1;
  Result.Y := var2 * cos(OB2000) - var3 * sin(OB2000);
  Result.Z := var2 * sin(OB2000) + var3 * cos(OB2000);
end;

{--------------------------------------------------------------------------}

function RADec(Planet, Sun : TStRectangularCord;
               ComputeElong : Boolean) : TStPlanetsRec;
var
  var1,
  var2,
  var3,
  var4,
  var5       : Double;
begin
  FillChar(Result, SizeOf(TStPlanetsRec), #0);

  var1 := Sun.X + Planet.X;
  var2 := Sun.Y + Planet.Y;
  var3 := Sun.Z + Planet.Z;

  var4 := arctan(var2/var1);
  var4 := RealAngle(var2, var1, var4) * radcor;

  var5 := sqrt(sqr(var1) + sqr(var2) + sqr(var3));
  var3 := StInvsin(var3/var5) * radcor;

  Result.RA := var4;
  Result.DC := var3;

  var4 := Result.RA / radcor;
  var3 := Result.DC / radcor;

  if (ComputeElong) then begin
    var1 := sin(SunEQ.DC/radcor) * sin(var3);
    var2 := cos(SunEQ.DC/radcor) * cos(var3) * cos(SunEQ.RA/radcor - var4);
    Result.Elong := StInvcos(var1+var2) * radcor;
  end;
end;

{--------------------------------------------------------------------------}

function MercuryPosition(JD : Double) : TStPlanetsRec;
begin
  PlanEC := ComputeMercury(JD);
  PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
  Result := RADec(PlanRC, SunRC, True);
end;

{--------------------------------------------------------------------------}

function VenusPosition(JD : Double) : TStPlanetsRec;
begin
  PlanEC := ComputeVenus(JD);
  PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
  Result := RADec(PlanRC, SunRC, True);
end;

{--------------------------------------------------------------------------}

function MarsPosition(JD : Double) : TStPlanetsRec;
begin
  PlanEC := ComputeMars(JD);
  PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
  Result := RADec(PlanRC, SunRC, True);
end;

{--------------------------------------------------------------------------}

function JupiterPosition(JD : Double) : TStPlanetsRec;
begin
  PlanEC := ComputeJupiter(JD);
  PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
  Result := RADec(PlanRC, SunRC, True);
end;

{--------------------------------------------------------------------------}

function SaturnPosition(JD : Double) : TStPlanetsRec;
begin
  PlanEC := ComputeSaturn(JD);
  PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
  Result := RADec(PlanRC, SunRC, True);
end;

{--------------------------------------------------------------------------}

function UranusPosition(JD : Double) : TStPlanetsRec;
begin
  PlanEC := ComputeUranus(JD);
  PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
  Result := RADec(PlanRC, SunRC, True);
end;

{--------------------------------------------------------------------------}

function NeptunePosition(JD : Double) : TStPlanetsRec;
begin
  PlanEC := ComputeNeptune(JD);
  PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
  Result := RADec(PlanRC, SunRC, True);
end;

{--------------------------------------------------------------------------}

function PlutoPosition(JD : Double) : TStPlanetsRec;
begin
  PlanEC := ComputePluto(JD);
  PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
  Result := RADec(PlanRC, SunRC, True);
end;

{--------------------------------------------------------------------------}

procedure PlanetsPos(JD : Double; var PA : TStPlanetsArray);
var
  I   : Integer;
  Sun : TStRectangularCord;
begin
  {find Sun's Rectangular Coordinates}
  SunRC := SunofDate(JD);

  FillChar(SunEQ, SizeOf(TStPlanetsRec), #0);
  FillChar(Sun, SizeOf(TStRectangularCord), #0);

  {find Sun's RA/Dec}
  SunEQ := RADec(SunRC, Sun, False);
  PA[1] := PlutoPosition(JD);

  {find RA/Dec of each planet}
  for I := 1 to 8 do begin
    case I of
      1 : PA[I] := MercuryPosition(JD);
      2 : PA[I] := VenusPosition(JD);
      3 : PA[I] := MarsPosition(JD);
      4 : PA[I] := JupiterPosition(JD);
      5 : PA[I] := SaturnPosition(JD);
      6 : PA[I] := UranusPosition(JD);
      7 : PA[I] := NeptunePosition(JD);
      8 : PA[I] := PlutoPosition(JD);
    end;
  end;
end;


end.
1	// Upgraded to Delphi 2009: Sebastian Zierer
2
3	(* *** BEGIN LICENSE BLOCK ***
4	* Version: MPL 1.1
5	*
6	* The contents of this file are subject to the Mozilla Public License Version
7	* 1.1 (the "License"); you may not use this file except in compliance with
8	* the License. You may obtain a copy of the License at
9	* http://www.mozilla.org/MPL/
10	*
11	* Software distributed under the License is distributed on an "AS IS" basis,
12	* WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
13	* for the specific language governing rights and limitations under the
14	* License.
15	*
16	* The Original Code is TurboPower SysTools
17	*
18	* The Initial Developer of the Original Code is
19	* TurboPower Software
20	*
21	* Portions created by the Initial Developer are Copyright (C) 1996-2002
22	* the Initial Developer. All Rights Reserved.
23	*
24	* Contributor(s):
25	*
26	* *** END LICENSE BLOCK *** *)
27
28	{*********************************************************}
29	{* SysTools: StAstroP.pas 4.04 *}
30	{*********************************************************}
31	{* SysTools: Astronomical Routines (general Planetary) *}
32	{*********************************************************}
33
34	{$I StDefine.inc}
35
36	{ ************************************************************** }
37	{ Sources: }
38	{ 1. Astronomical Algorithms, Jean Meeus, Willmann-Bell, 1991. }
39	{ }
40	{ 2. Planetary and Lunar Coordinates (1984-2000), U.S. Govt, }
41	{ 1983. }
42	{ }
43	{ 3. Supplement to the American Ephemeris and Nautical Almanac,}
44	{ U.S. Govt, 1964. }
45	{ }
46	{ 4. MPO96-MPO98 source files, Brian D. Warner, 1995-2000. }
47	{ }
48	{ ************************************************************** }
49
50	unit StAstroP;
51
52	interface
53
54	const
55	StdDate = 2451545.0; {Ast. Julian Date for J2000 Epoch}
56	OB2000 = 0.409092804; {J2000 obliquity of the ecliptic (radians)}
57
58	type
59	TStEclipticalCord = packed record
60	L0,
61	B0,
62	R0 : Double;
63	end;
64
65	TStRectangularCord = packed record
66	X,
67	Y,
68	Z : Double;
69	end;
70
71	TStPlanetsRec = packed record
72	RA,
73	DC,
74	Elong : Double;
75	end;
76	TStPlanetsArray = array[1..8] of TStPlanetsRec;
77
78
79	procedure PlanetsPos(JD : Double; var PA : TStPlanetsArray);
80
81
82	implementation
83
84	uses
85	Windows,
86	StDate, StMerc, StVenus, StMars, StJup, StSaturn, StUranus, StNeptun,
87	StPluto, StMath;
88
89	var
90	PlanEC : TStEclipticalCord;
91	PlanRC,
92	SunRC : TStRectangularCord;
93	SunEQ : TStPlanetsRec;
94
95
96	{--------------------------------------------------------------------------}
97
98	function RealAngle(Value2, Value1, Start : Double) : Double;
99	begin
100	Result := Start;
101
102	if (Value1 = 0) then begin
103	if Value2 > 0 then
104	Result := Pi / 2.0
105	else
106	Result := 3.0 * Pi / 2.0;
107	end else begin
108	if (Value2 > 0.0) then begin
109	if (Value1 < 0.0) then
110	Result := Start + Pi
111	else
112	Result := Start;
113	end else begin
114	if (Value2 = 0) then begin
115	if Value1 > 0 then
116	Result := 0
117	else
118	Result := Pi;
119	end else begin
120	if (Value2 < 0) then begin
121	if (Value1 < 0) then
122	Result := Start + Pi
123	else
124	Result := Start + (2.0 * Pi)
125	end;
126	end;
127	end;
128	end;
129	end;
130
131	{--------------------------------------------------------------------------}
132
133	function SunOfDate(JD : Double) : TStRectangularCord;
134	{-compute J2000 XYZ coordinates of the Sun}
135	var
136	T0,
137	A,
138	L,
139	B,
140	RV,
141	TX,
142	TY,
143	TZ : Double;
144
145	begin
146	T0 := (JD - StdDate) / 365250;
147
148	{solar longitude}
149	L := 175347046
150	+ 3341656 * cos(4.6692568 + 6283.07585*T0)
151	+ 34894 * cos(4.6261000 + 12566.1517*T0)
152	+ 3497 * cos(2.7441000 + 5753.3849*T0)
153	+ 3418 * cos(2.8289000 + 3.5231*T0)
154	+ 3136 * cos(3.6277000 + 77713.7715*T0)
155	+ 2676 * cos(4.4181000 + 7860.4194*T0)
156	+ 2343 * cos(6.1352000 + 3930.2097*T0)
157	+ 1324 * cos(0.7425000 + 11506.7698*T0)
158	+ 1273 * cos(2.0371000 + 529.6910*T0)
159	+ 1199 * cos(1.1096000 + 1577.3435*T0)
160	+ 990 * cos(5.2330000 + 5884.9270*T0)
161	+ 902 * cos(2.0450000 + 26.1490*T0)
162	+ 857 * cos(3.5080000 + 398.149*T0)
163	+ 780 * cos(1.1790000 + 5223.694*T0)
164	+ 753 * cos(2.5330000 + 5507.553*T0)
165	+ 505 * cos(4.5830000 + 18849.228*T0)
166	+ 492 * cos(4.2050000 + 775.523*T0)
167	+ 357 * cos(2.9200000 + 0.067*T0)
168	+ 317 * cos(5.8490000 + 11790.626*T0)
169	+ 284 * cos(1.8990000 + 796.298*T0)
170	+ 271 * cos(0.3150000 + 10977.079*T0)
171	+ 243 * cos(0.3450000 + 5486.778*T0)
172	+ 206 * cos(4.8060000 + 2544.314*T0)
173	+ 205 * cos(1.8690000 + 5573.143*T0)
174	+ 202 * cos(2.4580000 + 6069.777*T0)
175	+ 156 * cos(0.8330000 + 213.299*T0)
176	+ 132 * cos(3.4110000 + 2942.463*T0)
177	+ 126 * cos(1.0830000 + 20.775*T0)
178	+ 115 * cos(0.6450000 + 0.980*T0)
179	+ 103 * cos(0.6360000 + 4694.003*T0)
180	+ 102 * cos(0.9760000 + 15720.839*T0)
181	+ 102 * cos(4.2670000 + 7.114*T0)
182	+ 99 * cos(6.2100000 + 2146.170*T0)
183	+ 98 * cos(0.6800000 + 155.420*T0)
184	+ 86 * cos(5.9800000 +161000.690*T0)
185	+ 85 * cos(1.3000000 + 6275.960*T0)
186	+ 85 * cos(3.6700000 + 71430.700*T0)
187	+ 80 * cos(1.8100000 + 17260.150*T0);
188
189	A := 628307584999.0
190	+ 206059 * cos(2.678235 + 6283.07585*T0)
191	+ 4303 * cos(2.635100 + 12566.1517*T0)
192	+ 425 * cos(1.590000 + 3.523*T0)
193	+ 119 * cos(5.796000 + 26.298*T0)
194	+ 109 * cos(2.966000 + 1577.344*T0)
195	+ 93 * cos(2.590000 + 18849.23*T0)
196	+ 72 * cos(1.140000 + 529.69*T0)
197	+ 68 * cos(1.870000 + 398.15*T0)
198	+ 67 * cos(4.410000 + 5507.55*T0)
199	+ 59 * cos(2.890000 + 5223.69*T0)
200	+ 56 * cos(2.170000 + 155.42*T0)
201	+ 45 * cos(0.400000 + 796.30*T0)
202	+ 36 * cos(0.470000 + 775.52*T0)
203	+ 29 * cos(2.650000 + 7.11*T0)
204	+ 21 * cos(5.340000 + 0.98*T0)
205	+ 19 * cos(1.850000 + 5486.78*T0)
206	+ 19 * cos(4.970000 + 213.30*T0)
207	+ 17 * cos(2.990000 + 6275.96*T0)
208	+ 16 * cos(0.030000 + 2544.31*T0);
209	L := L + (A * T0);
210
211	A := 8722 * cos(1.0725 + 6283.0758*T0)
212	+ 991 * cos(3.1416)
213	+ 295 * cos(0.437 + 12566.1520*T0)
214	+ 27 * cos(0.050 + 3.52*T0)
215	+ 16 * cos(5.190 + 26.30*T0)
216	+ 16 * cos(3.69 + 155.42*T0)
217	+ 9 * cos(0.30 + 18849.23*T0)
218	+ 9 * cos(2.06 + 77713.77*T0);
219	L := L + (A * sqr(T0));
220
221	A := 289 * cos(5.842 + 6283.076*T0)
222	+ 21 * cos(6.05 + 12566.15*T0)
223	+ 3 * cos(5.20 + 155.42*T0)
224	+ 3 * cos(3.14);
225	L := L + (A * sqr(T0) * T0);
226	L := L / 1.0E+8;
227
228
229	{solar latitude}
230	B := 280 * cos(3.199 + 84334.662*T0)
231	+ 102 * cos(5.422 + 5507.553*T0)
232	+ 80 * cos(3.88 + 5223.69*T0)
233	+ 44 * cos(3.70 + 2352.87*T0)
234	+ 32 * cos(4.00 + 1577.34*T0);
235	B := B / 1.0E+8;
236
237	A := 227778 * cos(3.413766 + 6283.07585*T0)
238	+ 3806 * cos(3.3706 + 12566.1517*T0)
239	+ 3620
240	+ 72 * cos(3.33 + 18849.23*T0)
241	+ 8 * cos(3.89 + 5507.55*T0)
242	+ 8 * cos(1.79 + 5223.69*T0)
243	+ 6 * cos(5.20 + 2352.87*T0);
244	B := B + (A * T0 / 1.0E+8);
245
246	A := 9721 * cos(5.1519 + 6283.07585*T0)
247	+ 233 * cos(3.1416)
248	+ 134 * cos(0.644 + 12566.152*T0)
249	+ 7 * cos(1.07 + 18849.23*T0);
250	B := B + (A * sqr(T0) / 1.0E+8);
251
252	A := 276 * cos(0.595 + 6283.076*T0)
253	+ 17 * cos(3.14)
254	+ 4 * cos(0.12 + 12566.15*T0);
255	B := B + (A * sqr(T0) * T0 / 1.0E+8);
256
257
258	{solar radius vector (astronomical units)}
259	RV := 100013989
260	+ 1670700 * cos(3.0984635 + 6283.07585*T0)
261	+ 13956 * cos(3.05525 + 12566.15170*T0)
262	+ 3084 * cos(5.1985 + 77713.7715*T0)
263	+ 1628 * cos(1.1739 + 5753.3849*T0)
264	+ 1576 * cos(2.8649 + 7860.4194*T0)
265	+ 925 * cos(5.453 + 11506.770*T0)
266	+ 542 * cos(4.564 + 3930.210*T0)
267	+ 472 * cos(3.661 + 5884.927*T0)
268	+ 346 * cos(0.964 + 5507.553*T0)
269	+ 329 * cos(5.900 + 5223.694*T0)
270	+ 307 * cos(0.299 + 5573.143*T0)
271	+ 243 * cos(4.273 + 11790.629*T0)
272	+ 212 * cos(5.847 + 1577.344*T0)
273	+ 186 * cos(5.022 + 10977.079*T0)
274	+ 175 * cos(3.012 + 18849.228*T0)
275	+ 110 * cos(5.055 + 5486.778*T0)
276	+ 98 * cos(0.89 + 6069.78*T0)
277	+ 86 * cos(5.69 + 15720.84*T0)
278	+ 86 * cos(1.27 +161000.69*T0)
279	+ 65 * cos(0.27 + 17260.15*T0)
280	+ 63 * cos(0.92 + 529.69*T0)
281	+ 57 * cos(2.01 + 83996.85*T0)
282	+ 56 * cos(5.24 + 71430.70*T0)
283	+ 49 * cos(3.25 + 2544.31*T0)
284	+ 47 * cos(2.58 + 775.52*T0)
285	+ 45 * cos(5.54 + 9437.76*T0)
286	+ 43 * cos(6.01 + 6275.96*T0)
287	+ 39 * cos(5.36 + 4694.00*T0)
288	+ 38 * cos(2.39 + 8827.39*T0)
289	+ 37 * cos(0.83 + 19651.05*T0)
290	+ 37 * cos(4.90 + 12139.55*T0)
291	+ 36 * cos(1.67 + 12036.46*T0)
292	+ 35 * cos(1.84 + 2942.46*T0)
293	+ 33 * cos(0.24 + 7084.90*T0)
294	+ 32 * cos(0.18 + 5088.63*T0)
295	+ 32 * cos(1.78 + 398.15*T0)
296	+ 28 * cos(1.21 + 6286.60*T0)
297	+ 28 * cos(1.90 + 6279.55*T0)
298	+ 26 * cos(4.59 + 10447.39*T0);
299	RV := RV / 1.0E+8;
300
301	A := 103019 * cos(1.107490 + 6283.075850*T0)
302	+ 1721 * cos(1.0644 + 12566.1517*T0)
303	+ 702 * cos(3.142)
304	+ 32 * cos(1.02 + 18849.23*T0)
305	+ 31 * cos(2.84 + 5507.55*T0)
306	+ 25 * cos(1.32 + 5223.69*T0)
307	+ 18 * cos(1.42 + 1577.34*T0)
308	+ 10 * cos(5.91 + 10977.08*T0)
309	+ 9 * cos(1.42 + 6275.96*T0)
310	+ 9 * cos(0.27 + 5486.78*T0);
311	RV := RV + (A * T0 / 1.0E+8);
312
313	A := 4359 * cos(5.7846 + 6283.0758*T0)
314	+ 124 * cos(5.579 + 12566.152*T0)
315	+ 12 * cos(3.14)
316	+ 9 * cos(3.63 + 77713.77*T0)
317	+ 6 * cos(1.87 + 5573.14*T0)
318	+ 3 * cos(5.47 + 18849.23*T0);
319	RV := RV + (A * sqr(T0) / 1.0E+8);
320
321	L := (L + PI);
322	L := Frac(L / 2.0 / PI) * 2.0 * Pi;
323	if L < 0 then
324	L := L + (2.0*PI);
325	B := -B;
326
327	TX := RV * cos(B) * cos(L);
328	TY := RV * cos(B) * sin(L);
329	TZ := RV * sin(B);
330
331	Result.X := TX + 4.40360E-7 * TY - 1.90919E-7 * TZ;
332	Result.Y := -4.79966E-7 * TX + 0.917482137087 * TY - 0.397776982902 * TZ;
333	Result.Z := 0.397776982902 * TY + 0.917482137087 * TZ;
334	end;
335
336	{--------------------------------------------------------------------------}
337
338	function EclipticToRectangular(Longitude, Latitude,
339	RadiusVector : Double) : TStRectangularCord;
340	var
341	var1,
342	var2,
343	var3 : Double;
344	begin
345	var1 := RadiusVector * cos(Longitude) * cos(Latitude);
346	var2 := RadiusVector * sin(Longitude) * cos(Latitude);
347	var3 := RadiusVector * sin(Latitude);
348
349	Result.X := var1;
350	Result.Y := var2 * cos(OB2000) - var3 * sin(OB2000);
351	Result.Z := var2 * sin(OB2000) + var3 * cos(OB2000);
352	end;
353
354	{--------------------------------------------------------------------------}
355
356	function RADec(Planet, Sun : TStRectangularCord;
357	ComputeElong : Boolean) : TStPlanetsRec;
358	var
359	var1,
360	var2,
361	var3,
362	var4,
363	var5 : Double;
364	begin
365	FillChar(Result, SizeOf(TStPlanetsRec), #0);
366
367	var1 := Sun.X + Planet.X;
368	var2 := Sun.Y + Planet.Y;
369	var3 := Sun.Z + Planet.Z;
370
371	var4 := arctan(var2/var1);
372	var4 := RealAngle(var2, var1, var4) * radcor;
373
374	var5 := sqrt(sqr(var1) + sqr(var2) + sqr(var3));
375	var3 := StInvsin(var3/var5) * radcor;
376
377	Result.RA := var4;
378	Result.DC := var3;
379
380	var4 := Result.RA / radcor;
381	var3 := Result.DC / radcor;
382
383	if (ComputeElong) then begin
384	var1 := sin(SunEQ.DC/radcor) * sin(var3);
385	var2 := cos(SunEQ.DC/radcor) * cos(var3) * cos(SunEQ.RA/radcor - var4);
386	Result.Elong := StInvcos(var1+var2) * radcor;
387	end;
388	end;
389
390	{--------------------------------------------------------------------------}
391
392	function MercuryPosition(JD : Double) : TStPlanetsRec;
393	begin
394	PlanEC := ComputeMercury(JD);
395	PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
396	Result := RADec(PlanRC, SunRC, True);
397	end;
398
399	{--------------------------------------------------------------------------}
400
401	function VenusPosition(JD : Double) : TStPlanetsRec;
402	begin
403	PlanEC := ComputeVenus(JD);
404	PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
405	Result := RADec(PlanRC, SunRC, True);
406	end;
407
408	{--------------------------------------------------------------------------}
409
410	function MarsPosition(JD : Double) : TStPlanetsRec;
411	begin
412	PlanEC := ComputeMars(JD);
413	PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
414	Result := RADec(PlanRC, SunRC, True);
415	end;
416
417	{--------------------------------------------------------------------------}
418
419	function JupiterPosition(JD : Double) : TStPlanetsRec;
420	begin
421	PlanEC := ComputeJupiter(JD);
422	PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
423	Result := RADec(PlanRC, SunRC, True);
424	end;
425
426	{--------------------------------------------------------------------------}
427
428	function SaturnPosition(JD : Double) : TStPlanetsRec;
429	begin
430	PlanEC := ComputeSaturn(JD);
431	PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
432	Result := RADec(PlanRC, SunRC, True);
433	end;
434
435	{--------------------------------------------------------------------------}
436
437	function UranusPosition(JD : Double) : TStPlanetsRec;
438	begin
439	PlanEC := ComputeUranus(JD);
440	PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
441	Result := RADec(PlanRC, SunRC, True);
442	end;
443
444	{--------------------------------------------------------------------------}
445
446	function NeptunePosition(JD : Double) : TStPlanetsRec;
447	begin
448	PlanEC := ComputeNeptune(JD);
449	PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
450	Result := RADec(PlanRC, SunRC, True);
451	end;
452
453	{--------------------------------------------------------------------------}
454
455	function PlutoPosition(JD : Double) : TStPlanetsRec;
456	begin
457	PlanEC := ComputePluto(JD);
458	PlanRC := EclipticToRectangular(PlanEC.L0, PlanEC.B0, PlanEC.R0);
459	Result := RADec(PlanRC, SunRC, True);
460	end;
461
462	{--------------------------------------------------------------------------}
463
464	procedure PlanetsPos(JD : Double; var PA : TStPlanetsArray);
465	var
466	I : Integer;
467	Sun : TStRectangularCord;
468	begin
469	{find Sun's Rectangular Coordinates}
470	SunRC := SunofDate(JD);
471
472	FillChar(SunEQ, SizeOf(TStPlanetsRec), #0);
473	FillChar(Sun, SizeOf(TStRectangularCord), #0);
474
475	{find Sun's RA/Dec}
476	SunEQ := RADec(SunRC, Sun, False);
477	PA[1] := PlutoPosition(JD);
478
479	{find RA/Dec of each planet}
480	for I := 1 to 8 do begin
481	case I of
482	1 : PA[I] := MercuryPosition(JD);
483	2 : PA[I] := VenusPosition(JD);
484	3 : PA[I] := MarsPosition(JD);
485	4 : PA[I] := JupiterPosition(JD);
486	5 : PA[I] := SaturnPosition(JD);
487	6 : PA[I] := UranusPosition(JD);
488	7 : PA[I] := NeptunePosition(JD);
489	8 : PA[I] := PlutoPosition(JD);
490	end;
491	end;
492	end;
493
494
495	end.