1 /*
2 SPDX-FileCopyrightText: 2007 Jason Harris <kstars@30doradus.org>
3
4 SPDX-License-Identifier: GPL-2.0-or-later
5 */
6
7 #include "modcalcvizequinox.h"
8
9 #include "dms.h"
10 #include "kstarsdata.h"
11 #include "ksnotification.h"
12 #include "skyobjects/kssun.h"
13 #include "widgets/dmsbox.h"
14
15 #include <KLineEdit>
16 #include <KPlotAxis>
17 #include <KPlotObject>
18 #include <KPlotPoint>
19
20 #include <cmath>
21
modCalcEquinox(QWidget * parentSplit)22 modCalcEquinox::modCalcEquinox(QWidget *parentSplit) : QFrame(parentSplit), dSpring(), dSummer(), dAutumn(), dWinter()
23 {
24 setupUi(this);
25
26 connect(Year, SIGNAL(valueChanged(int)), this, SLOT(slotCompute()));
27 connect(InputFileBatch, SIGNAL(urlSelected(QUrl)), this, SLOT(slotCheckFiles()));
28 connect(OutputFileBatch, SIGNAL(urlSelected(QUrl)), this, SLOT(slotCheckFiles()));
29 connect(RunButtonBatch, SIGNAL(clicked()), this, SLOT(slotRunBatch()));
30 connect(ViewButtonBatch, SIGNAL(clicked()), this, SLOT(slotViewBatch()));
31
32 Plot->axis(KPlotWidget::LeftAxis)->setLabel(i18n("Sun's Declination"));
33 Plot->setTopPadding(40);
34 //Don't draw Top & Bottom axes; custom axes drawn as plot objects
35 Plot->axis(KPlotWidget::BottomAxis)->setVisible(false);
36 Plot->axis(KPlotWidget::TopAxis)->setVisible(false);
37
38 //This will call slotCompute():
39 Year->setValue(KStarsData::Instance()->lt().date().year());
40
41 RunButtonBatch->setEnabled(false);
42 ViewButtonBatch->setEnabled(false);
43
44 show();
45 }
46
dmonth(int i)47 double modCalcEquinox::dmonth(int i)
48 {
49 Q_ASSERT(i >= 0 && i < 12 && "Month must be in 0 .. 11 range");
50 return DMonth[i];
51 }
52
slotCheckFiles()53 void modCalcEquinox::slotCheckFiles()
54 {
55 RunButtonBatch->setEnabled(!InputFileBatch->lineEdit()->text().isEmpty() &&
56 !OutputFileBatch->lineEdit()->text().isEmpty());
57 }
58
slotRunBatch()59 void modCalcEquinox::slotRunBatch()
60 {
61 QString inputFileName = InputFileBatch->url().toLocalFile();
62
63 if (QFile::exists(inputFileName))
64 {
65 QFile f(inputFileName);
66 if (!f.open(QIODevice::ReadOnly))
67 {
68 KSNotification::sorry(i18n("Could not open file %1.", f.fileName()), i18n("Could Not Open File"));
69 inputFileName.clear();
70 return;
71 }
72
73 QTextStream istream(&f);
74 processLines(istream);
75
76 ViewButtonBatch->setEnabled(true);
77
78 f.close();
79 }
80 else
81 {
82 KSNotification::sorry(i18n("Invalid file: %1", inputFileName), i18n("Invalid file"));
83 inputFileName.clear();
84 return;
85 }
86 }
87
processLines(QTextStream & istream)88 void modCalcEquinox::processLines(QTextStream &istream)
89 {
90 QFile fOut(OutputFileBatch->url().toLocalFile());
91 fOut.open(QIODevice::WriteOnly);
92 QTextStream ostream(&fOut);
93 int originalYear = Year->value();
94
95 //Write header to output file
96 ostream << i18n("# Timing of Equinoxes and Solstices\n") << i18n("# computed by KStars\n#\n")
97 << i18n("# Vernal Equinox\t\tSummer Solstice\t\t\tAutumnal Equinox\t\tWinter Solstice\n#\n");
98
99 while (!istream.atEnd())
100 {
101 QString line = istream.readLine();
102 bool ok = false;
103 int year = line.toInt(&ok);
104
105 //for now I will simply change the value of the Year widget to trigger
106 //computation of the Equinoxes and Solstices.
107 if (ok)
108 {
109 //triggers slotCompute(), which sets values of dSpring et al.:
110 Year->setValue(year);
111
112 //Write to output file
113 ostream << QLocale().toString(dSpring.date(), QLocale::LongFormat) << "\t"
114 << QLocale().toString(dSummer.date(), QLocale::LongFormat) << "\t"
115 << QLocale().toString(dAutumn.date(), QLocale::LongFormat) << "\t"
116 << QLocale().toString(dWinter.date(), QLocale::LongFormat) << '\n';
117 }
118 }
119
120 if (Year->value() != originalYear)
121 Year->setValue(originalYear);
122 }
123
slotViewBatch()124 void modCalcEquinox::slotViewBatch()
125 {
126 QFile fOut(OutputFileBatch->url().toLocalFile());
127 fOut.open(QIODevice::ReadOnly);
128 QTextStream istream(&fOut);
129 QStringList text;
130
131 while (!istream.atEnd())
132 text.append(istream.readLine());
133
134 fOut.close();
135
136 KMessageBox::informationList(nullptr, i18n("Results of Sidereal time calculation"), text,
137 OutputFileBatch->url().toLocalFile());
138 }
139
slotCompute()140 void modCalcEquinox::slotCompute()
141 {
142 KStarsData *data = KStarsData::Instance();
143 KSSun Sun;
144 int year0 = Year->value();
145
146 KStarsDateTime dt(QDate(year0, 1, 1), QTime(0, 0, 0));
147 long double jd0 = dt.djd(); //save JD on Jan 1st
148 for (int imonth = 0; imonth < 12; imonth++)
149 {
150 KStarsDateTime kdt(QDate(year0, imonth + 1, 1), QTime(0, 0, 0));
151 DMonth[imonth] = kdt.djd() - jd0;
152 }
153
154 Plot->removeAllPlotObjects();
155
156 //Add the celestial equator, just a single line bisecting the plot horizontally
157 KPlotObject *ce = new KPlotObject(data->colorScheme()->colorNamed("EqColor"), KPlotObject::Lines, 2.0);
158 ce->addPoint(0.0, 0.0);
159 ce->addPoint(366.0, 0.0);
160 Plot->addPlotObject(ce);
161
162 // Tropic of cancer
163 KPlotObject *tcan = new KPlotObject(data->colorScheme()->colorNamed("EqColor"), KPlotObject::Lines, 2.0);
164 tcan->addPoint(0.0, 23.5);
165 tcan->addPoint(366.0, 23.5);
166 Plot->addPlotObject(tcan);
167
168 // Tropic of capricorn
169 KPlotObject *tcap = new KPlotObject(data->colorScheme()->colorNamed("EqColor"), KPlotObject::Lines, 2.0);
170 tcap->addPoint(0.0, -23.5);
171 tcap->addPoint(366.0, -23.5);
172 Plot->addPlotObject(tcap);
173
174 //Add Ecliptic. This is more complicated than simply incrementing the
175 //ecliptic longitude, because we want the x-axis to be time, not RA.
176 //For each day in the year, compute the Sun's position.
177 KPlotObject *ecl = new KPlotObject(data->colorScheme()->colorNamed("EclColor"), KPlotObject::Lines, 2);
178 ecl->setLinePen(QPen(ecl->pen().color(), 4));
179
180 Plot->setLimits(1.0, double(dt.date().daysInYear()), -30.0, 30.0);
181
182 //Add top and bottom axis lines, and custom tickmarks at each month
183 addDateAxes();
184
185 for (int i = 1; i <= dt.date().daysInYear(); i++)
186 {
187 KSNumbers num(dt.djd());
188 Sun.findPosition(&num);
189 ecl->addPoint(double(i), Sun.dec().Degrees());
190
191 dt = dt.addDays(1);
192 }
193 Plot->addPlotObject(ecl);
194
195 // Calculate dates for all solstices and equinoxes
196 findSolsticeAndEquinox(Year->value());
197
198 //Display the Date/Time of each event in the text fields
199 VEquinox->setText(QLocale().toString(dSpring, QLocale::LongFormat));
200 SSolstice->setText(QLocale().toString(dSummer, QLocale::LongFormat));
201 AEquinox->setText(QLocale().toString(dAutumn, QLocale::LongFormat));
202 WSolstice->setText(QLocale().toString(dWinter, QLocale::LongFormat));
203
204 //Add vertical dotted lines at times of the equinoxes and solstices
205 KPlotObject *poSpring = new KPlotObject(Qt::white, KPlotObject::Lines, 1);
206 poSpring->setLinePen(QPen(Qt::white, 1.0, Qt::DotLine));
207 poSpring->addPoint(dSpring.djd() - jd0, Plot->dataRect().top());
208 poSpring->addPoint(dSpring.djd() - jd0, Plot->dataRect().bottom());
209 Plot->addPlotObject(poSpring);
210 KPlotObject *poSummer = new KPlotObject(Qt::white, KPlotObject::Lines, 1);
211 poSummer->setLinePen(QPen(Qt::white, 1.0, Qt::DotLine));
212 poSummer->addPoint(dSummer.djd() - jd0, Plot->dataRect().top());
213 poSummer->addPoint(dSummer.djd() - jd0, Plot->dataRect().bottom());
214 Plot->addPlotObject(poSummer);
215 KPlotObject *poAutumn = new KPlotObject(Qt::white, KPlotObject::Lines, 1);
216 poAutumn->setLinePen(QPen(Qt::white, 1.0, Qt::DotLine));
217 poAutumn->addPoint(dAutumn.djd() - jd0, Plot->dataRect().top());
218 poAutumn->addPoint(dAutumn.djd() - jd0, Plot->dataRect().bottom());
219 Plot->addPlotObject(poAutumn);
220 KPlotObject *poWinter = new KPlotObject(Qt::white, KPlotObject::Lines, 1);
221 poWinter->setLinePen(QPen(Qt::white, 1.0, Qt::DotLine));
222 poWinter->addPoint(dWinter.djd() - jd0, Plot->dataRect().top());
223 poWinter->addPoint(dWinter.djd() - jd0, Plot->dataRect().bottom());
224 Plot->addPlotObject(poWinter);
225 }
226
227 //Add custom top/bottom axes with tickmarks for each month
addDateAxes()228 void modCalcEquinox::addDateAxes()
229 {
230 KPlotObject *poTopAxis = new KPlotObject(Qt::white, KPlotObject::Lines, 1);
231 poTopAxis->addPoint(0.0, Plot->dataRect().bottom()); //y-axis is reversed!
232 poTopAxis->addPoint(366.0, Plot->dataRect().bottom());
233 Plot->addPlotObject(poTopAxis);
234
235 KPlotObject *poBottomAxis = new KPlotObject(Qt::white, KPlotObject::Lines, 1);
236 poBottomAxis->addPoint(0.0, Plot->dataRect().top() + 0.02);
237 poBottomAxis->addPoint(366.0, Plot->dataRect().top() + 0.02);
238 Plot->addPlotObject(poBottomAxis);
239
240 //Tick mark for each month
241 for (int imonth = 0; imonth < 12; imonth++)
242 {
243 KPlotObject *poMonth = new KPlotObject(Qt::white, KPlotObject::Lines, 1);
244 poMonth->addPoint(dmonth(imonth), Plot->dataRect().top());
245 poMonth->addPoint(dmonth(imonth), Plot->dataRect().top() + 1.4);
246 Plot->addPlotObject(poMonth);
247 poMonth = new KPlotObject(Qt::white, KPlotObject::Lines, 1);
248 poMonth->addPoint(dmonth(imonth), Plot->dataRect().bottom());
249 poMonth->addPoint(dmonth(imonth), Plot->dataRect().bottom() - 1.4);
250 Plot->addPlotObject(poMonth);
251 }
252 }
253
254 // Calculate and store dates for all solstices and equinoxes
findSolsticeAndEquinox(uint32_t year)255 void modCalcEquinox::findSolsticeAndEquinox(uint32_t year)
256 {
257 // All the magic numbers are taken from Meeus - 1991 chapter 27
258 qreal m, jdme[4];
259 if(year > 1000)
260 {
261 m = (year-2000.0) / 1000.0;
262 jdme[0] = (-0.00057 * qPow(m, 4)) + (-0.00411 * qPow(m, 3)) + (0.05169 * qPow(m, 2)) + (365242.37404 * m) + 2451623.80984;
263 jdme[1] = (-0.0003 * qPow(m, 4)) + (0.00888 * qPow(m, 3)) + (0.00325 * qPow(m, 2)) + (365241.62603 * m) + 2451716.56767;
264 jdme[2] = (0.00078 * qPow(m, 4)) + (0.00337 * qPow(m, 3)) + (-0.11575 * qPow(m, 2)) + (365242.01767 * m) + 2451810.21715;
265 jdme[3] = (0.00032 * qPow(m, 4)) + (-0.00823 * qPow(m, 3)) + (-0.06223 * qPow(m, 2)) + (365242.74049 * m) + 2451900.05952;
266 }
267 else
268 {
269 m = year / 1000.0;
270 jdme[0] = (-0.00071 * qPow(m, 4)) + (0.00111 * qPow(m, 3)) + (0.06134 * qPow(m, 2)) + (365242.13740 * m) + 1721139.29189;
271 jdme[1] = (0.00025 * qPow(m, 4)) + (0.00907 * qPow(m, 3)) + (-0.05323 * qPow(m, 2)) + (365241.72562 * m) + 1721233.25401;
272 jdme[2] = (0.00074 * qPow(m, 4)) + (-0.00297 * qPow(m, 3)) + (-0.11677 * qPow(m, 2)) + (365242.49558 * m) + 1721325.70455;
273 jdme[3] = (-0.00006 * qPow(m, 4)) + (-0.00933 * qPow(m, 3)) + (-0.00769 * qPow(m, 2)) + (365242.88257 * m) + 1721414.39987;
274 }
275
276 qreal t[4];
277 for(int i = 0; i < 4; i++)
278 {
279 t[i] = (jdme[i] - 2451545)/36525;
280 }
281
282 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}};
283
284 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}};
285
286 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};
287
288 qreal d[4][24];
289
290 for (int i = 0; i < 4; i++)
291 {
292 for (int j = 0; j < 24; j++)
293 {
294 d[i][j] = c[j] * t[i];
295 }
296 }
297
298 for (int i = 0; i < 4; i++)
299 {
300 for (int j = 0; j < 24; j++)
301 {
302 d[i][j] = qCos(qDegreesToRadians(b[i][j] + d[i][j]));
303 }
304 }
305
306 for (int i = 0; i < 4; i++)
307 {
308 for (int j = 0; j < 24; j++)
309 {
310 d[i][j] = d[i][j] * a[i][j];
311 }
312 }
313
314 qreal e[4], f[4], g[4], jd[4];
315
316 for (int i = 0; i < 4; i++)
317 {
318 e[i] = 0;
319 for (int j = 0; j < 24; j++)
320 {
321 e[i] += d[i][j];
322 }
323 }
324
325 for (int i = 0; i < 4; i++)
326 {
327 f[i] = (35999.373*t[i]-2.47);
328 }
329
330 for (int i = 0; i < 4; i++)
331 {
332 g[i] = 1 + (0.0334 * qCos(qDegreesToRadians(f[i]))) + (0.0007 * qCos(qDegreesToRadians(2*f[i])));
333 }
334
335 for (int i = 0; i < 4; i++)
336 {
337 jd[i] = jdme[i] + (0.00001*e[i]/g[i]);
338 }
339
340 // Get correction
341 qreal correction = FindCorrection(year);
342
343 for (int i = 0; i < 4; i++)
344 {
345 KStarsDateTime dt(jd[i]);
346
347 // Apply correction
348 dt = dt.addSecs(correction);
349
350 switch(i)
351 {
352 case 0 :
353 dSpring = dt;
354 break;
355 case 1 :
356 dSummer = dt;
357 break;
358 case 2 :
359 dAutumn = dt;
360 break;
361 case 3 :
362 dWinter = dt;
363 break;
364 }
365 }
366 }
367
FindCorrection(uint32_t year)368 qreal modCalcEquinox::FindCorrection(uint32_t year)
369 {
370 uint32_t tblFirst = 1620, tblLast = 2002;
371
372 // Corrections taken from Meeus -1991 chapter 10
373 qreal tbl[] = {/*1620*/ 121,112,103, 95, 88, 82, 77, 72, 68, 63, 60, 56, 53, 51, 48, 46, 44, 42, 40, 38,
374 /*1660*/ 35, 33, 31, 29, 26, 24, 22, 20, 18, 16, 14, 12, 11, 10, 9, 8, 7, 7, 7, 7,
375 /*1700*/ 7, 7, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 11, 11, 11,
376 /*1740*/ 11, 11, 12, 12, 12, 12, 13, 13, 13, 14, 14, 14, 14, 15, 15, 15, 15, 15, 16, 16,
377 /*1780*/ 16, 16, 16, 16, 16, 16, 15, 15, 14, 13,
378 /*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,
379 /*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,
380 /*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,
381 /*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,
382 /*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,
383 /*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,
384 /*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};
385
386 qreal deltaT = 0;
387 qreal t = (year - 2000.0)/100.0;
388
389 if(year >= tblFirst && year <= tblLast)
390 {
391 if(year % 2)
392 {
393 deltaT = (tbl[(year-tblFirst-1)/2] + tbl[(year-tblFirst+1)/2]) / 2;
394 }
395 else
396 {
397 deltaT = tbl[(year-tblFirst)/2];
398 }
399 }
400 else if(year < 948)
401 {
402 deltaT = 2177 + 497*t + 44.1*qPow(t, 2);
403 }
404 else if(year >= 948)
405 {
406 deltaT = 102 + 102*t + 25.3*qPow(t, 2);
407 if (year>=2000 && year <=2100)
408 {
409 // Special correction to avoid discontinuity in 2000
410 deltaT += 0.37 * ( year - 2100.0 );
411 }
412 }
413 return -deltaT;
414 }
415