Kstars

calibration.cpp
1/*
2 SPDX-FileCopyrightText: 2020 Hy Murveit <hy@murveit.com>
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
14Calibration::Calibration()
15{
16 ROT_Z = GuiderUtils::Matrix(0);
17}
18
19void Calibration::setAngle(double rotationAngle)
20{
21 angle = rotationAngle;
22 ROT_Z = GuiderUtils::RotateZ(-M_PI * angle / 180.0);
23}
24
25void 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
43void Calibration::setBinningUsed(int x, int y)
44{
45 subBinXused = x;
46 subBinYused = y;
47}
48
49void Calibration::setRaPulseMsPerArcsecond(double rate)
50{
51 raPulseMsPerArcsecond = rate;
52}
53
54void Calibration::setDecPulseMsPerArcsecond(double rate)
55{
56 decPulseMsPerArcsecond = rate;
57}
58
59GuiderUtils::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
67GuiderUtils::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
75void Calibration::convertToPixels(double xArcseconds, double yArcseconds,
76 double *xPixel, double *yPixel) const
77{
78 *xPixel = xArcseconds / xArcsecondsPerPixel();
79 *yPixel = yArcseconds / yArcsecondsPerPixel();
80}
81
82GuiderUtils::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
90void 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
101double Calibration::binFactor() const
102{
103 return static_cast<double>(subBinXused) / static_cast<double>(subBinX);
104}
105
106double Calibration::inverseBinFactor() const
107{
108 return 1.0 / binFactor();
109}
110
111double 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
117double Calibration::yArcsecondsPerPixel() const
118{
119 return binFactor() * (206264.806 * ccd_pixel_height * subBinY) / focalMm;
120}
121
122double Calibration::xPixelsPerArcsecond() const
123{
124 return inverseBinFactor() * (focalMm / (206264.806 * ccd_pixel_width * subBinX));
125}
126
127double Calibration::yPixelsPerArcsecond() const
128{
129 return inverseBinFactor() * (focalMm / (206264.806 * ccd_pixel_height * subBinY));
130}
131
132double Calibration::raPulseMillisecondsPerArcsecond() const
133{
134 return raPulseMsPerArcsecond;
135}
136
137double Calibration::decPulseMillisecondsPerArcsecond() const
138{
139 return decPulseMsPerArcsecond;
140}
141
142double 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
168bool 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
174bool 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
211bool 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}
222bool 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
323void 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.
333void 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}
341QString Calibration::serialize() const
342{
343 QDateTime now = QDateTime::currentDateTime();
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
358namespace
359{
360bool 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
367bool 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
375bool 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
384bool 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
444void 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
451bool 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
459double 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
491bool 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. Encoding:" << encoding;
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. Encoding: %5")
541 .arg(angle).arg(decSwap ? "T" : "F").arg(raPulseMsPerArcsecond).arg(decPulseMsPerArcsecond).arg(encoding);
542 return true;
543}
dms
An angle, stored as degrees, but expressible in many ways.
Definition dms.h:38
dms::fromString
static dms fromString(const QString &s, bool deg)
Static function to create a DMS object from a QString.
Definition dms.cpp:429
dms::toDMSString
const QString toDMSString(const bool forceSign=false, const bool machineReadable=false, const bool highPrecision=false) const
Definition dms.cpp:287
dms::toHMSString
const QString toHMSString(const bool machineReadable=false, const bool highPrecision=false) const
Definition dms.cpp:378
dms::Degrees
const double & Degrees() const
Definition dms.h:141
QDateTime::currentDateTime
QDateTime currentDateTime()
QDateTime::toString
QString toString(QStringView format, QCalendar cal) const const
QList::size
qsizetype size() const const
QString::arg
QString arg(Args &&... args) const const
QString::mid
QString mid(qsizetype position, qsizetype n) const const
QString::number
QString number(double n, char format, int precision)
QString::size
qsizetype size() const const
Qt::dec
QTextStream & dec(QTextStream &stream)
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Fri Sep 13 2024 11:53:47 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.