{-------------------} {$N+}program Sun; constP1 = 3.14159265358;P2 = 2*P1;DR = P1/180;K1 = 15*DR*1.0027379; type TTime = recordHour, Min, Sec: Extended;end;TDate = recordYear, Month, Day: Extended;end; type RequestBlock = recordLatitude,Longitude: Extended;HourZone: Word;Date: TDate;end;ReplyBlock = recordSunRise: TTime;SunRiseAzimuth: Extended;SunSet: TTime;SunSetAzimuth: Extended;end; function Sign(val: Extended): Integer;beginif 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)}varL, G, V, U, W: Extended;beginL := 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;beginwith Rq dobeginL5 := 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) thenbeginS := 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 thenbeginf := 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 dobeginp := (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 thenelse 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) thenbeginV1 := 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 thenbeginD := 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) thenbeginRp.SunRise.Hour := Trunc(H3);Rp.SunRise.Min := Trunc(M3);Rp.SunRiseAzimuth := AZ;M8 := 1;end;if (V0 > 0) and (V2 < 0) thenbeginRp.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) thenbeginif (V2 < 0) then ComputeSunTime := $A3;if (V2 > 0) then ComputeSunTime := $A4;endelsebeginif (M8 = 0) then ComputeSunTime := $A1;if (W8 = 0) then ComputeSunTime := $A2;end;end;end; constSMessage = 'Заход солнца в ';RMessage = 'Восход солнца в ';Message1 = 'В этот день солнце не восходит ';Message2 = 'В этот день солнце не заходит ';Message3 = 'Солнце заходит весь день ';Message4 = 'Солнце восходит весь день '; varRq: RequestBlock;Rp: ReplyBlock; begin with Rq dobeginWrite(' Широта........:'); 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);elseWrite(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. |