// Astronomicke vypocty pro urceni aktualni polohy Slunce a Mesice //================================================================== //------------------------------------------------------------- // posun z lokalniho casu do UTC (pro SEC je to 1, pro SELC je to 2) // vstupem jsou globalni promenne "LOC_xxx" // vysledky jsou v globalnich promennych "astro_UTC_xxx" void z_LOC_na_Astro_UTC(int offset) { // pocet dni v mesici pro ruzne pocitani dni od zacatku roku (unor se pak pri prestupnem roce jeste upravuje) // x , Led, Uno, Bre, Dub, Kve, Cer, Cec, Srp, Zar, Rij, Lis, Pro int pomdny_d[] = { 0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 }; if ((LOC_rok % 4) == 0) pomdny_d[2] = 29; // pri prestupnem roku se nastavuje unor na 29 dni if (LOC_rok == 2100) pomdny_d[2] = 28; // rok 2100 bude specialita. I kdyz je deliteny 4, tak NEBUDE prestupny. astro_UTC_min = LOC_min; // na zacatku se UTC promenne nastavi na stejnou hodnotu jako mistni cas astro_UTC_rok = LOC_rok; astro_UTC_mes = LOC_mes; astro_UTC_den = LOC_den; astro_UTC_hod = LOC_hod - offset; // ... akorat od hodin se odecte casovy posun (1 hodina pro SEC, 2 hodiny pro SELC ...) if (astro_UTC_hod < 0) // kontrola, jestli po odecteni casove zony, nedoslo k podteceni casu do predchozicho dne { astro_UTC_hod = astro_UTC_hod + 24; astro_UTC_den = astro_UTC_den - 1; if (astro_UTC_den < 1) // kontrola podteceni do predchoziho mesice { astro_UTC_mes = astro_UTC_mes - 1; if (astro_UTC_mes < 1) // kontrola podteceni do predchoziho roku { astro_UTC_rok = astro_UTC_rok - 1; astro_UTC_mes = 12; } astro_UTC_den = pomdny_d[astro_UTC_mes]; // cislo dne se nastavi na nejvyssi cislo dne v predchozim mesici } } } //--------------------------------------------------------- //--------------------------------------------------------- // vypocty pro polohy Mesice a Slunce (opsano nekde z internetu) void astro_vypocty(void) { if (leto == true) casova_zona = letni_posun; else casova_zona = zimni_posun; z_LOC_na_Astro_UTC(casova_zona); // z globalnich promennych pro casove udaje v mistni casove zone (LOC_xxx) vypocte UTC datum a cas double days = day2000(astro_UTC_rok, astro_UTC_mes, astro_UTC_den, astro_UTC_hod + astro_UTC_min/60.0); double t = days / 36525.0; // Data Slunce double L1 = rozsah(280.466 + 36000.8 * t); double M1 = rozsah(357.529+35999*t - 0.0001536* t*t + t*t*t/24490000.0); double C1 = (1.915 - 0.004817* t - 0.000014* t * t)* dsin(M1); C1 = C1 + (0.01999 - 0.000101 * t)* dsin(2*M1); C1 = C1 + 0.00029 * dsin(3*M1); double V1 = M1 + C1; double EC1 = 0.01671 - 0.00004204 * t - 0.0000001236 * t*t; double R1 = 0.99972 / (1 + EC1 * dcos(V1)); double Th1 = L1 + C1; double OM1 = rozsah(125.04 - 1934.1 * t); double LaM1 = Th1 - 0.00569 - 0.00478 * dsin(OM1); double Obl = (84381.448 - 46.815 * t)/3600; double Ra1 = datan2(dsin(Th1) * dcos(Obl) - dtan(0)* dsin(Obl), dcos(Th1)); double DEC1 = dasin(dsin(0)* dcos(Obl) + dcos(0)*dsin(Obl)*dsin(Th1)); // Data Mesic double F = rozsah(93.2721 + 483202 * t - 0.003403 * t* t - t * t * t/3526000.0); double L2 = rozsah(218.316 + 481268 * t); double Om2 = rozsah(125.045 - 1934.14 * t + 0.002071 * t * t + t * t * t/450000.0); double M2 = rozsah(134.963 + 477199 * t + 0.008997 * t * t + t * t * t/69700.0); double D = rozsah(297.85 + 445267 * t - 0.00163 * t * t + t * t * t/545900.0); double D2 = 2*D; double R2 = 1 + (-20954 * dcos(M2) - 3699 * dcos(D2 - M2) - 2956 * dcos(D2)) / 385000.0; double R3 = (R2 / R1) / 379.168831168831; double Bm = 5.128 * dsin(F) + 0.2806 * dsin(M2 + F); Bm = Bm + 0.2777 * dsin(M2 - F) + 0.1732 * dsin(D2 - F); double Lm = 6.289 * dsin(M2) + 1.274 * dsin(D2 -M2) + 0.6583 * dsin(D2); Lm = Lm + 0.2136 * dsin(2*M2) - 0.1851 * dsin(M1) - 0.1143 * dsin(2 * F); Lm = Lm +0.0588 * dsin(D2 - 2*M2) ; Lm = Lm + 0.0572* dsin(D2 - M1 - M2) + 0.0533* dsin(D2 + M2); Lm = Lm + L2; double Ra2 = (datan2(dsin(Lm) * dcos(Obl) - dtan(Bm)* dsin(Obl), dcos(Lm))) / 15.0; double Dec2 = dasin(dsin(Bm)* dcos(Obl) + dcos(Bm)*dsin(Obl)*dsin(Lm)); double HLm = rozsah(LaM1 + 180 + (180/PI) * R3 * dcos(Bm) * dsin(LaM1 - Lm)); double HBm = R3 * Bm; double I = 1.54242; double W = Lm - Om2; double qY = dcos(W) * dcos(Bm); double qX = dsin(W) * dcos(Bm) * dcos(I) - dsin(Bm) * dsin(I); double A = datan2(qX, qY); // double EL = A - F; // double EB = dasin(-dsin(W) * dcos(Bm) * dsin(I) - dsin(Bm) * dcos(I)); W = rozsah(HLm - Om2); qY = dcos(W) * dcos(HBm); qX = dsin(W) * dcos(HBm) * dcos(I) - dsin(HBm) * dsin(I); A = datan2(qX, qY); // double SL = rozsah(A - F); // double SB = dasin(-dsin(W) * dcos(HBm) * dsin(I) - dsin(HBm) * dcos(I)); // double Co = 0; // if (SL < 90) Co = 90 - SL; // else Co = 450 - SL; A = dcos(Bm) * dcos(Lm - LaM1); double Psi = 90 - datan(A / sqrt(1- A*A)); qX = R1 * dsin(Psi); qY = R3 - R1* A; double Il = datan2(qX, qY); double K = (1 + dcos(Il))/2; RaToAzim ( 1 , Ra1 / 15.0, DEC1 , LOC_rok , LOC_mes , LOC_den , LOC_hod , LOC_min , casova_zona , GeoLon , GeoLat ); // Slunce RaToAzim ( 2 , Ra2 , Dec2 , LOC_rok , LOC_mes , LOC_den , LOC_hod , LOC_min , casova_zona , GeoLon , GeoLat ); // Mesic Mes_osvit = zaokrouhlit(K*100); Mes_elevace = zaokrouhlit(Mes_elevace) + 90; // do pameti se ukladaji jen kladna cisla, proto se elevace posouvaji o 90 stupnu vyse Slu_elevace = zaokrouhlit(Slu_elevace) + 90; // do pameti se ukladaji jen kladna cisla, proto se elevace posouvaji o 90 stupnu vyse Mes_azimut = zaokrouhlit(Mes_azimut); Slu_azimut = zaokrouhlit(Slu_azimut); } //--------------------------------------------------------- //------------------------------------------------------------- int zaokrouhlit(double n) { if (n >= 0) n = n + 0.5; // kladna cisla se zaokrouhluji prictenim 0.5 a oseknutim desetinne casti else n = n - 0.5; // zaporna cisla se zaokrouhluji odectenim 0.5 a oseknutim desetinne casti return (int)n; } //------------------------------------------------------------- //--------------------------------------------------------- // Ruzna aritmetika - opsano nekde z internetu (prevedeno z BASICu do PHP a pak sem do Wiringu) double rozsah(double x) { double b = x / 360.0; double a = 360 * (b - ipart(b)); if (a < 0 ) { a = a + 360; } return a; } double ipart(double x) { double a; if (x > 0) { a = floor(x); } else { a = ceil(x); } return a; } double dsin(double x) { return sin(PI / 180.0 * x); } double dcos(double x) { return cos(PI / 180.0 * x); } double dtan(double x) { return tan(PI / 180.0 * x); } // pocet dnu od J2000 double day2000(int y, int m, int d, double h) { // unsigned long greg = y*10000 + m*100 + d; if (m == 1 || m == 2) { y = y - 1; m = m + 12; } double a = floor(y/100); double b = 2 - a + floor(a/4); int c = floor(365.25 * y); int d1 = floor(30.6001 * (m + 1)); return ( b + c + d1 -730550.5 + d + h/24.0); } double datan2(double y, double x) { double a; if ((x == 0) && (y == 0)) { return 0; } else { a = datan(y / x); if (x < 0) { a = a + 180; } if (y < 0 && x > 0) { a = a + 360; } return a; } } double dasin(double x) { return 180/ PI * asin(x); } double dacos(double x) { return 180/ PI * acos(x); } double datan(double x) { return 180/ PI * atan(x); } //------------------------------------------------------------- //------------------------------------------------------------- // prevede Rektascenzi a Deklinaci na Azimut a Elevaci (podle lokalniho casu, mistni casove zony a zemepisnych souradnic) void RaToAzim( byte teleso , double qRA, double qDECL, int qrok, byte qmes, byte qden, byte qhod, byte qmin , int qutcdif, double qlong, double qlat) { double qcas = qhod + (qmin / 60.0); double dj2000 = floor(day2000(qrok, qmes, qden, qhod +12)) + ((qcas - qutcdif) / 24.0) - 0.5; double lst = 100.46 + 0.985647 * dj2000 + qlong + 15.0 * (qcas - qutcdif); while (lst <= 0) lst = lst + 360; while (lst >= 360) lst = lst - 360; double ha = lst - qRA * 15.0; while (ha <= 0) ha = ha + 360; while (ha >= 360) ha = ha - 360; double alt = dasin(dsin(qDECL) * dsin(qlat) + dcos(qDECL) * dcos(qlat) * dcos(ha)); double azm = datan2(dsin(ha), dcos(ha) * dsin(qlat) - dtan(qDECL) * dcos(qlat)) + 180; if (azm >= 360) azm = azm - 360; if (teleso == 1) // teleso je Slunce { Slu_elevace = alt; Slu_azimut = azm; } if (teleso == 2) // teleso je Mesic { Mes_elevace = alt; Mes_azimut = azm; } if (teleso == 3) // teleso neni zadane (pri prepoctu Ra-Dec na Az-El z menu) { obecna_elevace = alt; obecny_azimut = azm; } } //------------------------------------------------------------- //------------------------------------------------------------- // ze dvou azimutu a elevaci vypocte uhlovou vzdalenost ve stupnich float uhlova_vzdalenost(int azimut1, int elevace1, int azimut2, int elevace2) { float u_vzdal; u_vzdal = dacos((dcos(elevace1) * dcos(elevace2) * dcos(azimut1 - azimut2)) + (dsin(elevace1) * dsin(elevace2))); return u_vzdal; }