diff --git a/kstars/tools/modcalcvizequinox.h b/kstars/tools/modcalcvizequinox.h --- a/kstars/tools/modcalcvizequinox.h +++ b/kstars/tools/modcalcvizequinox.h @@ -48,8 +48,7 @@ private: void processLines(QTextStream &istream); void addDateAxes(); - KStarsDateTime findEquinox(int year, bool Spring, KPlotObject *po); - KStarsDateTime findSolstice(int year, bool Summer); + void findSolsticeAndEquinox(uint32_t year); public: KStarsDateTime dSpring, dSummer, dAutumn, dWinter; diff --git a/kstars/tools/modcalcvizequinox.cpp b/kstars/tools/modcalcvizequinox.cpp --- a/kstars/tools/modcalcvizequinox.cpp +++ b/kstars/tools/modcalcvizequinox.cpp @@ -203,10 +203,8 @@ } Plot->addPlotObject(ecl); - dSpring = findEquinox(Year->value(), true, ecl); - dSummer = findSolstice(Year->value(), true); - dAutumn = findEquinox(Year->value(), false, ecl); - dWinter = findSolstice(Year->value(), false); + // Calculate dates for all solstices and equinoxes + findSolsticeAndEquinox(Year->value()); //Display the Date/Time of each event in the text fields VEquinox->setText(QLocale().toString(dSpring, QLocale::LongFormat)); @@ -264,117 +262,110 @@ } } -KStarsDateTime modCalcEquinox::findEquinox(int year, bool Spring, KPlotObject *ecl) +// Calculate and store dates for all solstices and equinoxes +void modCalcEquinox::findSolsticeAndEquinox(uint32_t year) { - // Interpolate to find the moment when the Sun crosses the equator - // Set initial guess in February or August to be sure that this - // point is before equinox. - const int month = Spring ? 2 : 8; - int i = QDate(year, month, 1).dayOfYear(); - double dec1, dec2; - dec2 = ecl->points().at(i)->y(); - do + double m, jdme[4]; + if(year >= 1000) { - ++i; - dec1 = dec2; - dec2 = ecl->points().at(i)->y(); + m = (year-2000.0) / 1000.0; + jdme[0] = (-0.00057 * pow(m, 4)) + (-0.00411 * pow(m, 3)) + (0.05169 * pow(m, 2)) + (365242.37404 * m) + 2451623.80984; + jdme[1] = (-0.0003 * pow(m, 4)) + (0.00888 * pow(m, 3)) + (0.00325 * pow(m, 2)) + (365241.62603 * m) + 2451716.56767; + jdme[2] = (0.00078 * pow(m, 4)) + (0.00337 * pow(m, 3)) + (-0.11575 * pow(m, 2)) + (365242.01767 * m) + 2451810.21715; + jdme[3] = (0.00032 * pow(m, 4)) + (-0.00823 * pow(m, 3)) + (-0.06223 * pow(m, 2)) + (365242.74049 * m) + 2451900.05952; + } + else + { + m = year / 1000.0; + jdme[0] = (-0.00071 * pow(m, 4)) + (0.00111 * pow(m, 3)) + (0.06134 * pow(m, 2)) + (365242.13740 * m) + 1721139.29189; + jdme[1] = (0.00025 * pow(m, 4)) + (0.00907 * pow(m, 3)) + (-0.05323 * pow(m, 2)) + (365241.72562 * m) + 1721233.25401; + jdme[2] = (0.00074 * pow(m, 4)) + (-0.00297 * pow(m, 3)) + (-0.11677 * pow(m, 2)) + (365242.49558 * m) + 1721325.70455; + jdme[3] = (-0.00006 * pow(m, 4)) + (-0.00933 * pow(m, 3)) + (-0.00769 * pow(m, 2)) + (365242.88257 * m) + 1721414.39987; + } + + double t[4]; + for(int i = 0; i < 4; i++) + { + t[i] = (jdme[i] - 2451545)/36525; } - while (dec1 * dec2 > 0.0); //when dec1*dec2<0.0, we bracket the zero - double x1 = ecl->points().at(i - 1)->x(); - double x2 = ecl->points().at(i)->x(); - double d = fabs(dec2 - dec1); - double f = 1.0 - fabs(dec2) / d; //fractional distance of the zero, from point1 to point2 + double a[4][24] = {{485, 203, 199, 182, 156, 136, 77, 74, 70, 58, 52, 50, 45, 44, 29, 18, 17, 16, 14, 12, 12, 12, 9, 8}, {485, 203, 199, 182, 156, 136, 77, 74, 70, 58, 52, 50, 45, 44, 29, 18, 17, 16, 14, 12, 12, 12, 9, 8}, {485, 203, 199, 182, 156, 136, 77, 74, 70, 58, 52, 50, 45, 44, 29, 18, 17, 16, 14, 12, 12, 12, 9, 8}, {485, 203, 199, 182, 156, 136, 77, 74, 70, 58, 52, 50, 45, 44, 29, 18, 17, 16, 14, 12, 12, 12, 9, 8}}; - KStarsDateTime dt0(QDate(year, 1, 1), QTime(0, 0, 0)); - KStarsDateTime dt = dt0.addSecs(86400.0 * (x1 - 1 + f * (x2 - x1))); - return dt; -} + double b[4][24] = {{324.96, 337.23, 342.08, 27.85, 73.14, 171.52, 222.54, 296.72, 243.58, 119.81, 297.17, 21.02, 247.54, 325.15, 60.93, 155.12, 288.79, 198.04, 199.76, 95.39, 287.11, 320.81, 227.73, 15.45}, {324.96, 337.23, 342.08, 27.85, 73.14, 171.52, 222.54, 296.72, 243.58, 119.81, 297.17, 21.02, 247.54, 325.15, 60.93, 155.12, 288.79, 198.04, 199.76, 95.39, 287.11, 320.81, 227.73, 15.45}, {324.96, 337.23, 342.08, 27.85, 73.14, 171.52, 222.54, 296.72, 243.58, 119.81, 297.17, 21.02, 247.54, 325.15, 60.93, 155.12, 288.79, 198.04, 199.76, 95.39, 287.11, 320.81, 227.73, 15.45}, {324.96, 337.23, 342.08, 27.85, 73.14, 171.52, 222.54, 296.72, 243.58, 119.81, 297.17, 21.02, 247.54, 325.15, 60.93, 155.12, 288.79, 198.04, 199.76, 95.39, 287.11, 320.81, 227.73, 15.45}}; -KStarsDateTime modCalcEquinox::findSolstice(int year, bool Summer) -{ - //Find the moment when the Sun reaches maximum declination - //First find three points which bracket the maximum (i.e., x2 > x1,x3) - //Start at June 16th, which will always be approaching the solstice - - long double jd1, jd2, jd3, jd4; - double y2(0.0), y3(0.0), y4(0.0); - int month = 6; - if (!Summer) - month = 12; - - jd3 = KStarsDateTime(QDate(year, month, 16), QTime(0, 0, 0)).djd(); - KSNumbers num(jd3); - KSSun Sun; - Sun.findPosition(&num); - y3 = Sun.dec().Degrees(); + double c[] = {1934.136, 32964.467, 20.186, 445267.112, 45036.886, 22518.443, 65928.934, 3034.906, 9037.513, 33718.147, 150.678, 2281.226, 29929.562, 31555.956, 4443.417, 67555.328, 4562.452, 62894.029, 31436.921, 14577.848, 31931.756, 34777.259, 1222.114, 16859.074}; - int sgn = 1; - if (!Summer) - sgn = -1; //find minimum if the winter solstice is sought + double d[4][24]; - do + for (int i = 0; i < 4; i++) { - jd3 += 1.0; - num.updateValues(jd3); - Sun.findPosition(&num); - y2 = y3; - Sun.findPosition(&num); - y3 = Sun.dec().Degrees(); + for (int j = 0; j < 24; j++) + { + d[i][j] = c[j] * t[i]; + } } - while (y3 * sgn > y2 * sgn); - //Ok, now y2 is larger(smaller) than both y3 and y1. - // Never read back - // jd2 = jd3 - 1.0; - jd1 = jd3 - 2.0; - - //Choose a new starting jd2 that follows the golden ratio: - // a/b = 1.618; a+b = 2...a = 0.76394 - jd2 = jd1 + 0.76394; - num.updateValues(jd2); - Sun.findPosition(&num); - y2 = Sun.dec().Degrees(); + for (int i = 0; i < 4; i++) + { + for (int j = 0; j < 24; j++) + { + d[i][j] = cos(qDegreesToRadians(b[i][j] + d[i][j])); + } + } - while (jd3 - jd1 > 0.0005) //sub-minute precision + for (int i = 0; i < 4; i++) { - jd4 = jd1 + jd3 - jd2; + for (int j = 0; j < 24; j++) + { + d[i][j] = d[i][j] * a[i][j]; + } + } - num.updateValues(jd4); - Sun.findPosition(&num); - y4 = Sun.dec().Degrees(); + double e[4], f[4], g[4], jd[4]; - if (y4 * sgn > y2 * sgn) //make jd4 the new center + for (int i = 0; i < 4; i++) + { + e[i] = 0; + for (int j = 0; j < 24; j++) { - if (jd4 > jd2) - { - jd1 = jd2; - jd2 = jd4; - y2 = y4; - } - else - { - jd3 = jd2; - // Never read back - // y3 = y2; - jd2 = jd4; - y2 = y4; - } + e[i] += d[i][j]; } - else //make jd4 a new endpoint + } + + for (int i = 0; i < 4; i++) + { + f[i] = (35999.373*t[i]-2.7); + } + + for (int i = 0; i < 4; i++) + { + g[i] = 1 + 0.0334 * cos(qDegreesToRadians(f[i])) + 0.00007 * cos(qDegreesToRadians(2*f[i])); + } + + for (int i = 0; i < 4; i++) + { + jd[i] = jdme[i] + 0.00001*(e[i]/g[i]); + } + + KStarsDateTime dt(QDate(year, 1, 1), QTime(0, 0, 0)); + + for (int i = 0; i < 4; i++) + { + dt.setDJD(jd[i]); + switch(i) { - if (jd4 > jd2) - { - jd3 = jd4; - // Never read back - // y3 = y4; - } - else - { - jd1 = jd4; - } + case 0 : + dSpring = dt; + break; + case 1 : + dSummer = dt; + break; + case 2 : + dAutumn = dt; + break; + case 3 : + dWinter = dt; + break; } } - - return KStarsDateTime(jd2); -} +} \ No newline at end of file