// 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) *}
{*********************************************************}
{$IFDEF FPC}
{$mode DELPHI}
{$ENDIF}
//{$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.