Differential D27413 Diff 75889 kstars/tools/modcalcvizequinox.cpp

Changeset View

Standalone View

View Options

kstars/tools/modcalcvizequinox.cpp

Context not available.
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));
Context not available.
264 } 262​ }
265} 263​}
266 264​
267KStarsDateTime 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; ​
273 int i = QDate(year, month, 1).dayOfYear(); ​
274 double dec1, dec2; ​
275 dec2 = ecl->points().at(i)->y(); ​
276 do ​
277 { 271​ {
278 ++i; 272​ m = (year-2000.0) / 1000.0;
279 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;
280 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;
275​ jdme[2] = (0.00078 * qPow(m, 4)) + (0.00337 * qPow(m, 3)) + (-0.11575 * qPow(m, 2)) + (365242.01767 * m) + 2451810.21715;
276​ jdme[3] = (0.00032 * qPow(m, 4)) + (-0.00823 * qPow(m, 3)) + (-0.06223 * qPow(m, 2)) + (365242.74049 * m) + 2451900.05952;
277​ }
278​ else
279​ {
280​ m = year / 1000.0;
281​ jdme[0] = (-0.00071 * qPow(m, 4)) + (0.00111 * qPow(m, 3)) + (0.06134 * qPow(m, 2)) + (365242.13740 * m) + 1721139.29189;
282​ jdme[1] = (0.00025 * qPow(m, 4)) + (0.00907 * qPow(m, 3)) + (-0.05323 * qPow(m, 2)) + (365241.72562 * m) + 1721233.25401;
283​ jdme[2] = (0.00074 * qPow(m, 4)) + (-0.00297 * qPow(m, 3)) + (-0.11677 * qPow(m, 2)) + (365242.49558 * m) + 1721325.70455;
284​ jdme[3] = (-0.00006 * qPow(m, 4)) + (-0.00933 * qPow(m, 3)) + (-0.00769 * qPow(m, 2)) + (365242.88257 * m) + 1721414.39987;
281 } 285​ }
282 while (dec1 * dec2 > 0.0); //when dec1*dec2<0.0, we bracket the zero ​
283 286​
284 double x1 = ecl->points().at(i - 1)->x(); 287​ qreal t[4];
285 double x2 = ecl->points().at(i)->x(); 288​ for(int i = 0; i < 4; i++)
286 double d = fabs(dec2 - dec1); 289​ {
287 double f = 1.0 - fabs(dec2) / d; //fractional distance of the zero, from point1 to point2 290​ t[i] = (jdme[i] - 2451545)/36525;
291​ }
288 292​
289 KStarsDateTime dt0(QDate(year, 1, 1), QTime(0, 0, 0)); 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}};
290 KStarsDateTime dt = dt0.addSecs(86400.0 * (x1 - 1 + f * (x2 - x1))); ​
291 return dt; ​
292} ​
293 294​
294KStarsDateTime modCalcEquinox::findSolstice(int year, bool Summer) 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}};
295{ 296​
296 //Find the moment when the Sun reaches maximum declination 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};
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 ​
306 jd3 = KStarsDateTime(QDate(year, month, 16), QTime(0, 0, 0)).djd(); ​
307 KSNumbers num(jd3); ​
308 KSSun Sun; ​
309 Sun.findPosition(&num); ​
310 y3 = Sun.dec().Degrees(); ​
311 298​
312 int sgn = 1; 299​ qreal d[4][24];
313 if (!Summer) ​
314 sgn = -1; //find minimum if the winter solstice is sought ​
315 300​
316 do 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 ​
327 //Ok, now y2 is larger(smaller) than both y3 and y1. ​
328 // Never read back ​
329 // jd2 = jd3 - 1.0; ​
330 jd1 = jd3 - 2.0; ​
331 308​
332 //Choose a new starting jd2 that follows the golden ratio: 309​ for (int i = 0; i < 4; i++)
333 // a/b = 1.618; a+b = 2...a = 0.76394 310​ {
334 jd2 = jd1 + 0.76394; 311​ for (int j = 0; j < 24; j++)
335 num.updateValues(jd2); 312​ {
336 Sun.findPosition(&num); 313​ d[i][j] = qCos(qDegreesToRadians(b[i][j] + d[i][j]));
337 y2 = Sun.dec().Degrees(); 314​ }
315​ }
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++)
328​ {
329​ e[i] = 0;
330​ for (int j = 0; j < 24; j++)
348 { 331​ {
349 if (jd4 > jd2) 332​ e[i] += d[i][j];
350 { ​
351 jd1 = jd2; ​
352 jd2 = jd4; ​
353 y2 = y4; ​
354 } ​
355 else ​
356 { ​
357 jd3 = jd2; ​
358 // Never read back ​
359 // y3 = y2; ​
360 jd2 = jd4; ​
361 y2 = y4; ​
362 } ​
363 } 333​ }
364 else //make jd4 a new endpoint 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)
365 { 362​ {
366 if (jd4 > jd2) 363​ case 0 :
367 { 364​ dSpring = dt;
368 jd3 = jd4; 365​ break;
369 // Never read back 366​ case 1 :
370 // y3 = y4; 367​ dSummer = dt;
371 } 368​ break;
372 else 369​ case 2 :
373 { 370​ dAutumn = dt;
374 jd1 = jd4; 371​ break;
375 } 372​ case 3 :
373​ dWinter = dt;
374​ break;
376 } 375​ }
377 } 376​ }
378 ​
379 return KStarsDateTime(jd2); ​
380} 377​}
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)
401​ {
402​ if(year % 2)
403​ {
404​ deltaT = (tbl[(year-tblFirst-1)/2] + tbl[(year-tblFirst+1)/2]) / 2;
405​ }
406​ else
407​ {
408​ deltaT = tbl[(year-tblFirst)/2];
409​ }
410​ }
411​ else if(year < 948)
412​ {
413​ deltaT = 2177 + 497*t + 44.1*qPow(t, 2);
414​ }
415​ else if(year >= 948)
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;
425​}
426​No newline at end of file
Context not available.