Аглоритм (уравнение) для определения восхода/захода солнца и луны (BASIC) - Алгоритмы - Не стандартные примеры на Delphi - Каталог статей

10 ' Восход-заход солнца
20 GOSUB 300
30 INPUT "Долгота (град)";B5,L5
40 INPUT "Часовая зона (час)";H
50 L5=L5/360: Z0=H/24
60 GOSUB 1170: T=(J-2451545)+F
70 TT=T/36525+1: ' TT = столетия,
80 ' начиная с 1900.0
90 GOSUB 410: T=T+Z0
100 '
110 ' Получаем положение солнца
120 GOSUB 910: A(1)=A5: D(1)=D5
130 T=T+1
140 GOSUB 910: A(2)=A5: D(2)=D5
150 IF A(2)<A(1) THEN A(2)=A(2)+P2
160 Z1=DR*90.833: ' Вычисление зенита
170 S=SIN(B5*DR): C=COS(B5*DR)
180 Z=COS(Z1): M8=0: W8=0: PRINT
190 A0=A(1): D0=D(1)
200 DA=A(2)-A(1): DD=D(2)-D(1)
210 FOR C0=0 TO 23
220 P=(C0+1)/24
230 A2=A(1)+P*DA: D2=D(1)+P*DD
240 GOSUB 490
250 A0=A2: D0=D2: V0=V2
260 NEXT
270 GOSUB 820: ' Вывод информации?
280 END
290 '
300 ' Константы
310 DIM A(2),D(2)
320 P1=3.14159265: P2=2*P1
330 DR=P1/180: K1=15*DR*1.0027379
340 S$="Заход солнца в "
350 R$="Восход солнца в "
360 M1$="В этот день солнце не восходит"
370 M2$="В этот день солнце не заходит"
380 M3$="Солнце заходит весь день"
390 M4$="Солнце восходит весь день"
400 RETURN
410 ' Получение часового пояса
420 T0=T/36525
430 S=24110.5+8640184.813*T0
440 S=S+86636.6*Z0+86400*L5
450 S=S/86400: S=S-INT(S)
460 T0=S*360*DR
470 RETURN
480 '
490 ' Просматриваем возможные события на полученный час
500 L0=T0+C0*K1: L2=L0+K1
510 H0=L0-A0: H2=L2-A2
520 H1=(H2+H0)/2: ' Часовой угол,
530 D1=(D2+D0)/2: ' наклон в
540 ' получасе
550 IF C0>0 THEN 570
560 V0=S*SIN(D0)+C*COS(D0)*COS(H0)-Z
570 V2=S*SIN(D2)+C*COS(D2)*COS(H2)-Z
580 IF SGN(V0)=SGN(V2) THEN 800 
590 V1=S*SIN(D1)+C*COS(D1)*COS(H1)-Z
600 A=2*V2-4*V1+2*V0: B=4*V1-3*V0-V2
610 D=B*B-4*A*V0: IF D<0 THEN 800
620 D=SQR(D)
630 IF V0<0 AND V2>0 THEN PRINT R$;
640 IF V0<0 AND V2>0 THEN M8=1
650 IF V0>0 AND V2<0 THEN PRINT S$;
660 IF V0>0 AND V2<0 THEN W8=1
670 E=(-B+D)/(2*A)
680 IF E>1 OR E<0 THEN E=(-B-D)/(2*A)
690 T3=C0+E+1/120: ' Округление
700 H3=INT(T3): M3=INT((T3-H3)*60)
710 PRINT USING "##:##";H3;M3;
720 H7=H0+E*(H2-H0)
730 N7=-COS(D1)*SIN(H7)
740 D7=C*SIN(D1)-S*COS(D1)*COS(H7)
750 AZ=ATN(N7/D7)/DR
760 IF D7<0 THEN AZ=AZ+180
770 IF AZ<0 THEN AZ=AZ+360
780 IF AZ>360 THEN AZ=AZ-360
790 PRINT USING ", азимут ###.#";AZ
800 RETURN
810 '
820 ' Процедура вывода информации
830 IF M8=0 AND W8=0 THEN 870
840 IF M8=0 THEN PRINT M1$
850 IF W8=0 THEN PRINT M2$
860 GOTO 890
870 IF V2<0 THEN PRINT M3$
880 IF V2>0 THEN PRINT M4$
890 RETURN
900 '
910 ' Фундаментальные константы
920 ' (Van Flandern &
930 ' Pulkkinen, 1979)
940 L=.779072+.00273790931*T
950 G=.993126+.0027377785*T
960 L=L-INT(L): G=G-INT(G)
970 L=L*P2: G=G*P2
980 V=.39785*SIN(L)
990 V=V-.01000*SIN(L-G)
1000 V=V+.00333*SIN(L+G)
1010 V=V-.00021*TT*SIN(L)
1020 U=1-.03349*COS(G)
1030 U=U-.00014*COS(2*L)
1040 U=U+.00008*COS(L)
1050 W=-.00010-.04129*SIN(2*L)
1060 W=W+.03211*SIN(G)
1070 W=W+.00104*SIN(2*L-G)
1080 W=W-.00035*SIN(2*L+G)
1090 W=W-.00008*TT*SIN(G)
1100 '
1110 ' Вычисление солнечных координат
1120 S=W/SQR(U-V*V)
1130 A5=L+ATN(S/SQR(1-S*S))
1140 S=V/SQR(U):D5=ATN(S/SQR(1-S*S))
1150 R5=1.00021*SQR(U)
1160 RETURN
1165 '
1170 ' Календарь --> JD
1180 INPUT "Год, Месяц, День";Y,M,D
1190 G=1: IF Y<1583 THEN G=0
1200 D1=INT(D): F=D-D1-.5
1210 J=-INT(7*(INT((M+9)/12)+Y)/4)
1220 IF G=0 THEN 1260
1230 S=SGN(M-9): A=ABS(M-9)
1240 J3=INT(Y+S*INT(A/7))
1250 J3=-INT((INT(J3/100)+1)*3/4)
1260 J=J+INT(275*M/9)+D1+G*J3
1270 J=J+1721027+2*G+367*Y
1280 IF F>=0 THEN 1300
1290 F=F+1: J=J-1
1300 RETURN
1310 '
1320 ' Программа вычисляет время восхода и захода
1330 ' солнца по дате (с точностью до минуты) в пределах
1340 ' нескольких текущих столетий. Производит корректировку, если географическая
1350 ' точка находится в арктическом или антарктическом регионе, где заход или восход солнца
1360 ' на текущую дату может не состояться. Вводимые данные: положительная северная широта и
1370 ' отрицательная западная долгота. Часовой пояс указывается относительно Гринвича 
1380 ' (например, 5 для EST и 4 для EDT). Алгоритм обсуждался в
1390 ' "Sky & Telescope" за август 1994, страница 84.

Pascal

Великолепно! Теперь данный алгоритм портирован и в Delphi. Полностью привожу письмо и код модуля Александра Ермолаева, приславшего мне решение этой головоломки:

Аглоритм (уравнение) для определения восхода/захода солнца и луны (Pascal).

Написал я в Делфи. Проект: sunproject.dpr, main.dfm, main.pas. Создайте пустой проект и скопируйте тексты main.dfm, main.pas в соответствующие окошки. Ну сами наверно все знаете, на всякий случай напоминаю.

Александр Ермолаев, swift@glazov.udm.net

sunproject.dpr

program sunproject;

uses

Forms,

main in 'main.pas' {Sun};

{$R *.RES}

begin

Application.Initialize;

Application.Title := 'Sun';

Application.CreateForm(TSun, Sun);

Application.Run;

end.

main.dfm

object Sun: TSun

Left = 210

Top = 106

BorderIcons = [biSystemMenu, biMinimize]

BorderStyle = bsSingle

Caption = 'Sun'

ClientHeight = 257

ClientWidth = 299

Color = clBtnFace

Font.Charset = DEFAULT_CHARSET

Font.Color = clWindowText

Font.Height = -11

Font.Name = 'MS Sans Serif'

Font.Style = []

OldCreateOrder = False

Position = poDesktopCenter

OnCreate = CreateForm

PixelsPerInch = 96

TextHeight = 13

object GroupBoxInput: TGroupBox

Left = 4

Top = 4

Width = 173

Height = 93

Caption = ' Ввод '

TabOrder = 0

object LabelLongitude: TLabel

Left = 35

Top = 44

Width = 78

Height = 13

Alignment = taRightJustify

Caption = 'Долгота (град):'

end

object LabelTimeZone: TLabel

Left = 13

Top = 68

Width = 100

Height = 13

Alignment = taRightJustify

Caption = 'Часовая зона (час):'

end

object LabelAtitude: TLabel

Left = 40

Top = 20

Width = 73

Height = 13

Alignment = taRightJustify

Caption = 'Широта (град):'

end

object EditB5: TEdit

Tag = 1

Left = 120

Top = 16

Width = 37

Height = 21

TabOrder = 0

Text = '0'

end

object EditL5: TEdit

Tag = 2

Left = 120

Top = 40

Width = 37

Height = 21

TabOrder = 1

Text = '0'

end

object EditH: TEdit

Tag = 3

Left = 120

Top = 64

Width = 37

Height = 21

TabOrder = 2

Text = '0'

end

object GroupBoxCalendar: TGroupBox

Left = 184

Top = 4

Width = 109

Height = 93

Caption = ' Календарь '

TabOrder = 1

object LabelD: TLabel

Left = 19

Top = 20

Width = 30

Height = 13

Alignment = taRightJustify

Caption = 'День:'

end

object LabelM: TLabel

Left = 13

Top = 44

Width = 36

Height = 13

Alignment = taRightJustify

Caption = 'Месяц:'

end

object LabelY: TLabel

Left = 28

Top = 68

Width = 21

Height = 13

Alignment = taRightJustify

Caption = 'Год:'

end

object EditD: TEdit

Tag = 1

Left = 56

Top = 16

Width = 37

Height = 21

TabOrder = 0

Text = '0'

end

object EditM: TEdit

Tag = 2

Left = 56

Top = 40

Width = 37

Height = 21

TabOrder = 1

Text = '0'

end

object EditY: TEdit

Tag = 3

Left = 56

Top = 64

Width = 37

Height = 21

TabOrder = 2

Text = '0'

end

object ButtonCalc: TButton

Left = 12

Top = 227

Width = 169

Height = 25

Caption = '&Вычислить'

TabOrder = 2

OnClick = ButtonCalcClick

end

object ListBox: TListBox

Left = 4

Top = 104

Width = 289

Height = 117

ItemHeight = 13

TabOrder = 3

end

object ButtonClear: TButton

Left = 192

Top = 227

Width = 91

Height = 25

Caption = '&Очистить'

TabOrder = 4

OnClick = ButtonClearClick

end

main.pas

{

Программа вычисляет время восхода и захода

солнца по дате (с точностью до минуты) в пределах

нескольких текущих столетий. Производит корректировку, если

географическая

точка находится в арктическом или антарктическом регионе, где заход

или восход солнца

на текущую дату может не состояться. Вводимые данные: положительная

северная широта и

отрицательная западная долгота. Часовой пояс указывается относительно

Гринвича

(например, 5 для EST и 4 для EDT). Алгоритм обсуждался в

"Sky & Telescope" за август 1994, страница 84.

}

unit main;

interface

uses

Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms,

Dialogs,

StdCtrls;

type

TSun = class(TForm)

GroupBoxInput: TGroupBox;

LabelLongitude: TLabel;

EditB5: TEdit;

EditL5: TEdit;

LabelTimeZone: TLabel;

EditH: TEdit;

GroupBoxCalendar: TGroupBox;

LabelD: TLabel;

LabelM: TLabel;

LabelY: TLabel;

EditD: TEdit;

EditM: TEdit;

EditY: TEdit;

ButtonCalc: TButton;

ListBox: TListBox;

ButtonClear: TButton;

LabelAtitude: TLabel;

procedure Calendar; // Календарь

procedure GetTimeZone; // Получение часового пояса

procedure PosOfSun; // Получаем положение солнца

procedure OutInform; // Процедура вывода информации

procedure PossibleEvents(Hour: integer); // Возможные события на

полученный час

procedure GetDate; //Получить значения даты

procedure GetInput; //Получить значения широты,...

procedure ButtonCalcClick(Sender: TObject);

procedure CreateForm(Sender: TObject);

procedure ButtonClearClick(Sender: TObject);

private

function Sgn(Value: Double): integer; // Сигнум

public

{ Public declarations }

end;

var

Sun: TSun;

st: string;

aA, aD: array [1 .. 2] of double;

B5: integer;

L5: double;

H: integer;

Z, Z0, Z1: double;

D: double;

M, Y: integer;

A5, D5, R5: double;

J3: integer;

T, T0, TT, T3: double;

L0, L2: double;

H0, H1, H2, H7, N7, D7: double;

H3, M3: integer;

M8, W8: double;

A, B, A0, D0, A2, D1, D2, DA, DD: double;

E, F, J, S, C, P, L, G, V, U, W: double;

V0, V1, V2: double;

C0: integer;

AZ: double;

const

P2 = Pi * 2; // 2 * Pi

DR = Pi / 180; // Радиан на градус

K1 = 15 * DR * 1.0027379;

implementation

{$R *.DFM}

function TSun.Sgn(Value: Double): integer;
begin

{if Value = 0 then} Result := 0;

if Value > 0 then Result := 1;

if Value < 0 then Result := -1;

end;

procedure TSun.Calendar;
begin

G := 1;

if Y < 1583 then G := 0;

D1 := Trunc(D);

F := D - D1 - 0.5;

J := -Trunc(7 * (Trunc((M + 9) / 12) + Y) / 4);

if G = 1 then

begin

S := Sgn(M - 9);

A := Abs(M - 9);

J3 := Trunc(Y + S * Trunc(A / 7));

J3 := -Trunc((Trunc(J3 / 100) + 1) * 3 / 4);

end;

J := J + Trunc(275 * M / 9) + D1 + G * J3;

J := J + 1721027 + 2 * G + 367 * Y;

if F >= 0 then Exit;

F := F + 1;

J := J - 1;

end;

procedure TSun.GetTimeZone;
begin

T0 := T / 36525;

S := 24110.5 + 8640184.813 * T0;

S := S + 86636.6 * Z0 + 86400 * L5;

S := S / 86400;

S := S - Trunc(S);

T0 := S * 360 * DR;

end;

procedure TSun.PosOfSun;
begin

// Фундаментальные константы

// (Van Flandern & Pulkkinen, 1979)

L := 0.779072 + 0.00273790931 * T;

G := 0.993126 + 0.0027377785 * T;

L := L - Trunc(L);

G := G - Trunc(G);

L := L * P2;

G := G * P2;

V := 0.39785 * Sin(L);

V := V - 0.01000 * 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.00010 - 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;

procedure TSun.PossibleEvents(Hour: integer);
var num: string;
begin

st := '';

L0 := T0 + Hour * K1;

L2 := L0 + K1;

H0 := L0 - A0;

H2 := L2 - A2;

H1 := (H2 + H0) / 2; // Часовой угол,

D1 := (D2 + D0) / 2; // наклон в получасе

if Hour <= 0 then

V0 := S * Sin(D0) + C * Cos(D0) * Cos(H0) - Z;

V2 := S * Sin(D2) + C * Cos(D2) * Cos(H2) - Z;

if Sgn(V0) = Sgn(V2) then Exit;

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 Exit;

D := Sqrt(D);

if (V0 < 0) and (V2 > 0) then st := st + 'Восход солнца в ';

if (V0 < 0) and (V2 > 0) then M8 := 1;

if (V0 > 0) and (V2 < 0) then st := st + 'Заход солнца в ';

if (V0 > 0) and (V2 < 0) then W8 := 1;

E := (-B + D) / (2 * A);

if (E > 1) or (E < 0) then E := (-B - D) / (2 * A);

T3 := Hour + E + 1 / 120; // Округление

H3 := Trunc(T3);

M3 := Trunc((T3 - H3) * 60);

Str(H3:2, num);

st := st + num + ':';

Str(M3:2, num);

st := st + num;

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;

Str(AZ:4:1, num);

st := st + ', азимут ' + num;

end;

procedure TSun.OutInform;
begin

if (M8 = 0) and (W8 = 0) then

begin

if V2 < 0 then ListBox.Items.Add('Солнце заходит весь день ');

if V2 > 0 then ListBox.Items.Add('Солнце восходит весь день ');

end

else

begin

if M8 = 0 then ListBox.Items.Add('В этот день солнце не восходит ');

if W8 = 0 then ListBox.Items.Add('В этот день солнце не заходит ');

end;

end;

procedure TSun.GetDate;
begin

D := StrToInt(EditD.text);

M := StrToInt(EditM.text);

Y := StrToInt(EditY.text);

end;

procedure TSun.GetInput;
begin

B5 := StrToInt(EditB5.Text);

L5 := -StrToInt(EditL5.Text);

H := StrToInt(EditH.Text);

end;

procedure TSun.ButtonCalcClick(Sender: TObject);
var C0: integer;
begin

GetDate;

GetInput;

ListBox.Items.Add('Широта: ' + EditB5.Text +

' Долгота: ' + EditL5.Text +

' Зона: ' + EditH.Text +

' Дата: ' + EditD.Text +

'/' + EditM.Text +

'/' + EditY.Text);

L5 := L5 / 360;

Z0 := H / 24;

Calendar;

T := (J - 2451545) + F;

TT := T / 36525 + 1; // TT - столетия, начиная с 1900.0

GetTimeZone; // Получение часового пояса

T := T + Z0;

PosOfSun; // Получаем положение солнца

aA[1] := A5;

aD[1] := D5;

T := T + 1;

PosOfSun;

aA[2] := A5;

aD[2] := D5;

if aA[2] < aA[1] then aA[2] := aA[2] + P2;

Z1 := DR * 90.833; // Вычисление зенита

S := Sin(B5 * DR);

C := Cos(B5 * DR);

Z := Cos(Z1);

M8 := 0;

W8 := 0;

A0 := aA[1];

D0 := aD[1];

DA := aA[2] - aA[1];

DD := aD[2] - aD[1];

for C0 := 0 to 23 do

begin

P := (C0 + 1) / 24;

A2 := aA[1] + P * DA;

D2 := aD[1] + P * DD;

PossibleEvents(C0);

if st <> '' then ListBox.Items.Add(st);

A0 := A2;

D0 := D2;

V0 := V2;

end;

OutInform;

ListBox.Items.Add(''); // Разделяем данные

end;

procedure TSun.CreateForm(Sender: TObject);
begin

EditD.Text := FormatDateTime('d', Date);

EditM.Text := FormatDateTime('m', Date);

EditY.Text := FormatDateTime('yyyy', Date);

end;

procedure TSun.ButtonClearClick(Sender: TObject);
begin

ListBox.Clear;

end;
end.

Вот еще одно решение, от Yuri:

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

{$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;

{ } 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;

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.

Категория: Алгоритмы | Добавил: DelphiAiX (28.04.2012)

Просмотров: 2450 | Рейтинг: 0.0/0

Всего комментариев: 0

Добавлять комментарии могут только зарегистрированные пользователи.
[ Регистрация | Вход ]