📄 stjupsat.pas
字号:
ArgJup := 172.74 + (0.00111588 * DateDif);
ArgJup := Frac(ArgJup/360.0) * 360.0;
if (ArgJup < 0) then
ArgJup := 360.0 + ArgJup;
ArgJup := ArgJup / radcor;
AnomE := 357.529 + (0.9856003 * DateDif);
AnomE := Frac(AnomE/360.0) * 360.0;
if (AnomE < 0) then
AnomE := 360.0 + AnomE;
AnomE := AnomE / radcor;
AnomJ := 20.020 + (0.0830853 * DateDif + (0.329 * sin(ArgJup)));
AnomJ := Frac(AnomJ/360.0) * 360.0;
if (AnomJ < 0) then
AnomJ := 360.0 + AnomJ;
AnomJ := AnomJ / radcor;
DeltaLong := 66.115 + (0.9025179 * DateDif - (0.329 * sin(ArgJup)));
DeltaLong := Frac(DeltaLong/360.0) * 360.0;
if (DeltaLong < 0) then
DeltaLong := 360.0 + DeltaLong;
DeltaLong := DeltaLong / radcor;
ECenterE := 1.915 * sin(AnomE) + 0.020 * sin(2*AnomE);
ECenterE := ECenterE / radcor;
ECenterJ := 5.555 * sin(AnomJ) + 0.168 * sin(2*AnomJ);
ECenterJ := ECenterJ / radcor;
K := (DeltaLong + ECenterE - ECenterJ);
RVE := 1.00014 - (0.01671 * cos(AnomE)) - (0.00014 * cos(2*AnomE));
RVJ := 5.20872 - (0.25208 * cos(AnomJ)) - (0.00611 * cos(2*AnomJ));
EJDist := sqrt(sqr(RVJ) + sqr(RVE) - (2 * RVJ * RVE * cos(K)));
Phase := RVE/EJDist * sin(K);
Phase := StInvSin(Phase);
if ((sin(K) < 0) and (Phase > 0)) or
((sin(K) > 0) and (Phase < 0)) then
Phase := -Phase;
Lambda := 34.35 + (0.083091 * DateDif) + (0.329 * sin(ArgJup));
Lambda := Lambda / radcor + ECenterJ;
DS := 3.12 * sin(Lambda + 42.8 / radcor);
DE := DS - 2.22 * sin(Phase) * cos(Lambda + 22/radcor)
- 1.30 * ((RVJ - EJDist) / EJDist) * sin(Lambda - 100.5/radcor);
DE := DE / radcor;
Mu1 := 163.8067 + 203.4058643 * (DateDif - (EJDist / 173));
Mu1 := Frac(Mu1/360.0) * 360.0;
if (Mu1 < 0) then
Mu1 := 360.0 + Mu1;
Mu1 := Mu1 / radcor + Phase - ECenterJ;
Mu2 := 358.4108 + 101.2916334 * (DateDif - (EJDist / 173));
Mu2 := Frac(Mu2/360.0) * 360.0;
if (Mu2 < 0) then
Mu2 := 360.0 + Mu2;
Mu2 := Mu2 / radcor + Phase - ECenterJ;
Mu3 := 5.7129 + 50.2345179 * (DateDif - (EJDist / 173));
Mu3 := Frac(Mu3/360.0) * 360.0;
if (Mu3 < 0) then
Mu3 := 360.0 + Mu3;
Mu3 := Mu3 / radcor + Phase - ECenterJ;
Mu4 := 224.8151 + 21.4879801 * (DateDif - (EJDist / 173));
Mu4 := Frac(Mu4/360.0) * 360.0;
if (Mu4 < 0) then
Mu4 := 360.0 + Mu4;
Mu4 := Mu4 / radcor + Phase - ECenterJ;
G := 331.18 + 50.310482 * (DateDif - (EJDist / 173));
G := Frac(G/360.0) * 360.0;
if (G < 0) then
G := 360.0 + G;
G := G / radcor;
H := 87.40 + 21.569231 * (DateDif - (EJDist / 173));
H := Frac(H/360.0) * 360.0;
if (H < 0) then
H := 360.0 + H;
H := H / radcor;
TmpDbl1 := 0.473 * sin(2 * (Mu1 - Mu2)) / radcor;
TmpDbl2 := 1.065 * sin(2 * (Mu2 - Mu3)) / radcor;
R1 := 5.9073 - 0.0244 * cos(2 * (Mu1 - Mu2));
R2 := 9.3991 - 0.0882 * cos(2 * (Mu2 - Mu3));
R3 := 14.9924 - 0.0216 * cos(G);
R4 := 26.3699 - 0.1935 * cos(H);
Mu1 := Mu1 + TmpDbl1;
Mu2 := Mu2 + TmpDbl2;
Mu3 := Mu3 + (0.165 * sin(G)) / radcor;
Mu4 := Mu4 + (0.841 * sin(H)) / radcor;
Result.Io.X := R1 * sin(Mu1);
Result.Io.Y := -R1 * cos(Mu1) * sin(DE);
Result.Europa.X := R2 * sin(Mu2);
Result.Europa.Y := -R2 * cos(Mu2) * sin(DE);
Result.Ganymede.X := R3 * sin(Mu3);
Result.Ganymede.Y := -R3 * cos(Mu3) * sin(DE);
Result.Callisto.X := R4 * sin(Mu4);
Result.Callisto.Y := -R4 * cos(Mu4) * sin(DE);
end;
{-------------------------------------------------------------------------}
function JupSatsHi(AJD : Double; Shadows : Boolean) : TStJupSats;
var
I : longint;
SunPos : SunCoordsRec;
STUT : TStDateTimeRec;
JupPos : TStEclipticalCord;
SatX : array[1..5] of Double;
SatY : array[1..5] of Double;
SatZ : array[1..5] of Double;
TD1,
TD2,
Angle, {Temporary Double values}
LTime, {Tau}
AJDT, {AJD adjusted for light time (Tau)}
JupX,
JupY,
JupZ, {Jupiter's geocentric rectangular coordinates}
EJDist, {Delta}
Jup1,
Jup2, {/\, Alpha}
DateDif, {t}
L1, L2,
L3, L4, {script L1-4}
Pi1, Pi2,
Pi3, Pi4, {Pi1-4}
W1, W2,
W3, W4, {Omega1-4}
Inequality, {upside down L}
PhiLambda,
NodeJup, {Psi}
AnomJup, {G}
AnomSat, {G'}
LongPerJ,
S1, S2,
S3, S4, {Sum1-4}
TL1, TL2,
TL3, TL4, {Capital L1-4}
B1, B2,
B3, B4, {tangent of latitude}
R1, R2,
R3, R4, {radius vector}
T0, {Julian Centuries}
Precession, {P}
Inclination {I}
: Double;
Transforms : array[1..5] of TranformRec;
begin
FillChar(Result, SizeOf(TStJupSats), #0);
AJD := DateTimeToAJD(AJD);
SunPos := SunCoords(AJD);
if not Shadows then begin
TD1 := 5;
AJDT := AJD - 0.0057755183 * TD1; {first guess}
repeat
JupPos := ComputeJupiter(AJDT);
JupX := JupPos.R0 * cos(JupPos.B0) * cos(JupPos.L0)
+ SunPos.R * cos(SunPos.L);
JupY := JupPos.R0 * cos(JupPos.B0) * sin(JupPos.L0)
+ SunPos.R * sin(SunPos.L);
JupZ := JupPos.R0 * sin(JupPos.B0);
EJDist := sqrt(sqr(JupX) + sqr(JupY) + sqr(JupZ));
TD2 := abs(EJDist - TD1);
if abs(TD2) > 0.0005 then begin
AJDT := AJD - 0.0057755183 * ((EJDist + TD1) / 2);
TD1 := EJDist;
end;
until (TD2 <= 0.0005);
end else begin
JupPos := ComputeJupiter(AJD);
JupX := JupPos.R0 * cos(JupPos.B0) * cos(JupPos.L0);
JupY := JupPos.R0 * cos(JupPos.B0) * sin(JupPos.L0);
JupZ := JupPos.R0 * sin(JupPos.B0);
EJDist := sqrt(sqr(JupX+SunPos.X) +
sqr(JupY+SunPos.Y) + sqr(JupZ+SunPos.Z));
end;
Jup1 := StInvTan2(JupX, JupY);
Jup2 := ArcTan(JupZ / sqrt(sqr(JupX) + sqr(JupY)));
DateDif := AJD - 2443000.5 - (0.0057755183 * EJDist);
L1 := 106.07947 + 203.488955432 * DateDif;
L1 := Frac(L1/360.0) * 360.0;
if (L1 < 0) then
L1 := 360.0 + L1;
L1 := L1 / radcor;
L2 := 175.72938 + 101.374724550 * DateDif;
L2 := Frac(L2/360.0) * 360.0;
if (L2 < 0) then
L2 := 360.0 + L2;
L2 := L2 / radcor;
L3 := 120.55434 + 50.317609110 * DateDif;
L3 := Frac(L3/360.0) * 360.0;
if (L3 < 0) then
L3 := 360.0 + L3;
L3 := L3 / radcor;
L4 := 84.44868 + 21.571071314 * DateDif;
L4 := Frac(L4/360.0) * 360.0;
if (L4 < 0) then
L4 := 360.0 + L4;
L4 := L4 / radcor;
Pi1 := 58.3329 + 0.16103936 * DateDif;
Pi1 := Frac(Pi1/360.0) * 360.0;
if (Pi1 < 0) then
Pi1 := 360.0 + Pi1;
Pi1 := Pi1 / radcor;
Pi2 := 132.8959 + 0.04647985 * DateDif;
Pi2 := Frac(Pi2/360.0) * 360.0;
if (Pi2 < 0) then
Pi2 := 360.0 + Pi2;
Pi2 := Pi2 / radcor;
Pi3 := 187.2887 + 0.00712740 * DateDif;
Pi3 := Frac(Pi3/360.0) * 360.0;
if (Pi3 < 0) then
Pi3 := 360.0 + Pi3;
Pi3 := Pi3 / radcor;
Pi4 := 335.3418 + 0.00183998 * DateDif;
Pi4 := Frac(Pi4/360.0) * 360.0;
if (Pi4 < 0) then
Pi4 := 360.0 + Pi4;
Pi4 := Pi4 / radcor;
W1 := 311.0793 - 0.13279430 * DateDif;
W1 := Frac(W1/360.0) * 360.0;
if (W1 < 0) then
W1 := 360.0 + W1;
W1 := W1 / radcor;
W2 := 100.5099 - 0.03263047 * DateDif;
W2 := Frac(W2/360.0) * 360.0;
if (W2 < 0) then
W2 := 360.0 + W2;
W2 := W2 / radcor;
W3 := 119.1688 - 0.00717704 * DateDif;
W3 := Frac(W3/360.0) * 360.0;
if (W3 < 0) then
W3 := 360.0 + W3;
W3 := W3 / radcor;
W4 := 322.5729 - 0.00175934 * DateDif;
W4 := Frac(W4/360.0) * 360.0;
if (W4 < 0) then
W4 := 360.0 + W4;
W4 := W4 / radcor;
Inequality := 0.33033 * sin((163.679 + 0.0010512*DateDif) / radcor)
+ 0.03439 * sin((34.486 - 0.0161731*DateDif) / radcor);
Inequality := Inequality / radcor;
PhiLambda := 191.8132 + 0.17390023 * DateDif;
PhiLambda := Frac(PhiLambda / 360.0) * 360.0;
if (PhiLambda < 0) then
PhiLambda := 360.0 + PhiLambda;
PhiLambda := PhiLambda / radcor;
NodeJup := 316.5182 - 0.00000208 * DateDif;
NodeJup := Frac(NodeJup / 360.0) * 360.0;
if (NodeJup < 0) then
NodeJup := 360.0 + NodeJup;
NodeJup := NodeJup / radcor;
AnomJup := 30.23756 + 0.0830925701 * DateDif;
AnomJup := Frac(AnomJup / 360.0) * 360.0;
if (AnomJup < 0) then
AnomJup := 360.0 + AnomJup;
AnomJup := AnomJup / radcor + Inequality;
AnomSat := 31.97853 + 0.0334597339 * DateDif;
AnomSat := Frac(AnomSat / 360.0) * 360.0;
if (AnomSat < 0) then
AnomSat := 360.0 + AnomSat;
AnomSat := AnomSat / radcor;
LongPerJ := 13.469942 / radcor;
S1 := 0.47259 * sin(2*(L1-L2))
- 0.03480 * sin(Pi3-Pi4)
- 0.01756 * sin(Pi1 + Pi3 - 2*LongPerJ - 2*AnomJup)
+ 0.01080 * sin(L2 - 2*L3 + Pi3)
+ 0.00757 * sin(PhiLambda)
+ 0.00663 * sin(L2 - 2*L3 + Pi4)
+ 0.00453 * sin(L1 - Pi3)
+ 0.00453 * sin(L2 - 2*L3 + Pi2)
- 0.00354 * sin(L1-L2)
- 0.00317 * sin(2*NodeJup - 2*LongPerJ)
- 0.00269 * sin(L2 - 2*L3 + Pi1)
+ 0.00263 * sin(L1 - Pi4)
+ 0.00186 * sin(L1 - Pi1)
- 0.00186 * sin(AnomJup)
+ 0.00167 * sin(Pi2 - Pi3)
+ 0.00158 * sin(4*(L1-L2))
- 0.00155 * sin(L1 - L3)
- 0.00142 * sin(NodeJup + W3 - 2*LongPerJ - 2*AnomJup)
- 0.00115 * sin(2*(L1 - 2*L2 + W2))
+ 0.00089 * sin(Pi2 - Pi4)
+ 0.00084 * sin(W2 - W3)
+ 0.00084 * sin(L1 + Pi3 - 2*LongPerJ - 2*AnomJup)
+ 0.00053 * sin(NodeJup - W2);
S2 := 1.06476 * sin(2*(L2-L3))
+ 0.04253 * sin(L1 - 2*L2 + Pi3)
+ 0.03579 * sin(L2 - Pi3)
+ 0.02383 * sin(L1 - 2*L2 + Pi4)
+ 0.01977 * sin(L2 - Pi4)
- 0.01843 * sin(PhiLambda)
+ 0.01299 * sin(Pi3 - Pi4)
- 0.01142 * sin(L2 - L3)
+ 0.01078 * sin(L2 - Pi2)
- 0.01058 * sin(AnomJup)
+ 0.00870 * sin(L2 - 2*L3 + Pi2)
- 0.00775 * sin(2*(NodeJup - LongPerJ))
+ 0.00524 * sin(2*(L1-L2))
- 0.00460 * sin(L1-L3)
+ 0.00450 * sin(L2 - 2*L3 + Pi1)
+ 0.00327 * sin(NodeJup - 2*AnomJup + W3 - 2*LongPerJ)
- 0.00296 * sin(Pi1 + Pi3 - 2*LongPerJ - 2*AnomJup)
- 0.00151 * sin(2*AnomJup)
+ 0.00146 * sin(NodeJup - W3)
+ 0.00125 * sin(NodeJup - W4)
- 0.00117 * sin(L1 -2*L3 + Pi3)
- 0.00095 * sin(2*(L2-W2))
+ 0.00086 * sin(2*(L1-2*L2 +W2))
- 0.00086 * sin(5*AnomSat - 2*AnomJup + 52.225/radcor)
- 0.00078 * sin(L2-L4)
- 0.00064 * sin(L1 - 2*L3 + Pi4)
- 0.00063 * sin(3*L3 - 7*L4 + 4*Pi4)
+ 0.00061 * sin(Pi1 - Pi4)
+ 0.00058 * sin(2*(NodeJup - LongPerJ - AnomJup))
+ 0.00058 * sin(W3 - W4)
+ 0.00056 * sin(2*(L2-L4))
+ 0.00055 * sin(2*(L1-L3))
+ 0.00052 * sin(3*L3 - 7*L4 + Pi3 + 3*Pi4)
- 0.00043 * sin(L1 - Pi3)
+ 0.00042 * sin(Pi3 - Pi2)
+ 0.00041 * sin(5*(L2-L3))
+ 0.00041 * sin(Pi4 - LongPerJ)
+ 0.00038 * sin(L2 - Pi1)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -