📄 stastro.pas
字号:
end else if (ICount = 5) then begin
TV2 := TRise;
A2 := (Alt-H0) * radcor;
TRise := TV1 + (A1 / A2) * (TV1 - TV2);
break;
end;
end;
until (abs(DeltaR) < 0.0005); {0.0005d = 0.72 min}
{refine set time by interpolation/iteration}
ICount := 0;
TV1 := 0;
A1 := 0;
repeat
NST := ST + 360.985647 * TSet;
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 + TSet/2 * (N1 + N2 + TSet*N3);
N1 := PR[2].DC - PR[1].DC;
N2 := PR[3].DC - PR[2].DC;
N3 := N2 - N1;
DC := PR[2].DC + TSet/2 * (N1 + N2 + TSet*N3);
DC := DC/radcor;
HA := (NST + Longitude - RA) / radcor;
Alt := StInvSin(sin(LatR) * sin(DC) + cos(LatR) * cos(DC) * cos(HA));
DeltaS := ((Alt - H0) * radcor) / (360 * cos(DC) * cos(LatR) * sin(HA));
TSet := TSet + DeltaS;
Inc(ICount);
if (ICount > 3) and (abs(DeltaS) >= 0.0005) then begin
if (ICount = 4) then begin
TV1 := TSet;
A1 := (Alt-H0) * radcor;
end else if (ICount = 5) then begin
TV2 := TSet;
A2 := (Alt-H0) * radcor;
TSet := TV1 + (A1 / A2) * (TV1 - TV2);
break;
end;
end;
until (abs(DeltaS) < 0.0005); {0.0005d = 0.72 min}
end;
if (TRise >= 0) and (TRise < 1) then
Result.ORise := Trunc(TRise * SecondsInDay)
else begin
if TRise < 0 then
Result.ORise := Trunc((TRise + 1) * SecondsInDay)
else
Result.ORise := Trunc(Frac(TRise) * SecondsInDay);
end;
if Result.ORise < 0 then
Inc(Result.ORise, SecondsInDay);
if Result.ORise >= SecondsInDay then
Dec(Result.ORise, SecondsInDay);
if (TSet >= 0) and (TSet < 1) then
Result.OSet := Trunc(TSet * SecondsInDay)
else begin
if TSet < 0 then
Result.OSet := Trunc((TSet + 1) * SecondsInDay)
else
Result.OSet := Trunc(Frac(TSet) * SecondsInDay);
end;
if Result.OSet < 0 then
Inc(Result.OSet, SecondsInDay);
if Result.OSet >= SecondsInDay then
Dec(Result.OSet, SecondsInDay);
end;
function SunRiseSet(LD : TStDate; Longitude, Latitude : Double) : TStRiseSetRec;
{-compute the Sun rise or set time}
{the value for H0 accounts for approximate refraction of 0.5667 deg. and}
{that rise or set is based on the upper limb instead of the center of the solar disc}
var
I : Integer;
H0 : Double;
UT : TStDateTimeRec;
RP : TStPosRecArray;
begin
Dec(LD);
H0 := -0.8333;
UT.T := 0;
UT.D := LD-1;
if CheckDate(UT) then begin
UT.D := UT.D + 2;
if (not CheckDate(UT)) then begin
Result.ORise := -4;
Result.OSet := -4;
Exit;
end else
UT.D := UT.D-2;
end else begin
Result.ORise := -4;
Result.OSet := -4;
Exit;
end;
for I := 1 to 3 do begin
RP[I] := SunPos(UT);
if I >= 2 then begin
if RP[I].RA < RP[I-1].RA then
RP[I].RA := RP[I].RA + 360;
end;
Inc(UT.D);
end;
Result := RiseSetPrim(LD, Longitude, Latitude, H0, RP, False);
end;
function Twilight(LD : TStDate; Longitude, Latitude : Double;
TwiType : TStTwilight) : TStRiseSetRec;
{-compute the beginning or end of twilight}
{twilight computations are based on the zenith distance of the center }
{of the solar disc.}
{Civil = 6 deg. below the horizon}
{Nautical = 12 deg. below the horizon}
{Astronomical = 18 deg. below the horizon}
var
I : Integer;
H0 : Double;
UT : TStDateTimeRec;
RP : TStPosRecArray;
begin
UT.D := LD-1;
UT.T := 0;
if CheckDate(UT) then begin
UT.D := UT.D + 2;
if (not CheckDate(UT)) then begin
Result.ORise := -4;
Result.OSet := -4;
Exit;
end else
UT.D := UT.D-2;
end else begin
Result.ORise := -4;
Result.OSet := -4;
Exit;
end;
case TwiType of
ttCivil : H0 := -6.0;
ttNautical : H0 := -12.0;
ttAstronomical : H0 := -18.0;
else
H0 := -18.0;
end;
for I := 1 to 3 do begin
UT.D := LD + I-1;
RP[I] := SunPos(UT);
if (I > 1) then begin
if RP[I].RA < RP[I-1].RA then
RP[I].RA := RP[I].RA + 360.0;
end;
end;
Result := RiseSetPrim(LD, Longitude, Latitude, H0, RP, False);
end;
function FixedRiseSet(LD : TStDate;
RA, DC, Longitude, Latitude : Double) : TStRiseSetRec;
{-compute the rise/set time for a fixed object, e.g., star}
{the value for H0 accounts for approximate refraction of 0.5667 deg.}
{this routine does not refine the intial estimate and so may be off by five}
{minutes or so}
var
H0 : Double;
UT : TStDateTimeRec;
RP : TStPosRecArray;
begin
H0 := -0.5667;
UT.T := 0;
UT.D := LD;
if not CheckDate(UT) then begin
Result.ORise := -4;
Result.OSet := -4;
Exit;
end;
RP[2].RA := RA;
RP[2].DC := DC;
Result := RiseSetPrim(LD, Longitude, Latitude, H0, RP, True);
end;
function MoonPos(UT : TStDateTimeRec) : TStMoonPosRec;
{-compute the J2000 RA/Declination of the moon}
begin
if not CheckDate(UT) then begin
Result.RA := -1;
Result.DC := -1;
Exit;
end;
Result := MoonPosPrim(UT);
end;
function MoonRiseSet(LD : TStDate; Longitude, Latitude : Double) : TStRiseSetRec;
{compute the Moon rise and set time}
{the value for H0 accounts for approximate refraction of 0.5667 deg., }
{that rise or set is based on the upper limb instead of the center of the}
{lunar disc, and the lunar parallax. In accordance with American Ephemeris }
{practice, the phase of the moon is not taken into account, i.e., the time}
{is based on the upper limb whether it is lit or not}
var
I : Integer;
H0 : Double;
UT : TStDateTimeRec;
RP : TStPosRecArray;
MPR : TStMoonPosRec;
begin
H0 := 0.125; { default value }
Dec(LD);
UT.T := 0;
UT.D := LD;
if CheckDate(UT) then begin
UT.D := UT.D + 2;
if (not CheckDate(UT)) then begin
Result.ORise := -4;
Result.OSet := -4;
Exit;
end else
UT.D := UT.D-2;
end else begin
Result.ORise := -4;
Result.OSet := -4;
Exit;
end;
for I := 1 to 3 do begin
MPR := MoonPos(UT);
RP[I].RA := MPR.RA;
RP[I].DC := MPR.DC;
if I >= 2 then begin
if I = 2 then
H0 := 0.7275*MPR.Plx - 0.5667;
if RP[I].RA < RP[I-1].RA then
RP[I].RA := RP[I].RA + 360;
end;
Inc(UT.D);
end;
Result := RiseSetPrim(LD, Longitude, Latitude, H0, RP, False);
end;
function LunarPhase(UT : TStDateTimeRec) : Double;
{-compute the phase of the moon}
{The value is positive if between New and Full Moon}
{ negative if between Full and New Moon}
var
MPR : TStMoonPosRec;
begin
MPR := MoonPosPrim(UT);
Result := MPR.Phase;
end;
procedure GetPhases(K : Double; var PR : TStPhaseRecord);
{primitive routine to find the date/time of phases in a lunar cycle}
var
JD,
NK,
TD,
J1,
J2,
J3 : Double;
step : Integer;
E,
FP,
S1,
M1,
M2,
M3 : Double;
function AddCor : Double;
begin
AddCor := 0.000325 * sin((299.77 + 0.107408*K - 0.009173*J2)/radcor)
+ 0.000165 * sin((251.88 + 0.016321*K)/radcor)
+ 0.000164 * sin((251.83 + 26.651886*K)/radcor)
+ 0.000126 * sin((349.42 + 36.412478*K)/radcor)
+ 0.000110 * sin((84.660 + 18.206239*K)/radcor)
+ 0.000062 * sin((141.74 + 53.303771*K)/radcor)
+ 0.000060 * sin((207.14 + 2.453732*K)/radcor)
+ 0.000056 * sin((154.84 + 7.306860*K)/radcor)
+ 0.000047 * sin((34.520 + 27.261239*K)/radcor)
+ 0.000042 * sin((207.19 + 0.121824*K)/radcor)
+ 0.000040 * sin((291.34 + 1.844379*K)/radcor)
+ 0.000037 * sin((161.72 + 24.198154*K)/radcor)
+ 0.000035 * sin((239.56 + 25.513099*K)/radcor)
+ 0.000023 * sin((331.55 + 3.592518*K)/radcor);
end;
begin
NK := K;
FillChar(PR, SizeOf(TStPhaseRecord), #0);
for step := 0 to 3 do begin
K := NK + (step*0.25);
FP := Frac(K);
if FP < 0 then
FP := FP + 1;
{compute Julian Centuries}
J1 := K / 1236.85;
J2 := Sqr(J1);
J3 := J2 * J1;
{solar mean anomaly}
S1 := 2.5534
+ 29.1053569 * K
- 0.0000218 * J2
- 0.00000011 * J3;
S1 := Frac(S1 / 360.0) * 360;
if S1 < 0 then
S1 := S1 + 360.0;
{lunar mean anomaly}
M1 := 201.5643
+ 385.81693528 * K
+ 0.0107438 * J2
+ 0.00001239 * J3
- 0.000000058 * J2 * J2;
M1 := Frac(M1 / 360.0) * 360;
if M1 < 0 then
M1 := M1 + 360.0;
{lunar argument of latitude}
M2 := 160.7108
+ 390.67050274 * K
- 0.0016341 * J2
- 0.00000227 * J3
+ 0.000000011 * J2 * J2;
M2 := Frac(M2 / 360.0) * 360;
if M2 < 0 then
M2 := M2 + 360.0;
{lunar ascending node}
M3 := 124.7746
- 1.56375580 * K
+ 0.0020691 * J2
+ 0.00000215 * J3;
M3 := Frac(M3 / 360.0) * 360;
if M3 < 0 then
M3 := M3 + 360.0;
{convert to radians}
S1 := S1/radcor;
M1 := M1/radcor;
M2 := M2/radcor;
M3 := M3/radcor;
{mean Julian Date for phase}
JD := 2451550.09765
+ 29.530588853 * K
+ 0.0001337 * J2
- 0.000000150 * J3
+ 0.00000000073 * J2 * J2;
{earth's orbital eccentricity}
E := 1.0 - 0.002516 * J1 - 0.0000074 * J2;
{New Moon date time}
if FP < 0.01 then begin
TD := - 0.40720 * sin(M1)
+ 0.17241 * E * sin(S1)
+ 0.01608 * sin(2*M1)
+ 0.01039 * sin(2*M2)
+ 0.00739 * E * sin(M1-S1)
- 0.00514 * E * sin(M1 + S1)
+ 0.00208 * E * E * sin(2 * S1)
- 0.00111 * sin(M1 - 2*M2)
- 0.00057 * sin(M1 + 2*M2)
+ 0.00056 * E * sin(2*M1 + S1)
- 0.00042 * sin(3 * M1)
+ 0.00042 * E * sin(S1 + 2*M2)
+ 0.00038 * E * sin(S1 - 2*M2)
- 0.00024 * E * sin(2*(M1-S1))
- 0.00017 * sin(M3)
- 0.00007 * sin(M1 + 2*S1);
JD := JD + TD + AddCor;
PR.NMDate := JD;
end;
{Full Moon date/time}
if Abs(FP - 0.5) < 0.01 then begin
TD := - 0.40614 * sin(M1)
+ 0.17302 * E * sin(S1)
+ 0.01614 * sin(2*M1)
+ 0.01043 * sin(2*M2)
+ 0.00734 * E * sin(M1-S1)
- 0.00515 * E * sin(M1 + S1)
+ 0.00209 * E * E * sin(2 * S1)
- 0.00111 * sin(M1 - 2*M2)
- 0.00057 * sin(M1 + 2*M2)
+ 0.00056 * E * sin(2*M1 + S1)
- 0.00042 * sin(3 * M1)
+ 0.00042 * E * sin(S1 + 2*M2)
+ 0.00038 * E * sin(S1 - 2*M2)
- 0.00024 * E * sin(2*(M1-S1))
- 0.00017 * sin(M3)
- 0.00007 * sin(M1 + 2*S1);
JD := JD + TD + AddCor;
PR.FMDate := JD;
end;
{Quarters date/time}
if (abs(FP - 0.25) < 0.01) or (abs(FP - 0.75) < 0.01) then begin
TD := - 0.62801 * sin(M1)
+ 0.17172 * sin(S1) * E
- 0.01183 * sin(M1+S1) * E
+ 0.00862 * sin(2*M1)
+ 0.00804 * sin(2*M2)
+ 0.00454 * sin(M1-S1) * E
+ 0.00204 * sin(2*S1) * E * E
- 0.00180 * sin(M1 - 2*M2)
- 0.00070 * sin(M1 + 2*M2)
- 0.00040 * sin(3*M1)
- 0.00034 * sin(2*M1-S1) * E
+ 0.00032 * sin(S1 + 2*M2) * E
+ 0.00032 * sin(S1 - 2*M2) * E
- 0.00028 * sin(M1 + 2*S1) * E * E
+ 0.00027 * sin(2*M1 + S1) * E
- 0.00017 * sin(M3)
- 0.00005 * sin(M1 - S1 - 2*M2);
JD := JD + TD + AddCor;
{adjustment to computed Julian Date}
TD := 0.00306
- 0.00038 * E * cos(S1)
+ 0.00026 * cos(M1)
- 0.00002 * cos(M1-S1)
+ 0.00002 * cos(M1+S1)
+ 0.00002 * cos(2*M2);
if Abs(FP - 0.25) < 0.01 then
PR.FQDate := JD + TD
else
PR.LQDate := JD - TD;
end;
end;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -