ВНИМАНИЕ! На форуме началось голосование в конкурсе - астрофотография месяца АПРЕЛЬ!
0 Пользователей и 1 Гость просматривают эту тему.
равно -1.47127914 радиан.
азимут (A) и высоту (h) любого светила, который занесен в её родной каталог, в любое время и в любом месте на Земле может указать Winstars http://winstars.net/ Есть еще в первой версии. Два клика правой кнопкой мышки - и там указан не только азимут...
По какому принципу вы будете добавлять 180 градусов?
так как круговые функции не однозначны
анализировать знаки синуса и косинуса...
program sun;{$APPTYPE CONSOLE}uses Windows;type Time = type Double; TimeB = record Day : Byte; Mon : Byte; Yr : WORD; Hr : Byte; Min : Byte; Sec : Byte; MS : WORD; end;type PDayTable = ^TDayTable; TDayTable = array[1..12] of Word; const MonthDays: array [Boolean] of TDayTable = ((31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31), (31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31)); HoursPerDay = 24; MSecsPerSec = 1000; MinsPerDay = HoursPerDay * 60; SecsPerDay = MinsPerDay * 60; SecPerHour = 3600; MSecsPerDay = SecsPerDay * MSecsPerSec; MSecsPerMin = MSecsPerSec * 60; MSecsPerHour= MSecsPerMin * 60; JulMonth = 30.6001; JulYear = 365.25; GST_A = 0.0657098; GST_C = 1.002738; GST_D = 0.997270; e_g = 278.833540; w_g = 282.596403; e_ = 0.016718; r_0 = 1.495985*100000000; t_0 = 0.533128; rad = Pi/180; int360 = 360;//тригонометрические функции взяты из модуля Mathfunction ArcTan2(const Y, X: Extended): Extended;asm FLD Y FLD X FPATAN FWAITend;function Tan(const X: Extended): Extended;{ Tan := Sin(X) / Cos(X) }asm FLD X FPTAN FSTP ST(0) { FPTAN pushes 1.0 after result } FWAITend;function ArcCos(const X: Extended): Extended;begin Result := ArcTan2(Sqrt(1 - X * X), X);end;function ArcSin(const X: Extended): Extended;begin Result := ArcTan2(X, Sqrt(1 - X * X))end;function DegToRad(const Degrees: Double): Double;begin Result := Degrees * rad;end;function RadToDeg(const Radians: Double): Double;begin Result := Radians / rad;end;function IsLeapYear(const Year: Word): Boolean;beginResult := (Year mod 4 = 0) and ((Year mod 100 <> 0) or (Year mod 400 = 0));end;//Взято из стандартного модуля SysUtilsfunction DateToInt(T:TimeB): Time;var // $3 I: Integer; DayTable: PDayTable;begin DayTable := @MonthDays[IsLeapYear(T.Yr)]; if (T.Yr >= 1) and (T.Yr <= 9999) and (T.Mon >= 1) and (T.Mon <= 12) and (T.Day >= 1) and (T.Day <= DayTable^[T.Mon]) then begin Result := T.Day; for I := 1 to T.Mon - 1 do Result := Result+DayTable^[I]; end else Result := -1;end;function TimeToDec(const T:TimeB): Time;begin // $7 Result := (((((T.MS/1000)+T.Sec)/60)+T.Min)/60)+T.Hr; {Так как функция применяется не только для времени а и для перевода градусов то выключил за ненадобностью if T.Hr < 24... if (T.Hr < 24) and (T.Min < 60) and (T.Sec < 60) and (T.MS < 1000) then Result := (((((T.MS/1000)+T.Sec)/60)+T.Min)/60)+T.Hr else Result := -1;}end;function DecToTime(const T:Time): TimeB;var // $8 m,s:Real;begin with Result do begin Hr := Trunc(T); m := (T-Hr)*60; Min := Trunc(m); s := (m-Min)*60; Sec := Trunc(s); MS := Trunc((s-Sec)*1000); end;end;function TimeToStr(Time:Time):String;function ziro(S:String):String;begin if Length(S)<2then Result:='0'+S else Result := S;end;var T:TimeB; S:String;begin T := DecToTime(Time); Str(T.Hr,S); Result := ziro(S)+':'; Str(T.Min,S); Result := Result+ziro(S)+':'; Str(T.Sec,S); Result := Result+ziro(S)+'.'; Str(T.MS,S); Result := Result+ziro(S);end;function Jul(const Day,Month:Byte;const Year:Word): Time;begin Result := (1461 * (Year + 4800 + (Month - 14) div 12)) div 4 + (367 * (Month - 2 - 12 * ((Month - 14) div 12))) div 12 - (3 * ((Year + 4900 + (Month - 14) div 12) div 100)) div 4 + Day - 32075.5;//Взято из стандартного модуля DateUtilsend;function BackTo360(const Val:Time):Time;begin Result := Val; while Result >= int360 do Result := Result-int360; while Result < 0 do Result := Result+int360;end;function BackTo24(const Val:Time):Time;begin Result := Val; while Result >= HoursPerDay do Result := Result-HoursPerDay; while Result < 0 do Result := Result+HoursPerDay;end;function GST2GMT(const T:TimeB;const TimeInDec:Time):Time;function GetGST_B(const Year:WORD):Double;var T,R:Double;begin T := (Jul(0,1,Year) - 2415020)/36525; R := 6.6460656+T*(2400.051262+0.00002581*T); Result := 24-(R-24*(Year-1900));end;begin // $13 Гринвическое время! Result := BackTo24(DateToInt(T)*GST_A-GetGST_B(T.Yr)); Result := BackTo24(TimeInDec-Result); Result := Result*GST_D;end;function LST2GST(const Angle:Time;const East:Boolean;T:Time):Time;var // $15 R : Time;begin R := Angle / 15; if East then Result := BackTo24(T-R) else Result := BackTo24(T+R);end;procedure EclipToEquCoords(const YearNow:WORD; Lam,Ba:Time; var Alif,Dal:Time);function Get_T(const Year:WORD):Time; begin Result := (Jul(Year,1,0)-2415020)/36525;end;function Get_E(const Year:WORD):Time;var T : Time; const c : array [1..4] of Double = (23.452294,46.845,0.0059,0.001811); begin T := Get_T(Year);Result := c[1]-(T*(c[2]-T*(c[3]+c[4]*T))/3600);end;var // $27 sinD,x,y,E:Double;begin Lam := DegToRad(Lam); Ba := DegToRad(Ba); E := DegToRad(23.441884);//Get_E(YearNow); sinD:= Sin(Ba) * Cos(E) + Cos(Ba) * Sin(E) * Sin(Lam); Dal := RadToDeg(ArcSin(SinD)); y := Sin(Lam)* Cos(E) - Tan(Ba) * Sin(E); x := Cos(Lam); alif:= RadToDeg(ArcTan2{SmartATan}(Y,X))/15;end;procedure VoshodZahod__(Alif,Dal, NS:Time; var VosAz,ZahAz,Voshod,Zahod:Time);var H:Time;begin VosAz := RadToDeg(arccos(sin(dal)/cos(NS))); ZahAz := 360-VosAz; H := RadToDeg(arccos(-tan(NS)*tan(Dal)))/15; Voshod:= BackTo24(HoursPerDay+Alif-H); Zahod := BackTo24(Alif+H);end;procedure VoshodZahod(const T:TimeB; Alif,Dal, NS,WE:Time; var VosAz,ZahAz,Voshod,Zahod:Time);begin// $32 Dal := DegToRad(Dal); NS := DegToRad(NS); VoshodZahod__(Alif,Dal,NS,VosAz,ZahAz,Voshod,Zahod); if WE >= 0 then begin Voshod:= LST2GST(WE,True,Voshod); Zahod := LST2GST(WE,True,Zahod); end else begin Voshod:= LST2GST(abs(WE),False,Voshod); Zahod := LST2GST(abs(WE),False,Zahod); end; Voshod:= GST2GMT(T,Voshod); Zahod := GST2GMT(T,Zahod);end;procedure VoshodZahodPlus(const T:TimeB; Ref,Alif,Dal, NS,WE:Time; var VosAz,ZahAz,Voshod,Zahod:Time);var // $32 psi,DeltaA, A_r,A_s, y,DeltaT, LSTr,LSTs : Time;begin //Alif := DegToRad(Alif); Dal := DegToRad(Dal); Ref := DegToRad(Ref); NS := DegToRad(NS); VoshodZahod__(Alif,Dal,NS,A_r,A_s,LSTr,LSTs); psi := arccos(sin(NS)/cos(dal)); DeltaA:= arcsin(tan(Ref)/tan(psi)); VosAz := A_r-DeltaA; ZahAz := A_s+DeltaA; y := arcsin(sin(Ref)/sin(psi)); DeltaT:= RadToDeg(240*y/cos(dal))/SecPerHour; LSTr := LSTr - DeltaT; LSTs := LSTs + DeltaT; if WE >= 0 then begin Voshod:= LST2GST(WE,True,LSTr); Zahod := LST2GST(WE,True,LSTs); end else begin Voshod:= LST2GST(abs(WE),False,LSTr); Zahod := LST2GST(abs(WE),False,LSTs); end; Voshod:= GST2GMT(T,Voshod); Zahod := GST2GMT(T,Zahod);end;procedure SunCentr(T:TimeB;var Lam,Alif,Dal:Time);function R2(const M:Time):Time;const eps = 0.0000001;var E,D,DE:Extended;begin E := M; D := E-e_*sin(E)-M; while abs(D) > eps do begin DE := D/(1-(e_*cos(E))); E := E-DE; D := E-e_*sin(E)-M; end; Result := E;end;var D,E,M,N,Tg05N,Atn05N:Time;begin D := Jul(T.Day,T.Mon,T.Yr)-Jul(0,1,1980); N := BackTo360(((360/365.2422)*(D))); M := BackTo360(N + e_g - w_g); E := R2(DegToRad(M)); tg05N := Sqrt((1+e_)/(1-e_))*Tan(E/2); Atn05N := arctan(tg05N); N := RadToDeg(2*Atn05N); Lam := BackTo360(N+w_g); EclipToEquCoords(T.Yr,Lam,0,Alif,Dal);end;procedure SunSet( const T:TimeB; const Visota, NS,WE,TimeZone, Dim,Par,Ref : Time; var Voshod, Zahod:Time);var Lam, Alif,Dal, VosAz,ZahAz, STr,STs : Array [1..2] of Time; Tr,Ts, DalU,psi,x,y,DeltaT: Time;begin SunCentr(T,Lam[1],Alif[1],Dal[1]); //косяк 1 //результаты функции расходятся с книгой на сотые доли //поэтому я взял данные полученные авторами из книги Lam[2] := Lam[1]+0.985647; EclipToEquCoords(T.Yr,Lam[2],0,Alif[2],Dal[2]); VoshodZahod(T,Alif[1],Dal[1],NS,WE,VosAz[1],ZahAz[1],STr[1],STs[1]); VoshodZahod(T,Alif[2],Dal[2],NS,WE,VosAz[2],ZahAz[2],STr[2],STs[2]); //косяк 2 //результаты функции по черному расходятся с книгой Tr := (24.07*STr[1]) / (24.07+STr[1]-STr[2]); TS := (24.07*STs[1]) / (24.07+STs[1]-STs[2]); DalU := DegToRad((Dal[1]+Dal[2])/2); psi := arccos(sin(DegToRad(NS))/cos(DalU)); x := (Dim/2)+Par+Ref; y := RadToDeg(arcsin(sin(DegToRad(x))/sin(psi))); DeltaT:= (240*y/cos(DalU))/SecPerHour; Voshod:= Tr-DeltaT+TimeZone; Zahod := Ts+DeltaT+TimeZone; WriteLn(DeltaT:6:6,#10, Tr:6:6,' ',Ts:6:6,#10, TimeToStr(Voshod),' ',TimeToStr(Zahod));end;var T:TimeB; Visota, NS,WE,V,Z, Dim,Par,Ref:Time;begin T.Day := 26; T.Mon := 7; T.Yr := 2009; T.Hr:=51;T.Min:=30;T.Sec:=00; NS := 44.95;//TimeToDec(T);//широта WE := 34.1;//долгота T.Hr:=00;T.Min:=34;T.Sec:=00; Ref := TimeToDec(T);//Рефракция T.Hr:=00;T.Min:=00;T.Sec:=8;T.MS:=790; Par := TimeToDec(T);//Паралакс T.Hr:=00;T.Min:=16;T.Sec:=0;T.MS:=0; Dim := 0.533;//TimeToDec(T);//диаметр солнца Visota := 0; //V [Z] - время восхода [захода] SunSet(T,Visota,NS,WE,3,Dim,Par,Ref,V,Z);//смотрите процедуру ReadLn;end.
//косяк 1//результаты функции расходятся с книгой на сотые доли//поэтому я взял данные полученные авторами из книги
//результаты функции по черному расходятся с книгой
Формулы есть в Астрономическом календаре, постоянная часть, 1981 год
IMHO, не тратьте время на этот код. Требуемую точность вы не получите. К тому же, вы сами-то заметили, что подменили одну задачу другой? Вам вроде бы нужно было высоту Солнца вычислить в зависимости от времени, а этот код в результате дает время восхода/захода Солнца и азимутальные координаты восхода/захода.