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	torben	2671	// 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.