1
0
Files
applications
bindings
components
Comba_Animation
aboutcomponent
acs
beepfp
callite
chelper
chemtext
cmdline
cmdlinecfg
colorpalette
cryptini
csvdocument
epiktimer
everettrandom
examplecomponent
extrasyn
fpexif
fpsound
fpspreadsheet
fractions
freetypepascal
geckoport
gradcontrols
grid_semaphor
industrialstuff
iosdesigner
iphonelazext
jujiboutils
jvcllaz
kcontrols
lazautoupdate
lazbarcodes
lazmapviewer
lclextensions
longtimer
manualdock
mbColorLib
mplayer
multithreadprocs
nvidia-widgets
onguard
orpheus
playsoundpackage
poweredby
powerpdf
rgbgraphics
richmemo
richview
rtfview
rx
scrolltext
smnetgradient
spktoolbar
svn
systools
examples
images
source
db
general
design
run
st2dbarc.pas
stastro.pas
stastrop.pas
stbarc.pas
stbarpn.pas
stbase.pas
stbcd.pas
stbits.pas
stccy.dat
stccycnv.dat
stcoll.pas
stconst.pas
stcrc.pas
stdate.pas
stdatest.pas
stdecmth.pas
stdict.pas
stdque.pas
steclpse.pas
stexpr.pas
stexpr.txt
stfin.pas
sthash.pas
stinistm.pas
stjup.pas
stjupsat.pas
stlarr.pas
stlist.pas
stmars.pas
stmath.pas
stmerc.pas
stmerge.pas
stmoney.pas
stneptun.pas
stnvbits.pas
stnvcoll.pas
stnvcont.pas
stnvdict.pas
stnvdq.pas
stnvlary.pas
stnvlist.pas
stnvlmat.pas
stnvscol.pas
stnvtree.pas
stpluto.pas
stpqueue.pas
stptrns.pas
strandom.pas
stregex.pas
stsaturn.pas
ststat.pas
ststrl.pas
ststrms.pas
ststrs.pas
sttext.pas
sttohtml.pas
sttree.pas
sttxtdat.pas
sturanus.pas
stutils.pas
stvarr.pas
stvenus.pas
include
windows_only
laz_systools.lpk
laz_systools.pas
laz_systools_all.lpg
laz_systools_design.lpk
laz_systools_design.pas
laz_systoolsdb.lpk
laz_systoolsdb.pas
laz_systoolsdb_design.lpk
laz_systoolsdb_design.pas
laz_systoolswin.lpk
laz_systoolswin.pas
laz_systoolswin_design.lpk
laz_systoolswin_design.pas
readme-orig.txt
readme.txt
readme404pre.txt
tdi
thtmlport
tparadoxdataset
tvplanit
xdev_toolkit
zlibar
zmsql
examples
image_sources
lclbindings
wst
lazarus-ccr/components/systools/source/general/run/stastrop.pas

500 lines
15 KiB
ObjectPascal
Raw Normal View History

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