Kstars

modcalcvizequinox.cpp
1 /*
2  SPDX-FileCopyrightText: 2007 Jason Harris <[email protected]>
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 
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 
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 
53 void modCalcEquinox::slotCheckFiles()
54 {
55  RunButtonBatch->setEnabled(!InputFileBatch->lineEdit()->text().isEmpty() &&
56  !OutputFileBatch->lineEdit()->text().isEmpty());
57 }
58 
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 
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 
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 
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
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);
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);
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);
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
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  {
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
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 
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 }
void append(const T &value)
Extension of QDateTime for KStars KStarsDateTime can represent the date/time as a Julian Day,...
void findPosition(const KSNumbers *num, const CachingDms *lat=nullptr, const CachingDms *LST=nullptr, const KSPlanetBase *Earth=nullptr)
Find position, including correction for Figure-of-the-Earth.
void clear()
bool exists() const const
QMetaObject::Connection connect(const QObject *sender, const char *signal, const QObject *receiver, const char *method, Qt::ConnectionType type)
void setLinePen(const QPen &p)
Provides necessary information about the Sun.
Definition: kssun.h:23
QString i18n(const char *text, const TYPE &arg...)
Store several time-dependent astronomical quantities.
Definition: ksnumbers.h:42
const CachingDms & dec() const
Definition: skypoint.h:269
bool atEnd() const const
char * toString(const T &value)
ColorScheme * colorScheme()
Definition: kstarsdata.h:171
QString readLine(qint64 maxlen)
QString toString(qlonglong i) const const
int toInt(bool *ok, int base) const const
void setupUi(QWidget *widget)
void show()
void addPoint(const QPointF &p, const QString &label=QString(), double barWidth=0.0)
const double & Degrees() const
Definition: dms.h:141
QColor color() const const
void informationList(QWidget *parent, const QString &text, const QStringList &strlist, const QString &title=QString(), const QString &dontShowAgainName=QString(), Options options=Notify)
const QPen & pen() const
QColor colorNamed(const QString &name) const
Retrieve a color by name.
Definition: colorscheme.cpp:86
This file is part of the KDE documentation.
Documentation copyright © 1996-2022 The KDE developers.
Generated on Fri Aug 19 2022 03:57:52 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.