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

KDE's Doxygen guidelines are available online.