Instrument Neutral Distributed Interface INDI  2.0.2
scopesim_helper.cpp
Go to the documentation of this file.
1 /*
2  * scopesim_helper.cpp
3  *
4  * Copyright 2020 Chris Rowland <chris.rowland@cherryfield.me.uk>
5  *
6  * implementation of telescope simulator helper classes Angle, Axis and Alignment
7  *
8  * This program is free software; you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation; either version 2 of the License, or
11  * (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with this program; if not, write to the Free Software
20  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
21  * MA 02110-1301, USA.
22  */
23 
24 #include "scopesim_helper.h"
25 
26 #include "indilogger.h"
27 
29 
30 // Angle implementation
31 
33 {
34  switch(type)
35  {
36  case DEGREES:
37  angle = range(value);
38  break;
39  case HOURS:
40  angle = range(value * 15.0);
41  break;
42  case RADIANS:
43  angle = range(value * 180.0 / M_PI);
44  break;
45  }
46 }
47 
49 {
50  return this->angle * M_PI / 180.0;
51 }
52 
53 bool Angle::operator== (const Angle &a)
54 {
55  return std::abs(difference(a)) < 10E-6;
56 }
57 
58 bool Angle::operator!= (const Angle &a)
59 {
60  return std::abs(difference(a)) >= 10E-6;
61 }
62 
64 
65 // Axis Implementation
66 
67 void Axis::setDegrees(double degrees)
68 {
69  this->position = degrees;
70  this->target = degrees;
71 }
72 
73 void Axis::setHours(double hours)
74 {
75  this->position = hours * 15.0;
76  this->target = hours * 15.0;
77 }
78 
80 {
81  LOGF_DEBUG("%s StartSlew to %f", axisName, angle.Degrees());
82  target = angle;
83  isSlewing = true;
84 }
85 
86 void Axis::Tracking(bool enabled)
87 {
88  tracking = enabled;
89  LOGF_EXTRA1("%s Teacking enabled %s", axisName, enabled ? "true" : "false");
90 }
91 
93 {
94  trackingRate = rate;
95  switch (trackingRate)
96  {
99  break;
100  case AXIS_TRACK_RATE::SIDEREAL:
101  TrackingRateDegSec = siderealRate;
102  break;
103  case AXIS_TRACK_RATE::SOLAR:
104  TrackingRateDegSec = solarRate;
105  break;
106  case AXIS_TRACK_RATE::LUNAR:
107  TrackingRateDegSec = lunarRate;
108  break;
109  }
110  LOGF_EXTRA1("%s: TrackRate %i, trackingRateDegSec %f arcsec", axisName, trackingRate, TrackingRateDegSec.Degrees() * 3600);
111 }
112 
114 {
115  return trackingRate;
116 }
117 
118 void Axis::StartGuide(double rate, uint32_t durationMs)
119 {
120  // rate is fraction of sidereal, signed to give the direction
121 
122  guideRateDegSec = (360.0 / 86400) * rate;
123  guideDuration = durationMs / 1000.0;
124  LOGF_DEBUG("%s StartGuide rate %f=>%f, dur %d =>%f", axisName, rate, guideRateDegSec.Degrees(), durationMs, guideDuration);
125 }
126 
127 void Axis::update() // called about once a second to update the position and mode
128 {
129  struct timeval currentTime
130  {
131  0, 0
132  };
133 
134  /* update elapsed time since last poll, don't presume exactly POLLMS */
135  gettimeofday(&currentTime, nullptr);
136 
137  if (lastTime.tv_sec == 0 && lastTime.tv_usec == 0)
138  lastTime = currentTime;
139 
140  // Time diff in seconds
141  double interval = currentTime.tv_sec - lastTime.tv_sec + (currentTime.tv_usec - lastTime.tv_usec) / 1e6;
142  lastTime = currentTime;
143  double change = 0;
144 
145  //LOGF_DEBUG("%s: position %f, target %f, interval %f", axisName, position.Degrees(), target.Degrees(), interval);
146 
147  // tracking
148  if (isTracking())
149  {
150  position += TrackingRateDegSec * interval;
151  target += TrackingRateDegSec * interval;
152  //LOGF_EXTRA1("%s: tracking, rate %f, position %f, target %f", axisName, TrackingRateDegSec.Degrees(), position.Degrees(), target.Degrees());
153  }
154 
155  // handle the slew
156  if (isSlewing)
157  {
158  // get positions relative to the rotate centre
159  Angle trc = target - rotateCentre;
160  Angle prc = position - rotateCentre;
161  // get the change, don't use Angle so the change goes through the rotate centre
162  double delta = trc.Degrees() - prc.Degrees();
163  double fastChange = mcRates[4].Degrees() * interval;
164  double slowChange = fastChange / 5;
165  //LOGF_DEBUG("slew: trc %f prc %f, delta %f", trc.Degrees(), prc.Degrees(), delta);
166  // apply the change to the relative position
167  if (delta < -fastChange)
168  {
169  change = -fastChange;
170  }
171  else if (delta < -slowChange)
172  {
173  change = -slowChange;
174  }
175  else if (delta > fastChange)
176  {
177  change = fastChange;
178  }
179  else if (delta > slowChange)
180  {
181  change = slowChange;
182  }
183  else
184  {
185  position = target;
186  isSlewing = false;
187  //OnMoveFinished();
188  }
189  position += change;
190  //LOGF_DEBUG("move %s: change %f, position %f", b, change, position.Degrees());
191  }
192 
193  // handle the motion control
194  if (mcRate < 0)
195  {
196  change = -mcRates[-mcRate].Degrees() * interval;
197  //LOGF_DEBUG("mcRate %d, rate %f, change %f", mcRate, mcRates[-mcRate].Degrees(), change);
198  position += change;
199  }
200  else if (mcRate > 0)
201  {
202  change = mcRates[mcRate].Degrees() * interval;
203  //LOGF_DEBUG("mcRate %d, rate %f, change %f", mcRate, mcRates[mcRate].Degrees(), change);
204  position += change;
205  }
206 
207  // handle guiding
208  if (guideDuration > 0)
209  {
210  change = guideRateDegSec.Degrees() * (guideDuration > interval ? interval : guideDuration);
211  guideDuration -= interval;
212  //LOGF_DEBUG("guide rate %f, remaining duration %f, change %f", guideRateDegSec.Degrees(), guideDuration, change);
213  position += change;
214  }
215 }
216 
218 
219 // Alignment methods
220 
222 {
224 }
225 
226 void Alignment::mountToApparentHaDec(Angle primary, Angle secondary, Angle * apparentHa, Angle* apparentDec)
227 {
228  Angle prio, seco;
229  // get instrument place
230  switch (mountType)
231  {
232  case MOUNT_TYPE::ALTAZ:
233  case MOUNT_TYPE::EQ_FORK:
234  seco = (latitude >= 0) ? secondary : -secondary;
235  prio = primary;
236  break;
237  case MOUNT_TYPE::EQ_GEM:
238  seco = (latitude >= 0) ? secondary : -secondary;
239  prio = primary;
240  if (seco > 90 || seco < -90)
241  {
242  // pointing state inverted
243  seco = Angle(180.0 - seco.Degrees());
244  prio += 180.0;
245  }
246  break;
247  }
248  // instrument to observed, ignore apparent
249  instrumentToObserved(prio, seco, apparentHa, apparentDec);
250  // finally rotate an AltAzm mount to the Ha/Dec coordinates
251  if (mountType == MOUNT_TYPE::ALTAZ)
252  {
253  Angle rot = latitude - Angle(90);
254  Vector haDec = Vector(prio, seco).rotateY(rot);
255  *apparentHa = haDec.primary();
256  *apparentDec = haDec.secondary();
257  LOGF_EXTRA1("m2a Azm Alt %f, %f Ha Dec %f, %f rot %f", prio.Degrees(), seco.Degrees(), apparentHa->Degrees(),
258  apparentDec->Degrees(), rot.Degrees());
259  }
260 }
261 
262 void Alignment::mountToApparentRaDec(Angle primary, Angle secondary, Angle * apparentRa, Angle* apparentDec)
263 {
264  Angle ha;
265  mountToApparentHaDec(primary, secondary, &ha, apparentDec);
266  *apparentRa = lst() - ha;
267  LOGF_EXTRA1("mountToApparentRaDec %f, %f to ha %f, ra %f, %f", primary.Degrees(), secondary.Degrees(), ha.Degrees(),
268  apparentRa->Degrees(), apparentDec->Degrees());
269 }
270 
271 void Alignment::apparentHaDecToMount(Angle apparentHa, Angle apparentDec, Angle* primary, Angle* secondary)
272 {
273  // convert to Alt Azm first
274  if (mountType == MOUNT_TYPE::ALTAZ)
275  {
276  // rotate the apparent HaDec vector to the vertical
277  // TODO sort out Southern Hemisphere
278  Vector altAzm = Vector(apparentHa, apparentDec).rotateY(Angle(90) - latitude);
279  // for now we are making no mount corrections
280  // this all leaves me wondering if the GEM corrections should be done before the mount model
281  *primary = altAzm.primary();
282  *secondary = altAzm.secondary();
283  LOGF_EXTRA1("a2M haDec %f, %f Azm Alt %f, %f", apparentHa.Degrees(), apparentDec.Degrees(), primary->Degrees(),
284  secondary->Degrees() );
285  }
286  // ignore diurnal abberations and refractions to get observed ha, dec
287  // apply telescope pointing to get instrument
288  Angle instrumentHa, instrumentDec;
289  observedToInstrument(apparentHa, apparentDec, &instrumentHa, &instrumentDec);
290 
291  switch (mountType)
292  {
293  case MOUNT_TYPE::ALTAZ:
294  break;
295  case MOUNT_TYPE::EQ_FORK:
296  *secondary = (latitude >= 0) ? instrumentDec : -instrumentDec;
297  *primary = instrumentHa;
298  break;
299  case MOUNT_TYPE::EQ_GEM:
300  *secondary = instrumentDec;
301  *primary = instrumentHa;
302  // use the instrument Ha to select the pointing state
303  if (instrumentHa < flipHourAngle)
304  {
305  // pointing state inverted
306  *primary += Angle(180);
307  *secondary = Angle(180) - instrumentDec;
308  }
309  if (latitude < 0)
310  *secondary = -*secondary;
311  break;
312  }
313 }
314 
315 void Alignment::apparentRaDecToMount(Angle apparentRa, Angle apparentDec, Angle* primary, Angle* secondary)
316 {
317  Angle ha = lst() - apparentRa;
318  apparentHaDecToMount(ha, apparentDec, primary, secondary);
319  LOGF_EXTRA1("apparentRaDecToMount ra %f, ha %f, %f to %f, %f", apparentRa.Degrees(), ha.Degrees(), apparentDec.Degrees(),
320  primary->Degrees(), secondary->Degrees());
321 }
322 
323 #define NEW
324 
325 void Alignment::instrumentToObserved(Angle instrumentHa, Angle instrumentDec, Angle * observedHa, Angle* observedDec)
326 {
327  static Angle maxDec = Angle(89);
328 
329  // do the corrections consecutively
330  // apply Ha and Dec zero offsets
331  *observedHa = instrumentHa + IH;
332  *observedDec = instrumentDec + ID;
333 
334  // apply collimation (cone) error, limited to CH * 10
335  double cosDec = *observedDec < maxDec ? std::cos(observedDec->radians()) : std::cos(maxDec.radians());
336  double tanDec = *observedDec < maxDec ? std::tan(observedDec->radians()) : std::tan(maxDec.radians());
337 
338  *observedHa += (CH / cosDec);
339  // apply Ha and Dec axis non perpendiculary, limited to NP * 10
340  *observedHa += (NP * tanDec);
341 
342  // Use rotations so the corrections work at the pole
343  // apply polar axis Azimuth error using rotation in the East, West, Pole plane (X)
344  Vector vMa = Vector(*observedHa, *observedDec).rotateX(MA);
345 
346  // apply polar axis elevation error using rotation in the North, South, Pole plane (Y)
347  Vector vMe = vMa.rotateY(Angle(ME));
348 
349  *observedHa = vMe.primary();
350  *observedDec = vMe.secondary();
351 }
352 
353 void Alignment::observedToInstrument(Angle observedHa, Angle observedDec, Angle * instrumentHa, Angle *instrumentDec)
354 {
355 #ifdef NEW
356 
357  // avoid errors near dec 90 by limiting sec and tan dec to 57
358  static Angle maxDec = Angle(89);
359 
360  // do the corrections consecutively
361  // use vector rotations for MA and ME so they work close to the pole
362  // rotate about the EW axis (Y)
363  Vector vMe = Vector(observedHa, observedDec).rotateY(Angle(-ME));
364 
365  // apply polar axis Azimuth error
366  Vector vMa = vMe.rotateX(-MA);
367 
368  *instrumentHa = vMa.primary();
369  *instrumentDec = vMa.secondary();
370 
371  // apply Ha and Dec axis non perpendiculary, limited to maxDec
372  double tanDec = *instrumentDec < maxDec ? std::tan(instrumentDec->radians()) : std::tan(maxDec.radians());
373  *instrumentHa -= (NP * tanDec);
374  // apply collimation (cone) error, limited to maxDec
375  double cosDec = *instrumentDec < maxDec ? std::cos(instrumentDec->radians()) : std::cos(maxDec.radians());
376  *instrumentHa -= (CH / cosDec);
377 
378  // apply Ha and Dec zero offsets
379  *instrumentHa -= IH;
380  *instrumentDec -= ID;
381 
382 #else
383 
384 
385  Angle correctionHa, correctionDec;
386  correction(observedHa, observedDec, &correctionHa, &correctionDec);
387  // subtract the correction to get a first pass at the instrument
388  *instrumentHa = observedHa - correctionHa;
389  *instrumentDec = observedDec - correctionDec;
390 
391  static const Angle minDelta(1.0 / 3600.0); // 1.0 arc seconds
392  Angle dHa(180);
393  Angle dDec(180);
394  int num = 0;
395  while ( num < 10)
396  {
397  // use the previously calculated instrument to get a new position
398  Angle newHa, newDec;
399  correction(*instrumentHa, *instrumentDec, &correctionHa, &correctionDec);
400  newHa = *instrumentHa + correctionHa;
401  newDec = *instrumentDec + correctionDec;
402  // get the deltas
403  dHa = observedHa - newHa;
404  dDec = observedDec - newDec;
405  // LOGF_DEBUG("observedToInstrument %i: observed %f, %f, new %f, %f, delta %f, %f", num,
406  // observedHa.Degrees(), observedDec.Degrees(),
407  // newHa.Degrees(), newDec.Degrees(), dHa.Degrees(), dDec.Degrees());
408  if (Angle(fabs(dHa.Degrees())) < minDelta && Angle(fabs(dDec.Degrees())) < minDelta)
409  {
410  return;
411  }
412  // apply the delta
413  *instrumentHa += dHa;
414  *instrumentDec += dDec;
415  num++;
416  }
417  LOGF_DEBUG("mountToObserved iterations %i, delta %f, %f", num, dHa.Degrees(), dDec.Degrees());
418 #endif
419 }
420 
421 // corrections based on the instrument position, add to instrument to get observed
422 // see Patrick Wallace's white paper for details
423 void Alignment::correction(Angle instrumentHa, Angle instrumentDec, Angle * correctionHa, Angle * correctionDec)
424 {
425  // apply Ha and Dec zero offsets
426  *correctionHa = IH;
427  *correctionDec = ID;
428 
429  // avoid errors near dec 90 by limiting sec and tan dec to 100
430  static const double minCos = 0.01;
431  static const double maxTan = 100.0;
432 
433  double cosDec = std::cos(instrumentDec.radians());
434  if (cosDec >= 0 && cosDec < minCos) cosDec = minCos;
435  else if (cosDec <= 0 && cosDec > -minCos) cosDec = -minCos;
436 
437  double tanDec = std::tan(instrumentDec.radians());
438  if (tanDec > maxTan) tanDec = maxTan;
439  if (tanDec < -maxTan) tanDec = -maxTan;
440 
441  double cosHa = std::cos(instrumentHa.radians());
442  double sinHa = std::sin(instrumentHa.radians());
443 
444  // apply collimation (cone) error, limited to CH * 10
445  *correctionHa += (CH / cosDec);
446  // apply Ha and Dec axis non perpendiculary, limited to NP * 10
447  *correctionHa += (NP * tanDec);
448 
449  // apply polar axis Azimuth error
450  *correctionHa += (-MA * cosHa * tanDec);
451  *correctionDec += (MA * sinHa);
452 
453  // apply polar axis elevation error
454  *correctionHa += (ME * sinHa * tanDec);
455  *correctionDec += (ME * cosHa);
456 
457  LOGF_EXTRA1("correction %f, %f", correctionHa->Degrees(), correctionDec->Degrees());
458 }
459 
460 #ifdef FALSE
462 {
463  // rotate by IH and ID
464  Vector vIh = HaDec.rotateZ(Angle(IH));
465  Vector vId = vIh.rotateX(Angle(ID));
466 }
467 #endif
468 
469 void Alignment::setCorrections(double ih, double id, double ch, double np, double ma, double me)
470 {
471  IH = ih;
472  ID = id;
473  CH = ch;
474  NP = np;
475  MA = ma;
476  ME = me;
477  LOGF_DEBUG("setCorrections IH %f, ID %f, CH %f, NP %f, MA %f, ME %f", IH, ID, CH, NP, MA, ME);
478 }
479 
481 
482 // Vector methods
483 
484 Vector::Vector(double x, double y, double z)
485 {
486  double len = std::sqrt(x * x + y * y + z * z);
487  L = x / len;
488  M = y / len;
489  N = z / len;
490 }
491 
492 Vector::Vector(Angle primary, Angle secondary)
493 {
494  double sp = std::sin(primary.radians());
495  double cp = std::cos(primary.radians());
496  double ss = std::sin(secondary.radians());
497  double cs = std::cos(secondary.radians());
498 
499  L = cs * cp;
500  M = cs * sp;
501  N = ss;
502 }
503 
505 {
506  double len = length();
507  L /= len;
508  M /= len;
509  N /= len;
510 }
511 
513 {
514  Angle a = Angle(atan2(M, L), Angle::RADIANS);
515  //IDLog("primary %f", a.Degrees());
516  return a;
517 }
518 
520 {
521  Angle a = Angle(std::asin(N), Angle::RADIANS);
522  //IDLog("secondary %f", a.Degrees());
523  return a;
524 }
525 
527 {
528  double ca = std::cos(angle.radians());
529  double sa = std::sin(angle.radians());
530  return Vector(L, M * ca + N * sa, N * ca - M * sa);
531 }
532 
534 {
535  double ca = std::cos(angle.radians());
536  double sa = std::sin(angle.radians());
537  return Vector(L * ca - N * sa, M, L * sa + N * ca);
538 }
539 
541 {
542  double ca = std::cos(angle.radians());
543  double sa = std::sin(angle.radians());
544  return Vector(L * ca + M * sa, M * ca - L * sa, N);
545 }
double ch()
void setCorrections(double ih, double id, double ch, double np, double ma, double me)
setCorrections set the values of the six corrections
void apparentRaDecToMount(Angle apparentRa, Angle apparentDec, Angle *primary, Angle *secondary)
void mountToApparentHaDec(Angle primary, Angle secondary, Angle *apparentHa, Angle *apparentDec)
mountToApparentHaDec: convert mount position to apparent Ha, Dec
double ih()
double np()
double id()
void apparentHaDecToMount(Angle apparentHa, Angle apparentDec, Angle *primary, Angle *secondary)
apparentHaDecToMount
double ma()
void observedToInstrument(Angle observedHa, Angle observedDec, Angle *instrumentHa, Angle *instrumentDec)
void mountToApparentRaDec(Angle primary, Angle secondary, Angle *apparentRa, Angle *apparentDec)
mountToApparentRaDec: convert mount position to apparent Ra, Dec
void instrumentToObserved(Angle instrumentHa, Angle instrumentDec, Angle *observedHa, Angle *observedDec)
double me()
MOUNT_TYPE mountType
The Angle class This class implements an angle type. This holds an angle that is always in the range ...
double Degrees360()
Degrees360.
double difference(Angle a)
double radians()
radians
bool operator==(const Angle &a)
bool operator!=(const Angle &a)
double Degrees()
Degrees.
Angle TrackingRateDegSec
TrackingRateDegSec.
bool isSlewing
const char * axisName
void setDegrees(double degrees)
void StartGuide(double rate, uint32_t durationMs)
StartGuide start guiding.
void Tracking(bool enabled)
AXIS_TRACK_RATE TrackRate()
TrackRate returns the AXIS_TRACK_RATE property.
void StartSlew(Angle angle)
void update()
Angle position
bool isTracking()
AXIS_TRACK_RATE
The AXIS_TRACK_RATE enum defines the common track rates.
int mcRate
void setHours(double hours)
The Vector class This implements the Directional Cosine used by Taki in his Matrix method....
Vector rotateY(Angle angle)
rotateY rotates this vector through angle about the Y axis
double M
double L
double N
Vector rotateX(Angle angle)
rotateX rotates this vector through angle about the X axis
double length()
Vector()
Vector creates an empty vector.
Vector rotateZ(Angle angle)
rotateZ rotates this vector through angle about the Z axis
void normalise()
Angle secondary()
secondary returns the secondary angle (dec, altitude) from this vector
Angle primary()
primary returns the primary angle (Ra, Ha, Azimuth) from this vector
double get_local_sidereal_time(double longitude)
get_local_sidereal_time Returns local sideral time given longitude and system clock.
#define LOGF_DEBUG(fmt,...)
Definition: indilogger.h:83
#define LOGF_EXTRA1(fmt,...)
Definition: indilogger.h:84
const char * OFF
__le16 type
Definition: pwc-ioctl.h:0