7#include "calibration.h"
9#include "ekos_guide_debug.h"
12#include "indi/indimount.h"
14Calibration::Calibration()
16 ROT_Z = GuiderUtils::Matrix(0);
19void Calibration::setAngle(
double rotationAngle)
21 angle = rotationAngle;
22 ROT_Z = GuiderUtils::RotateZ(-M_PI * angle / 180.0);
25void Calibration::setParameters(
double ccd_pix_width,
double ccd_pix_height,
28 ISD::Mount::PierSide currentPierSide,
29 const dms &mountRa,
const dms &mountDec)
31 focalMm = focalLengthMm;
32 ccd_pixel_width = ccd_pix_width;
33 ccd_pixel_height = ccd_pix_height;
34 calibrationRA = mountRa;
35 calibrationDEC = mountDec;
38 subBinXused = subBinX;
39 subBinYused = subBinY;
40 calibrationPierSide = currentPierSide;
43void Calibration::setBinningUsed(
int x,
int y)
49void Calibration::setRaPulseMsPerArcsecond(
double rate)
51 raPulseMsPerArcsecond = rate;
54void Calibration::setDecPulseMsPerArcsecond(
double rate)
56 decPulseMsPerArcsecond = rate;
59GuiderUtils::Vector Calibration::convertToArcseconds(
const GuiderUtils::Vector &input)
const
61 GuiderUtils::Vector arcseconds;
62 arcseconds.x = input.x * xArcsecondsPerPixel();
63 arcseconds.y = input.y * yArcsecondsPerPixel();
67GuiderUtils::Vector Calibration::convertToPixels(
const GuiderUtils::Vector &input)
const
69 GuiderUtils::Vector arcseconds;
70 arcseconds.x = input.x / xArcsecondsPerPixel();
71 arcseconds.y = input.y / yArcsecondsPerPixel();
75void Calibration::convertToPixels(
double xArcseconds,
double yArcseconds,
76 double *xPixel,
double *yPixel)
const
78 *xPixel = xArcseconds / xArcsecondsPerPixel();
79 *yPixel = yArcseconds / yArcsecondsPerPixel();
82GuiderUtils::Vector Calibration::rotateToRaDec(
const GuiderUtils::Vector &input)
const
84 GuiderUtils::Vector in;
90void Calibration::rotateToRaDec(
double dx,
double dy,
91 double *ra,
double *dec)
const
93 GuiderUtils::Vector input;
96 GuiderUtils::Vector out = rotateToRaDec(input);
101double Calibration::binFactor()
const
103 return static_cast<double>(subBinXused) /
static_cast<double>(subBinX);
106double Calibration::inverseBinFactor()
const
108 return 1.0 / binFactor();
111double Calibration::xArcsecondsPerPixel()
const
114 return binFactor() * (206264.806 * ccd_pixel_width * subBinX) / focalMm;
117double Calibration::yArcsecondsPerPixel()
const
119 return binFactor() * (206264.806 * ccd_pixel_height * subBinY) / focalMm;
122double Calibration::xPixelsPerArcsecond()
const
124 return inverseBinFactor() * (focalMm / (206264.806 * ccd_pixel_width * subBinX));
127double Calibration::yPixelsPerArcsecond()
const
129 return inverseBinFactor() * (focalMm / (206264.806 * ccd_pixel_height * subBinY));
132double Calibration::raPulseMillisecondsPerArcsecond()
const
134 return raPulseMsPerArcsecond;
137double Calibration::decPulseMillisecondsPerArcsecond()
const
139 return decPulseMsPerArcsecond;
142double Calibration::calculateRotation(
double x,
double y)
150 if ((!GuiderUtils::Vector(x, y, 0)) < 1)
154 if (fabs(x) < fabs(y) / 1000000.0)
156 phi = y > 0 ? 90.0 : 270;
160 phi = 180.0 / M_PI * atan2(y, x);
168bool Calibration::calculate1D(
double start_x,
double start_y,
169 double end_x,
double end_y,
int RATotalPulse)
171 return calculate1D(end_x - start_x, end_y - start_y, RATotalPulse);
174bool Calibration::calculate1D(
double dx,
double dy,
int RATotalPulse)
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)
181 qCDebug(KSTARS_EKOS_GUIDE)
182 <<
QString(
"Bad input to calculate1D: ra %1 %2 total pulse %3")
187 const double phi = calculateRotation(arcSecondsX, arcSecondsY);
192 calibrationAngle = phi;
193 calibrationAngleRA = phi;
194 calibrationAngleDEC = -1;
195 decSwap = calibrationDecSwap =
false;
197 if (RATotalPulse > 0)
198 setRaPulseMsPerArcsecond(RATotalPulse / arcSeconds);
200 if (raPulseMillisecondsPerArcsecond() > 10000)
202 qCDebug(KSTARS_EKOS_GUIDE)
203 <<
"Calibration computed unreasonable pulse-milliseconds-per-arcsecond: "
204 << raPulseMillisecondsPerArcsecond() <<
" & " << decPulseMillisecondsPerArcsecond();
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)
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);
222bool Calibration::calculate2D(
223 double ra_dx,
double ra_dy,
double dec_dx,
double dec_dy,
224 bool *swap_dec,
int RATotalPulse,
int DETotalPulse)
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)
234 qCDebug(KSTARS_EKOS_GUIDE)
235 <<
QString(
"Bad input to calculate2D: ra %1 %2 dec %3 %4 total pulses %5 %6")
244 GuiderUtils::Vector ra_vect = GuiderUtils::Normalize(GuiderUtils::Vector(raArcsecondsX, -raArcsecondsY, 0));
245 GuiderUtils::Vector dec_vect = GuiderUtils::Normalize(GuiderUtils::Vector(decArcsecondsX, -decArcsecondsY, 0));
247 GuiderUtils::Vector try_increase = dec_vect * GuiderUtils::RotateZ(M_PI / 2);
248 GuiderUtils::Vector try_decrease = dec_vect * GuiderUtils::RotateZ(-M_PI / 2);
250 double cos_increase = try_increase & ra_vect;
251 double cos_decrease = try_decrease & ra_vect;
253 bool do_increase = cos_increase > cos_decrease ? true :
false;
255 phi_ra = calculateRotation(raArcsecondsX, raArcsecondsY);
259 phi_dec = calculateRotation(decArcsecondsX, decArcsecondsY);
264 calibrationAngleRA = phi_ra;
265 calibrationAngleDEC = phi_dec;
277 if (fabs(phi_dec - phi_ra) > 180)
279 if (phi_ra > phi_dec)
286 phi = (phi_ra + phi_dec) / 2;
294 calibrationAngle = phi;
298 *swap_dec = decSwap = do_increase ? false :
true;
299 calibrationDecSwap = decSwap;
301 if (RATotalPulse > 0)
302 setRaPulseMsPerArcsecond(RATotalPulse / raArcseconds);
304 if (DETotalPulse > 0)
305 setDecPulseMsPerArcsecond(DETotalPulse / decArcseconds);
308 if (raPulseMillisecondsPerArcsecond() > 10000 || decPulseMillisecondsPerArcsecond() > 10000)
310 qCDebug(KSTARS_EKOS_GUIDE)
311 <<
"Calibration computed unreasonable pulse-milliseconds-per-arcsecond: "
312 << raPulseMillisecondsPerArcsecond() <<
" & " << decPulseMillisecondsPerArcsecond();
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());
323void Calibration::computeDrift(
const GuiderUtils::Vector &detection,
const GuiderUtils::Vector &reference,
324 double *raDrift,
double *decDrift)
const
326 GuiderUtils::Vector drift = detection - reference;
327 drift = rotateToRaDec(drift);
333void Calibration::setDeclinationSwapEnabled(
bool value)
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");
341QString Calibration::serialize()
const
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);
360bool parseString(
const QStringRef &ref,
const QString &
id,
QString *result)
362 if (!ref.startsWith(
id))
return false;
363 *result = ref.
mid(
id.size()).toString();
367bool parseDouble(
const QStringRef &ref,
const QString &
id,
double *result)
370 if (!ref.startsWith(
id))
return false;
371 *result = ref.mid(
id.size()).toDouble(&ok);
375bool parseInt(
const QStringRef &ref,
const QString &
id,
int *result)
378 if (!ref.startsWith(
id))
return false;
379 *result = ref.mid(
id.size()).toInt(&ok);
384bool Calibration::restore(
const QString &encoding)
387 if (items.
size() != 17)
return false;
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;
402 if (!parseDouble(items[++i],
"ramspas=", &raPulseMsPerArcsecond))
405 double raPulseMsPerPixel;
406 if (parseDouble(items[i],
"ramspp=", &raPulseMsPerPixel) && (xArcsecondsPerPixel() > 0))
409 raPulseMsPerArcsecond = raPulseMsPerPixel / xArcsecondsPerPixel();
413 if (!parseDouble(items[++i],
"decmspas=", &decPulseMsPerArcsecond))
416 double decPulseMsPerPixel;
417 if (parseDouble(items[i],
"decmspp=", &decPulseMsPerPixel) && (yArcsecondsPerPixel() > 0))
420 decPulseMsPerArcsecond = decPulseMsPerPixel / yArcsecondsPerPixel();
426 if (!parseInt(items[++i],
"swap=", &tempInt))
return false;
427 decSwap =
static_cast<bool>(tempInt);
428 calibrationDecSwap = decSwap;
430 if (!parseString(items[++i],
"ra=", &tempStr))
return false;
433 if (!parseString(items[++i],
"dec=", &tempStr))
return 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;
439 if (items[++i] !=
"calEnd")
return false;
444void Calibration::save()
const
446 QString encoding = serialize();
447 Options::setSerializedCalibration(encoding);
448 qCDebug(KSTARS_EKOS_GUIDE) <<
QString(
"Saved calibration: %1").
arg(encoding);
451bool Calibration::restore(ISD::Mount::PierSide currentPierSide,
452 bool reverseDecOnPierChange,
int currentBinX,
int currentBinY,
453 const dms *declination)
455 return restore(Options::serializedCalibration(), currentPierSide,
456 reverseDecOnPierChange, currentBinX, currentBinY, declination);
459double Calibration::correctRA(
double raMsPerArcsec,
const dms &calibrationDec,
const dms ¤tDec)
461 constexpr double MAX_DEC = 60.0;
463 if (std::isnan(calibrationDec.
Degrees()) || std::isnan(currentDec.
Degrees()))
465 qCDebug(KSTARS_EKOS_GUIDE) <<
QString(
"Did not convert calibration RA rate. Input DEC invalid");
466 return raMsPerArcsec;
468 if ((calibrationDec.
Degrees() > MAX_DEC) || (calibrationDec.
Degrees() < -MAX_DEC))
470 qCDebug(KSTARS_EKOS_GUIDE) <<
QString(
"Didn't correct calibration RA rate. Calibration DEC %1 too close to pole")
472 return raMsPerArcsec;
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.")
488 return adjustedMsPerArcsecond;
491bool Calibration::restore(
const QString &encoding, ISD::Mount::PierSide currentPierSide,
492 bool reverseDecOnPierChange,
int currentBinX,
int currentBinY,
493 const dms *currentDeclination)
496 if (!restore(encoding))
498 qCDebug(KSTARS_EKOS_GUIDE) <<
"Could not restore calibration--couldn't read.";
503 if (calibrationPierSide == ISD::Mount::PIER_UNKNOWN ||
504 currentPierSide == ISD::Mount::PIER_UNKNOWN)
506 qCDebug(KSTARS_EKOS_GUIDE) <<
"Could not restore calibration--pier side unknown.";
510 if (currentDeclination !=
nullptr)
511 raPulseMsPerArcsecond = correctRA(raPulseMsPerArcsecond, calibrationDEC, *currentDeclination);
513 subBinXused = currentBinX;
514 subBinYused = currentBinY;
517 if (currentPierSide == calibrationPierSide)
519 qCDebug(KSTARS_EKOS_GUIDE) <<
"Restored calibration--same pier side. Encoding:" << encoding;
524 angle = angle + 180.0;
525 while (angle >= 360.0)
526 angle = angle - 360.0;
528 angle = angle + 360.0;
536 if (!reverseDecOnPierChange)
537 decSwap = !calibrationDecSwap;
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);
An angle, stored as degrees, but expressible in many ways.
static dms fromString(const QString &s, bool deg)
Static function to create a DMS object from a QString.
const QString toDMSString(const bool forceSign=false, const bool machineReadable=false, const bool highPrecision=false) const
const QString toHMSString(const bool machineReadable=false, const bool highPrecision=false) const
const double & Degrees() const
QDateTime currentDateTime()
QString toString(QStringView format, QCalendar cal) const const
qsizetype size() const const
QString arg(Args &&... args) const const
QString mid(qsizetype position, qsizetype n) const const
QString number(double n, char format, int precision)
qsizetype size() const const
QTextStream & dec(QTextStream &stream)