paweld/lazarus-ccr
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
b10c38f5d7
lazarus-ccr/components/systools/source/general/run/steclpse.pas
wp_xxyyzz 543cdf06d9 systools: Rearrange units and packages 
git-svn-id: https://svn.code.sf.net/p/lazarus-ccr/svn@6159 8e941d3f-bd1b-0410-a28a-d453659cc2b4
2018-01-30 16:17:37 +00:00

738 lines
22 KiB
ObjectPascal
Raw Blame 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: StEclpse.pas 4.04 *}
{*********************************************************}
{* SysTools: Lunar/Solar Eclipses *}
{*********************************************************}
{$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-98 source files, Brian D. Warner, 1995-98. }
{ }
{ ************************************************************** }
unit StEclpse;
interface
uses
{$IFDEF UseMathUnit} Math, {$ENDIF}
StBase, StList, StDate, StAstro, StMath;
type
TStEclipseType = (etLunarPenumbral, etLunarPartial, etLunarTotal,
etSolarPartial, etSolarAnnular, etSolarTotal,
etSolarAnnularTotal);
TStHemisphereType = (htNone, htNorthern, htSouthern);
TStContactTimes = packed record
UT1, {start of lunar penumbral phase}
UT2, {end of lunar penumbral phase}
FirstContact, {start of partial eclipse}
SecondContact, {start of totality}
MidEclipse, {mid-eclipse}
ThirdContact, {end of totality}
FourthContact : TDateTime; {end of partial phase}
end;
TStLongLat = packed record
JD : TDateTime;
Longitude,
Latitude,
Duration : Double;
end;
PStLongLat = ^TStLongLat;
TStEclipseRecord = packed record
EType : TStEclipseType; {type of Eclipse}
Magnitude : Double; {magnitude of eclipse}
Hemisphere : TStHemisphereType; {favored hemisphere - solar}
LContacts : TStContactTimes; {Universal Times of critical points}
{ in lunar eclipses}
Path : TStList; {longitude/latitude of moon's shadow}
end; { during solar eclipse}
PStEclipseRecord = ^TStEclipseRecord;
TStBesselianRecord = packed record
JD : TDateTime;
Delta,
Angle,
XAxis,
YAxis,
L1,
L2 : Double;
end;
TStEclipses = class(TStList)
{.Z+}
protected {private}
FBesselianElements : array[1..25] of TStBesselianRecord;
F0,
FUPrime,
FDPrime : Double;
function GetEclipse(Idx : longint) : PStEclipseRecord;
procedure CentralEclipseTime(JD, K, J2,
SunAnom, MoonAnom,
ArgLat, AscNode, EFac : Double;
var F1, A1, CentralTime : Double);
procedure CheckForEclipse(K : Double);
procedure TotalLunarEclipse(CentralJD, MoonAnom, Mu,
PMag, UMag, Gamma : Double);
procedure PartialLunarEclipse(CentralJD, MoonAnom, Mu,
PMag, UMag, Gamma : Double);
procedure PenumbralLunarEclipse(CentralJD, MoonAnom, Mu,
PMag, UMag, Gamma : Double);
procedure GetBesselianElements(CentralJD : Double);
procedure GetShadowPath(I1, I2 : Integer; Path : TStList);
procedure NonPartialSolarEclipse(CentralJD, Mu, Gamma : Double);
procedure PartialSolarEclipse(CentralJD, Mu, Gamma : Double);
{.Z-}
public
constructor Create(NodeClass : TStNodeClass);
override;
procedure FindEclipses(Year : integer);
property Eclipses[Idx : longint] : PStEclipseRecord
read GetEclipse;
end;
implementation
procedure DisposePathData(Data : Pointer); far;
begin
Dispose(PStLongLat(Data));
end;
procedure DisposeEclipseRecord(Data : Pointer); far;
var
EcData : TStEclipseRecord;
begin
EcData := TStEclipseRecord(Data^);
if (Assigned(EcData.Path)) then
EcData.Path.Free;
Dispose(PStEclipseRecord(Data));
end;
constructor TStEclipses.Create(NodeClass : TStNodeClass);
begin
inherited Create(NodeClass);
DisposeData := DisposeEclipseRecord;
end;
function TStEclipses.GetEclipse(Idx : longint) : PStEclipseRecord;
begin
if (Idx < 0) or (Idx > pred(Count)) then
Result := nil
else
Result := PStEclipseRecord(Items[Idx].Data);
end;
procedure TStEclipses.FindEclipses(Year : integer);
var
K,
MeanJD,
JDofFirst,
JDofLast : Double;
begin
JDofFirst := AstJulianDatePrim(Year, 1, 1, 0);
JDofLast := AstJulianDatePrim(Year, 12, 31, pred(SecondsInDay));
K := Int((Year - 2000) * 12.3685 - 1);
repeat
MeanJD := 2451550.09765 + 29.530588853 * K;
if (MeanJD < JDofFirst) then
K := K + 0.5;
until (MeanJD >= JDofFirst);
while (MeanJD < JDofLast) do begin
CheckForEclipse(K);
K := K + 0.5;
MeanJD := 2451550.09765 + 29.530588853 * K;
end;
end;
procedure TStEclipses.CentralEclipseTime(JD, K, J2,
SunAnom, MoonAnom,
ArgLat, AscNode, EFac : Double;
var F1, A1, CentralTime : Double);
{the mean error of this routine is 0.36 minutes in a test between}
{1951 through 2050 with a maximum of 1.1 - Meeus}
begin
F1 := ArgLat - (0.02665/radcor) * sin(AscNode);
A1 := (299.77/radcor)
+ (0.107408/radcor) * K
- (0.009173/radcor) * J2;
if (Frac(K) > 0.1) then
{correction at Full Moon - Lunar eclipse}
CentralTime := JD
- 0.4065 * sin(MoonAnom)
+ 0.1727 * EFac * sin(SunAnom)
else
{correction at New Moon - solar eclipse}
CentralTime := JD
- 0.4075 * sin(MoonAnom)
+ 0.1721 * EFac * sin(SunAnom);
CentralTime := CentralTime
+ 0.0161 * sin(2 * MoonAnom)
- 0.0097 * sin(2 * F1)
+ 0.0073 * sin(MoonAnom - SunAnom) * EFac
- 0.0050 * sin(MoonAnom + SunAnom) * EFac
- 0.0023 * sin(MoonAnom - 2*F1)
+ 0.0021 * sin(2*SunAnom) * EFac
+ 0.0012 * sin(MoonAnom + 2*F1)
+ 0.0006 * sin(2*MoonAnom + SunAnom) * EFac
- 0.0004 * sin(3*MoonAnom)
- 0.0003 * sin(SunAnom + 2*F1) * EFac
+ 0.0003 * sin(A1)
- 0.0002 * sin(SunAnom - 2*F1) * EFac
- 0.0002 * sin(2*MoonAnom - SunAnom) * EFac
- 0.0002 * sin(AscNode);
end;
procedure TStEclipses.CheckForEclipse(K : Double);
var
MeanJD,
J1, J2, J3,
PMag, UMag,
CentralJD,
SunAnom,
MoonAnom,
ArgLat,
AscNode,
EFac,
DeltaT,
F1, A1,
P, Q, W,
Gamma, Mu : Double;
begin
{compute Julian Centuries}
J1 := K / 1236.85;
J2 := Sqr(J1);
J3 := J2 * J1;
{mean Julian Date for phase}
MeanJD := 2451550.09765
+ 29.530588853 * K
+ 0.0001337 * J2
- 0.000000150 * J3
+ 0.00000000073 * J2 * J2;
{solar mean anomaly}
SunAnom := 2.5534
+ 29.1053569 * K
- 0.0000218 * J2
- 0.00000011 * J3;
SunAnom := Frac(SunAnom / 360.0) * 360;
if (SunAnom < 0) then
SunAnom := SunAnom + 360.0;
{lunar mean anomaly}
MoonAnom := 201.5643
+ 385.81693528 * K
+ 0.0107438 * J2
+ 0.00001239 * J3
- 0.000000058 * J2 * J2;
MoonAnom := Frac(MoonAnom / 360.0) * 360;
if (MoonAnom < 0) then
MoonAnom := MoonAnom + 360.0;
{lunar argument of latitude}
ArgLat := 160.7108
+ 390.67050274 * K
- 0.0016341 * J2
- 0.00000227 * J3
+ 0.000000011 * J2 * J2;
ArgLat := Frac(ArgLat / 360.0) * 360;
if (ArgLat < 0) then
ArgLat := ArgLat + 360.0;
{lunar ascending node}
AscNode := 124.7746
- 1.56375580 * K
+ 0.0020691 * J2
+ 0.00000215 * J3;
AscNode := Frac(AscNode / 360.0) * 360;
if (AscNode < 0) then
AscNode := AscNode + 360.0;
{convert to radians}
SunAnom := SunAnom/radcor;
MoonAnom := MoonAnom/radcor;
ArgLat := ArgLat/radcor;
AscNode := AscNode/radcor;
{correction factor}
EFac := 1.0 - 0.002516 * J1 - 0.0000074 * J2;
{if AscNode > 21 deg. from 0 or 180 then no eclipse}
if (abs(sin(ArgLat)) > (sin(21.0/radcor))) then Exit;
{there is probably an eclipse - what kind? when?}
CentralEclipseTime(MeanJD, K, J2, SunAnom, MoonAnom,
ArgLat, AscNode, EFac,
F1, A1, CentralJD);
{Central JD is in Dynamical Time. Sun/Moon Positions are based on UT}
{An APPROXIMATE conversion is made to UT. This has limited accuracy}
DeltaT := (-15 + (sqr(CentralJD - 2382148) / 41048480)) / 86400;
CentralJD := CentralJD - DeltaT;
P := 0.2070 * sin(SunAnom) * EFac
+ 0.0024 * sin(2*SunAnom) * EFac
- 0.0392 * sin(MoonAnom)
+ 0.0116 * sin(2*MoonAnom)
- 0.0073 * sin(SunAnom + MoonAnom) * EFac
+ 0.0067 * sin(MoonAnom - SunAnom) * EFac
+ 0.0118 * sin(2*F1);
Q := 5.2207
- 0.0048 * cos(SunAnom) * EFac
+ 0.0020 * cos(2*SunAnom) * EFac
- 0.3299 * cos(MoonAnom)
- 0.0060 * cos(SunAnom + MoonAnom) * EFac
+ 0.0041 * cos(MoonAnom - SunAnom) * EFac;
W := abs(cos(F1));
Gamma := (P * cos(F1) + Q * sin(F1)) * (1 - 0.0048 * W);
Mu := 0.0059
+ 0.0046 * cos(SunAnom) * EFac
- 0.0182 * cos(MoonAnom)
+ 0.0004 * cos(2*MoonAnom)
- 0.0005 * cos(SunAnom + MoonAnom);
if (Frac(abs(K)) > 0.1) then begin
{Check for Lunar Eclipse possibilities}
PMag := (1.5573 + Mu - abs(Gamma)) / 0.5450;
UMag := (1.0128 - Mu - abs(Gamma)) / 0.5450;
if (UMag >= 1.0) then
TotalLunarEclipse(CentralJD, MoonAnom, Mu, PMag, UMag, Gamma)
else if (UMag > 0) then
PartialLunarEclipse(CentralJD, MoonAnom, Mu, PMag, UMag, Gamma)
else if ((UMag < 0) and (PMag > 0)) then
PenumbralLunarEclipse(CentralJD, MoonAnom, Mu, PMag, UMag, Gamma);
end else begin
{Check for Solar Eclipse possibilities}
{
Non-partial eclipses
--------------------
central Axis of moon's umbral shadow touches earth's surface
Can be total, annular, or both
non-central Axis of moon's umbral shadow does not touch earth's surface
Eclipse is usually partial but can be one of possibilities
for central eclipse if very near one of the earth's poles
Partial eclipses
----------------
partial None of the moon's umbral shadow touches the earth's surface
}
if (abs(Gamma) <= (0.9972 + abs(Mu))) then
NonPartialSolarEclipse(CentralJD, Mu, Gamma)
else begin
if (abs(Gamma) < 1.5433 + Mu) then
PartialSolarEclipse(CentralJD, Mu, Gamma);
end;
end;
end;
procedure TStEclipses.TotalLunarEclipse(CentralJD, MoonAnom, Mu,
PMag, UMag, Gamma : Double);
var
TLE : PStEclipseRecord;
PartialSemiDur,
TotalSemiDur,
Dbl1 : Double;
begin
New(TLE);
TLE^.Magnitude := UMag;
TLE^.Hemisphere := htNone;
TLE^.EType := etLunarTotal;
TLE^.Path := nil;
CentralJD := AJDToDateTime(CentralJD);
PartialSemiDur := 1.0128 - Mu;
TotalSemiDur := 0.4678 - Mu;
Dbl1 := 0.5458 + 0.0400 * cos(MoonAnom);
PartialSemiDur := 60/Dbl1 * sqrt(sqr(PartialSemiDur) - sqr(Gamma)) / 1440;
TotalSemiDur := 60/Dbl1 * sqrt(sqr(TotalSemiDur) - sqr(Gamma)) / 1440;
TLE^.LContacts.FirstContact := CentralJD - PartialSemiDur;
TLE^.LContacts.SecondContact := CentralJD - TotalSemiDur;
TLE^.LContacts.MidEclipse := CentralJD;
TLE^.LContacts.ThirdContact := CentralJD + TotalSemiDur;
TLE^.LContacts.FourthContact := CentralJD + PartialSemiDur;
PartialSemiDur := 60/Dbl1 * sqrt(sqr(1.5573 + Mu) - sqr(Gamma)) / 1440;
TLE^.LContacts.UT1 := CentralJD - PartialSemiDur;
TLE^.LContacts.UT2 := CentralJD + PartialSemiDur;
Self.Append(TLE);
end;
procedure TStEclipses.PartialLunarEclipse(CentralJD, MoonAnom, Mu,
PMag, UMag, Gamma : Double);
var
TLE : PStEclipseRecord;
PartialSemiDur,
Dbl1 : Double;
begin
New(TLE);
TLE^.Magnitude := UMag;
TLE^.Hemisphere := htNone;
TLE^.EType := etLunarPartial;
TLE^.Path := nil;
CentralJD := AJDToDateTime(CentralJD);
PartialSemiDur := 1.0128 - Mu;
Dbl1 := 0.5458 + 0.0400 * cos(MoonAnom);
PartialSemiDur := 60/Dbl1 * sqrt(sqr(PartialSemiDur) - sqr(Gamma)) / 1440;
TLE^.LContacts.FirstContact := CentralJD - PartialSemiDur;
TLE^.LContacts.SecondContact := 0;
TLE^.LContacts.MidEclipse := CentralJD;
TLE^.LContacts.ThirdContact := 0;
TLE^.LContacts.FourthContact := CentralJD + PartialSemiDur;
PartialSemiDur := 60/Dbl1 * sqrt(sqr(1.5573 + Mu) - sqr(Gamma)) / 1440;
TLE^.LContacts.UT1 := CentralJD - PartialSemiDur;
TLE^.LContacts.UT2 := CentralJD + PartialSemiDur;
Self.Append(TLE);
end;
procedure TStEclipses.PenumbralLunarEclipse(CentralJD, MoonAnom, Mu,
PMag, UMag, Gamma : Double);
var
TLE : PStEclipseRecord;
PartialSemiDur,
Dbl1 : Double;
begin
New(TLE);
TLE^.Magnitude := PMag;
TLE^.Hemisphere := htNone;
TLE^.EType := etLunarPenumbral;
TLE^.Path := nil;
CentralJD := AJDToDateTime(CentralJD);
TLE^.LContacts.FirstContact := 0;
TLE^.LContacts.SecondContact := 0;
TLE^.LContacts.MidEclipse := CentralJD;
TLE^.LContacts.ThirdContact := 0;
TLE^.LContacts.FourthContact := 0;
Dbl1 := 0.5458 + 0.0400 * cos(MoonAnom);
PartialSemiDur := 60/Dbl1 * sqrt(sqr(1.5573 + Mu) - sqr(Gamma)) / 1440;
TLE^.LContacts.UT1 := CentralJD - PartialSemiDur;
TLE^.LContacts.UT2 := CentralJD + PartialSemiDur;
Self.Append(TLE);
end;
procedure TStEclipses.GetBesselianElements(CentralJD : Double);
var
I,
Mins : LongInt;
CurJD,
SidTime,
SunDist,
MoonDist,
DistRatio,
Alpha,
Theta,
Sun1,
Sun2,
Moon1,
Moon2,
Dbl3,
F1, F2 : Double;
DTRec : TStDateTimeRec;
SunXYZ : TStSunXYZRec;
Sun : TStPosRec;
Moon : TStMoonPosRec;
begin
{compute BesselianElements every 10 minutes starting 2 hours prior to CentralJD}
{but forcing positions to be at multiple of ten minutes}
for I := 1 to 25 do begin
CurJD := CentralJD + ((I-13) * (10/1440));
DTRec.D := AstJulianDateToStDate(CurJD, True);
if ((Frac(CurJD) + 0.5) >= 1) then
Mins := Trunc(((Frac(CurJD) + 0.5)-1) * 1440)
else
Mins := Trunc((Frac(CurJD) + 0.5) * 1440);
{changed because, for example, both 25 and 35 rounded to 30}
Mins := ((Mins + 5) div 10) * 10;
if (Mins = 1440) then
Mins := 0;
DTRec.T := Mins * 60;
SidTime := SiderealTime(DTRec) / radcor;
SunXYZ := SunPosPrim(DTRec);
Sun := SunPos(DTRec);
Moon := MoonPos(DTRec);
Sun1 := Sun.RA/radcor;
Sun2 := Sun.DC/radcor;
Moon1 := Moon.RA/radcor;
Moon2 := Moon.DC/radcor;
FBesselianElements[I].JD := StDateToDateTime(DTRec.D)
+ StTimeToDateTime(DTRec.T);
SunDist := 1.0 / sin(8.794/SunXYZ.RV/3600.0/radcor);
MoonDist := 1.0 / sin(Moon.Plx/radcor);
DistRatio := MoonDist / SunDist;
Theta := DistRatio/cos(Sun2) * cos(Moon2) * (Moon1 - Sun1);
Theta := Theta/(1.0-DistRatio);
Alpha := Sun1 - Theta;
Theta := DistRatio/(1.0 - DistRatio) * (Moon2 - Sun2);
FBesselianElements[I].Delta := Sun2 - Theta;
FBesselianElements[I].Angle := SidTime - Alpha;
FBesselianElements[I].XAxis := MoonDist * cos(Moon2) * (sin(Moon1 - Alpha));
Dbl3 := FBesselianElements[I].Delta;
FBesselianElements[I].YAxis := MoonDist * (sin(Moon2) * cos(Dbl3)
- cos(Moon2) * sin(Dbl3) * cos(Moon1 - Alpha));
Theta := DistRatio * SunXYZ.RV;
Theta := SunXYZ.RV - Theta;
F1 := StInvSin(0.004664012/Theta);
F2 := StInvSin(0.004640787/Theta);
Theta := MoonDist * (sin(Moon2) * sin(Dbl3) + cos(Moon2)
* cos(Dbl3) * cos(Moon1 - Alpha));
FBesselianElements[I].L1 := (Theta + 0.272453/sin(F1)) * StTan(F1);
FBesselianElements[I].L2 := (Theta - 0.272453/sin(F2)) * StTan(F2);
if (I = 13) then
F0 := StTan(F2);
if (I = 16) then begin
FUPrime := FBesselianElements[16].Angle - FBesselianElements[10].Angle;
FDPrime := FBesselianElements[16].Delta - FBesselianElements[10].Delta;
end;
end;
end;
procedure TStEclipses.GetShadowPath(I1, I2 : Integer; Path : TStList);
var
J,
I3, I4,
I5, I6 : integer;
Delta,
Dbl1,
Dbl2,
P1,
XAxis,
YAxis,
Eta,
R1, R2,
D1, D2,
D3, D4,
V3, V4,
V5, V6, V7,
XPrime,
YPrime : Double;
Position : PStLongLat;
begin
for J := I1 to I2 do begin
Eta := 0.00669454;
Delta := FBesselianElements[J].Delta;
XAxis := FBesselianElements[J].XAxis;
YAxis := FBesselianElements[J].YAxis;
R1 := sqrt(1.0 - Eta * sqr(cos(Delta)));
R2 := sqrt(1.0 - Eta * sqr(sin(Delta)));
D1 := sin(Delta) / R1;
D2 := sqrt(1.0 - Eta) * cos(Delta) / R1;
D3 := Eta * sin(Delta) * cos(Delta) / R1 / R2;
D4 := sqrt(1.0 - Eta) / R1 / R2;
V3 := YAxis / R1;
V4 := sqrt(1.0 - sqr(XAxis) - sqr(V3));
V5 := R2 * (V4 * D4 - V3 * D3);
V6 := FUPrime * (-YAxis * sin(Delta) + V5 * cos(Delta));
V7 := FUPrime * XAxis * sin(Delta) - FDPrime * V5;
if ((I2-I1) div 2) >= 4 then begin
I3 := (I2-I1) div 2;
I4 := I1 + I3;
I5 := I4 - 3;
I6 := I4 + 3;
XPrime := FBesselianElements[I6].XAxis
- FBesselianElements[I5].XAxis;
YPrime := FBesselianElements[I6].YAxis
- FBesselianElements[I5].YAxis;
end else begin
XPrime := (FBesselianElements[J+1].XAxis -
FBesselianElements[J-1].XAxis) * 3;
YPrime := (FBesselianElements[J+1].YAxis -
FBesselianElements[J-1].YAxis) * 3;
end;
New(Position);
Dbl1 := sqrt(sqr(XPrime - V6) + sqr(YPrime - V7));
Position^.JD := FBesselianElements[J].JD;
Dbl2 := (FBesselianElements[J].L2 - V5 * F0) / Dbl1;
Dbl2 := abs(Dbl2 * 120);
Position^.Duration := int(Dbl2) + frac(Dbl2) * 0.6;
Dbl1 := -V3 * D1 + V4 * D2;
P1 := StInvTan2(Dbl1, XAxis);
Dbl2 := (FBesselianElements[J].Angle - P1) * radcor;
Dbl2 := frac(Dbl2 / 360.0) * 360;
if (Dbl2 < 0) then
Dbl2 := Dbl2 + 360.0;
if (Dbl2 > 180.0) and (Dbl2 < 360.0) then
Dbl2 := Dbl2 - 360.0;
Position^.Longitude := Dbl2;
Dbl1 := StInvSin(V3 * D2 + V4 * D1);
Dbl2 := ArcTan(1.003364 * StTan(Dbl1)) * radcor;
Position^.Latitude := Dbl2;
Path.Append(Position);
end;
end;
procedure TStEclipses.NonPartialSolarEclipse(CentralJD, Mu, Gamma : Double);
var
SolEc : PStEclipseRecord;
I1, I2 : Integer;
begin
New(SolEc);
if (Mu < 0) then
SolEc^.EType := etSolarTotal
else if (Mu > 0.0047) then
SolEc^.EType := etSolarAnnular
else begin
if (Mu < (0.00464 * sqrt(1 - sqr(Gamma)))) then
SolEc^.EType := etSolarAnnularTotal
else
SolEc^.EType := etSolarTotal;
end;
SolEc^.Magnitude := -1;
if (Gamma > 0) then
SolEc^.Hemisphere := htNorthern
else
SolEc^.Hemisphere := htSouthern;
FillChar(SolEc^.LContacts, SizeOf(TStContactTimes), #0);
SolEc^.LContacts.MidEclipse := AJDtoDateTime(CentralJD);
GetBesselianElements(CentralJD);
{find limits - then go one step inside at each end}
I1 := 1;
while (sqr(FBesselianElements[I1].XAxis) +
sqr(FBesselianElements[I1].YAxis) >= 1.0) and (I1 < 25) do
Inc(I1);
Inc(I1);
I2 := I1;
repeat
if (I2 < 25) then begin
if (sqr(FBesselianElements[I2+1].XAxis) +
sqr(FBesselianElements[I2+1].YAxis) < 1) then
Inc(I2)
else
break;
end;
until (sqr(FBesselianElements[I2].XAxis) +
sqr(FBesselianElements[I2].YAxis) >= 1) or (I2 >= 25);
Dec(I2);
{this test accounts for non-central eclipses, i.e., those that are total}
{and/or annular but the axis of the moon's shadow does not touch the earth}
if (I1 <> I2) and (I1 < 26) and (I2 < 26) then begin
SolEc^.Path := TStList.Create(TStListNode);
SolEc^.Path.DisposeData := DisposePathData;
GetShadowPath(I1, I2, SolEc^.Path);
end else
SolEc^.Path := nil;
Self.Append(SolEc);
end;
procedure TStEclipses.PartialSolarEclipse(CentralJD, Mu, Gamma : Double);
var
SolEc : PStEclipseRecord;
begin
New(SolEc);
SolEc^.EType := etSolarPartial;
SolEc^.Magnitude := (1.5433 + Mu - abs(Gamma)) / (0.5461 + 2*Mu);
if (Gamma > 0) then
SolEc^.Hemisphere := htNorthern
else
SolEc^.Hemisphere := htSouthern;
FillChar(SolEc^.LContacts, SizeOf(TStContactTimes), #0);
SolEc^.LContacts.MidEclipse := AJDToDateTime(CentralJD);
SolEc^.Path := nil;
Self.Append(SolEc);
end;
end.
Reference in New Issue View Git Blame Copy Permalink