Алгоритм для определения восхода/захода солнца и луны - Алгоритм для определения восхода/захода солнца и луны. Вариант 2

Written on . Posted in [Delphi] Алгоритмы

ОГЛАВЛЕНИЕ

Страница 2 из 2
 

Вариант 2.

 

{-------------------}

{$N+}
program Sun;

const
P1       = 3.14159265358;
P2       = 2*P1;
DR       = P1/180;
K1       = 15*DR*1.0027379;

type
TTime = record
Hour, Min, Sec: Extended;
end;
TDate = record
Year, Month, Day: Extended;
end;

type
RequestBlock = record
Latitude,
Longitude: Extended;
HourZone: Word;
Date: TDate;
end;
ReplyBlock = record
SunRise: TTime;
SunRiseAzimuth: Extended;
SunSet: TTime;
SunSetAzimuth: Extended;
end;

function Sign(val: Extended): Integer;
begin
if val < 0 then Sign := -1;
if val > 0 then Sign := 1;
if val = 0 then Sign := 0;
end;

var
A_MATR, D_MATR: array [1..2] of Extended;
B5, L5: Extended; {Широта и долгота}
H: Extended; {Часовая зона}
z0, z1, z: Extended;
g: Extended;
D1, f, J, J3, S, C, A, B, D, E: Extended;
T, TT, T0: Extended;
A5, D5, R5: Extended;
M8, W8: Extended;
A0, D0: Extended;
dA, dD: Extended;
c0: Integer;
p: Extended;
A2, D2: Extended;
L0, L2, H0, H2, H1, t3: Extended;
V0, V1, V2: Extended;
H3, M3: Extended;
H7, N7, D7, AZ: Extended;


procedure ComputeVars;
{
Фундаментальные константы
(Van Flandern & Pulkkinen, 1979)
}
var
L, G, V, U, W: Extended;
begin
L := 0.779072 + 0.00273790931 * T;
G := 0.993126 + 0.0027377785 * T;
L := L - Int(L);
G := G - Int(G);
L := L * P2;
G := G * P2;
V := 0.39785 * Sin(L);
V := V - 0.01 * Sin(L - G);
V := V + 0.00333 * Sin(L + G);
V := V - 0.00021 * TT * Sin(L);
U := 1 - 0.03349 * Cos(G);
U := U - 0.00014 * Cos(2 * L);
U := U + 0.00008 * Cos(L);
W := -0.0001 - 0.04129 * Sin(2 * L);
W := W + 0.03211 * Sin(G);
W := W + 0.00104 * Sin(2 * L - G);
W := W - 0.00035 * Sin(2 * L + G);
W := W - 0.00008 * TT * Sin(G);
{ Вычисление солнечных координат }
S := W / Sqrt(U - V * V);
A5 := L + ArcTan(S / Sqrt(1 - S * S));
S  := V / Sqrt(U);
D5 := ArcTan(S / Sqrt(1 - S * S));
R5 := 1.00021 * Sqrt(U);
end;


function ComputeSunTime(Rq: RequestBlock; var Rp: ReplyBlock): Byte;
begin
with Rq do
begin
L5 := Longitude/360;
z0 := HourZone /24;
G := 1;
if Date.Year < 1583 then G := 0;
D1 := Int(Date.Day);
f := Date.Day - D1 - 0.5;
J := -Int(7*(Int((Date.Month + 9)/12) + Date.Year)/4);

if (g <> 0) then
begin
S := Sign(Rq.Date.Month - 9);
A := Abs(Date.Month - 9);
J3 := Int(Date.Year + S*Int(A/7));
J3 := -Int((Int(J3/100) + 1)*3/4);
end;
J := J + Int(275*Date.Month/9) + D1 + G*J3;
J := J + 1721027 + 2*G + 367*Rq.Date.Year;
if f < 0 then
begin
f := f+1;
J := J-1;
end;

T := (J - 2451545)+f;
TT := T/36525+1;          {TT = столетия, начиная с 1900.0}

{ Получение часового пояса }
T0 := T/36525;
S := 24110.5 + 8640184.813*T0;
S := S + 86636.6*z0 + 86400*L5;
S := S/86400;
S := S - Int(S);
T0 := S*360*DR;

T := T + z0;
{ Получаем положение Солнца }
ComputeVars;
A_MATR[1] := A5;
D_MATR[1] := D5;
T := T + 1;
ComputeVars;
A_MATR[2] := A5;
D_MATR[2] := D5;
if A_MATR[2] < A_MATR[1] then A_MATR[2] := A_MATR[2] + P2;

{ Вычисление зенита }
z1 := DR*90.833;
S := Sin(Latitude*DR);
C := Cos(Latitude*DR);
z := Cos(z1);
M8 := 0;
W8 := 0;
A0 := A_MATR[1];
D0 := D_MATR[1];
dA := A_MATR[2] - A_MATR[1];
dD := D_MATR[2] - D_MATR[1];

for c0 := 0 to 23 do
begin
p := (c0 + 1)/24;
A2 := A_MATR[1] + p*dA;
D2 := D_MATR[1] + p*dD;
{ Просматриваем возможные события на полученный час }
L0 := T0 + c0*K1;
L2 := L0 + K1;
H0 := L0 - A0;
H2 := L2 - A2;
H1 := (H2 + H0)/2; { Часовой угол, }
D1 := (D2 + D0)/2; { наклон в получасе }

if c0 > 0 then
else V0 := S*Sin(D0) + C*Cos(D0)*Cos(H0) - z;
V2 := S*Sin(D2) + C*Cos(D2)*Cos(H2) - z;

if Sign(V0) <> Sign(V2) then
begin
V1 := S*Sin(D1) + C*Cos(D1)*Cos(H1) - z;
A  := 2*V2 - 4*V1 + 2*V0;
B  := 4*V1 - 3*V0 - V2;
D  := B*B - 4*A*V0;
if D >= 0 then
begin
D := Sqrt(D);
E := (-B + D)/(2*A);
if (E > 1) or (E < 0) then E := (-B - D)/(2*A);
t3 := c0 + E + 1/120; { округление }
H3 := Int(t3);
M3 := Int((t3 - H3)*60);

H7 := H0 + E*(H2 - H0);
N7 := -Cos(D1)*Sin(H7);
D7 := C*Sin(D1) - S*Cos(D1)*Cos(H7);
AZ := ArcTan(N7/D7)/DR;
if D7 < 0 then AZ := AZ + 180;
if AZ < 0 then AZ := AZ + 360;
if AZ > 360 then AZ := AZ - 360;

if (V0 < 0) and (V2 > 0) then
begin
Rp.SunRise.Hour := Trunc(H3);
Rp.SunRise.Min := Trunc(M3);
Rp.SunRiseAzimuth := AZ;
M8 := 1;
end;
if (V0 > 0) and (V2 < 0) then
begin
Rp.SunSet.Hour := Trunc(H3);
Rp.SunSet.Min := Trunc(M3);
Rp.SunSetAzimuth := AZ;
W8 := 1;
end;
end;
end;
{ }         A0 := A2;
D0 := D2;
V0 := V2;
end;

{ Вывод информации? }
if (M8 = 0) and (W8 = 0) then
begin
if (V2 < 0) then ComputeSunTime := $A3;
if (V2 > 0) then ComputeSunTime := $A4;
end
else
begin
if (M8 = 0) then ComputeSunTime := $A1;
if (W8 = 0) then ComputeSunTime := $A2;
end;
end;
end;

const
SMessage = 'Заход солнца в ';
RMessage = 'Восход солнца в ';
Message1 = 'В этот день солнце не восходит ';
Message2 = 'В этот день солнце не заходит ';
Message3 = 'Солнце заходит весь день ';
Message4 = 'Солнце восходит весь день ';

var
Rq: RequestBlock;
Rp: ReplyBlock;

begin

with Rq do
begin
Write(' Широта........:'); ReadLn(Latitude);
Write(' Долгота.......:'); ReadLn(Longitude);
Write(' Часовой пояс..:'); ReadLn(HourZone);
Write(' Год...........:'); ReadLn(Date.Year);
Write(' Месяц.........:'); ReadLn(Date.Month);
Write(' День..........:'); ReadLn(Date.Day);
end;
WriteLn;

case ComputeSunTime(Rq, Rp) of
$A1: WriteLn(Message1);
$A2: WriteLn(Message2);
$A3: WriteLn(Message3);
$A4: WriteLn(Message4);
else
Write(RMessage, Rp.SunRise.Hour:0:0, ':', Rp.SunRise.Min:0:0,
';      ');
WriteLn('Азимут: ', Rp.SunRiseAzimuth:2:2);
Write(SMessage, Rp.SunSet.Hour:0:0, ':', Rp.SunSet.Min:0:0, ';
');
WriteLn('Азимут: ', Rp.SunSetAzimuth:2:2);
end;
WriteLn;
end.
delphi форум