Kstars

calibration.cpp
1 /*
2  SPDX-FileCopyrightText: 2020 Hy Murveit <[email protected]>
3 
4  SPDX-License-Identifier: GPL-2.0-or-later
5 */
6 
7 #include "calibration.h"
8 #include "Options.h"
9 #include "ekos_guide_debug.h"
10 #include <QStringRef>
11 #include <QDateTime>
12 #include "indi/indimount.h"
13 
14 Calibration::Calibration()
15 {
16  ROT_Z = GuiderUtils::Matrix(0);
17 }
18 
19 void Calibration::setAngle(double rotationAngle)
20 {
21  angle = rotationAngle;
22  ROT_Z = GuiderUtils::RotateZ(-M_PI * angle / 180.0);
23 }
24 
25 void Calibration::setParameters(double ccd_pix_width, double ccd_pix_height,
26  double focalLengthMm,
27  int binX, int binY,
28  ISD::Mount::PierSide currentPierSide,
29  const dms &mountRa, const dms &mountDec)
30 {
31  focalMm = focalLengthMm;
32  ccd_pixel_width = ccd_pix_width;
33  ccd_pixel_height = ccd_pix_height;
34  calibrationRA = mountRa;
35  calibrationDEC = mountDec;
36  subBinX = binX;
37  subBinY = binY;
38  subBinXused = subBinX;
39  subBinYused = subBinY;
40  calibrationPierSide = currentPierSide;
41 }
42 
43 void Calibration::setBinningUsed(int x, int y)
44 {
45  subBinXused = x;
46  subBinYused = y;
47 }
48 
49 void Calibration::setRaPulseMsPerArcsecond(double rate)
50 {
51  raPulseMsPerArcsecond = rate;
52 }
53 
54 void Calibration::setDecPulseMsPerArcsecond(double rate)
55 {
56  decPulseMsPerArcsecond = rate;
57 }
58 
59 GuiderUtils::Vector Calibration::convertToArcseconds(const GuiderUtils::Vector &input) const
60 {
61  GuiderUtils::Vector arcseconds;
62  arcseconds.x = input.x * xArcsecondsPerPixel();
63  arcseconds.y = input.y * yArcsecondsPerPixel();
64  return arcseconds;
65 }
66 
67 GuiderUtils::Vector Calibration::convertToPixels(const GuiderUtils::Vector &input) const
68 {
69  GuiderUtils::Vector arcseconds;
70  arcseconds.x = input.x / xArcsecondsPerPixel();
71  arcseconds.y = input.y / yArcsecondsPerPixel();
72  return arcseconds;
73 }
74 
75 void Calibration::convertToPixels(double xArcseconds, double yArcseconds,
76  double *xPixel, double *yPixel) const
77 {
78  *xPixel = xArcseconds / xArcsecondsPerPixel();
79  *yPixel = yArcseconds / yArcsecondsPerPixel();
80 }
81 
82 GuiderUtils::Vector Calibration::rotateToRaDec(const GuiderUtils::Vector &input) const
83 {
84  GuiderUtils::Vector in;
85  in.x = input.x;
86  in.y = -input.y;
87  return (in * ROT_Z);
88 }
89 
90 void Calibration::rotateToRaDec(double dx, double dy,
91  double *ra, double *dec) const
92 {
93  GuiderUtils::Vector input;
94  input.x = dx;
95  input.y = dy;
96  GuiderUtils::Vector out = rotateToRaDec(input);
97  *ra = out.x;
98  *dec = out.y;
99 }
100 
101 double Calibration::binFactor() const
102 {
103  return static_cast<double>(subBinXused) / static_cast<double>(subBinX);
104 }
105 
106 double Calibration::inverseBinFactor() const
107 {
108  return 1.0 / binFactor();
109 }
110 
111 double Calibration::xArcsecondsPerPixel() const
112 {
113  // arcs = 3600*180/pi * (pix*ccd_pix_sz) / focal_len
114  return binFactor() * (206264.806 * ccd_pixel_width * subBinX) / focalMm;
115 }
116 
117 double Calibration::yArcsecondsPerPixel() const
118 {
119  return binFactor() * (206264.806 * ccd_pixel_height * subBinY) / focalMm;
120 }
121 
122 double Calibration::xPixelsPerArcsecond() const
123 {
124  return inverseBinFactor() * (focalMm / (206264.806 * ccd_pixel_width * subBinX));
125 }
126 
127 double Calibration::yPixelsPerArcsecond() const
128 {
129  return inverseBinFactor() * (focalMm / (206264.806 * ccd_pixel_height * subBinY));
130 }
131 
132 double Calibration::raPulseMillisecondsPerArcsecond() const
133 {
134  return raPulseMsPerArcsecond;
135 }
136 
137 double Calibration::decPulseMillisecondsPerArcsecond() const
138 {
139  return decPulseMsPerArcsecond;
140 }
141 
142 double Calibration::calculateRotation(double x, double y)
143 {
144  double phi;
145 
146  y = -y;
147 
148  //if( (!GuiderUtils::Vector(delta_x, delta_y, 0)) < 2.5 )
149  // JM 2015-12-10: Lower threshold to 1 pixel
150  if ((!GuiderUtils::Vector(x, y, 0)) < 1)
151  return -1;
152 
153  // 90 or 270 degrees
154  if (fabs(x) < fabs(y) / 1000000.0)
155  {
156  phi = y > 0 ? 90.0 : 270;
157  }
158  else
159  {
160  phi = 180.0 / M_PI * atan2(y, x);
161  if (phi < 0)
162  phi += 360.0;
163  }
164 
165  return phi;
166 }
167 
168 bool Calibration::calculate1D(double start_x, double start_y,
169  double end_x, double end_y, int RATotalPulse)
170 {
171  return calculate1D(end_x - start_x, end_y - start_y, RATotalPulse);
172 }
173 
174 bool Calibration::calculate1D(double dx, double dy, int RATotalPulse)
175 {
176  const double arcSecondsX = dx * xArcsecondsPerPixel();
177  const double arcSecondsY = dy * yArcsecondsPerPixel();
178  const double arcSeconds = std::hypot(arcSecondsX, arcSecondsY);
179  if (arcSeconds < .1 || RATotalPulse <= 0)
180  {
181  qCDebug(KSTARS_EKOS_GUIDE)
182  << QString("Bad input to calculate1D: ra %1 %2 total pulse %3")
183  .arg(dx).arg(dy).arg(RATotalPulse);
184  return false;
185  }
186 
187  const double phi = calculateRotation(arcSecondsX, arcSecondsY);
188  if (phi < 0)
189  return false;
190 
191  setAngle(phi);
192  calibrationAngle = phi;
193  calibrationAngleRA = phi;
194  calibrationAngleDEC = -1;
195  decSwap = calibrationDecSwap = false;
196 
197  if (RATotalPulse > 0)
198  setRaPulseMsPerArcsecond(RATotalPulse / arcSeconds);
199 
200  if (raPulseMillisecondsPerArcsecond() > 10000)
201  {
202  qCDebug(KSTARS_EKOS_GUIDE)
203  << "Calibration computed unreasonable pulse-milliseconds-per-arcsecond: "
204  << raPulseMillisecondsPerArcsecond() << " & " << decPulseMillisecondsPerArcsecond();
205  }
206 
207  initialized = true;
208  return true;
209 }
210 
211 bool Calibration::calculate2D(
212  double start_ra_x, double start_ra_y, double end_ra_x, double end_ra_y,
213  double start_dec_x, double start_dec_y, double end_dec_x, double end_dec_y,
214  bool *swap_dec, int RATotalPulse, int DETotalPulse)
215 {
216  return calculate2D((end_ra_x - start_ra_x),
217  (end_ra_y - start_ra_y),
218  (end_dec_x - start_dec_x),
219  (end_dec_y - start_dec_y),
220  swap_dec, RATotalPulse, DETotalPulse);
221 }
222 bool Calibration::calculate2D(
223  double ra_dx, double ra_dy, double dec_dx, double dec_dy,
224  bool *swap_dec, int RATotalPulse, int DETotalPulse)
225 {
226  const double raArcsecondsX = ra_dx * xArcsecondsPerPixel();
227  const double raArcsecondsY = ra_dy * yArcsecondsPerPixel();
228  const double decArcsecondsX = dec_dx * xArcsecondsPerPixel();
229  const double decArcsecondsY = dec_dy * yArcsecondsPerPixel();
230  const double raArcseconds = std::hypot(raArcsecondsX, raArcsecondsY);
231  const double decArcseconds = std::hypot(decArcsecondsX, decArcsecondsY);
232  if (raArcseconds < .1 || decArcseconds < .1 || RATotalPulse <= 0 || DETotalPulse <= 0)
233  {
234  qCDebug(KSTARS_EKOS_GUIDE)
235  << QString("Bad input to calculate2D: ra %1 %2 dec %3 %4 total pulses %5 %6")
236  .arg(ra_dx).arg(ra_dy).arg(dec_dx).arg(dec_dy).arg(RATotalPulse).arg(DETotalPulse);
237  return false;
238  }
239 
240  double phi_ra = 0; // angle calculated by GUIDE_RA drift
241  double phi_dec = 0; // angle calculated by GUIDE_DEC drift
242  double phi = 0;
243 
244  GuiderUtils::Vector ra_vect = GuiderUtils::Normalize(GuiderUtils::Vector(raArcsecondsX, -raArcsecondsY, 0));
245  GuiderUtils::Vector dec_vect = GuiderUtils::Normalize(GuiderUtils::Vector(decArcsecondsX, -decArcsecondsY, 0));
246 
247  GuiderUtils::Vector try_increase = dec_vect * GuiderUtils::RotateZ(M_PI / 2);
248  GuiderUtils::Vector try_decrease = dec_vect * GuiderUtils::RotateZ(-M_PI / 2);
249 
250  double cos_increase = try_increase & ra_vect;
251  double cos_decrease = try_decrease & ra_vect;
252 
253  bool do_increase = cos_increase > cos_decrease ? true : false;
254 
255  phi_ra = calculateRotation(raArcsecondsX, raArcsecondsY);
256  if (phi_ra < 0)
257  return false;
258 
259  phi_dec = calculateRotation(decArcsecondsX, decArcsecondsY);
260  if (phi_dec < 0)
261  return false;
262 
263  // Store the calibration angles.
264  calibrationAngleRA = phi_ra;
265  calibrationAngleDEC = phi_dec;
266 
267  if (do_increase)
268  phi_dec += 90;
269  else
270  phi_dec -= 90;
271 
272  if (phi_dec > 360)
273  phi_dec -= 360.0;
274  if (phi_dec < 0)
275  phi_dec += 360.0;
276 
277  if (fabs(phi_dec - phi_ra) > 180)
278  {
279  if (phi_ra > phi_dec)
280  phi_ra -= 360;
281  else
282  phi_dec -= 360;
283  }
284 
285  // average angles
286  phi = (phi_ra + phi_dec) / 2;
287  if (phi < 0)
288  phi += 360.0;
289 
290  // setAngle sets the angle we'll use when guiding.
291  // calibrationAngle is the saved angle to be stored.
292  // They're the same now, but if we flip pier sides, angle may change.
293  setAngle(phi);
294  calibrationAngle = phi;
295 
296  // check DEC
297  if (swap_dec)
298  *swap_dec = decSwap = do_increase ? false : true;
299  calibrationDecSwap = decSwap;
300 
301  if (RATotalPulse > 0)
302  setRaPulseMsPerArcsecond(RATotalPulse / raArcseconds);
303 
304  if (DETotalPulse > 0)
305  setDecPulseMsPerArcsecond(DETotalPulse / decArcseconds);
306 
307  // Check for unreasonable values.
308  if (raPulseMillisecondsPerArcsecond() > 10000 || decPulseMillisecondsPerArcsecond() > 10000)
309  {
310  qCDebug(KSTARS_EKOS_GUIDE)
311  << "Calibration computed unreasonable pulse-milliseconds-per-arcsecond: "
312  << raPulseMillisecondsPerArcsecond() << " & " << decPulseMillisecondsPerArcsecond();
313  return false;
314  }
315 
316  qCDebug(KSTARS_EKOS_GUIDE) << QString("Set RA ms/as = %1ms / %2as = %3. DEC: %4ms / %5px = %6.")
317  .arg(RATotalPulse).arg(raArcseconds).arg(raPulseMillisecondsPerArcsecond())
318  .arg(DETotalPulse).arg(decArcseconds).arg(decPulseMillisecondsPerArcsecond());
319  initialized = true;
320  return true;
321 }
322 
323 void Calibration::computeDrift(const GuiderUtils::Vector &detection, const GuiderUtils::Vector &reference,
324  double *raDrift, double *decDrift) const
325 {
326  GuiderUtils::Vector drift = detection - reference;
327  drift = rotateToRaDec(drift);
328  *raDrift = drift.x;
329  *decDrift = drift.y;
330 }
331 
332 // Not sure why we allow this.
333 void Calibration::setDeclinationSwapEnabled(bool value)
334 {
335  qCDebug(KSTARS_EKOS_GUIDE) << QString("decSwap set to %1, was %2, cal %3")
336  .arg(value ? "T" : "F")
337  .arg(decSwap ? "T" : "F")
338  .arg(calibrationDecSwap ? "T" : "F");
339  decSwap = value;
340 }
341 QString Calibration::serialize() const
342 {
344  QString dateStr = now.toString("yyyy-MM-dd hh:mm:ss");
345  QString raStr = std::isnan(calibrationRA.Degrees()) ? "" : calibrationRA.toDMSString(false, true, true);
346  QString decStr = std::isnan(calibrationDEC.Degrees()) ? "" : calibrationDEC.toHMSString(true, true);
347  QString s =
348  QString("Cal v1.0,bx=%1,by=%2,pw=%3,ph=%4,fl=%5,ang=%6,angR=%7,angD=%8,"
349  "ramspas=%9,decmspas=%10,swap=%11,ra=%12,dec=%13,side=%14,when=%15,calEnd")
350  .arg(subBinX).arg(subBinY).arg(ccd_pixel_width).arg(ccd_pixel_height)
351  .arg(focalMm).arg(calibrationAngle).arg(calibrationAngleRA)
352  .arg(calibrationAngleDEC).arg(raPulseMsPerArcsecond)
353  .arg(decPulseMsPerArcsecond).arg(calibrationDecSwap ? 1 : 0)
354  .arg(raStr).arg(decStr).arg(static_cast<int>(calibrationPierSide)).arg(dateStr);
355  return s;
356 }
357 
358 namespace
359 {
360 bool parseString(const QStringRef &ref, const QString &id, QString *result)
361 {
362  if (!ref.startsWith(id)) return false;
363  *result = ref.mid(id.size()).toString();
364  return true;
365 }
366 
367 bool parseDouble(const QStringRef &ref, const QString &id, double *result)
368 {
369  bool ok;
370  if (!ref.startsWith(id)) return false;
371  *result = ref.mid(id.size()).toDouble(&ok);
372  return ok;
373 }
374 
375 bool parseInt(const QStringRef &ref, const QString &id, int *result)
376 {
377  bool ok;
378  if (!ref.startsWith(id)) return false;
379  *result = ref.mid(id.size()).toInt(&ok);
380  return ok;
381 }
382 } // namespace
383 
384 bool Calibration::restore(const QString &encoding)
385 {
386  QVector<QStringRef> items = encoding.splitRef(',');
387  if (items.size() != 17) return false;
388  int i = 0;
389  if (items[i] != "Cal v1.0") return false;
390  if (!parseInt(items[++i], "bx=", &subBinX)) return false;
391  if (!parseInt(items[++i], "by=", &subBinY)) return false;
392  if (!parseDouble(items[++i], "pw=", &ccd_pixel_width)) return false;
393  if (!parseDouble(items[++i], "ph=", &ccd_pixel_height)) return false;
394  if (!parseDouble(items[++i], "fl=", &focalMm)) return false;
395  if (!parseDouble(items[++i], "ang=", &calibrationAngle)) return false;
396  setAngle(calibrationAngle);
397  if (!parseDouble(items[++i], "angR=", &calibrationAngleRA)) return false;
398  if (!parseDouble(items[++i], "angD=", &calibrationAngleDEC)) return false;
399 
400  // Switched from storing raPulseMsPerPixel to ...PerArcsecond and similar for DEC.
401  // The below allows back-compatibility for older stored calibrations.
402  if (!parseDouble(items[++i], "ramspas=", &raPulseMsPerArcsecond))
403  {
404  // Try the old version
405  double raPulseMsPerPixel;
406  if (parseDouble(items[i], "ramspp=", &raPulseMsPerPixel) && (xArcsecondsPerPixel() > 0))
407  {
408  // The previous calibration only worked on square pixels.
409  raPulseMsPerArcsecond = raPulseMsPerPixel / xArcsecondsPerPixel();
410  }
411  else return false;
412  }
413  if (!parseDouble(items[++i], "decmspas=", &decPulseMsPerArcsecond))
414  {
415  // Try the old version
416  double decPulseMsPerPixel;
417  if (parseDouble(items[i], "decmspp=", &decPulseMsPerPixel) && (yArcsecondsPerPixel() > 0))
418  {
419  // The previous calibration only worked on square pixels.
420  decPulseMsPerArcsecond = decPulseMsPerPixel / yArcsecondsPerPixel();
421  }
422  else return false;
423  }
424 
425  int tempInt;
426  if (!parseInt(items[++i], "swap=", &tempInt)) return false;
427  decSwap = static_cast<bool>(tempInt);
428  calibrationDecSwap = decSwap;
429  QString tempStr;
430  if (!parseString(items[++i], "ra=", &tempStr)) return false;
431  dms nullDms;
432  calibrationRA = tempStr.size() == 0 ? nullDms : dms::fromString(tempStr, true);
433  if (!parseString(items[++i], "dec=", &tempStr)) return false;
434  calibrationDEC = tempStr.size() == 0 ? nullDms : dms::fromString(tempStr, false);
435  if (!parseInt(items[++i], "side=", &tempInt)) return false;
436  calibrationPierSide = static_cast<ISD::Mount::PierSide>(tempInt);
437  if (!parseString(items[++i], "when=", &tempStr)) return false;
438  // Don't keep the time. It's just for reference.
439  if (items[++i] != "calEnd") return false;
440  initialized = true;
441  return true;
442 }
443 
444 void Calibration::save() const
445 {
446  QString encoding = serialize();
447  Options::setSerializedCalibration(encoding);
448  qCDebug(KSTARS_EKOS_GUIDE) << QString("Saved calibration: %1").arg(encoding);
449 }
450 
451 bool Calibration::restore(ISD::Mount::PierSide currentPierSide,
452  bool reverseDecOnPierChange, int currentBinX, int currentBinY,
453  const dms *declination)
454 {
455  return restore(Options::serializedCalibration(), currentPierSide,
456  reverseDecOnPierChange, currentBinX, currentBinY, declination);
457 }
458 
459 double Calibration::correctRA(double raMsPerArcsec, const dms &calibrationDec, const dms &currentDec)
460 {
461  constexpr double MAX_DEC = 60.0;
462  // Don't use uninitialized dms values.
463  if (std::isnan(calibrationDec.Degrees()) || std::isnan(currentDec.Degrees()))
464  {
465  qCDebug(KSTARS_EKOS_GUIDE) << QString("Did not convert calibration RA rate. Input DEC invalid");
466  return raMsPerArcsec;
467  }
468  if ((calibrationDec.Degrees() > MAX_DEC) || (calibrationDec.Degrees() < -MAX_DEC))
469  {
470  qCDebug(KSTARS_EKOS_GUIDE) << QString("Didn't correct calibration RA rate. Calibration DEC %1 too close to pole")
471  .arg(QString::number(calibrationDec.Degrees(), 'f', 1));
472  return raMsPerArcsec;
473  }
474  const double toRadians = M_PI / 180.0;
475  const double numer = std::cos(currentDec.Degrees() * toRadians);
476  const double denom = std::cos(calibrationDec.Degrees() * toRadians);
477  if (raMsPerArcsec == 0) return raMsPerArcsec;
478  const double rate = 1.0 / raMsPerArcsec;
479  if (denom == 0) return raMsPerArcsec;
480  const double adjustedRate = numer * rate / denom;
481  if (adjustedRate == 0) return raMsPerArcsec;
482  const double adjustedMsPerArcsecond = 1.0 / adjustedRate;
483  qCDebug(KSTARS_EKOS_GUIDE) << QString("Corrected calibration RA rate. %1 --> %2. Calibration DEC %3 current DEC %4.")
484  .arg(QString::number(raMsPerArcsec, 'f', 1))
485  .arg(QString::number(adjustedMsPerArcsecond, 'f', 1))
486  .arg(QString::number(calibrationDec.Degrees(), 'f', 1))
487  .arg(QString::number(currentDec.Degrees(), 'f', 1));
488  return adjustedMsPerArcsecond;
489 }
490 
491 bool Calibration::restore(const QString &encoding, ISD::Mount::PierSide currentPierSide,
492  bool reverseDecOnPierChange, int currentBinX, int currentBinY,
493  const dms *currentDeclination)
494 {
495  // Fail if we couldn't read the calibration.
496  if (!restore(encoding))
497  {
498  qCDebug(KSTARS_EKOS_GUIDE) << "Could not restore calibration--couldn't read.";
499  return false;
500  }
501  // We don't restore calibrations where either the calibration or the current pier side
502  // is unknown.
503  if (calibrationPierSide == ISD::Mount::PIER_UNKNOWN ||
504  currentPierSide == ISD::Mount::PIER_UNKNOWN)
505  {
506  qCDebug(KSTARS_EKOS_GUIDE) << "Could not restore calibration--pier side unknown.";
507  return false;
508  }
509 
510  if (currentDeclination != nullptr)
511  raPulseMsPerArcsecond = correctRA(raPulseMsPerArcsecond, calibrationDEC, *currentDeclination);
512 
513  subBinXused = currentBinX;
514  subBinYused = currentBinY;
515 
516  // Succeed if the calibration was on the same side of the pier as we're currently using.
517  if (currentPierSide == calibrationPierSide)
518  {
519  qCDebug(KSTARS_EKOS_GUIDE) << "Restored calibration--same pier side.";
520  return true;
521  }
522 
523  // Otherwise, swap the angles and succeed.
524  angle = angle + 180.0;
525  while (angle >= 360.0)
526  angle = angle - 360.0;
527  while (angle < 0.0)
528  angle = angle + 360.0;
529  // Although angle is already set, we do this to set the rotation matrix.
530  setAngle(angle);
531 
532  // Note the negation in testing the option.
533  // Since above the angle rotated by 180 degrees, both RA and DEC have already reversed.
534  // If the user says not to reverse DEC, then we re-reverse by setting decSwap.
535  // Doing it this way makes the option consistent with other programs the user may be familiar with (PHD2).
536  if (!reverseDecOnPierChange)
537  decSwap = !calibrationDecSwap;
538 
539  qCDebug(KSTARS_EKOS_GUIDE)
540  << QString("Restored calibration--flipped angles. Angle %1, swap %2 ms/as: %3 %4")
541  .arg(angle).arg(decSwap ? "T" : "F").arg(raPulseMsPerArcsecond).arg(decPulseMsPerArcsecond);
542  return true;
543 }
QString number(int n, int base)
int size() const const
QDateTime currentDateTime()
void ref()
QVector< QStringRef > splitRef(const QString &sep, QString::SplitBehavior behavior, Qt::CaseSensitivity cs) const const
QTextStream & dec(QTextStream &stream)
An angle, stored as degrees, but expressible in many ways.
Definition: dms.h:37
QString arg(qlonglong a, int fieldWidth, int base, QChar fillChar) const const
const double & Degrees() const
Definition: dms.h:141
int size() const const
QString toString(Qt::DateFormat format) const const
static dms fromString(const QString &s, bool deg)
Static function to create a DMS object from a QString.
Definition: dms.cpp:421
This file is part of the KDE documentation.
Documentation copyright © 1996-2022 The KDE developers.
Generated on Fri Aug 12 2022 04:00:52 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.