Changeset View
Changeset View
Standalone View
Standalone View
kstars/tools/modcalcvizequinox.cpp
Show First 20 Lines • Show All 197 Lines • ▼ Show 20 Line(s) | 197 | { | |||
---|---|---|---|---|---|
198 | KSNumbers num(dt.djd()); | 198 | KSNumbers num(dt.djd()); | ||
199 | Sun.findPosition(&num); | 199 | Sun.findPosition(&num); | ||
200 | ecl->addPoint(double(i), Sun.dec().Degrees()); | 200 | ecl->addPoint(double(i), Sun.dec().Degrees()); | ||
201 | 201 | | |||
202 | dt = dt.addDays(1); | 202 | dt = dt.addDays(1); | ||
203 | } | 203 | } | ||
204 | Plot->addPlotObject(ecl); | 204 | Plot->addPlotObject(ecl); | ||
205 | 205 | | |||
206 | dSpring = findEquinox(Year->value(), true, ecl); | 206 | // Calculate dates for all solstices and equinoxes | ||
207 | dSummer = findSolstice(Year->value(), true); | 207 | findSolsticeAndEquinox(Year->value()); | ||
208 | dAutumn = findEquinox(Year->value(), false, ecl); | | |||
209 | dWinter = findSolstice(Year->value(), false); | | |||
210 | 208 | | |||
211 | //Display the Date/Time of each event in the text fields | 209 | //Display the Date/Time of each event in the text fields | ||
212 | VEquinox->setText(QLocale().toString(dSpring, QLocale::LongFormat)); | 210 | VEquinox->setText(QLocale().toString(dSpring, QLocale::LongFormat)); | ||
213 | SSolstice->setText(QLocale().toString(dSummer, QLocale::LongFormat)); | 211 | SSolstice->setText(QLocale().toString(dSummer, QLocale::LongFormat)); | ||
214 | AEquinox->setText(QLocale().toString(dAutumn, QLocale::LongFormat)); | 212 | AEquinox->setText(QLocale().toString(dAutumn, QLocale::LongFormat)); | ||
215 | WSolstice->setText(QLocale().toString(dWinter, QLocale::LongFormat)); | 213 | WSolstice->setText(QLocale().toString(dWinter, QLocale::LongFormat)); | ||
216 | 214 | | |||
217 | //Add vertical dotted lines at times of the equinoxes and solstices | 215 | //Add vertical dotted lines at times of the equinoxes and solstices | ||
▲ Show 20 Lines • Show All 41 Lines • ▼ Show 20 Line(s) | 253 | { | |||
259 | Plot->addPlotObject(poMonth); | 257 | Plot->addPlotObject(poMonth); | ||
260 | poMonth = new KPlotObject(Qt::white, KPlotObject::Lines, 1); | 258 | poMonth = new KPlotObject(Qt::white, KPlotObject::Lines, 1); | ||
261 | poMonth->addPoint(dmonth(imonth), Plot->dataRect().bottom()); | 259 | poMonth->addPoint(dmonth(imonth), Plot->dataRect().bottom()); | ||
262 | poMonth->addPoint(dmonth(imonth), Plot->dataRect().bottom() - 1.4); | 260 | poMonth->addPoint(dmonth(imonth), Plot->dataRect().bottom() - 1.4); | ||
263 | Plot->addPlotObject(poMonth); | 261 | Plot->addPlotObject(poMonth); | ||
264 | } | 262 | } | ||
265 | } | 263 | } | ||
266 | 264 | | |||
267 | KStarsDateTime modCalcEquinox::findEquinox(int year, bool Spring, KPlotObject *ecl) | 265 | // Calculate and store dates for all solstices and equinoxes | ||
266 | void modCalcEquinox::findSolsticeAndEquinox(uint32_t year) | ||||
268 | { | 267 | { | ||
269 | // Interpolate to find the moment when the Sun crosses the equator | 268 | // All the magic numbers are taken from Meeus - 1991 chapter 27 | ||
270 | // Set initial guess in February or August to be sure that this | 269 | qreal m, jdme[4]; | ||
271 | // point is before equinox. | 270 | if(year > 1000) | ||
272 | const int month = Spring ? 2 : 8; | 271 | { | ||
273 | int i = QDate(year, month, 1).dayOfYear(); | 272 | m = (year-2000.0) / 1000.0; | ||
274 | double dec1, dec2; | 273 | jdme[0] = (-0.00057 * qPow(m, 4)) + (-0.00411 * qPow(m, 3)) + (0.05169 * qPow(m, 2)) + (365242.37404 * m) + 2451623.80984; | ||
275 | dec2 = ecl->points().at(i)->y(); | 274 | jdme[1] = (-0.0003 * qPow(m, 4)) + (0.00888 * qPow(m, 3)) + (0.00325 * qPow(m, 2)) + (365241.62603 * m) + 2451716.56767; | ||
276 | do | 275 | jdme[2] = (0.00078 * qPow(m, 4)) + (0.00337 * qPow(m, 3)) + (-0.11575 * qPow(m, 2)) + (365242.01767 * m) + 2451810.21715; | ||
277 | { | 276 | jdme[3] = (0.00032 * qPow(m, 4)) + (-0.00823 * qPow(m, 3)) + (-0.06223 * qPow(m, 2)) + (365242.74049 * m) + 2451900.05952; | ||
278 | ++i; | 277 | } | ||
279 | dec1 = dec2; | 278 | else | ||
280 | dec2 = ecl->points().at(i)->y(); | 279 | { | ||
281 | } | 280 | m = year / 1000.0; | ||
282 | while (dec1 * dec2 > 0.0); //when dec1*dec2<0.0, we bracket the zero | 281 | jdme[0] = (-0.00071 * qPow(m, 4)) + (0.00111 * qPow(m, 3)) + (0.06134 * qPow(m, 2)) + (365242.13740 * m) + 1721139.29189; | ||
283 | 282 | jdme[1] = (0.00025 * qPow(m, 4)) + (0.00907 * qPow(m, 3)) + (-0.05323 * qPow(m, 2)) + (365241.72562 * m) + 1721233.25401; | |||
284 | double x1 = ecl->points().at(i - 1)->x(); | 283 | jdme[2] = (0.00074 * qPow(m, 4)) + (-0.00297 * qPow(m, 3)) + (-0.11677 * qPow(m, 2)) + (365242.49558 * m) + 1721325.70455; | ||
285 | double x2 = ecl->points().at(i)->x(); | 284 | jdme[3] = (-0.00006 * qPow(m, 4)) + (-0.00933 * qPow(m, 3)) + (-0.00769 * qPow(m, 2)) + (365242.88257 * m) + 1721414.39987; | ||
286 | double d = fabs(dec2 - dec1); | 285 | } | ||
287 | double f = 1.0 - fabs(dec2) / d; //fractional distance of the zero, from point1 to point2 | | |||
288 | | ||||
289 | KStarsDateTime dt0(QDate(year, 1, 1), QTime(0, 0, 0)); | | |||
290 | KStarsDateTime dt = dt0.addSecs(86400.0 * (x1 - 1 + f * (x2 - x1))); | | |||
291 | return dt; | | |||
292 | } | | |||
293 | | ||||
294 | KStarsDateTime modCalcEquinox::findSolstice(int year, bool Summer) | | |||
295 | { | | |||
296 | //Find the moment when the Sun reaches maximum declination | | |||
297 | //First find three points which bracket the maximum (i.e., x2 > x1,x3) | | |||
298 | //Start at June 16th, which will always be approaching the solstice | | |||
299 | | ||||
300 | long double jd1, jd2, jd3, jd4; | | |||
301 | double y2(0.0), y3(0.0), y4(0.0); | | |||
302 | int month = 6; | | |||
303 | if (!Summer) | | |||
304 | month = 12; | | |||
305 | 286 | | |||
306 | jd3 = KStarsDateTime(QDate(year, month, 16), QTime(0, 0, 0)).djd(); | 287 | qreal t[4]; | ||
307 | KSNumbers num(jd3); | 288 | for(int i = 0; i < 4; i++) | ||
308 | KSSun Sun; | 289 | { | ||
309 | Sun.findPosition(&num); | 290 | t[i] = (jdme[i] - 2451545)/36525; | ||
310 | y3 = Sun.dec().Degrees(); | 291 | } | ||
292 | | ||||
293 | qreal 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}}; | ||||
294 | | ||||
295 | qreal 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}}; | ||||
311 | 296 | | |||
312 | int sgn = 1; | 297 | qreal 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}; | ||
313 | if (!Summer) | | |||
314 | sgn = -1; //find minimum if the winter solstice is sought | | |||
315 | 298 | | |||
316 | do | 299 | qreal d[4][24]; | ||
300 | | ||||
301 | for (int i = 0; i < 4; i++) | ||||
317 | { | 302 | { | ||
318 | jd3 += 1.0; | 303 | for (int j = 0; j < 24; j++) | ||
319 | num.updateValues(jd3); | 304 | { | ||
320 | Sun.findPosition(&num); | 305 | d[i][j] = c[j] * t[i]; | ||
321 | y2 = y3; | 306 | } | ||
322 | Sun.findPosition(&num); | | |||
323 | y3 = Sun.dec().Degrees(); | | |||
324 | } | 307 | } | ||
325 | while (y3 * sgn > y2 * sgn); | | |||
326 | 308 | | |||
327 | //Ok, now y2 is larger(smaller) than both y3 and y1. | 309 | for (int i = 0; i < 4; i++) | ||
328 | // Never read back | 310 | { | ||
329 | // jd2 = jd3 - 1.0; | 311 | for (int j = 0; j < 24; j++) | ||
330 | jd1 = jd3 - 2.0; | 312 | { | ||
331 | 313 | d[i][j] = qCos(qDegreesToRadians(b[i][j] + d[i][j])); | |||
332 | //Choose a new starting jd2 that follows the golden ratio: | 314 | } | ||
333 | // a/b = 1.618; a+b = 2...a = 0.76394 | 315 | } | ||
334 | jd2 = jd1 + 0.76394; | | |||
335 | num.updateValues(jd2); | | |||
336 | Sun.findPosition(&num); | | |||
337 | y2 = Sun.dec().Degrees(); | | |||
338 | 316 | | |||
339 | while (jd3 - jd1 > 0.0005) //sub-minute precision | 317 | for (int i = 0; i < 4; i++) | ||
340 | { | 318 | { | ||
341 | jd4 = jd1 + jd3 - jd2; | 319 | for (int j = 0; j < 24; j++) | ||
320 | { | ||||
321 | d[i][j] = d[i][j] * a[i][j]; | ||||
322 | } | ||||
323 | } | ||||
342 | 324 | | |||
343 | num.updateValues(jd4); | 325 | qreal e[4], f[4], g[4], jd[4]; | ||
344 | Sun.findPosition(&num); | | |||
345 | y4 = Sun.dec().Degrees(); | | |||
346 | 326 | | |||
347 | if (y4 * sgn > y2 * sgn) //make jd4 the new center | 327 | for (int i = 0; i < 4; i++) | ||
348 | { | 328 | { | ||
349 | if (jd4 > jd2) | 329 | e[i] = 0; | ||
330 | for (int j = 0; j < 24; j++) | ||||
350 | { | 331 | { | ||
351 | jd1 = jd2; | 332 | e[i] += d[i][j]; | ||
352 | jd2 = jd4; | | |||
353 | y2 = y4; | | |||
354 | } | 333 | } | ||
355 | else | 334 | } | ||
335 | | ||||
336 | for (int i = 0; i < 4; i++) | ||||
337 | { | ||||
338 | f[i] = (35999.373*t[i]-2.47); | ||||
339 | } | ||||
340 | | ||||
341 | for (int i = 0; i < 4; i++) | ||||
342 | { | ||||
343 | g[i] = 1 + (0.0334 * qCos(qDegreesToRadians(f[i]))) + (0.0007 * qCos(qDegreesToRadians(2*f[i]))); | ||||
344 | } | ||||
345 | | ||||
346 | for (int i = 0; i < 4; i++) | ||||
347 | { | ||||
348 | jd[i] = jdme[i] + (0.00001*e[i]/g[i]); | ||||
349 | } | ||||
350 | | ||||
351 | // Get correction | ||||
352 | qreal correction = FindCorrection(year); | ||||
353 | | ||||
354 | for (int i = 0; i < 4; i++) | ||||
355 | { | ||||
356 | KStarsDateTime dt(jd[i]); | ||||
357 | | ||||
358 | // Apply correction | ||||
359 | dt = dt.addSecs(correction); | ||||
360 | | ||||
361 | switch(i) | ||||
356 | { | 362 | { | ||
357 | jd3 = jd2; | 363 | case 0 : | ||
358 | // Never read back | 364 | dSpring = dt; | ||
359 | // y3 = y2; | 365 | break; | ||
360 | jd2 = jd4; | 366 | case 1 : | ||
361 | y2 = y4; | 367 | dSummer = dt; | ||
368 | break; | ||||
369 | case 2 : | ||||
370 | dAutumn = dt; | ||||
371 | break; | ||||
372 | case 3 : | ||||
373 | dWinter = dt; | ||||
374 | break; | ||||
375 | } | ||||
362 | } | 376 | } | ||
363 | } | 377 | } | ||
364 | else //make jd4 a new endpoint | 378 | | ||
379 | qreal modCalcEquinox::FindCorrection(uint32_t year) | ||||
380 | { | ||||
381 | int tblFirst = 1620, tblLast = 2002; | ||||
382 | | ||||
383 | // Corrections taken from Meeus -1991 chapter 10 | ||||
384 | qreal tbl[] = {/*1620*/ 121,112,103, 95, 88, 82, 77, 72, 68, 63, 60, 56, 53, 51, 48, 46, 44, 42, 40, 38, | ||||
385 | /*1660*/ 35, 33, 31, 29, 26, 24, 22, 20, 18, 16, 14, 12, 11, 10, 9, 8, 7, 7, 7, 7, | ||||
386 | /*1700*/ 7, 7, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 11, 11, 11, | ||||
387 | /*1740*/ 11, 11, 12, 12, 12, 12, 13, 13, 13, 14, 14, 14, 14, 15, 15, 15, 15, 15, 16, 16, | ||||
388 | /*1780*/ 16, 16, 16, 16, 16, 16, 15, 15, 14, 13, | ||||
389 | /*1800*/ 13.1, 12.5, 12.2, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 11.9, 11.6, 11.0, 10.2, 9.2, 8.2, | ||||
390 | /*1830*/ 7.1, 6.2, 5.6, 5.4, 5.3, 5.4, 5.6, 5.9, 6.2, 6.5, 6.8, 7.1, 7.3, 7.5, 7.6, | ||||
391 | /*1860*/ 7.7, 7.3, 6.2, 5.2, 2.7, 1.4, -1.2, -2.8, -3.8, -4.8, -5.5, -5.3, -5.6, -5.7, -5.9, | ||||
392 | /*1890*/ -6.0, -6.3, -6.5, -6.2, -4.7, -2.8, -0.1, 2.6, 5.3, 7.7, 10.4, 13.3, 16.0, 18.2, 20.2, | ||||
393 | /*1920*/ 21.1, 22.4, 23.5, 23.8, 24.3, 24.0, 23.9, 23.9, 23.7, 24.0, 24.3, 25.3, 26.2, 27.3, 28.2, | ||||
394 | /*1950*/ 29.1, 30.0, 30.7, 31.4, 32.2, 33.1, 34.0, 35.0, 36.5, 38.3, 40.2, 42.2, 44.5, 46.5, 48.5, | ||||
395 | /*1980*/ 50.5, 52.5, 53.8, 54.9, 55.8, 56.9, 58.3, 60.0, 61.6, 63.0, 63.8, 64.3}; | ||||
396 | | ||||
397 | qreal deltaT = 0; | ||||
398 | qreal t = (year - 2000.0)/100.0; | ||||
399 | | ||||
400 | if(year >= tblFirst && year <= tblLast) | ||||
365 | { | 401 | { | ||
366 | if (jd4 > jd2) | 402 | if(year % 2) | ||
367 | { | 403 | { | ||
368 | jd3 = jd4; | 404 | deltaT = (tbl[(year-tblFirst-1)/2] + tbl[(year-tblFirst+1)/2]) / 2; | ||
369 | // Never read back | | |||
370 | // y3 = y4; | | |||
371 | } | 405 | } | ||
372 | else | 406 | else | ||
373 | { | 407 | { | ||
374 | jd1 = jd4; | 408 | deltaT = tbl[(year-tblFirst)/2]; | ||
375 | } | 409 | } | ||
376 | } | 410 | } | ||
411 | else if(year < 948) | ||||
412 | { | ||||
413 | deltaT = 2177 + 497*t + 44.1*qPow(t, 2); | ||||
377 | } | 414 | } | ||
378 | 415 | else if(year >= 948) | |||
379 | return KStarsDateTime(jd2); | 416 | { | ||
417 | deltaT = 102 + 102*t + 25.3*qPow(t, 2); | ||||
418 | if (year>=2000 && year <=2100) | ||||
419 | { | ||||
420 | // Special correction to avoid discontinuity in 2000 | ||||
421 | deltaT += 0.37 * ( year - 2100.0 ); | ||||
422 | } | ||||
423 | } | ||||
424 | return -deltaT; | ||||
380 | } | 425 | } | ||
426 | No newline at end of file |