/[projects]/dao/DelphiScanner/Components/tpsystools_4.04/source/StJupSat.pas
ViewVC logotype

Annotation of /dao/DelphiScanner/Components/tpsystools_4.04/source/StJupSat.pas

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2671 - (hide annotations) (download)
Tue Aug 25 18:15:15 2015 UTC (8 years, 10 months ago) by torben
File size: 39121 byte(s)
Added tpsystools component
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: StJupsat.pas 4.04 *}
30     {*********************************************************}
31     {* SysTools: Astronomical Routines *}
32     {* (for the four "Gallilean" moons of Jupiter *}
33     {* Callisto, Europa, Ganymede, and Io) *}
34     {*********************************************************}
35    
36     {$I StDefine.inc}
37    
38     { ************************************************************** }
39     { Sources: }
40     { 1. Astronomical Algorithms, Jean Meeus, Willmann-Bell, 1991. }
41     { }
42     { 2. Planetary and Lunar Coordinates (1984-2000), U.S. Govt, }
43     { 1983. }
44     { }
45     { 3. Supplement to the American Ephemeris and Nautical Almanac,}
46     { U.S. Govt, 1964. }
47     { }
48     { 4. MPO96-98 source files, Brian D. Warner, 1995-98. }
49     { }
50     { ************************************************************** }
51    
52     (* **************************************************************
53    
54     The formulae in this unit are based on DYNAMICAL time, which is based
55     on the actual rotation of the Earth and is gradually slowing. Universal
56     Time is based on a fixed rotation rate. To directly compare results in
57     the astronomical literature (and Meeus' examples), you must use a Universal
58     Time that is adjusted by the value Delta T. This value is approximately 1 min
59     in the latter part of the 20th century and will be about 80 seconds in the
60     year 2020. As an example, to compare the High Precision positions for
61     1992 December 16 (Meeus' example), you must use a Universal Time of
62     1992 December 16 at 00:00:59 which equals 00:00:00 Dynamical Time
63    
64    
65     The Shadows parameter is used for high precision calculations only. If True,
66     the positions are calculated as seen from the SUN, not the Earth. For eclipses
67     of the satellites by Jupiter, in effect, the position is in reference to the
68     SHADOW of Jupiter in space, not the planet itself. For shadow transits, where
69     the shadow of the satellite is projected onto the surface of the planet, the
70     position is that of the satellite's shadow in reference to Jupiter and not
71     the satellite itself.
72    
73     The purpose of the Shadows parameter is to complete the prediction for
74     satellite phenomenon. For example, using Shadow := False, the result may
75     indicate that a satellite goes behind Jupiter at a given instant but not
76     if the satellite is visible because it is in Jupiter's shadow. Setting
77     Shadow := True for the same time will indicate if the planet is in or out
78     of Jupiter's shadow.
79    
80     (Shadow := FALSE) and (abs(satellite X-coordinate) = 1)
81     -------------------------------------------------------
82     If the X-value is negative and heading towards 0, the satellite is entering
83     the front of the planet.
84     If the X-value is negative and increasing, the satellite is coming from
85     behind the planet.
86     If the X-value is positive and heading towards 0, the satellite is going
87     behind the planet.
88     If the X-value is positive and increasing, the satellite is leaving
89     the front of the planet.
90    
91     (Shadow := TRUE) and (abs(satellite X-coordinate) = 1)
92     -------------------------------------------------------
93     If the X-value is negative and heading towards 0, the satellite's shadow is
94     entering the planet's disc.
95     If the X-value is negative and increasing, the satellite is leaving the
96     planet's shadow.
97     If the X-value is positive and heading towards 0, the satellite entering
98     the planet's shadow.
99     If the X-value is positive and increasing, the satellite's shadow is
100     leaving the planet.
101    
102     The X and Y coordinates are based on the equatorial radius of Jupiter. Because
103     the planet is considerably flattened by its rapid rotation, the polar diameter
104     is less than 1. To avoid dealing with an elliptical disc for Jupiter, multiply
105     the Y values only by 1.071374. This creates a "circular Jupiter" and so makes
106     determining if the satellite is above or below Jupiter easier
107     (abs(Y-coordinate) > 1).
108    
109     ****************************************************************** *)
110    
111     unit StJupsat;
112    
113     interface
114    
115     type
116     TStJupSatPos = packed record
117     X : Double;
118     Y : Double;
119     end;
120    
121     TStJupSats = packed record
122     Io : TStJupSatPos;
123     Europa : TStJupSatPos;
124     Ganymede : TStJupSatPos;
125     Callisto : TStJupSatPos;
126     end;
127    
128     function GetJupSats(JD : TDateTime; HighPrecision, Shadows : Boolean) : TStJupSats;
129    
130     implementation
131     uses
132     StDate, StAstro, StAstroP, StJup, StMath;
133    
134     type
135     SunCoordsRec = packed record
136     X, Y, Z : Double;
137     L, B, R : Double;
138     end;
139    
140     TranformRec = packed record
141     A, B, C : array[1..6] of Double;
142     end;
143    
144    
145     function SunCoords(JD : Double) : SunCoordsRec;
146     var
147     L, B, R,
148     T0, TM,
149     RS,
150     OB, A : Double;
151     begin
152     T0 := (JD - StdDate) / 365250;
153     TM := T0/100;
154     RS := radcor * 3600;
155     OB := 0.4090928042223
156     - 4680.93/RS * TM
157     - 1.55/RS * sqr(TM)
158     + 1999.25/RS * sqr(TM) * TM
159     - 51.38/RS * sqr(sqr(TM))
160     - 249.67/RS * sqr(sqr(TM)) * TM
161     - 39.05/RS * sqr(sqr(TM)) * sqr(TM)
162     + 7.12/RS * sqr(sqr(TM)) * sqr(TM) * TM
163     + 27.87/RS * sqr(sqr(sqr(TM)));
164    
165     L := 175347046
166     + 3341656 * cos(4.6692568 + 6283.07585*T0)
167     + 34894 * cos(4.6261000 + 12566.1517*T0)
168     + 3497 * cos(2.7441000 + 5753.3849*T0)
169     + 3418 * cos(2.8289000 + 3.5231*T0)
170     + 3136 * cos(3.6277000 + 77713.7715*T0)
171     + 2676 * cos(4.4181000 + 7860.4194*T0)
172     + 2343 * cos(6.1352000 + 3930.2097*T0)
173     + 1324 * cos(0.7425000 + 11506.7698*T0)
174     + 1273 * cos(2.0371000 + 529.6910*T0)
175     + 1199 * cos(1.1096000 + 1577.3435*T0)
176     + 990 * cos(5.2330000 + 5884.9270*T0)
177     + 902 * cos(2.0450000 + 26.1490*T0)
178     + 857 * cos(3.5080000 + 398.149*T0)
179     + 780 * cos(1.1790000 + 5223.694*T0)
180     + 753 * cos(2.5330000 + 5507.553*T0)
181     + 505 * cos(4.5830000 + 18849.228*T0)
182     + 492 * cos(4.2050000 + 775.523*T0)
183     + 357 * cos(2.9200000 + 0.067*T0)
184     + 317 * cos(5.8490000 + 11790.626*T0)
185     + 284 * cos(1.8990000 + 796.298*T0)
186     + 271 * cos(0.3150000 + 10977.079*T0)
187     + 243 * cos(0.3450000 + 5486.778*T0)
188     + 206 * cos(4.8060000 + 2544.314*T0)
189     + 205 * cos(1.8690000 + 5573.143*T0)
190     + 202 * cos(2.4580000 + 6069.777*T0)
191     + 156 * cos(0.8330000 + 213.299*T0)
192     + 132 * cos(3.4110000 + 2942.463*T0)
193     + 126 * cos(1.0830000 + 20.775*T0)
194     + 115 * cos(0.6450000 + 0.980*T0)
195     + 103 * cos(0.6360000 + 4694.003*T0)
196     + 102 * cos(0.9760000 + 15720.839*T0)
197     + 102 * cos(4.2670000 + 7.114*T0)
198     + 99 * cos(6.2100000 + 2146.170*T0)
199     + 98 * cos(0.6800000 + 155.420*T0)
200     + 86 * cos(5.9800000 +161000.690*T0)
201     + 85 * cos(1.3000000 + 6275.960*T0)
202     + 85 * cos(3.6700000 + 71430.700*T0)
203     + 80 * cos(1.8100000 + 17260.150*T0);
204    
205     A := 628331966747.0
206     + 206059 * cos(2.678235 + 6283.07585*T0)
207     + 4303 * cos(2.635100 + 12566.1517*T0)
208     + 425 * cos(1.590000 + 3.523*T0)
209     + 119 * cos(5.796000 + 26.298*T0)
210     + 109 * cos(2.966000 + 1577.344*T0)
211     + 93 * cos(2.590000 + 18849.23*T0)
212     + 72 * cos(1.140000 + 529.69*T0)
213     + 68 * cos(1.870000 + 398.15*T0)
214     + 67 * cos(4.410000 + 5507.55*T0)
215     + 59 * cos(2.890000 + 5223.69*T0)
216     + 56 * cos(2.170000 + 155.42*T0)
217     + 45 * cos(0.400000 + 796.30*T0)
218     + 36 * cos(0.470000 + 775.52*T0)
219     + 29 * cos(2.650000 + 7.11*T0)
220     + 21 * cos(5.340000 + 0.98*T0)
221     + 19 * cos(1.850000 + 5486.78*T0)
222     + 19 * cos(4.970000 + 213.30*T0)
223     + 17 * cos(2.990000 + 6275.96*T0)
224     + 16 * cos(0.030000 + 2544.31*T0);
225     L := L + (A * T0);
226    
227     A := 52919
228     + 8720 * cos(1.0721 + 6283.0758*T0)
229     + 309 * cos(0.867 + 12566.152*T0)
230     + 27 * cos(0.050 + 3.52*T0)
231     + 16 * cos(5.190 + 26.30*T0)
232     + 16 * cos(3.68 + 155.42*T0)
233     + 10 * cos(0.76 + 18849.23*T0)
234     + 9 * cos(2.06 + 77713.77*T0)
235     + 7 * cos(0.83 + 775.52*T0)
236     + 5 * cos(4.66 + 1577.34*T0);
237     L := L + (A * sqr(T0));
238    
239     A := 289 * cos(5.844 + 6283.076*T0)
240     + 35
241     + 17 * cos(5.49 + 12566.15*T0)
242     + 3 * cos(5.20 + 155.42*T0)
243     + 1 * cos(4.72 + 3.52*T0);
244     L := L + (A * sqr(T0) * T0);
245    
246     A := 114 * cos(3.142);
247     L := L + (A * sqr(sqr(T0)));
248     L := L / 1.0E+8;
249    
250     {solar latitude}
251     B := 280 * cos(3.199 + 84334.662*T0)
252     + 102 * cos(5.422 + 5507.553*T0)
253     + 80 * cos(3.88 + 5223.69*T0)
254     + 44 * cos(3.70 + 2352.87*T0)
255     + 32 * cos(4.00 + 1577.34*T0);
256    
257     A := 9 * cos(3.90 + 5507.550*T0)
258     + 6 * cos(1.73 + 5223.690*T0);
259     B := B + (A * T0);
260     B := B / 1.0E+8;
261    
262    
263     {solar radius vector (astronomical units)}
264     R := 100013989
265     + 1670700 * cos(3.0984635 + 6283.07585*T0)
266     + 13956 * cos(3.05525 + 12566.15170*T0)
267     + 3084 * cos(5.1985 + 77713.7715*T0)
268     + 1628 * cos(1.1739 + 5753.3849*T0)
269     + 1576 * cos(2.8649 + 7860.4194*T0)
270     + 925 * cos(5.453 + 11506.770*T0)
271     + 542 * cos(4.564 + 3930.210*T0)
272     + 472 * cos(3.661 + 5884.927*T0)
273     + 346 * cos(0.964 + 5507.553*T0)
274     + 329 * cos(5.900 + 5223.694*T0)
275     + 307 * cos(0.299 + 5573.143*T0)
276     + 243 * cos(4.273 + 11790.629*T0)
277     + 212 * cos(5.847 + 1577.344*T0)
278     + 186 * cos(5.022 + 10977.079*T0)
279     + 175 * cos(3.012 + 18849.228*T0)
280     + 110 * cos(5.055 + 5486.778*T0)
281     + 98 * cos(0.89 + 6069.78*T0)
282     + 86 * cos(5.69 + 15720.84*T0)
283     + 86 * cos(1.27 +161000.69*T0)
284     + 65 * cos(0.27 + 17260.15*T0)
285     + 63 * cos(0.92 + 529.69*T0)
286     + 57 * cos(2.01 + 83996.85*T0)
287     + 56 * cos(5.24 + 71430.70*T0)
288     + 49 * cos(3.25 + 2544.31*T0)
289     + 47 * cos(2.58 + 775.52*T0)
290     + 45 * cos(5.54 + 9437.76*T0)
291     + 43 * cos(6.01 + 6275.96*T0)
292     + 39 * cos(5.36 + 4694.00*T0)
293     + 38 * cos(2.39 + 8827.39*T0)
294     + 37 * cos(0.83 + 19651.05*T0)
295     + 37 * cos(4.90 + 12139.55*T0)
296     + 36 * cos(1.67 + 12036.46*T0)
297     + 35 * cos(1.84 + 2942.46*T0)
298     + 33 * cos(0.24 + 7084.90*T0)
299     + 32 * cos(0.18 + 5088.63*T0)
300     + 32 * cos(1.78 + 398.15*T0)
301     + 28 * cos(1.21 + 6286.60*T0)
302     + 28 * cos(1.90 + 6279.55*T0)
303     + 26 * cos(4.59 + 10447.39*T0);
304     R := R / 1.0E+8;
305    
306     A := 103019 * cos(1.107490 + 6283.075850*T0)
307     + 1721 * cos(1.0644 + 12566.1517*T0)
308     + 702 * cos(3.142)
309     + 32 * cos(1.02 + 18849.23*T0)
310     + 31 * cos(2.84 + 5507.55*T0)
311     + 25 * cos(1.32 + 5223.69*T0)
312     + 18 * cos(1.42 + 1577.34*T0)
313     + 10 * cos(5.91 + 10977.08*T0)
314     + 9 * cos(1.42 + 6275.96*T0)
315     + 9 * cos(0.27 + 5486.78*T0);
316     R := R + (A * T0 / 1.0E+8);
317    
318     A := 4359 * cos(5.7846 + 6283.0758*T0)
319     + 124 * cos(5.579 + 12566.152*T0)
320     + 12 * cos(3.14)
321     + 9 * cos(3.63 + 77713.77*T0)
322     + 6 * cos(1.87 + 5573.14*T0)
323     + 3 * cos(5.47 + 18849.23*T0);
324     R := R + (A * sqr(T0) / 1.0E+8);
325    
326     L := (L + PI);
327     L := Frac(L / 2.0 / PI) * 2.0 * Pi;
328     if L < 0 then
329     L := L + (2.0*PI);
330     B := -B;
331    
332     Result.L := L;
333     Result.B := B;
334     Result.R := R;
335     Result.X := R * cos(B) * cos(L);
336     Result.Y := R * (cos(B) * sin(L) * cos(OB) - sin(B) * sin(OB));
337     Result.Z := R * (cos(B) * sin(L) * sin(OB) + sin(B) * cos(OB));
338     end;
339    
340    
341     {-------------------------------------------------------------------------}
342    
343     function JupSatsLo(AJD : Double) : TStJupSats;
344     var
345     DateDif, {d}
346     ArgJup, {V}
347     AnomE, {M}
348     AnomJ, {N}
349     DeltaLong, {J}
350     ECenterE, {A}
351     ECenterJ, {B}
352     K,
353     RVE, {R}
354     RVJ, {r}
355     EJDist, {Delta}
356     Phase, {Psi}
357     Lambda, {Lambda}
358     DS, DE, {DS, DE}
359     Mu1, Mu2,
360     Mu3, Mu4, {Mu1 - Mu4}
361     G, H, {G, H}
362     TmpDbl1,
363     TmpDbl2,
364     R1, R2,
365     R3, R4 {R1 - R4}
366     : Double;
367    
368     begin
369     AJD := DateTimeToAJD(AJD);
370     DateDif := AJD - 2451545.0;
371     ArgJup := 172.74 + (0.00111588 * DateDif);
372     ArgJup := Frac(ArgJup/360.0) * 360.0;
373     if (ArgJup < 0) then
374     ArgJup := 360.0 + ArgJup;
375     ArgJup := ArgJup / radcor;
376    
377     AnomE := 357.529 + (0.9856003 * DateDif);
378     AnomE := Frac(AnomE/360.0) * 360.0;
379     if (AnomE < 0) then
380     AnomE := 360.0 + AnomE;
381     AnomE := AnomE / radcor;
382    
383     AnomJ := 20.020 + (0.0830853 * DateDif + (0.329 * sin(ArgJup)));
384     AnomJ := Frac(AnomJ/360.0) * 360.0;
385     if (AnomJ < 0) then
386     AnomJ := 360.0 + AnomJ;
387     AnomJ := AnomJ / radcor;
388    
389     DeltaLong := 66.115 + (0.9025179 * DateDif - (0.329 * sin(ArgJup)));
390     DeltaLong := Frac(DeltaLong/360.0) * 360.0;
391     if (DeltaLong < 0) then
392     DeltaLong := 360.0 + DeltaLong;
393     DeltaLong := DeltaLong / radcor;
394    
395     ECenterE := 1.915 * sin(AnomE) + 0.020 * sin(2*AnomE);
396     ECenterE := ECenterE / radcor;
397    
398     ECenterJ := 5.555 * sin(AnomJ) + 0.168 * sin(2*AnomJ);
399     ECenterJ := ECenterJ / radcor;
400    
401     K := (DeltaLong + ECenterE - ECenterJ);
402    
403     RVE := 1.00014 - (0.01671 * cos(AnomE)) - (0.00014 * cos(2*AnomE));
404     RVJ := 5.20872 - (0.25208 * cos(AnomJ)) - (0.00611 * cos(2*AnomJ));
405    
406     EJDist := sqrt(sqr(RVJ) + sqr(RVE) - (2 * RVJ * RVE * cos(K)));
407    
408     Phase := RVE/EJDist * sin(K);
409     Phase := StInvSin(Phase);
410    
411     if ((sin(K) < 0) and (Phase > 0)) or
412     ((sin(K) > 0) and (Phase < 0)) then
413     Phase := -Phase;
414    
415     Lambda := 34.35 + (0.083091 * DateDif) + (0.329 * sin(ArgJup));
416     Lambda := Lambda / radcor + ECenterJ;
417    
418     DS := 3.12 * sin(Lambda + 42.8 / radcor);
419     DE := DS - 2.22 * sin(Phase) * cos(Lambda + 22/radcor)
420     - 1.30 * ((RVJ - EJDist) / EJDist) * sin(Lambda - 100.5/radcor);
421     DE := DE / radcor;
422    
423     Mu1 := 163.8067 + 203.4058643 * (DateDif - (EJDist / 173));
424     Mu1 := Frac(Mu1/360.0) * 360.0;
425     if (Mu1 < 0) then
426     Mu1 := 360.0 + Mu1;
427     Mu1 := Mu1 / radcor + Phase - ECenterJ;
428    
429     Mu2 := 358.4108 + 101.2916334 * (DateDif - (EJDist / 173));
430     Mu2 := Frac(Mu2/360.0) * 360.0;
431     if (Mu2 < 0) then
432     Mu2 := 360.0 + Mu2;
433     Mu2 := Mu2 / radcor + Phase - ECenterJ;
434    
435     Mu3 := 5.7129 + 50.2345179 * (DateDif - (EJDist / 173));
436     Mu3 := Frac(Mu3/360.0) * 360.0;
437     if (Mu3 < 0) then
438     Mu3 := 360.0 + Mu3;
439     Mu3 := Mu3 / radcor + Phase - ECenterJ;
440    
441     Mu4 := 224.8151 + 21.4879801 * (DateDif - (EJDist / 173));
442     Mu4 := Frac(Mu4/360.0) * 360.0;
443     if (Mu4 < 0) then
444     Mu4 := 360.0 + Mu4;
445     Mu4 := Mu4 / radcor + Phase - ECenterJ;
446    
447     G := 331.18 + 50.310482 * (DateDif - (EJDist / 173));
448     G := Frac(G/360.0) * 360.0;
449     if (G < 0) then
450     G := 360.0 + G;
451     G := G / radcor;
452     H := 87.40 + 21.569231 * (DateDif - (EJDist / 173));
453     H := Frac(H/360.0) * 360.0;
454     if (H < 0) then
455     H := 360.0 + H;
456     H := H / radcor;
457    
458     TmpDbl1 := 0.473 * sin(2 * (Mu1 - Mu2)) / radcor;
459     TmpDbl2 := 1.065 * sin(2 * (Mu2 - Mu3)) / radcor;
460    
461     R1 := 5.9073 - 0.0244 * cos(2 * (Mu1 - Mu2));
462     R2 := 9.3991 - 0.0882 * cos(2 * (Mu2 - Mu3));
463     R3 := 14.9924 - 0.0216 * cos(G);
464     R4 := 26.3699 - 0.1935 * cos(H);
465    
466     Mu1 := Mu1 + TmpDbl1;
467     Mu2 := Mu2 + TmpDbl2;
468     Mu3 := Mu3 + (0.165 * sin(G)) / radcor;
469     Mu4 := Mu4 + (0.841 * sin(H)) / radcor;
470    
471     Result.Io.X := R1 * sin(Mu1);
472     Result.Io.Y := -R1 * cos(Mu1) * sin(DE);
473    
474     Result.Europa.X := R2 * sin(Mu2);
475     Result.Europa.Y := -R2 * cos(Mu2) * sin(DE);
476    
477     Result.Ganymede.X := R3 * sin(Mu3);
478     Result.Ganymede.Y := -R3 * cos(Mu3) * sin(DE);
479    
480     Result.Callisto.X := R4 * sin(Mu4);
481     Result.Callisto.Y := -R4 * cos(Mu4) * sin(DE);
482     end;
483    
484     {-------------------------------------------------------------------------}
485    
486     function JupSatsHi(AJD : Double; Shadows : Boolean) : TStJupSats;
487     var
488     I : longint;
489     SunPos : SunCoordsRec;
490     STUT : TStDateTimeRec;
491     JupPos : TStEclipticalCord;
492    
493     SatX : array[1..5] of Double;
494     SatY : array[1..5] of Double;
495     SatZ : array[1..5] of Double;
496    
497     TD1,
498     TD2,
499     Angle, {Temporary Double values}
500     LTime, {Tau}
501     AJDT, {AJD adjusted for light time (Tau)}
502     JupX,
503     JupY,
504     JupZ, {Jupiter's geocentric rectangular coordinates}
505     EJDist, {Delta}
506     Jup1,
507     Jup2, {/\, Alpha}
508     DateDif, {t}
509     L1, L2,
510     L3, L4, {script L1-4}
511     Pi1, Pi2,
512     Pi3, Pi4, {Pi1-4}
513     W1, W2,
514     W3, W4, {Omega1-4}
515     Inequality, {upside down L}
516     PhiLambda,
517     NodeJup, {Psi}
518     AnomJup, {G}
519     AnomSat, {G'}
520     LongPerJ,
521     S1, S2,
522     S3, S4, {Sum1-4}
523     TL1, TL2,
524     TL3, TL4, {Capital L1-4}
525     B1, B2,
526     B3, B4, {tangent of latitude}
527     R1, R2,
528     R3, R4, {radius vector}
529     T0, {Julian Centuries}
530     Precession, {P}
531     Inclination {I}
532    
533     : Double;
534     Transforms : array[1..5] of TranformRec;
535    
536     begin
537     FillChar(Result, SizeOf(TStJupSats), #0);
538     AJD := DateTimeToAJD(AJD);
539     SunPos := SunCoords(AJD);
540    
541     if not Shadows then begin
542     TD1 := 5;
543     AJDT := AJD - 0.0057755183 * TD1; {first guess}
544     repeat
545     JupPos := ComputeJupiter(AJDT);
546    
547     JupX := JupPos.R0 * cos(JupPos.B0) * cos(JupPos.L0)
548     + SunPos.R * cos(SunPos.L);
549     JupY := JupPos.R0 * cos(JupPos.B0) * sin(JupPos.L0)
550     + SunPos.R * sin(SunPos.L);
551     JupZ := JupPos.R0 * sin(JupPos.B0);
552    
553     EJDist := sqrt(sqr(JupX) + sqr(JupY) + sqr(JupZ));
554     TD2 := abs(EJDist - TD1);
555     if abs(TD2) > 0.0005 then begin
556     AJDT := AJD - 0.0057755183 * ((EJDist + TD1) / 2);
557     TD1 := EJDist;
558     end;
559     until (TD2 <= 0.0005);
560     end else begin
561     JupPos := ComputeJupiter(AJD);
562    
563     JupX := JupPos.R0 * cos(JupPos.B0) * cos(JupPos.L0);
564     JupY := JupPos.R0 * cos(JupPos.B0) * sin(JupPos.L0);
565     JupZ := JupPos.R0 * sin(JupPos.B0);
566     EJDist := sqrt(sqr(JupX+SunPos.X) +
567     sqr(JupY+SunPos.Y) + sqr(JupZ+SunPos.Z));
568     end;
569    
570     Jup1 := StInvTan2(JupX, JupY);
571     Jup2 := ArcTan(JupZ / sqrt(sqr(JupX) + sqr(JupY)));
572    
573     DateDif := AJD - 2443000.5 - (0.0057755183 * EJDist);
574    
575     L1 := 106.07947 + 203.488955432 * DateDif;
576     L1 := Frac(L1/360.0) * 360.0;
577     if (L1 < 0) then
578     L1 := 360.0 + L1;
579     L1 := L1 / radcor;
580    
581     L2 := 175.72938 + 101.374724550 * DateDif;
582     L2 := Frac(L2/360.0) * 360.0;
583     if (L2 < 0) then
584     L2 := 360.0 + L2;
585     L2 := L2 / radcor;
586    
587     L3 := 120.55434 + 50.317609110 * DateDif;
588     L3 := Frac(L3/360.0) * 360.0;
589     if (L3 < 0) then
590     L3 := 360.0 + L3;
591     L3 := L3 / radcor;
592    
593     L4 := 84.44868 + 21.571071314 * DateDif;
594     L4 := Frac(L4/360.0) * 360.0;
595     if (L4 < 0) then
596     L4 := 360.0 + L4;
597     L4 := L4 / radcor;
598    
599     Pi1 := 58.3329 + 0.16103936 * DateDif;
600     Pi1 := Frac(Pi1/360.0) * 360.0;
601     if (Pi1 < 0) then
602     Pi1 := 360.0 + Pi1;
603     Pi1 := Pi1 / radcor;
604    
605     Pi2 := 132.8959 + 0.04647985 * DateDif;
606     Pi2 := Frac(Pi2/360.0) * 360.0;
607     if (Pi2 < 0) then
608     Pi2 := 360.0 + Pi2;
609     Pi2 := Pi2 / radcor;
610    
611     Pi3 := 187.2887 + 0.00712740 * DateDif;
612     Pi3 := Frac(Pi3/360.0) * 360.0;
613     if (Pi3 < 0) then
614     Pi3 := 360.0 + Pi3;
615     Pi3 := Pi3 / radcor;
616    
617     Pi4 := 335.3418 + 0.00183998 * DateDif;
618     Pi4 := Frac(Pi4/360.0) * 360.0;
619     if (Pi4 < 0) then
620     Pi4 := 360.0 + Pi4;
621     Pi4 := Pi4 / radcor;
622    
623     W1 := 311.0793 - 0.13279430 * DateDif;
624     W1 := Frac(W1/360.0) * 360.0;
625     if (W1 < 0) then
626     W1 := 360.0 + W1;
627     W1 := W1 / radcor;
628    
629     W2 := 100.5099 - 0.03263047 * DateDif;
630     W2 := Frac(W2/360.0) * 360.0;
631     if (W2 < 0) then
632     W2 := 360.0 + W2;
633     W2 := W2 / radcor;
634    
635     W3 := 119.1688 - 0.00717704 * DateDif;
636     W3 := Frac(W3/360.0) * 360.0;
637     if (W3 < 0) then
638     W3 := 360.0 + W3;
639     W3 := W3 / radcor;
640    
641     W4 := 322.5729 - 0.00175934 * DateDif;
642     W4 := Frac(W4/360.0) * 360.0;
643     if (W4 < 0) then
644     W4 := 360.0 + W4;
645     W4 := W4 / radcor;
646    
647     Inequality := 0.33033 * sin((163.679 + 0.0010512*DateDif) / radcor)
648     + 0.03439 * sin((34.486 - 0.0161731*DateDif) / radcor);
649     Inequality := Inequality / radcor;
650    
651     PhiLambda := 191.8132 + 0.17390023 * DateDif;
652     PhiLambda := Frac(PhiLambda / 360.0) * 360.0;
653     if (PhiLambda < 0) then
654     PhiLambda := 360.0 + PhiLambda;
655     PhiLambda := PhiLambda / radcor;
656    
657     NodeJup := 316.5182 - 0.00000208 * DateDif;
658     NodeJup := Frac(NodeJup / 360.0) * 360.0;
659     if (NodeJup < 0) then
660     NodeJup := 360.0 + NodeJup;
661     NodeJup := NodeJup / radcor;
662    
663     AnomJup := 30.23756 + 0.0830925701 * DateDif;
664     AnomJup := Frac(AnomJup / 360.0) * 360.0;
665     if (AnomJup < 0) then
666     AnomJup := 360.0 + AnomJup;
667     AnomJup := AnomJup / radcor + Inequality;
668    
669     AnomSat := 31.97853 + 0.0334597339 * DateDif;
670     AnomSat := Frac(AnomSat / 360.0) * 360.0;
671     if (AnomSat < 0) then
672     AnomSat := 360.0 + AnomSat;
673     AnomSat := AnomSat / radcor;
674    
675     LongPerJ := 13.469942 / radcor;
676    
677     S1 := 0.47259 * sin(2*(L1-L2))
678     - 0.03480 * sin(Pi3-Pi4)
679     - 0.01756 * sin(Pi1 + Pi3 - 2*LongPerJ - 2*AnomJup)
680     + 0.01080 * sin(L2 - 2*L3 + Pi3)
681     + 0.00757 * sin(PhiLambda)
682     + 0.00663 * sin(L2 - 2*L3 + Pi4)
683     + 0.00453 * sin(L1 - Pi3)
684     + 0.00453 * sin(L2 - 2*L3 + Pi2)
685     - 0.00354 * sin(L1-L2)
686     - 0.00317 * sin(2*NodeJup - 2*LongPerJ)
687     - 0.00269 * sin(L2 - 2*L3 + Pi1)
688     + 0.00263 * sin(L1 - Pi4)
689     + 0.00186 * sin(L1 - Pi1)
690     - 0.00186 * sin(AnomJup)
691     + 0.00167 * sin(Pi2 - Pi3)
692     + 0.00158 * sin(4*(L1-L2))
693     - 0.00155 * sin(L1 - L3)
694     - 0.00142 * sin(NodeJup + W3 - 2*LongPerJ - 2*AnomJup)
695     - 0.00115 * sin(2*(L1 - 2*L2 + W2))
696     + 0.00089 * sin(Pi2 - Pi4)
697     + 0.00084 * sin(W2 - W3)
698     + 0.00084 * sin(L1 + Pi3 - 2*LongPerJ - 2*AnomJup)
699     + 0.00053 * sin(NodeJup - W2);
700    
701     S2 := 1.06476 * sin(2*(L2-L3))
702     + 0.04253 * sin(L1 - 2*L2 + Pi3)
703     + 0.03579 * sin(L2 - Pi3)
704     + 0.02383 * sin(L1 - 2*L2 + Pi4)
705     + 0.01977 * sin(L2 - Pi4)
706     - 0.01843 * sin(PhiLambda)
707     + 0.01299 * sin(Pi3 - Pi4)
708     - 0.01142 * sin(L2 - L3)
709     + 0.01078 * sin(L2 - Pi2)
710     - 0.01058 * sin(AnomJup)
711     + 0.00870 * sin(L2 - 2*L3 + Pi2)
712     - 0.00775 * sin(2*(NodeJup - LongPerJ))
713     + 0.00524 * sin(2*(L1-L2))
714     - 0.00460 * sin(L1-L3)
715     + 0.00450 * sin(L2 - 2*L3 + Pi1)
716     + 0.00327 * sin(NodeJup - 2*AnomJup + W3 - 2*LongPerJ)
717     - 0.00296 * sin(Pi1 + Pi3 - 2*LongPerJ - 2*AnomJup)
718     - 0.00151 * sin(2*AnomJup)
719     + 0.00146 * sin(NodeJup - W3)
720     + 0.00125 * sin(NodeJup - W4)
721     - 0.00117 * sin(L1 -2*L3 + Pi3)
722     - 0.00095 * sin(2*(L2-W2))
723     + 0.00086 * sin(2*(L1-2*L2 +W2))
724     - 0.00086 * sin(5*AnomSat - 2*AnomJup + 52.225/radcor)
725     - 0.00078 * sin(L2-L4)
726     - 0.00064 * sin(L1 - 2*L3 + Pi4)
727     - 0.00063 * sin(3*L3 - 7*L4 + 4*Pi4)
728     + 0.00061 * sin(Pi1 - Pi4)
729     + 0.00058 * sin(2*(NodeJup - LongPerJ - AnomJup))
730     + 0.00058 * sin(W3 - W4)
731     + 0.00056 * sin(2*(L2-L4))
732     + 0.00055 * sin(2*(L1-L3))
733     + 0.00052 * sin(3*L3 - 7*L4 + Pi3 + 3*Pi4)
734     - 0.00043 * sin(L1 - Pi3)
735     + 0.00042 * sin(Pi3 - Pi2)
736     + 0.00041 * sin(5*(L2-L3))
737     + 0.00041 * sin(Pi4 - LongPerJ)
738     + 0.00038 * sin(L2 - Pi1)
739     + 0.00032 * sin(W2 - W3)
740     + 0.00032 * sin(2*(L3 - AnomJup - LongPerJ))
741     + 0.00029 * sin(Pi1 - Pi3);
742    
743     S3 := 0.16477 * sin(L3 - Pi3)
744     + 0.09062 * sin(L3 - Pi4)
745     - 0.06907 * sin(L2 - L3)
746     + 0.03786 * sin(Pi3 - Pi4)
747     + 0.01844 * sin(2*(L3-L4))
748     - 0.01340 * sin(AnomJup)
749     + 0.00703 * sin(L2 - 2*L3 + Pi3)
750     - 0.00670 * sin(2*(NodeJup - LongPerJ))
751     - 0.00540 * sin(L3-L4)
752     + 0.00481 * sin(Pi1 + Pi3 -2*LongPerJ - 2*AnomJup)
753     - 0.00409 * sin(L2 - 2*L3 + Pi2)
754     + 0.00379 * sin(L2 - 2*L3 + Pi4)
755     + 0.00235 * sin(NodeJup - W3)
756     + 0.00198 * sin(NodeJup - W4)
757     + 0.00180 * sin(PhiLambda)
758     + 0.00129 * sin(3*(L3-L4))
759     + 0.00124 * sin(L1-L3)
760     - 0.00119 * sin(5*AnomSat - 2*AnomJup + 52.225/radcor)
761     + 0.00109 * sin(L1-L2)
762     - 0.00099 * sin(3*L3 - 7*L4 + 4*Pi4)
763     + 0.00091 * sin(W3 - W4)
764     + 0.00081 * sin(3*L3 - 7*L4 + Pi3 + 3*Pi4)
765     - 0.00076 * sin(2*L2 - 3*L3 + Pi3)
766     + 0.00069 * sin(Pi4 - LongPerJ)
767     - 0.00058 * sin(2*L3 - 3*L4 + Pi4)
768     + 0.00057 * sin(L3 + Pi3 - 2*LongPerJ - 2*AnomJup)
769     - 0.00057 * sin(L3 - 2*L4 + Pi4)
770     - 0.00052 * sin(Pi2 - Pi3)
771     - 0.00052 * sin(L2 - 2*L3 + Pi1)
772     + 0.00048 * sin(L3 - 2*L4 + Pi3)
773     - 0.00045 * sin(2*L2 - 3*L3 + Pi4)
774     - 0.00041 * sin(Pi2 - Pi4)
775     - 0.00038 * sin(2*AnomJup)
776     - 0.00033 * sin(Pi3 - Pi4 + W3 - W4)
777     - 0.00032 * sin(3*L3 - 7*L4 + 2*Pi3 + 2*Pi4)
778     + 0.00030 * sin(4*(L3-L4))
779     - 0.00029 * sin(W3 + NodeJup - 2*LongPerJ - 2*AnomJup)
780     + 0.00029 * sin(L3 + Pi4 - 2*LongPerJ - 2*AnomJup)
781     + 0.00026 * sin(L3 - LongPerJ - AnomJup)
782     + 0.00024 * sin(L2 - 3*L3 + 2*L4)
783     + 0.00021 * sin(2*(L3 - LongPerJ - AnomJup))
784     - 0.00021 * sin(L3 - Pi2)
785     + 0.00017 * sin(2*(L3 - Pi3));
786    
787     S4 := 0.84109 * sin(L4 - Pi4)
788     + 0.03429 * sin(Pi4 - Pi3)
789     - 0.03305 * sin(2*(NodeJup - LongPerJ))
790     - 0.03211 * sin(AnomJup)
791     - 0.01860 * sin(L4 - Pi3)
792     + 0.01182 * sin(NodeJup - W4)
793     + 0.00622 * sin(L4 + Pi4 - 2*AnomJup - 2*LongPerJ)
794     + 0.00385 * sin(2*(L4 - Pi4))
795     - 0.00284 * sin(5*AnomSat - 2*AnomJup + 52.225/radcor)
796     - 0.00233 * sin(2*(NodeJup - Pi4))
797     - 0.00223 * sin(L3 - L4)
798     - 0.00208 * sin(L4 - LongPerJ)
799     + 0.00177 * sin(NodeJup + W4 - 2*Pi4)
800     + 0.00134 * sin(Pi4 - LongPerJ)
801     + 0.00125 * sin(2*(L4 - AnomJup - LongPerJ))
802     - 0.00117 * sin(2*AnomJup)
803     - 0.00112 * sin(2*(L3 - L4))
804     + 0.00106 * sin(3*L3 - 7*L4 + 4*Pi4)
805     + 0.00102 * sin(L4 - AnomJup - LongPerJ)
806     + 0.00096 * sin(2*L4 - NodeJup - W4)
807     + 0.00087 * sin(2*(NodeJup - W4))
808     - 0.00087 * sin(3*L3 - 7*L4 + Pi3 + 3*Pi4)
809     + 0.00085 * sin(L3 - 2*L4 + Pi4)
810     - 0.00081 * sin(2*(L4 - NodeJup))
811     + 0.00071 * sin(L4 + Pi4 -2*LongPerJ - 3*AnomJup)
812     + 0.00060 * sin(L1 - L4)
813     - 0.00056 * sin(NodeJup - W3)
814     - 0.00055 * sin(L3 - 2*L4 + Pi3)
815     + 0.00051 * sin(L2 - L4)
816     + 0.00042 * sin(2*(NodeJup - AnomJup - LongPerJ))
817     + 0.00039 * sin(2*(Pi4 - W4))
818     + 0.00036 * sin(NodeJup + LongPerJ - Pi4 - W4)
819     + 0.00035 * sin(2*AnomSat - AnomJup + 188.37/radcor)
820     - 0.00035 * sin(L4 - Pi4 + 2*LongPerJ - 2*NodeJup)
821     - 0.00032 * sin(L4 + Pi4 - 2*LongPerJ - AnomJup)
822     + 0.00030 * sin(3*L3 - 7*L4 + 2*Pi3 + 2*Pi4)
823     + 0.00030 * sin(2*AnomSat - 2*AnomJup + 149.15/radcor)
824     + 0.00028 * sin(L4 - Pi4 + 2*NodeJup - 2*LongPerJ)
825     - 0.00028 * sin(2*(L4 - W4))
826     - 0.00027 * sin(Pi3 - Pi4 + W3 - W4)
827     - 0.00026 * sin(5*AnomSat - 3*AnomJup + 188.37/radcor)
828     + 0.00025 * sin(W4 - W3)
829     - 0.00025 * sin(L2 - 3*L3 + 2*L4)
830     - 0.00023 * sin(3*(L3 - L4))
831     + 0.00021 * sin(2*L4 - 2*LongPerJ - 3*AnomJup)
832     - 0.00021 * sin(2*L3 - 3*L4 + Pi4)
833     + 0.00019 * sin(L4 - Pi4 - AnomJup)
834     - 0.00019 * sin(2*L4 - Pi3 - Pi4)
835     - 0.00018 * sin(L4 - Pi4 + AnomJup)
836     - 0.00016 * sin(L4 + Pi3 -2*LongPerJ - 2*AnomJup);
837    
838     S1 := S1/radcor;
839     S2 := S2/radcor;
840     S3 := S3/radcor;
841     S4 := S4/radcor;
842    
843     TL1 := L1 + S1;
844     TL2 := L2 + S2;
845     TL3 := L3 + S3;
846     TL4 := L4 + S4;
847    
848     B1 := 0.0006502 * sin(TL1 - W1)
849     + 0.0001835 * sin(TL1 - W2)
850     + 0.0000329 * sin(TL1 - W3)
851     - 0.0000311 * sin(TL1 - NodeJup)
852     + 0.0000093 * sin(TL1 - W4)
853     + 0.0000075 * sin(3*TL1 - 4*L2 - 1.9927/radcor * S1 + W2)
854     + 0.0000046 * sin(TL1 + NodeJup - 2*LongPerJ - 2*AnomJup);
855     B1 := ArcTan(B1);
856    
857     B2 := 0.0081275 * sin(TL2 - W2)
858     + 0.0004512 * sin(TL2 - W3)
859     - 0.0003286 * sin(TL2 - NodeJup)
860     + 0.0001164 * sin(TL2 - W4)
861     + 0.0000273 * sin(L1 - 2*L3 + 1.0146/radcor * S2 + W2)
862     + 0.0000143 * sin(TL2 + NodeJup - 2*LongPerJ - 2*AnomJup)
863     - 0.0000143 * sin(TL2 - W1)
864     + 0.0000035 * sin(TL2 - NodeJup + AnomJup)
865     - 0.0000028 * sin(L1 - 2*L3 + 1.0146/radcor * S2 + W3);
866     B2 := ArcTan(B2);
867    
868     B3 := 0.0032364 * sin(TL3 - W3)
869     - 0.0016911 * sin(TL3 - NodeJup)
870     + 0.0006849 * sin(TL3 - W4)
871     - 0.0002806 * sin(TL3 - W2)
872     + 0.0000321 * sin(TL3 + NodeJup - 2*LongPerJ - 2*AnomJup)
873     + 0.0000051 * sin(TL3 - NodeJup + AnomJup)
874     - 0.0000045 * sin(TL3 - NodeJup - AnomJup)
875     - 0.0000045 * sin(TL3 + NodeJup - 2*LongPerJ)
876     + 0.0000037 * sin(TL3 + NodeJup - 2*LongPerJ - 3*AnomJup)
877     + 0.0000030 * sin(2*L2 - 3*TL3 + 4.03/radcor * S3 + W2)
878     - 0.0000021 * sin(2*L2 - 3*TL3 + 4.03/radcor * S3 + W3);
879     B3 := ArcTan(B3);
880    
881     B4 := -0.0076579 * sin(TL4 - NodeJup)
882     + 0.0044148 * sin(TL4 - W4)
883     - 0.0005106 * sin(TL4 - W3)
884     + 0.0000773 * sin(TL4 + NodeJup - 2*LongPerJ - 2*AnomJup)
885     + 0.0000104 * sin(TL4 - NodeJup + AnomJup)
886     - 0.0000102 * sin(TL4 - NodeJup - AnomJup)
887     + 0.0000088 * sin(TL4 + NodeJup - 2*LongPerJ - 3*AnomJup)
888     - 0.0000038 * sin(TL4 + NodeJup - 2*LongPerJ - AnomJup);
889     B4 := ArcTan(B4);
890    
891     R1 := -0.0041339 * cos(2*(L1-L2))
892     - 0.0000395 * cos(L1 - Pi3)
893     - 0.0000214 * cos(L1 - Pi4)
894     + 0.0000170 * cos(L1 - L2)
895     - 0.0000162 * cos(L1 - Pi1)
896     - 0.0000130 * cos(4*(L1-L2))
897     + 0.0000106 * cos(L1 - L3)
898     - 0.0000063 * cos(L1 + Pi3 - 2*LongPerJ - 2*AnomJup);
899    
900     R2 := 0.0093847 * cos(L1-L2)
901     - 0.0003114 * cos(L2 - Pi3)
902     - 0.0001738 * cos(L2 - Pi4)
903     - 0.0000941 * cos(L2 - Pi2)
904     + 0.0000553 * cos(L2 - L3)
905     + 0.0000523 * cos(L1 - L3)
906     - 0.0000290 * cos(2*(L1-L2))
907     + 0.0000166 * cos(2*(L2-W2))
908     + 0.0000107 * cos(L1 - 2*L3 + Pi3)
909     - 0.0000102 * cos(L2 - Pi1)
910     - 0.0000091 * cos(2*(L1-L3));
911    
912     R3 := -0.0014377 * cos(L3 - Pi3)
913     - 0.0007904 * cos(L3 - Pi4)
914     + 0.0006342 * cos(L2 - L3)
915     - 0.0001758 * cos(2*(L3-L4))
916     + 0.0000294 * cos(L3 - L4)
917     - 0.0000156 * cos(3*(L3-L4))
918     + 0.0000155 * cos(L1 - L3)
919     - 0.0000153 * cos(L1 - L2)
920     + 0.0000070 * cos(2*L2 - 3*L3 + Pi3)
921     - 0.0000051 * cos(L3 + Pi3 - 2*LongPerJ - 2*AnomJup);
922    
923     R4 := -0.0073391 * cos(L4 - Pi4)
924     + 0.0001620 * cos(L4 - Pi3)
925     + 0.0000974 * cos(L3 - L4)
926     - 0.0000541 * cos(L4 + Pi4 - 2*LongPerJ - 2*AnomJup)
927     - 0.0000269 * cos(2*(L4-Pi4))
928     + 0.0000182 * cos(L4- LongPerJ)
929     + 0.0000177 * cos(2*(L3-L4))
930     - 0.0000167 * cos(2*L4 - NodeJup - W4)
931     + 0.0000167 * cos(NodeJup - W4)
932     - 0.0000155 * cos(2*(L4-LongPerj-AnomJup))
933     + 0.0000142 * cos(2*(L4-NodeJup))
934     + 0.0000104 * cos(L1 - L4)
935     + 0.0000092 * cos(L2 - L4)
936     - 0.0000089 * cos(L4 - LongPerJ - AnomJup)
937     - 0.0000062 * cos(L4 + Pi4 - 2*LongPerJ - 3*AnomJup)
938     + 0.0000048 * cos(2*(L4-W4));
939    
940     R1 := 5.90730 * (1 + R1);
941     R2 := 9.39912 * (1 + R2);
942     R3 := 14.99240 * (1 + R3);
943     R4 := 26.36990 * (1 + R4);
944    
945     T0 := (AJD - 2433282.423) / 36525;
946     Precession := (1.3966626*T0 + 0.0003088*sqr(T0)) / radcor;
947    
948     TL1 := TL1 + Precession;
949     TL2 := TL2 + Precession;
950     TL3 := TL3 + Precession;
951     TL4 := TL4 + Precession;
952     NodeJup := NodeJup + Precession;
953    
954     T0 := (AJD - AstJulianDatePrim(1900, 1, 1, 0)) / 36525;
955     Inclination := (3.120262 + 0.0006*T0) / radcor;
956    
957     SatX[1] := R1 * cos(TL1 - NodeJup) * cos(B1);
958     SatY[1] := R1 * sin(TL1 - NodeJup) * cos(B1);
959     SatZ[1] := R1 * sin(B1);
960    
961     SatX[2] := R2 * cos(TL2 - NodeJup) * cos(B2);
962     SatY[2] := R2 * sin(TL2 - NodeJup) * cos(B2);
963     SatZ[2] := R2 * sin(B2);
964    
965     SatX[3] := R3 * cos(TL3 - NodeJup) * cos(B3);
966     SatY[3] := R3 * sin(TL3 - NodeJup) * cos(B3);
967     SatZ[3] := R3 * sin(B3);
968    
969     SatX[4] := R4 * cos(TL4 - NodeJup) * cos(B4);
970     SatY[4] := R4 * sin(TL4 - NodeJup) * cos(B4);
971     SatZ[4] := R4 * sin(B4);
972    
973     SatX[5] := 0;
974     SatY[5] := 0;
975     SatZ[5] := 1;
976    
977     T0 := (AJD - 2451545.0) / 36525.0;
978     TD1 := 100.464441
979     + 1.0209550 * T0
980     + 0.00040117 * sqr(T0)
981     + 0.000000569 * sqr(T0) * T0;
982     TD1 := TD1 / radcor;
983    
984     TD2 := 1.303270
985     - 0.0054966 * T0
986     + 0.00000465 * sqr(T0)
987     - 0.000000004 * sqr(T0) * T0;
988     TD2 := TD2 / radcor;
989    
990     for I := 1 to 5 do begin
991     Transforms[I].A[1] := SatX[I];
992     Transforms[I].B[1] := SatY[I] * cos(Inclination)
993     - SatZ[I] * sin(Inclination);
994     Transforms[I].C[1] := SatY[I] * sin(Inclination)
995     + SatZ[I] * cos(Inclination);
996    
997     Transforms[I].A[2] := Transforms[I].A[1] * cos(NodeJup - TD1)
998     - Transforms[I].B[1] * sin(NodeJup - TD1);
999     Transforms[I].B[2] := Transforms[I].A[1] * sin(NodeJup - TD1)
1000     + Transforms[I].B[1] * cos(NodeJup - TD1);
1001     Transforms[I].C[2] := Transforms[I].C[1];
1002    
1003     Transforms[I].A[3] := Transforms[I].A[2];
1004     Transforms[I].B[3] := Transforms[I].B[2] * cos(TD2)
1005     - Transforms[I].C[2] * sin(TD2);
1006     Transforms[I].C[3] := Transforms[I].B[2] * sin(TD2)
1007     + Transforms[I].C[2] * cos(TD2);
1008    
1009     Transforms[I].A[4] := Transforms[I].A[3] * cos(TD1)
1010     - Transforms[I].B[3] * sin(TD1);
1011     Transforms[I].B[4] := Transforms[I].A[3] * sin(TD1)
1012     + Transforms[I].B[3] * cos(TD1);
1013     Transforms[I].C[4] := Transforms[I].C[3];
1014    
1015     Transforms[I].A[5] := Transforms[I].A[4] * sin(Jup1)
1016     - Transforms[I].B[4] * cos(Jup1);
1017     Transforms[I].B[5] := Transforms[I].A[4] * cos(Jup1)
1018     + Transforms[I].B[4] * sin(Jup1);
1019     Transforms[I].C[5] := Transforms[I].C[4];
1020    
1021     Transforms[I].A[6] := Transforms[I].A[5];
1022     Transforms[I].B[6] := Transforms[I].C[5] * sin(Jup2)
1023     + Transforms[I].B[5] * cos(Jup2);
1024     Transforms[I].C[6] := Transforms[I].C[5] * cos(Jup2)
1025     - Transforms[I].B[5] * sin(Jup2);
1026     end;
1027    
1028     Angle := StInvTan2(Transforms[5].C[6], Transforms[5].A[6]);
1029    
1030     {Io calculations}
1031     Result.Io.X := Transforms[1].A[6] * cos(Angle)
1032     - Transforms[1].C[6] * sin(Angle);
1033     Result.Io.Y := Transforms[1].A[6] * sin(Angle)
1034     + Transforms[1].C[6] * cos(Angle);
1035     TD1 := Transforms[1].B[6];
1036    
1037     {correct for light time}
1038     TD2 := abs(TD1) / 17295 * sqrt(1 - sqr(Result.Io.X/R1));
1039     Result.Io.X := Result.Io.X + TD2;
1040    
1041     {correct for perspective}
1042     TD2 := EJDist / (EJDist + TD1/2095);
1043     Result.Io.X := Result.Io.X * TD2;
1044     Result.Io.Y := Result.Io.Y * TD2;
1045    
1046     {Europa calculations}
1047     Result.Europa.X := Transforms[2].A[6] * cos(Angle)
1048     - Transforms[2].C[6] * sin(Angle);
1049     Result.Europa.Y := Transforms[2].A[6] * sin(Angle)
1050     + Transforms[2].C[6] * cos(Angle);
1051     TD1 := Transforms[2].B[6];
1052    
1053     {correct for light time}
1054     TD2 := abs(TD1) / 21819 * sqrt(1 - sqr(Result.Europa.X/R2));
1055     Result.Europa.X := Result.Europa.X + TD2;
1056    
1057     {correct for perspective}
1058     TD2 := EJDist / (EJDist + TD1/2095);
1059     Result.Europa.X := Result.Europa.X * TD2;
1060     Result.Europa.Y := Result.Europa.Y * TD2;
1061    
1062     {Ganymede calculations}
1063     Result.Ganymede.X := Transforms[3].A[6] * cos(Angle)
1064     - Transforms[3].C[6] * sin(Angle);
1065     Result.Ganymede.Y := Transforms[3].A[6] * sin(Angle)
1066     + Transforms[3].C[6] * cos(Angle);
1067     TD1 := Transforms[3].B[6];
1068    
1069     {correct for light time}
1070     TD2 := abs(TD1) / 27558 * sqrt(1 - sqr(Result.Ganymede.X/R3));
1071     Result.Ganymede.X := Result.Ganymede.X + TD2;
1072    
1073     {correct for perspective}
1074     TD2 := EJDist / (EJDist + TD1/2095);
1075     Result.Ganymede.X := Result.Ganymede.X * TD2;
1076     Result.Ganymede.Y := Result.Ganymede.Y * TD2;
1077    
1078     {Callisto calculations}
1079     Result.Callisto.X := Transforms[4].A[6] * cos(Angle)
1080     - Transforms[4].C[6] * sin(Angle);
1081     Result.Callisto.Y := Transforms[4].A[6] * sin(Angle)
1082     + Transforms[4].C[6] * cos(Angle);
1083     TD1 := Transforms[4].B[6];
1084    
1085     {correct for light time}
1086     TD2 := abs(TD1) / 36548 * sqrt(1 - sqr(Result.Callisto.X/R4));
1087     Result.Callisto.X := Result.Callisto.X + TD2;
1088    
1089     {correct for perspective}
1090     TD2 := EJDist / (EJDist + TD1/2095);
1091     Result.Callisto.X := Result.Callisto.X * TD2;
1092     Result.Callisto.Y := Result.Callisto.Y * TD2;
1093     end;
1094    
1095     {-------------------------------------------------------------------------}
1096    
1097     function GetJupSats(JD : TDateTime; HighPrecision, Shadows : Boolean) : TStJupSats;
1098     begin
1099     if not HighPrecision then
1100     Result := JupSatsLo(JD)
1101     else
1102     Result := JupSatsHi(JD, Shadows);
1103     end;
1104    
1105     end.

  ViewVC Help
Powered by ViewVC 1.1.20