📄 steclpse.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: StEclpse.pas 4.03 *}
{*********************************************************}
{* SysTools: Lunar/Solar Eclipses *}
{*********************************************************}
{$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
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -