📄 stjupsat.pas
字号:
(* ***** 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: StJupsat.pas 4.03 *}
{*********************************************************}
{* SysTools: Astronomical Routines *}
{* (for the four "Gallilean" moons of Jupiter *}
{* Callisto, Europa, Ganymede, and Io) *}
{*********************************************************}
{$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-98 source files, Brian D. Warner, 1995-98. }
{ }
{ ************************************************************** }
(* **************************************************************
The formulae in this unit are based on DYNAMICAL time, which is based
on the actual rotation of the Earth and is gradually slowing. Universal
Time is based on a fixed rotation rate. To directly compare results in
the astronomical literature (and Meeus' examples), you must use a Universal
Time that is adjusted by the value Delta T. This value is approximately 1 min
in the latter part of the 20th century and will be about 80 seconds in the
year 2020. As an example, to compare the High Precision positions for
1992 December 16 (Meeus' example), you must use a Universal Time of
1992 December 16 at 00:00:59 which equals 00:00:00 Dynamical Time
The Shadows parameter is used for high precision calculations only. If True,
the positions are calculated as seen from the SUN, not the Earth. For eclipses
of the satellites by Jupiter, in effect, the position is in reference to the
SHADOW of Jupiter in space, not the planet itself. For shadow transits, where
the shadow of the satellite is projected onto the surface of the planet, the
position is that of the satellite's shadow in reference to Jupiter and not
the satellite itself.
The purpose of the Shadows parameter is to complete the prediction for
satellite phenomenon. For example, using Shadow := False, the result may
indicate that a satellite goes behind Jupiter at a given instant but not
if the satellite is visible because it is in Jupiter's shadow. Setting
Shadow := True for the same time will indicate if the planet is in or out
of Jupiter's shadow.
(Shadow := FALSE) and (abs(satellite X-coordinate) = 1)
-------------------------------------------------------
If the X-value is negative and heading towards 0, the satellite is entering
the front of the planet.
If the X-value is negative and increasing, the satellite is coming from
behind the planet.
If the X-value is positive and heading towards 0, the satellite is going
behind the planet.
If the X-value is positive and increasing, the satellite is leaving
the front of the planet.
(Shadow := TRUE) and (abs(satellite X-coordinate) = 1)
-------------------------------------------------------
If the X-value is negative and heading towards 0, the satellite's shadow is
entering the planet's disc.
If the X-value is negative and increasing, the satellite is leaving the
planet's shadow.
If the X-value is positive and heading towards 0, the satellite entering
the planet's shadow.
If the X-value is positive and increasing, the satellite's shadow is
leaving the planet.
The X and Y coordinates are based on the equatorial radius of Jupiter. Because
the planet is considerably flattened by its rapid rotation, the polar diameter
is less than 1. To avoid dealing with an elliptical disc for Jupiter, multiply
the Y values only by 1.071374. This creates a "circular Jupiter" and so makes
determining if the satellite is above or below Jupiter easier
(abs(Y-coordinate) > 1).
****************************************************************** *)
unit StJupsat;
interface
type
TStJupSatPos = packed record
X : Double;
Y : Double;
end;
TStJupSats = packed record
Io : TStJupSatPos;
Europa : TStJupSatPos;
Ganymede : TStJupSatPos;
Callisto : TStJupSatPos;
end;
function GetJupSats(JD : TDateTime; HighPrecision, Shadows : Boolean) : TStJupSats;
implementation
uses
StDate, StAstro, StAstroP, StJup, StMath;
type
SunCoordsRec = packed record
X, Y, Z : Double;
L, B, R : Double;
end;
TranformRec = packed record
A, B, C : array[1..6] of Double;
end;
function SunCoords(JD : Double) : SunCoordsRec;
var
L, B, R,
T0, TM,
RS,
OB, A : Double;
begin
T0 := (JD - StdDate) / 365250;
TM := T0/100;
RS := radcor * 3600;
OB := 0.4090928042223
- 4680.93/RS * TM
- 1.55/RS * sqr(TM)
+ 1999.25/RS * sqr(TM) * TM
- 51.38/RS * sqr(sqr(TM))
- 249.67/RS * sqr(sqr(TM)) * TM
- 39.05/RS * sqr(sqr(TM)) * sqr(TM)
+ 7.12/RS * sqr(sqr(TM)) * sqr(TM) * TM
+ 27.87/RS * sqr(sqr(sqr(TM)));
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 := 628331966747.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 := 52919
+ 8720 * cos(1.0721 + 6283.0758*T0)
+ 309 * cos(0.867 + 12566.152*T0)
+ 27 * cos(0.050 + 3.52*T0)
+ 16 * cos(5.190 + 26.30*T0)
+ 16 * cos(3.68 + 155.42*T0)
+ 10 * cos(0.76 + 18849.23*T0)
+ 9 * cos(2.06 + 77713.77*T0)
+ 7 * cos(0.83 + 775.52*T0)
+ 5 * cos(4.66 + 1577.34*T0);
L := L + (A * sqr(T0));
A := 289 * cos(5.844 + 6283.076*T0)
+ 35
+ 17 * cos(5.49 + 12566.15*T0)
+ 3 * cos(5.20 + 155.42*T0)
+ 1 * cos(4.72 + 3.52*T0);
L := L + (A * sqr(T0) * T0);
A := 114 * cos(3.142);
L := L + (A * sqr(sqr(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);
A := 9 * cos(3.90 + 5507.550*T0)
+ 6 * cos(1.73 + 5223.690*T0);
B := B + (A * T0);
B := B / 1.0E+8;
{solar radius vector (astronomical units)}
R := 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);
R := R / 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);
R := R + (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);
R := R + (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;
Result.L := L;
Result.B := B;
Result.R := R;
Result.X := R * cos(B) * cos(L);
Result.Y := R * (cos(B) * sin(L) * cos(OB) - sin(B) * sin(OB));
Result.Z := R * (cos(B) * sin(L) * sin(OB) + sin(B) * cos(OB));
end;
{-------------------------------------------------------------------------}
function JupSatsLo(AJD : Double) : TStJupSats;
var
DateDif, {d}
ArgJup, {V}
AnomE, {M}
AnomJ, {N}
DeltaLong, {J}
ECenterE, {A}
ECenterJ, {B}
K,
RVE, {R}
RVJ, {r}
EJDist, {Delta}
Phase, {Psi}
Lambda, {Lambda}
DS, DE, {DS, DE}
Mu1, Mu2,
Mu3, Mu4, {Mu1 - Mu4}
G, H, {G, H}
TmpDbl1,
TmpDbl2,
R1, R2,
R3, R4 {R1 - R4}
: Double;
begin
AJD := DateTimeToAJD(AJD);
DateDif := AJD - 2451545.0;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -