1 // SPDX-License-Identifier: LGPL-2.1-or-later
2 //
3 // SPDX-FileCopyrightText: 2012 Gerhard HOLTKAMP
4 //
5
6 /***************************************************************************
7 * Calculate Spacecraft around other planets *
8 * *
9 * *
10 * Open Source Code. License: GNU LGPL Version 2+ *
11 * *
12 * Author: Gerhard HOLTKAMP, 26-AUG-2012 *
13 ***************************************************************************/
14
15 /*------------ include files and definitions -----------------------------*/
16 #include "planetarySats.h"
17 #include <fstream>
18 #include <iostream>
19 #include <cstdlib>
20 #include <cstring>
21 #include <cmath>
22 #include <ctime>
23 using namespace std;
24
25 #include "astrolib.h"
26
27 // ################ Planetary Sats Class ####################
28
PlanetarySats()29 PlanetarySats::PlanetarySats()
30 {
31 plsatinit();
32 }
33
~PlanetarySats()34 PlanetarySats::~PlanetarySats()
35 {
36
37 }
38
atan23(double y,double x)39 double PlanetarySats::atan23 (double y, double x)
40 {
41 // redefine atan2 so that it doesn't crash when both x and y are 0
42 double result;
43
44 if ((x == 0) && (y == 0)) result = 0;
45 else result = atan2 (y, x);
46
47 return result;
48 }
49
plsatinit()50 void PlanetarySats::plsatinit()
51 {
52 // initialize planetary sat data
53 pls_moonflg = false;
54 pls_day = 1;
55 pls_month = 1;
56 pls_year = 2012;
57 pls_hour = 0;
58 pls_minute = 0;
59 pls_second = 0;
60 pls_del_auto = 1;
61 pls_step = 60.0;
62 pls_delta_rt = 0.0;
63 getTime();
64 getMars();
65
66 strcpy (pls_satelmfl, "./planetarysats.txt");
67
68 }
69
getTime()70 void PlanetarySats::getTime () // Get System Time and Date
71 {
72 time_t tt;
73 int hh, mm, ss;
74 double jd, hr;
75
76 tt = time(NULL);
77 jd = 40587.0 + tt/86400.0; // seconds since 1-JAN-1970
78
79 jd = jd + pls_delta_rt / 24.0;
80 caldat(jd, hh, mm, ss, hr);
81 pls_year = ss;
82 pls_month = mm;
83 pls_day = hh;
84
85 dms(hr, hh, mm, jd);
86 pls_hour = hh;
87 pls_minute = mm;
88 pls_second = int(jd);
89 if (pls_del_auto) pls_del_tdut = DefTdUt(pls_year);
90 setMJD(pls_year, pls_month, pls_day, hh, mm, jd);
91 };
92
setStepWidth(double s)93 void PlanetarySats::setStepWidth(double s)
94 {
95 // set the step width (in seconds) used for calculations
96 if (s < 0.01) pls_step = 0.01;
97 else pls_step = s;
98 }
99
setDeltaRT(double drt)100 void PlanetarySats::setDeltaRT(double drt)
101 {
102 pls_delta_rt = drt;
103 }
104
setDeltaTAI_UTC(double d)105 void PlanetarySats::setDeltaTAI_UTC(double d)
106 {
107 // c is the difference between TAI and UTC according to the IERS
108 // we have to add 32.184 sec to get to the difference TT - UT
109 // which is used in the calculations here
110
111 pls_del_tdut = d + 32.184;
112 pls_del_auto = 0;
113 }
114
setAutoTAI_UTC()115 void PlanetarySats::setAutoTAI_UTC()
116 {
117 // set the difference between TAI and UTC according to the IERS
118 pls_del_auto = true;
119 pls_del_tdut = DefTdUt(pls_year);
120 }
121
setMJD(int year,int month,int day,int hour,int min,double sec)122 void PlanetarySats::setMJD(int year, int month, int day, int hour, int min, double sec)
123 {
124 // set the (MJD-) time currently used for calculations to year, month, day, hour, min, sec
125 double jd;
126
127 pls_year = year;
128 pls_month = month;
129 pls_day = day;
130 pls_hour = hour;
131 pls_minute = min;
132 pls_second = sec;
133
134 jd = ddd(hour, min, sec);
135 jd = mjd(day, month, year, jd);
136
137 pls_time = jd;
138
139 if (pls_del_auto) pls_del_tdut = DefTdUt(pls_year);
140
141 }
142
getDatefromMJD(double mjd,int & year,int & month,int & day,int & hour,int & min,double & sec)143 void PlanetarySats::getDatefromMJD(double mjd, int &year, int &month, int &day, int &hour, int &min, double &sec)
144 {
145 // convert times given in Modified Julian Date (MJD) into conventional date and time
146
147 double magn;
148
149 caldat((mjd), day, month, year, magn);
150 dms (magn,hour,min,sec);
151 if (sec>30.0) min++;
152 if (min>59)
153 {
154 hour++;
155 min=0;
156 };
157 }
158
setSatFile(char * fname)159 void PlanetarySats::setSatFile(char* fname)
160 {
161 strcpy (pls_satelmfl, fname);
162
163 }
164
setStateVector(double mjd,double x,double y,double z,double vx,double vy,double vz)165 void PlanetarySats::setStateVector(double mjd, double x, double y, double z, double vx, double vy, double vz)
166 {
167 pls_rep[0] = x;
168 pls_rep[1] = y;
169 pls_rep[2] = z;
170 pls_vep[0] = vx;
171 pls_vep[1] = vy;
172 pls_vep[2] = vz;
173
174 int year = 0;
175 int month = 0;
176 int day = 0;
177 int hour = 0;
178 int min = 0;
179 double sec = 0;
180 getDatefromMJD(mjd, year, month, day, hour, min, sec);
181 setMJD(year, month, day, hour, min, sec);
182 pls_tepoch = pls_time;
183 //pls_tepoch = pls_time + pls_del_tdut / 86400.0; // epoch in TT
184 }
185
getStateVector(int nsat)186 int PlanetarySats::getStateVector(int nsat)
187 {
188 // read the state vector from the planetary sat file
189 // nsat = number of eligible sat to select (1 if first or only sat)
190 // RETURN number of eligible sat selected, 0 if no suitable sats of file problems
191
192 int fsc, j, k, nst, utc;
193 int yr, month, day, hour, min;
194 double sec;
195 bool searching;
196 ifstream cfle;
197 char satname[40]; // name of satellite
198 char plntname[40]; // name of planet
199
200 strcpy(satname, "");
201 strcpy(plntname, "");
202 nst = 0;
203
204 searching = true;
205 fsc = 0;
206
207 cfle.open(pls_satelmfl, ios::in);
208 if (!cfle)
209 {
210 searching = false;
211 cfle.close();
212 };
213 if(searching)
214 {
215 while (searching)
216 {
217 fsc = 1;
218
219 if (!cfle.getline(satname,40)) fsc = 0;
220 else
221 {
222 k = strlen(satname);
223 if ((k>1) && (satname[0] == '#'))
224 {
225 for (j=1; j<k; ++j)
226 {
227 pls_satname[j-1] = satname[j];
228 if(pls_satname[j-1] == '\n') pls_satname[j-1] = '\0';
229 };
230 pls_satname[k-1] = '\0';
231 }
232 else fsc = 0;
233 };
234
235 if (cfle.eof())
236 {
237 fsc = 0;
238 searching = false;
239 };
240
241 if (fsc)
242 {
243 if (!cfle.getline(plntname, 40)) fsc = 0;
244 else
245 {
246 k = strlen(plntname);
247 if (k>0)
248 {
249 if(plntname[k-1] == '\n') plntname[k-1] = '\0';
250 if((k>1) && (plntname[k-2] == '\r')) plntname[k-2] = '\0';
251 };
252 };
253 };
254
255 if (fsc)
256 {
257 cfle >> yr >> month >> day >> hour >> min >> sec >> utc;
258 if (cfle.bad()) fsc = 0;
259 if (cfle.eof())
260 {
261 fsc = 0;
262 searching = false;
263 };
264 };
265
266 if (fsc)
267 {
268 cfle >> pls_rep[0] >> pls_rep[1] >> pls_rep[2];
269 if (cfle.bad()) fsc = 0;
270 if (cfle.eof())
271 {
272 fsc = 0;
273 searching = false;
274 };
275 };
276
277 if (fsc)
278 {
279 cfle >> pls_vep[0] >> pls_vep[1] >> pls_vep[2];
280 if (cfle.bad()) fsc = 0;
281 };
282
283 if (fsc)
284 {
285 if (strncmp(pls_plntname, plntname, 4) == 0) nst++;
286 if (nst == nsat)
287 {
288 searching = false;
289
290 setMJD(yr, month, day, hour, min, sec);
291 pls_tepoch = pls_time;
292 if (utc) pls_tepoch = pls_tepoch + pls_del_tdut/86400.0; // epoch in TT
293 };
294 };
295 };
296 cfle.close();
297 };
298
299 if (fsc == 0) nst = 0;
300
301 return nst;
302 }
303
setPlanet(char * pname)304 void PlanetarySats::setPlanet(char* pname)
305 {
306 pls_moonflg = false;
307 strcpy(pls_plntname, pname);
308 if (strncmp("Mars", pname, 4) == 0) getMars();
309 if (strncmp("Venus", pname, 4) == 0) getVenus();
310 if (strncmp("Mercury", pname, 4) == 0) getMercury();
311 if (strncmp("Moon", pname, 4) == 0) getMoon();
312 }
313
stateToKepler()314 void PlanetarySats::stateToKepler()
315 {
316 // convert state vector (mean equatorial J2000.0) into planetary Kepler elements
317
318 double t, dt, ag, gm, re, j2;
319 double n, c, w, a, ecc, inc;
320 Vec3 r1, v1;
321 Mat3 mx;
322
323 dt = (pls_tepoch - 51544.5) / 36525.0;
324 gm = pls_GM * 7.4649600000; // convert from m^3/s^2 into km^3/d^2
325 re = pls_R0;
326 j2 = pls_J2;
327
328 // convert into planet equatorial reference frame
329 if (pls_moonflg)
330 {
331 mx = mxidn();
332 r1 = mxvct (mx, pls_rep);
333 v1 = mxvct (mx, pls_vep);
334
335 }
336 else
337 {
338 ag = (pls_axl0 + pls_axl1 * dt) * M_PI / 180.0;
339 mx = zrot(ag + M_PI / 2.0);
340 r1 = mxvct (mx, pls_rep);
341 v1 = mxvct (mx, pls_vep);
342
343 ag = (pls_axb0 + pls_axb1 * dt) * M_PI / 180.0;
344 mx = xrot(M_PI / 2.0 - ag);
345 r1 = mxvct (mx, r1);
346 v1 = mxvct (mx, v1);
347 };
348
349 v1 *= 86400.0; // convert into km / day
350
351 oscelm(gm, pls_tepoch, r1, v1, t, pls_m0, pls_a, pls_ecc, pls_ra, pls_per, pls_inc);
352
353 // now the mean motion
354 a = pls_a;
355 ecc = pls_ecc;
356 inc = pls_inc;
357
358 // preliminary n
359 if (a == 0) a = 1.0; // just in case
360 if (a < 0) a = -a; // just in case
361 n = sqrt (gm / (a*a*a));
362
363 // correct for J2 term
364 w = 1.0 - ecc*ecc;
365 if (w > 0)
366 {
367 w = pow (w, -1.5);
368 c = sin (inc*M_PI/180.0);
369 n = n*(1.0 + 1.5*j2*re*re/(a*a) * w * (1.0 - 1.5*c*c));
370 }
371 else n = 1.0; // do something to avoid a domain error
372
373 n = n / (2.0*M_PI);
374 if (n > 1000.0) n = 1000.0; // avoid possible errors
375
376 pls_n0 = n;
377
378 }
379
getKeplerElements(double & perc,double & apoc,double & inc,double & ecc,double & ra,double & tano,double & m0,double & a,double & n0)380 void PlanetarySats::getKeplerElements(double &perc, double &apoc, double &inc, double &ecc, double &ra, double &tano, double &m0, double &a, double &n0)
381 {
382 // get Kepler elements of orbit with regard to the planetary equator and prime meridian.
383
384 double t, gm;
385 Vec3 r1, v1;
386 Mat3 mx;
387
388
389 if (pls_moonflg) // for the Moon we normally work in J2000. Now get it into planetary
390 {
391 gm = pls_GM * 7.4649600000; // convert from m^3/s^2 into km^3/d^2
392
393 mx = getSelenographic(pls_tepoch);
394 r1 = mxvct (mx, pls_rep);
395 v1 = mxvct (mx, pls_vep);
396 v1 *= 86400.0; // convert into km / day
397
398 oscelm(gm, pls_tepoch, r1, v1, t, m0, a, ecc, ra, tano, inc);
399
400 // now the mean motion
401 if (a == 0) a = 1.0; // just in case
402 if (a < 0) a = -a; // just in case
403 n0 = sqrt (gm / (a*a*a));
404 n0 = n0 / (2.0*M_PI);
405 }
406 else
407 {
408 a = pls_a;
409 n0 = pls_n0;
410 m0 = pls_m0;
411 tano = pls_per;
412 ra = pls_ra;
413 ecc = pls_ecc;
414 inc = pls_inc;
415 };
416
417 perc = pls_a * (1.0 - pls_ecc) - pls_R0;
418 apoc = pls_a * (1.0 + pls_ecc) - pls_R0;
419
420 }
421
selectSat(char * sname)422 int PlanetarySats::selectSat(char* sname)
423 {
424 // select specified satellite
425 // RETURN 1 if successful, 0 if no suitable satellite found
426
427 int nst, res, sl;
428 bool searching;
429
430 searching = true;
431 nst = 1;
432 sl = strlen(sname);
433
434 while (searching)
435 {
436 res = getStateVector(nst);
437 if (res)
438 {
439 if (strncmp(pls_satname, sname, sl) == 0) searching = false;
440 }
441 else searching = false;
442 nst++;
443 };
444
445 return res;
446 }
447
getSatName(char * sname) const448 void PlanetarySats::getSatName(char* sname) const
449 {
450 strcpy (sname, pls_satname);
451 }
452
currentPos()453 void PlanetarySats::currentPos()
454 {
455 getSatPos(pls_time);
456 }
457
nextStep()458 void PlanetarySats::nextStep()
459 {
460 pls_time = pls_time + pls_step / 86400.0;
461 getSatPos(pls_time);
462 }
463
getLastMJD() const464 double PlanetarySats::getLastMJD() const
465 {
466 return pls_time;
467 }
468
getPlanetographic(double & lng,double & lat,double & height) const469 void PlanetarySats::getPlanetographic(double &lng, double &lat, double &height) const
470 {
471 // planetographic coordinates from current state vector
472
473 lng = pls_lng;
474 lat = pls_lat;
475 height = pls_height;
476
477 }
478
getFixedFrame(double & x,double & y,double & z,double & vx,double & vy,double & vz)479 void PlanetarySats::getFixedFrame(double &x, double &y, double &z, double &vx, double &vy, double &vz)
480 {
481 // last state vector coordinates in planetary fixed frame
482
483 x = pls_r[0];
484 y = pls_r[1];
485 z = pls_r[2];
486 vx = pls_v[0];
487 vy = pls_v[1];
488 vz = pls_v[2];
489 }
490
getSatPos(double tutc)491 void PlanetarySats::getSatPos (double tutc)
492 {
493 // Get Position of Satellite at MJD-time t (UTC)
494
495 const double mp2 = 2.0*M_PI;
496
497 double t, dt, m0, ran, aper, inc, a, ecc, n0, re;
498 double f, e, c, s, k, sh, j2, gm, fac, b1, b2, b3;
499 int j;
500
501 Vec3 r1, v1, rg1, s2;
502 Mat3 m1, m2;
503
504 // prepare orbit calculation
505
506 t = tutc + pls_del_tdut/86400.0;
507 dt = t - pls_tepoch;
508
509 ecc = pls_ecc;
510 if (ecc >= 1.0) ecc = 0.999; // to avoid crashes
511 a = pls_a;
512 n0 = mp2 * pls_n0;
513 if (a < 1.0) a = 1.0; // avoid possible crashes later on
514 re = pls_R0;
515 f = pls_flat;
516 j2 = pls_J2;
517 gm = pls_GM * 7.4649600000; // convert from m^3/s^2 into km^3/d^2
518
519 // get current orbit elements ready
520 aper = (re / a) / (1.0 - ecc*ecc);
521 aper = 1.5 * j2 * aper*aper * n0;
522 m0 = pls_inc*M_PI/180.0;
523 ran = -aper * cos(m0) * dt;
524 m0 = sin (m0);
525 aper = aper * (2.0 - 2.5*m0*m0) * dt;
526 ran = pls_ra * M_PI/180.0 + ran;
527 aper = pls_per * M_PI/180.0 + aper;
528 m0 = pls_m0*M_PI/180.0 + n0*dt;
529 inc = pls_inc * M_PI/180.0;
530
531 // solve Kepler equation
532 if (a < 1.0) a = 1.0; // avoid possible crashes later on
533 e = eccanom(m0, ecc);
534 fac = sqrt(1.0 - ecc*ecc);
535 c = cos(e);
536 s = sin(e);
537 r1.assign (a*(c - ecc), a*fac*s, 0.0);
538
539 m0 = 1.0 - ecc*c;
540 k = sqrt(gm/a);
541 v1.assign (-k*s/m0, k*fac*c/m0, 0.0);
542
543 // convert into reference plane
544 m1 = zrot (-aper);
545 m2 = xrot (-inc);
546 m1 *= m2;
547 m2 = zrot (-ran);
548 m2 = m2 * m1;
549 r1 = mxvct (m2, r1);
550 v1 = mxvct (m2, v1);
551 v1 /= 86400.0;
552
553 // save state vector in planet-fixed frame
554
555 if (pls_moonflg) m1 = getSelenographic(t);
556 else m1 = zrot((pls_W + pls_Wd*(t-51544.5))*(M_PI/180.0));
557 pls_r = mxvct(m1,r1);
558 pls_v = mxvct(m1,v1);
559 pls_r *= 1000.0;
560 pls_v *= 1000.0;
561
562 // get the groundtrack coordinates
563
564 rg1 = mxvct(m1,r1);
565
566 s2 = carpol(rg1);
567 pls_lat = s2[2]; // just preliminary
568 pls_lng = s2[1]; // measured with the motion of rotation!
569 if (pls_lng > mp2) pls_lng -= mp2;
570 if (pls_lng < -M_PI) pls_lng += mp2;
571 if (pls_lng > M_PI) pls_lng -= mp2;
572
573 // get height above ellipsoid and geodetic latitude
574 if (abs(r1) > 0.1)
575 {
576 if (f == 0) pls_height = abs(r1) - re;
577 else
578 {
579 c = f*(2.0 - f);
580 s = r1[0]*r1[0] + r1[1]*r1[1];
581 sh = c*r1[2];
582
583 for (j=0; j<4; ++j)
584 {
585 b1 = r1[2] + sh;
586 b3 = sqrt(s+b1*b1);
587 if (b3 < 1e-5) b3 = sin(pls_lat); // just in case
588 else b3 = b1 / b3;
589 b2 = re / sqrt(1.0 - c*b3*b3);
590 sh = b2 * c * b3;
591 };
592
593 sh = r1[2] + sh;
594 pls_lat = atan20(sh,sqrt(s));
595 sh = sqrt(s+sh*sh) - b2;
596 pls_height = sh;
597 }
598 }
599 else pls_height = 0; // this should never happen
600
601 pls_lat = pls_lat * 180.0 / M_PI;
602 pls_lng = pls_lng * 180.0 / M_PI;
603
604 }
605
606
getMars()607 void PlanetarySats::getMars() // Mars planetary constants
608 {
609 pls_J2 = 1.964e-3;
610 pls_R0 = 3397.2;
611 pls_flat = 0.00647630;
612 pls_axl0 = 317.681;
613 pls_axl1 = -0.108;
614 pls_axb0 = 52.886;
615 pls_axb1 = -0.061;
616 pls_W = 176.868;
617 pls_Wd = 350.8919830;
618 pls_GM = 4.282828596416e+13; // 4.282837405582e+13
619 }
620
getVenus()621 void PlanetarySats::getVenus() // Venus planetary constants
622 {
623 pls_J2 = 0.027e-3;
624 pls_R0 = 6051.9;
625 pls_flat = 0.0;
626 pls_axl0 = 272.72;
627 pls_axl1 = 0.0;
628 pls_axb0 = 67.15;
629 pls_axb1 = 0.0;
630 pls_W = 160.26;
631 pls_Wd = -1.4813596;
632 pls_GM = 3.24858761e+14;
633 }
634
getMercury()635 void PlanetarySats::getMercury() // Mercury planetary constants
636 {
637 pls_J2 = 0.0;
638 pls_R0 = 2439.7;
639 pls_flat = 0.0;
640 pls_axl0 = 281.01;
641 pls_axl1 = -0.033;
642 pls_axb0 = 61.45;
643 pls_axb1 = -0.005;
644 pls_W = 329.71;
645 pls_Wd = 6.1385025;
646 pls_GM = 2.20320802e+13;
647 }
648
getMoon()649 void PlanetarySats::getMoon() // Moon planetary constants
650 {
651 pls_moonflg = true;
652 pls_J2 = 0.2027e-3;
653 pls_R0 = 1738.0;
654 pls_flat = 0.0;
655 pls_axl0 = 0.0;
656 pls_axl1 = 0.0;
657 pls_axb0 = 90.0;
658 pls_axb1 = 0.0;
659 pls_W = 0.0;
660 pls_Wd = 13.17635898;
661 pls_GM = 4.90279412e+12;
662 }
663
getSelenographic(double jd)664 Mat3 PlanetarySats::getSelenographic (double jd)
665 {
666 // Calculate the Matrix to transform from Mean of J2000 into selenographic
667 // coordinates at MJD time jd.
668
669 double t, gam, gmp, l, omg, mln;
670 double a, b, c, ic, gn, gp, omp;
671 double const degrad = M_PI / 180.0;
672 Mat3 m1, m2;
673
674 t = (jd - 15019.5) / 36525.0;
675 gam = 281.2208333 + ((0.33333333e-5*t + 0.45277778e-3)*t + 1.7191750)*t;
676 gmp = 334.3295556 + ((-0.125e-4*t - 0.10325e-1)*t + 4069.0340333)*t;
677 l = 279.6966778 + (0.3025e-3*t + 36000.768925)*t;
678 omg = 259.1832750 + ((0.22222222e-5*t + 0.20777778e-2)*t - 1934.1420083)*t;
679 b = 23.45229444 + ((0.50277778e-6*t -0.16388889e-5)*t - 0.130125e-1)*t;
680 mln = 270.4341639 + ((0.1888889e-5*t -0.11333e-2)*t + 481267.8831417)*t;
681 ic = 1.535*degrad;
682 gn = (l - gam)*degrad;
683 gp = (mln - gmp)*degrad;
684 omp = (gmp - omg)*degrad;
685 a = -107.0*cos(gp) + 37.0*cos(gp+2.0*omp) -11.0*cos(2.0*(gp+omp));
686 a = a*0.000277778*degrad + ic;
687 c = (-109.0*sin(gp) + 37.0*sin(gp+2.0*omp) - 11.0*sin(2.0*(gp+omp)))/sin(ic);
688 c = (c*0.000277778 + omg)*degrad;
689 gn = -12.0*sin(gp) + 59.0*sin(gn) + 18.0*sin(2.0*omp);
690 gn = gn*0.000277778*degrad; // tau
691
692 b *= degrad;
693 gam = cos(a)*cos(b) + sin(a)*sin(b)*cos(c);
694 gmp = gam*gam;
695 if(gmp > 1.0) gmp = 0;
696 else gmp = sqrt(1.0 - gmp);
697 gam = atan23(gmp,gam); // theta
698 t = cos(a)*sin(b) - sin(a)*sin(b)*cos(c);
699 l = -sin(a)*sin(c);
700
701 gmp = atan23(l,t); // phi
702 t = sin(a)*cos(b) - cos(a)*sin(b)*cos(c);
703 l = -sin(b)*sin(c);
704 a = atan23(l,t); // delta
705 c = a + mln*degrad + gn - c; // psi
706
707 // libration rotation matrix from Mean equator to true selenographic
708 m1 = zrot(gmp);
709 m2 = xrot(gam);
710 m1 = m2 * m1;
711 m2 = zrot(c);
712 m1 = m2 * m1;
713
714 t = julcent(jd);
715 m2 = pmatequ(0,t); // convert from mean of J2000 to mean of epoch
716 m1 = m1 * m2;
717
718 return m1;
719 }
720
721