📄 stastro.pas
字号:
var
JD,
TD,
JCent,
JC2, JC3, JC4,
LP, D,
M, MP,
F, I,
A1,A2,A3,
MoonLong,
MoonLat,
MoonDst,
S1,C1,
SunRA,
SunDC,
EE,EES : Double;
SP : TStSunXYZRec;
begin
JD := AstJulianDate(UT.D) + UT.T/86400;
JCent := (JD - 2451545) / 36525;
JC2 := sqr(JCent);
JC3 := JC2 * JCent;
JC4 := sqr(JC2);
SP := SunPosPrim(UT);
SunRA := StInvTan2(SP.SunX, SP.SunY) * radcor;
SunDC := StInvSin(SP.SunZ / SP.RV) * radcor;
{Lunar mean longitude}
LP := 218.3164591 + (481267.88134236 * JCent)
- (1.3268E-3 * JC2) + (JC3 / 538841) - (JC4 / 65194000);
LP := Frac(LP/360) * 360;
if LP < 0 then
LP := LP + 360;
LP := LP/radcor;
{Lunar mean elongation}
D := 297.8502042 + (445267.1115168 * JCent)
- (1.63E-3 * JC2) + (JC3 / 545868) - (JC4 / 113065000);
D := Frac(D/360) * 360;
if D < 0 then
D := D + 360;
D := D/radcor;
{Solar mean anomaly}
M := 357.5291092 + (35999.0502909 * JCent)
- (1.536E-4 * JC2) + (JC3 / 24490000);
M := Frac(M/360) * 360;
if M < 0 then
M := M + 360;
M := M/radcor;
{Lunar mean anomaly}
MP := 134.9634114 + (477198.8676313 * JCent)
+ (8.997E-3 * JC2) + (JC3 / 69699) - (JC4 / 14712000);
MP := Frac(MP/360) * 360;
if MP < 0 then
MP := MP + 360;
MP := MP/radcor;
{Lunar argument of latitude}
F := 93.2720993 + (483202.0175273 * JCent)
- (3.4029E-3 * JC2) - (JC3 / 3526000) + (JC4 / 863310000);
F := Frac(F/360) * 360;
if F < 0 then
F := F + 360;
F := F/radcor;
{Other required arguments}
A1 := 119.75 + (131.849 * JCent);
A1 := Frac(A1/360) * 360;
if A1 < 0 then
A1 := A1 + 360;
A1 := A1/radcor;
A2 := 53.09 + (479264.290 * JCent);
A2 := Frac(A2/360) * 360;
if A2 < 0 then
A2 := A2 + 360;
A2 := A2/radcor;
A3 := 313.45 + (481266.484 * JCent);
A3 := Frac(A3/360) * 360;
if A3 < 0 then
A3 := A3 + 360;
A3 := A3/radcor;
{Earth's orbital eccentricity}
EE := 1.0 - (2.516E-3 * JCent) - (7.4E-6 * JC2);
EES := sqr(EE);
MoonLong := 6288774 * sin(MP)
+ 1274027 * sin(2*D - MP)
+ 658314 * sin(2*D)
+ 213618 * sin(2*MP)
- 185116 * sin(M) * EE
- 114332 * sin(2*F)
+ 58793 * sin(2*(D-MP))
+ 57066 * sin(2*D-M-MP) * EE
+ 53322 * sin(2*D-MP)
+ 45758 * sin(2*D-M) * EE
- 40923 * sin(M-MP) * EE
- 34720 * sin(D)
- 30383 * sin(M+MP) * EE
+ 15327 * sin(2*(D-F))
- 12528 * sin(MP+2*F)
+ 10980 * sin(MP-2*F)
+ 10675 * sin(4*D-MP)
+ 10034 * sin(3*MP)
+ 8548 * sin(4*D-2*MP)
- 7888 * sin(2*D+M-MP) * EE
- 6766 * sin(2*D+M) * EE
- 5163 * sin(D-MP)
+ 4987 * sin(D+M) * EE
+ 4036 * sin(2*D-M+MP) * EE
+ 3994 * sin(2*(D+MP))
+ 3861 * sin(4*D)
+ 3665 * sin(2*D-3*MP)
- 2689 * sin(M-2*MP) * EE
- 2602 * sin(2*D-MP+2*F)
+ 2390 * sin(2*D-M-2*MP) * EE
- 2348 * sin(D-MP)
+ 2236 * sin(2*(D-M)) * EES
- 2120 * sin(M-2*MP) * EE
- 2069 * sin(2*M) * EE
+ 2048 * sin(2*D-2*M-MP) * EES
- 1773 * sin(2*D+MP-2*F)
- 1595 * sin(2*(D+F))
+ 1215 * sin(4*D-M-MP) * EE
- 1110 * sin(2*(MP+F))
- 892 * sin(3*D-MP)
- 810 * sin(2*D-M-MP) * EE
+ 759 * sin(4*D-M-2*MP) * EE
- 713 * sin(2*M-MP) * EE
- 700 * sin(2*D+2*M-MP) * EES
+ 691 * sin(2*D+M-2*MP) * EE
+ 596 * sin(2*D-M-2*F) * EE
+ 549 * sin(4*D+MP)
+ 537 * sin(4*MP)
+ 520 * sin(4*D-M) * EE;
MoonDst := - 20905355 * cos(MP)
- 3699111 * cos(2*D - MP)
- 2955968 * cos(2*D)
- 569925 * cos(2*MP)
+ 48888 * cos(M) * EE
- 3149 * cos(2*F)
+ 246158 * cos(2*(D-MP))
- 152138 * cos(2*D-M-MP) * EE
- 170733 * cos(2*D-MP)
- 204586 * cos(2*D-M) * EE
- 129620 * cos(M-MP) * EE
+ 108743 * cos(D)
+ 104755 * cos(M-MP) * EE
+ 10321 * cos(2*D-2*F)
+ 79661 * cos(MP-2*F)
- 34782 * cos(4*D-MP)
- 23210 * cos(3*MP)
- 21636 * cos(4*D-2*MP)
+ 24208 * cos(2*D+M-MP) * EE
+ 30824 * cos(2*D-M) * EE
- 8379 * cos(D-MP)
- 16675 * cos(D+M) * EE
- 12831 * cos(2*D-M+MP) * EE
- 10445 * cos(2*D+2*MP)
- 11650 * cos(4*D) * EE
+ 14403 * cos(2*D+3*MP)
- 7003 * cos(M-2*MP) * EE
+ 10056 * cos(2*D-M-2*MP) * EE
+ 6322 * cos(D+MP)
- 9884 * cos(2*D-2*M) * EES
+ 5751 * cos(M+2*MP) * EE
- 4950 * cos(2*D-2*M-MP) * EES
+ 4130 * cos(2*D+MP+2*F)
- 3958 * cos(4*D-M-MP) * EE
+ 3258 * cos(3*D-MP)
+ 2616 * cos(2*D+M+MP) * EE
- 1897 * cos(4*D-M-2*MP) * EE
- 2117 * cos(2*M-MP) * EES
+ 2354 * cos(2*D+2*M-MP) * EES
- 1423 * cos(4*D+MP)
- 1117 * cos(4*MP)
- 1571 * cos(4*D-M) * EE
- 1739 * cos(D-2*MP)
- 4421 * cos(2*MP-2*F)
+ 1165 * cos(2*M+MP)
+ 8752 * cos(2*D-MP-2*F);
MoonLat := 5128122 * sin(F)
+ 280602 * sin(MP+F)
+ 277693 * sin(MP-F)
+ 173237 * sin(2*D-F)
+ 55413 * sin(2*D-MP+F)
+ 46271 * sin(2*D-MP-F)
+ 32573 * sin(2*D+F)
+ 17198 * sin(2*MP+F)
+ 9266 * sin(2*D+MP-F)
+ 8822 * sin(2*MP-F)
+ 8216 * sin(2*D-M-F) * EE
+ 4324 * sin(2*D-2*MP-F)
+ 4200 * sin(2*D+MP+F)
- 3359 * sin(2*D+M-F) * EE
+ 2463 * sin(2*D-M-MP+F) * EE
+ 2211 * sin(2*D-M+F) * EE
+ 2065 * sin(2*D-M-MP-F) * EE
- 1870 * sin(M-MP-F) * EE
+ 1828 * sin(4*D-MP-F)
- 1794 * sin(M+F) * EE
- 1749 * sin(3*F)
- 1565 * sin(M-MP+F) * EE
- 1491 * sin(D+F)
- 1475 * sin(M+MP+F) * EE
- 1410 * sin(M+MP-F) * EE
- 1344 * sin(M-F) * EE
- 1335 * sin(D-F)
+ 1107 * sin(3*MP+F)
+ 1021 * sin(4*D-F)
+ 833 * sin(4*D-MP+F)
+ 777 * sin(MP-3*F)
+ 671 * sin(4*D-2*MP+F)
+ 607 * sin(2*D-3*F)
+ 596 * sin(2*D+2*MP-F)
+ 491 * sin(2*D-M+MP-F) * EE
- 451 * sin(2*D-2*MP+F)
+ 439 * sin(3*MP-F)
+ 422 * sin(2*D+2*MP+F)
+ 421 * sin(2*D-3*MP-F)
- 366 * sin(2*D+M-MP+F) * EE
- 351 * sin(2*D+M+F) * EE
+ 331 * sin(4*D+F)
+ 315 * sin(2*D-M+MP+F) * EE
+ 302 * sin(2*D-2*M-F) * EES
- 283 * sin(MP + 3*F)
- 229 * sin(2*D+M+MP-F) * EE
+ 223 * sin(D+M-F) * EE
+ 223 * sin(D+M+F) * EE;
MoonLong := MoonLong
+ 3958 * sin(A1)
+ 1962 * sin(LP-F)
+ 318 * sin(A2);
MoonLat := MoonLat
- 2235 * sin(LP)
+ 382 * sin(A3)
+ 175 * sin(A1-F)
+ 175 * sin(A1+F)
+ 127 * sin(LP-MP)
- 115 * sin(LP+MP);
MoonLong := LP + MoonLong/1000000/radcor;
MoonLat := MoonLat/1000000/radcor;
MoonDst := 385000.56 + MoonDst/1000;
Result.Plx := StInvSin(6378.14/MoonDst) * radcor;
Result.Dia := 358473400 / MoonDst * 2 / 3600;
S1 := sin(MoonLong) * cos(OB2000) - StTan(MoonLat) * sin(OB2000);
C1 := cos(MoonLong);
Result.RA := StInvTan2(C1, S1) * radcor;
TD := sin(MoonLat) * cos(OB2000)
+ cos(MoonLat) * sin(OB2000) * sin(MoonLong);
TD := StInvSin(TD);
Result.DC := TD * radcor;
I := sin(SunDC/radcor) * sin(TD)
+ cos(SunDC/radcor) * cos(TD) * cos((SunRA-Result.RA)/radcor);
Result.Phase := (1 - I) / 2;
I := StInvCos(I) * radcor;
Result.Elong := (Result.RA - SunRA);
if Result.Elong < 0 then
Result.Elong := 360 + Result.Elong;
if Result.Elong >= 180 then begin
Result.Phase := -Result.Phase; {waning moon}
Result.Elong := -I
end else
Result.Elong := I;
end;
function AmountOfSunlight(LD : TStDate; Longitude, Latitude : Double): TStTime;
{-compute the hours, min, sec of sunlight on a given date}
var
RS : TStRiseSetRec;
begin
RS := SunRiseSet(LD, Longitude, Latitude);
with RS do begin
if ORise = -3 then begin
{sun is always above the horizon}
Result := SecondsInDay;
Exit;
end;
if ORise = -2 then begin
{sun is always below horizon}
Result := 0;
Exit;
end;
if (ORise > -1) then begin
if (OSet > -1) then
Result := OSet - ORise
else
Result := SecondsInDay - ORise;
end else begin
if (OSet > -1) then
Result := OSet
else
Result := 0;
end;
end;
if (Result < 0) then
Result := Result + SecondsInDay
else if (Result >= SecondsInDay) then
Result := Result - SecondsInDay;
end;
function SunPos(UT : TStDateTimeRec) : TStPosRec;
{-compute the RA/Declination of the Sun}
var
SP : TStSunXYZRec;
begin
if not CheckDate(UT) then begin
Result.RA := -1;
Result.DC := -99;
Exit;
end;
SP := SunPosPrim(UT);
Result.RA := StInvTan2(SP.SunX, SP.SunY) * radcor;
Result.DC := StInvSin(SP.SunZ / SP.RV) * radcor;
end;
function RiseSetPrim(LD : TStDate;
Longitude, Latitude, H0 : Double;
PR : TStPosRecArray;
ApproxOnly : Boolean) : TStRiseSetRec;
{-primitive routine for finding rise/set times}
var
ST,
NST,
HA,
LatR,
N1,
N2,
N3,
TTran,
TRise,
TSet,
TV1,
TV2,
A1,
A2,
DeltaR,
DeltaS,
RA,
DC,
Alt : Double;
ICount : SmallInt;
UT : TStDateTimeRec;
begin
H0 := H0/radcor;
UT.D := LD;
UT.T := 0;
ST := SiderealTime(UT);
LatR := Latitude/radcor;
{check if object never rises/sets}
N1 := sin(H0) - sin(LatR) * sin(PR[2].DC/radcor);
N2 := cos(LatR) * cos(PR[2].DC/radcor);
HA := N1 / N2;
if (abs(HA) >= 1) then begin
{ if ((Latitude - 90) >= 90) then begin}
if (HA > 0) then begin
{object never rises}
Result.ORise := -2;
Result.OSet := -2;
end else begin
{object never sets, i.e., it is circumpolar}
Result.ORise := -3;
Result.OSet := -3;
end;
Exit;
end;
HA := StInvCos(HA) * radcor;
if HA > 180 then
HA := HA - 180;
if HA < 0 then
HA := HA + 180;
{compute approximate hour angle at transit}
TTran := (PR[2].RA - Longitude - ST) / 360;
if abs(TTran) >= 1 then
TTran:= Frac(TTran);
if TTran < 0 then
TTran := TTran + 1;
TRise := TTran - HA/360;
TSet := TTran + HA/360;
if abs(TRise) >= 1 then
TRise:= Frac(TRise);
if TRise < 0 then
TRise := TRise + 1;
if abs(TSet) >= 1 then
TSet := Frac(TSet);
if TSet < 0 then
TSet := TSet + 1;
if not ApproxOnly then begin
{refine rise time by interpolation/iteration}
ICount := 0;
TV1 := 0;
A1 := 0;
repeat
NST := ST + 360.985647 * TRise;
NST := Frac(NST / 360.0) * 360;
N1 := PR[2].RA - PR[1].RA;
N2 := PR[3].RA - PR[2].RA;
N3 := N2 - N1;
RA := PR[2].RA + TRise/2 * (N1 + N2 + TRise*N3);
N1 := PR[2].DC - PR[1].DC;
N2 := PR[3].DC - PR[2].DC;
N3 := N2 - N1;
DC := PR[2].DC + TRise/2 * (N1 + N2 + TRise*N3);
DC := DC/radcor;
HA := (NST + Longitude - RA) / radcor;
Alt := StInvSin(sin(LatR) * sin(DC) + cos(LatR) * cos(DC) * cos(HA));
DeltaR := ((Alt - H0) * radcor) / (360 * cos(DC) * cos(LatR) * sin(HA));
TRise := TRise + DeltaR;
Inc(ICount);
if (ICount > 3) and (abs(DeltaR) >= 0.0005) then begin
if (ICount = 4) then begin
TV1 := TRise;
A1 := (Alt-H0) * radcor;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -