Kstars

satellite.cpp
1 /*
2  SPDX-FileCopyrightText: 2009 Jerome SONRIER <[email protected]>
3 
4  SPDX-License-Identifier: GPL-2.0-or-later
5 */
6 
7 #include "satellite.h"
8 
9 #include "ksplanetbase.h"
10 #ifndef KSTARS_LITE
11 #include "kspopupmenu.h"
12 #endif
13 #include "kstarsdata.h"
14 #include "kssun.h"
15 #include "Options.h"
16 #include "skymapcomposite.h"
17 
18 #include <QDebug>
19 
20 #include <cmath>
21 #include <typeinfo>
22 
23 // Define some constants
24 // WGS-72 constants
25 #define RADIUSEARTHKM 6378.135 // Earth radius (km)
26 #define XKE 0.07436691613317 // 60.0 / sqrt(RADIUSEARTHKM^3/MU)
27 #define J2 0.001082616 // The second gravitational zonal harmonic of the Earth
28 #define J4 -0.00000165597 // The fourth gravitational zonal harmonic of the Earth
29 #define J3OJ2 -2.34506972e-3 // J3 / J2
30 
31 // Mathematical constants
32 #define TWOPI 6.2831853071795864769 // 2*PI
33 #define PIO2 1.5707963267948966192 // PI/2
34 #define X2O3 .66666666666666666667 // 2/3
35 #define DEG2RAD 1.745329251994330e-2 // Deg -> Rad
36 
37 // Other constants
38 #define MINPD 1440 // Minutes per day
39 #define MEANALT 0.84 // Mean altitude (km)
40 #define SR 6.96000e5 // Solar radius - km (IAU 76)
41 #define AU 1.49597870691e8 // Astronomical unit - km (IAU 76)
42 #define XPDOTP 229.1831180523293 // 1440.0 / (2.0 * pi)
43 #define SS 1.0122292801892716288 // Parameter for the SGP4 density function
44 #define QZMS2T 1.8802791590152706439e-9 // (( 120.0 - 78.0) / RADIUSEARTHKM )^4
45 #define F 3.35281066474748e-3 // Flattening factor
46 #define MFACTOR 7.292115e-5
47 
48 Satellite::Satellite(const QString &name, const QString &line1, const QString &line2)
49 {
50  //m_name = name;
51  m_number = line1.midRef(2, 5).toInt();
52  m_class = line1.at(7);
53  m_id = line1.mid(9, 8);
54  m_epoch = line1.midRef(18, 14).toDouble();
55  m_first_deriv = line1.midRef(33, 10).toDouble() / (XPDOTP * MINPD);
56  m_second_deriv =
57  line1.midRef(44, 6).toDouble() * (1.0e-5 / pow(10.0, line1.midRef(51, 1).toDouble())) / (XPDOTP * MINPD * MINPD);
58  m_bstar = line1.midRef(53, 6).toDouble() * 1.0e-5 / pow(10.0, line1.midRef(60, 1).toDouble());
59  m_ephem_type = line1.midRef(62, 1).toInt();
60  m_elem_number = line1.midRef(64, 4).toInt();
61  m_inclination = line2.midRef(8, 8).toDouble() * DEG2RAD;
62  m_ra = line2.midRef(17, 8).toDouble() * DEG2RAD;
63  m_eccentricity = line2.midRef(26, 7).toDouble() * 1.0e-7;
64  m_arg_perigee = line2.midRef(34, 8).toDouble() * DEG2RAD;
65  m_mean_anomaly = line2.midRef(43, 8).toDouble() * DEG2RAD;
66  m_mean_motion = line2.midRef(52, 11).toDouble() * TWOPI / MINPD;
67  m_nb_revolution = line2.midRef(63, 5).toInt();
68  m_tle = name + "\n" + line1 + "\n" + line2;
69 
70  setName(name);
71  setName2(name);
72  setLongName(name + " (" + m_id + ')');
73  setType(SkyObject::SATELLITE);
74  setMag(0.0);
75 
76  m_is_selected = Options::selectedSatellites().contains(name);
77 
78  // Convert TLE epoch to Julian date
79  double day = modf(m_epoch * 1.e-3, &m_epoch_year) * 1.e3;
80  if (m_epoch_year < 57.)
81  m_epoch_year += 2000.;
82  else
83  m_epoch_year += 1900.;
84  double year = m_epoch_year - 1.;
85  long i = year / 100;
86  long A = i;
87  i = A / 4;
88  long B = 2 - A + i;
89  i = 365.25 * year;
90  i += 30.6001 * 14;
91  m_tle_jd = i + 1720994.5 + B + day;
92 
93  init();
94 }
95 
97 {
98  Q_ASSERT(typeid(this) == typeid(static_cast<const Satellite *>(this))); // Ensure we are not slicing a derived class
99  return new Satellite(*this);
100 }
101 
102 void Satellite::init()
103 {
104  double ao, cosio, sinio, cosio2, omeosq, posq, rp, rteosq, eccsq, con42, cnodm, snodm, cosim, sinim, cosomm, sinomm,
105  cc1sq, cc2, cc3, coef, coef1, cosio4, day, em, emsq, eeta, etasq, gam, inclm, nm, perige, pinvsq, psisq,
106  qzms24, rtemsq, s1, s2, s3, s4, s5, s6, s7, sfour, ss1(0), ss2(0), ss3(0), ss4(0), ss5(0), ss6(0), ss7(0),
107  sz1(0), sz2(0), sz3(0), sz11(0), sz12(0), sz13(0), sz21(0), sz22(0), sz23(0), sz31(0), sz32(0), sz33(0), tc,
108  temp, temp1, temp2, temp3, tsi, xpidot, xhdot1, z1, z2, z3, z11, z12, z13, z21, z22, z23, z31, z32, z33, ak, d1,
109  del, adel, po, ds70, ts70, tfrac, c1, thgr70, fk5r, c1p2p;
110  // double dndt;
111 
112  // Init near earth variables
113  isimp = false;
114  aycof = 0.;
115  con41 = 0.;
116  cc1 = 0.;
117  cc4 = 0.;
118  cc5 = 0.;
119  d2 = 0.;
120  d3 = 0.;
121  d4 = 0.;
122  delmo = 0.;
123  eta = 0.;
124  argpdot = 0.;
125  omgcof = 0.;
126  sinmao = 0.;
127  t = 0.;
128  t2cof = 0.;
129  t3cof = 0.;
130  t4cof = 0.;
131  t5cof = 0.;
132  x1mth2 = 0.;
133  x7thm1 = 0.;
134  mdot = 0.;
135  nodedot = 0.;
136  xlcof = 0.;
137  xmcof = 0.;
138  nodecf = 0.;
139 
140  // Init deep space variables
141  irez = 0;
142  d2201 = 0.;
143  d2211 = 0.;
144  d3210 = 0.;
145  d3222 = 0.;
146  d4410 = 0.;
147  d4422 = 0.;
148  d5220 = 0.;
149  d5232 = 0.;
150  d5421 = 0.;
151  d5433 = 0.;
152  dedt = 0.;
153  del1 = 0.;
154  del2 = 0.;
155  del3 = 0.;
156  didt = 0.;
157  dmdt = 0.;
158  dnodt = 0.;
159  domdt = 0.;
160  e3 = 0.;
161  ee2 = 0.;
162  peo = 0.;
163  pgho = 0.;
164  pho = 0.;
165  pinco = 0.;
166  plo = 0.;
167  se2 = 0.;
168  se3 = 0.;
169  sgh2 = 0.;
170  sgh3 = 0.;
171  sgh4 = 0.;
172  sh2 = 0.;
173  sh3 = 0.;
174  si2 = 0.;
175  si3 = 0.;
176  sl2 = 0.;
177  sl3 = 0.;
178  sl4 = 0.;
179  gsto = 0.;
180  xfact = 0.;
181  xgh2 = 0.;
182  xgh3 = 0.;
183  xgh4 = 0.;
184  xh2 = 0.;
185  xh3 = 0.;
186  xi2 = 0.;
187  xi3 = 0.;
188  xl2 = 0.;
189  xl3 = 0.;
190  xl4 = 0.;
191  xlamo = 0.;
192  zmol = 0.;
193  zmos = 0.;
194  atime = 0.;
195  xli = 0.;
196  xni = 0.;
197 
198  method = 'n';
199 
200  m_is_visible = false;
201 
202  // Divisor for divide by zero check on inclination
203  const double temp4 = 1.5e-12;
204 
205  /*----- Initializes variables for sgp4 -----*/
206  // Calculate auxiliary epoch quantities
207  eccsq = m_eccentricity * m_eccentricity;
208  omeosq = 1.0 - eccsq;
209  rteosq = sqrt(omeosq);
210  cosio = cos(m_inclination);
211  cosio2 = cosio * cosio;
212 
213  // Un-kozai the mean motion
214  ak = pow(XKE / m_mean_motion, X2O3);
215  d1 = 0.75 * J2 * (3.0 * cosio2 - 1.0) / (rteosq * omeosq);
216  del = d1 / (ak * ak);
217  adel = ak * (1.0 - del * del - del * (1.0 / 3.0 + 134.0 * del * del / 81.0));
218  del = d1 / (adel * adel);
219  m_mean_motion = m_mean_motion / (1.0 + del);
220 
221  ao = pow(XKE / m_mean_motion, X2O3);
222  sinio = sin(m_inclination);
223  po = ao * omeosq;
224  con42 = 1.0 - 5.0 * cosio2;
225  con41 = -con42 - (2.0 * cosio2);
226  posq = po * po;
227  rp = ao * (1.0 - m_eccentricity);
228  method = 'n';
229 
230  // Find sidereal time
231  ts70 = m_tle_jd - 2433281.5 - 7305.0;
232  ds70 = floor(ts70 + 1.0e-8);
233  tfrac = ts70 - ds70;
234  // find greenwich location at epoch
235  c1 = 1.72027916940703639e-2;
236  thgr70 = 1.7321343856509374;
237  fk5r = 5.07551419432269442e-15;
238  c1p2p = c1 + TWOPI;
239  gsto = fmod(thgr70 + c1 * ds70 + c1p2p * tfrac + ts70 * ts70 * fk5r, TWOPI);
240  if (gsto < 0.0)
241  gsto = gsto + TWOPI;
242 
243  if ((omeosq >= 0.0) || (m_mean_motion >= 0.0))
244  {
245  if (rp < (220.0 / RADIUSEARTHKM + 1.0))
246  isimp = true;
247  sfour = SS;
248  qzms24 = QZMS2T;
249  perige = (rp - 1.0) * RADIUSEARTHKM;
250 
251  // For perigees below 156 km, s and qoms2t are altered
252  if (perige < 156.0)
253  {
254  sfour = perige - 78.0;
255  if (perige < 98.0)
256  sfour = 20.0;
257  qzms24 = pow(((120.0 - sfour) / RADIUSEARTHKM), 4.0);
258  sfour = sfour / RADIUSEARTHKM + 1.0;
259  }
260  pinvsq = 1.0 / posq;
261 
262  tsi = 1.0 / (ao - sfour);
263  eta = ao * m_eccentricity * tsi;
264  etasq = eta * eta;
265  eeta = m_eccentricity * eta;
266  psisq = fabs(1.0 - etasq);
267  coef = qzms24 * pow(tsi, 4.0);
268  coef1 = coef / pow(psisq, 3.5);
269  cc2 = coef1 * m_mean_motion *
270  (ao * (1.0 + 1.5 * etasq + eeta * (4.0 + etasq)) +
271  0.375 * J2 * tsi / psisq * con41 * (8.0 + 3.0 * etasq * (8.0 + etasq)));
272  cc1 = m_bstar * cc2;
273  cc3 = 0.0;
274  if (m_eccentricity > 1.0e-4)
275  cc3 = -2.0 * coef * tsi * J3OJ2 * m_mean_motion * sinio / m_eccentricity;
276  x1mth2 = 1.0 - cosio2;
277  cc4 = 2.0 * m_mean_motion * coef1 * ao * omeosq *
278  (eta * (2.0 + 0.5 * etasq) + m_eccentricity * (0.5 + 2.0 * etasq) -
279  J2 * tsi / (ao * psisq) *
280  (-3.0 * con41 * (1.0 - 2.0 * eeta + etasq * (1.5 - 0.5 * eeta)) +
281  0.75 * x1mth2 * (2.0 * etasq - eeta * (1.0 + etasq)) * cos(2.0 * m_arg_perigee)));
282  cc5 = 2.0 * coef1 * ao * omeosq * (1.0 + 2.75 * (etasq + eeta) + eeta * etasq);
283 
284  cosio4 = cosio2 * cosio2;
285  temp1 = 1.5 * J2 * pinvsq * m_mean_motion;
286  temp2 = 0.5 * temp1 * J2 * pinvsq;
287  temp3 = -0.46875 * J4 * pinvsq * pinvsq * m_mean_motion;
288  mdot = m_mean_motion + 0.5 * temp1 * rteosq * con41 +
289  0.0625 * temp2 * rteosq * (13.0 - 78.0 * cosio2 + 137.0 * cosio4);
290  argpdot = -0.5 * temp1 * con42 + 0.0625 * temp2 * (7.0 - 114.0 * cosio2 + 395.0 * cosio4) +
291  temp3 * (3.0 - 36.0 * cosio2 + 49.0 * cosio4);
292  xhdot1 = -temp1 * cosio;
293  nodedot = xhdot1 + (0.5 * temp2 * (4.0 - 19.0 * cosio2) + 2.0 * temp3 * (3.0 - 7.0 * cosio2)) * cosio;
294  xpidot = argpdot + nodedot;
295  omgcof = m_bstar * cc3 * cos(m_arg_perigee);
296  xmcof = 0.0;
297  if (m_eccentricity > 1.0e-4)
298  xmcof = -X2O3 * coef * m_bstar / eeta;
299  nodecf = 3.5 * omeosq * xhdot1 * cc1;
300  t2cof = 1.5 * cc1;
301  // Do not divide by zero
302  if (fabs(1.0 + cosio) > 1.5e-12)
303  xlcof = -0.25 * J3OJ2 * sinio * (3.0 + 5.0 * cosio) / (1.0 + cosio);
304  else
305  xlcof = -0.25 * J3OJ2 * sinio * (3.0 + 5.0 * cosio) / temp4;
306  aycof = -0.5 * J3OJ2 * sinio;
307  delmo = pow((1.0 + eta * cos(m_mean_anomaly)), 3);
308  sinmao = sin(m_mean_anomaly);
309  x7thm1 = 7.0 * cosio2 - 1.0;
310 
311  // Deep space initialization
312  if ((TWOPI / m_mean_motion) >= 225.0)
313  {
314  method = 'd';
315  isimp = true;
316  tc = 0.0;
317  inclm = m_inclination;
318 
319  // Init deep space common variables
320  // Define some constants
321  const double zes = 0.01675;
322  const double zel = 0.05490;
323  const double c1ss = 2.9864797e-6;
324  const double c1l = 4.7968065e-7;
325  const double zsinis = 0.39785416;
326  const double zcosis = 0.91744867;
327  const double zcosgs = 0.1945905;
328  const double zsings = -0.98088458;
329 
330  int lsflg;
331  double a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, betasq, cc, ctem, stem, x1, x2, x3, x4, x5, x6, x7, x8,
332  xnodce, xnoi, zcosg, zcosgl, zcosh, zcoshl, zcosi, zcosil, zsing, zsingl, zsinh, zsinhl, zsini, zsinil,
333  zx, zy;
334 
335  nm = m_mean_motion;
336  em = m_eccentricity;
337  snodm = sin(m_ra);
338  cnodm = cos(m_ra);
339  sinomm = sin(m_arg_perigee);
340  cosomm = cos(m_arg_perigee);
341  sinim = sin(m_inclination);
342  cosim = cos(m_inclination);
343  emsq = em * em;
344  betasq = 1.0 - emsq;
345  rtemsq = sqrt(betasq);
346 
347  // Initialize lunar solar terms
348  peo = 0.0;
349  pinco = 0.0;
350  plo = 0.0;
351  pgho = 0.0;
352  pho = 0.0;
353  day = m_tle_jd - 2433281.5 + 18261.5 + tc / 1440.0;
354  xnodce = fmod(4.5236020 - 9.2422029e-4 * day, TWOPI);
355  stem = sin(xnodce);
356  ctem = cos(xnodce);
357  zcosil = 0.91375164 - 0.03568096 * ctem;
358  zsinil = sqrt(1.0 - zcosil * zcosil);
359  zsinhl = 0.089683511 * stem / zsinil;
360  zcoshl = sqrt(1.0 - zsinhl * zsinhl);
361  gam = 5.8351514 + 0.0019443680 * day;
362  zx = 0.39785416 * stem / zsinil;
363  zy = zcoshl * ctem + 0.91744867 * zsinhl * stem;
364  zx = atan2(zx, zy);
365  zx = gam + zx - xnodce;
366  zcosgl = cos(zx);
367  zsingl = sin(zx);
368 
369  // Solar terms
370  zcosg = zcosgs;
371  zsing = zsings;
372  zcosi = zcosis;
373  zsini = zsinis;
374  zcosh = cnodm;
375  zsinh = snodm;
376  cc = c1ss;
377  xnoi = 1.0 / nm;
378 
379  for (lsflg = 1; lsflg <= 2; ++lsflg)
380  {
381  a1 = zcosg * zcosh + zsing * zcosi * zsinh;
382  a3 = -zsing * zcosh + zcosg * zcosi * zsinh;
383  a7 = -zcosg * zsinh + zsing * zcosi * zcosh;
384  a8 = zsing * zsini;
385  a9 = zsing * zsinh + zcosg * zcosi * zcosh;
386  a10 = zcosg * zsini;
387  a2 = cosim * a7 + sinim * a8;
388  a4 = cosim * a9 + sinim * a10;
389  a5 = -sinim * a7 + cosim * a8;
390  a6 = -sinim * a9 + cosim * a10;
391 
392  x1 = a1 * cosomm + a2 * sinomm;
393  x2 = a3 * cosomm + a4 * sinomm;
394  x3 = -a1 * sinomm + a2 * cosomm;
395  x4 = -a3 * sinomm + a4 * cosomm;
396  x5 = a5 * sinomm;
397  x6 = a6 * sinomm;
398  x7 = a5 * cosomm;
399  x8 = a6 * cosomm;
400 
401  z31 = 12.0 * x1 * x1 - 3.0 * x3 * x3;
402  z32 = 24.0 * x1 * x2 - 6.0 * x3 * x4;
403  z33 = 12.0 * x2 * x2 - 3.0 * x4 * x4;
404  z1 = 3.0 * (a1 * a1 + a2 * a2) + z31 * emsq;
405  z2 = 6.0 * (a1 * a3 + a2 * a4) + z32 * emsq;
406  z3 = 3.0 * (a3 * a3 + a4 * a4) + z33 * emsq;
407  z11 = -6.0 * a1 * a5 + emsq * (-24.0 * x1 * x7 - 6.0 * x3 * x5);
408  z12 = -6.0 * (a1 * a6 + a3 * a5) + emsq * (-24.0 * (x2 * x7 + x1 * x8) - 6.0 * (x3 * x6 + x4 * x5));
409  z13 = -6.0 * a3 * a6 + emsq * (-24.0 * x2 * x8 - 6.0 * x4 * x6);
410  z21 = 6.0 * a2 * a5 + emsq * (24.0 * x1 * x5 - 6.0 * x3 * x7);
411  z22 = 6.0 * (a4 * a5 + a2 * a6) + emsq * (24.0 * (x2 * x5 + x1 * x6) - 6.0 * (x4 * x7 + x3 * x8));
412  z23 = 6.0 * a4 * a6 + emsq * (24.0 * x2 * x6 - 6.0 * x4 * x8);
413  z1 = z1 + z1 + betasq * z31;
414  z2 = z2 + z2 + betasq * z32;
415  z3 = z3 + z3 + betasq * z33;
416  s3 = cc * xnoi;
417  s2 = -0.5 * s3 / rtemsq;
418  s4 = s3 * rtemsq;
419  s1 = -15.0 * em * s4;
420  s5 = x1 * x3 + x2 * x4;
421  s6 = x2 * x3 + x1 * x4;
422  s7 = x2 * x4 - x1 * x3;
423 
424  // Lunar terms
425  if (lsflg == 1)
426  {
427  ss1 = s1;
428  ss2 = s2;
429  ss3 = s3;
430  ss4 = s4;
431  ss5 = s5;
432  ss6 = s6;
433  ss7 = s7;
434  sz1 = z1;
435  sz2 = z2;
436  sz3 = z3;
437  sz11 = z11;
438  sz12 = z12;
439  sz13 = z13;
440  sz21 = z21;
441  sz22 = z22;
442  sz23 = z23;
443  sz31 = z31;
444  sz32 = z32;
445  sz33 = z33;
446  zcosg = zcosgl;
447  zsing = zsingl;
448  zcosi = zcosil;
449  zsini = zsinil;
450  zcosh = zcoshl * cnodm + zsinhl * snodm;
451  zsinh = snodm * zcoshl - cnodm * zsinhl;
452  cc = c1l;
453  }
454  }
455 
456  zmol = fmod(4.7199672 + 0.22997150 * day - gam, TWOPI);
457  zmos = fmod(6.2565837 + 0.017201977 * day, TWOPI);
458 
459  // Solar terms
460  se2 = 2.0 * ss1 * ss6;
461  se3 = 2.0 * ss1 * ss7;
462  si2 = 2.0 * ss2 * sz12;
463  si3 = 2.0 * ss2 * (sz13 - sz11);
464  sl2 = -2.0 * ss3 * sz2;
465  sl3 = -2.0 * ss3 * (sz3 - sz1);
466  sl4 = -2.0 * ss3 * (-21.0 - 9.0 * emsq) * zes;
467  sgh2 = 2.0 * ss4 * sz32;
468  sgh3 = 2.0 * ss4 * (sz33 - sz31);
469  sgh4 = -18.0 * ss4 * zes;
470  sh2 = -2.0 * ss2 * sz22;
471  sh3 = -2.0 * ss2 * (sz23 - sz21);
472 
473  // Lunar terms
474  ee2 = 2.0 * s1 * s6;
475  e3 = 2.0 * s1 * s7;
476  xi2 = 2.0 * s2 * z12;
477  xi3 = 2.0 * s2 * (z13 - z11);
478  xl2 = -2.0 * s3 * z2;
479  xl3 = -2.0 * s3 * (z3 - z1);
480  xl4 = -2.0 * s3 * (-21.0 - 9.0 * emsq) * zel;
481  xgh2 = 2.0 * s4 * z32;
482  xgh3 = 2.0 * s4 * (z33 - z31);
483  xgh4 = -18.0 * s4 * zel;
484  xh2 = -2.0 * s2 * z22;
485  xh3 = -2.0 * s2 * (z23 - z21);
486 
487  // Apply deep space long period periodic contributions to the mean elements
488  // double f2, f3, sinzf, zf, zm;
489  double ses, sghl, sghs, shll, shs, sis, sls;
490 
491  // Define some constants
492  const double zns = 1.19459e-5;
493  const double znl = 1.5835218e-4;
494 
495  // Calculate time varying periodics
496  // These results are never used, should we remove them?
497  // zm = zmos;
498  // zf = zm + 2.0 * zes * sin(zm);
499  // sinzf = sin(zf);
500  // f2 = 0.5 * sinzf * sinzf - 0.25;
501  // f3 = -0.5 * sinzf * cos(zf);
502  // ses = se2 * f2 + se3 * f3;
503  // sis = si2 * f2 + si3 * f3;
504  // sls = sl2 * f2 + sl3 * f3 + sl4 * sinzf;
505  // sghs = sgh2 * f2 + sgh3 * f3 + sgh4 * sinzf;
506  // shs = sh2 * f2 + sh3 * f3;
507  // zm = zmol;
508 
509  // zf = zm + 2.0 * zel * sin(zm);
510  // sinzf = sin(zf);
511  // f2 = 0.5 * sinzf * sinzf - 0.25;
512  // f3 = -0.5 * sinzf * cos(zf);
513  // sghl = xgh2 * f2 + xgh3 * f3 + xgh4 * sinzf;
514  // shll = xh2 * f2 + xh3 * f3;
515 
516  // Deep space contributions to mean motion dot due to geopotential resonance with half day and one day orbits
517  double ainv2, aonv = 0.0, cosisq, eoc, f220, f221, f311, f321, f322, f330, f441, f442, f522, f523, f542,
518  f543, g200, g201, g211, g300, g310, g322, g410, g422, g520, g521, g532, g533, sgs, sini2,
519  temp, temp1, theta, xno2, emsqo;
520  // double emo;
521 
522  // Define some constant
523  const double q22 = 1.7891679e-6;
524  const double q31 = 2.1460748e-6;
525  const double q33 = 2.2123015e-7;
526  const double root22 = 1.7891679e-6;
527  const double root44 = 7.3636953e-9;
528  const double root54 = 2.1765803e-9;
529  const double rptim = 4.37526908801129966e-3; // this equates to 7.29211514668855e-5 rad/sec
530  const double root32 = 3.7393792e-7;
531  const double root52 = 1.1428639e-7;
532 
533  // Deep space initialization
534  irez = 0;
535  if ((nm < 0.0052359877) && (nm > 0.0034906585))
536  irez = 1;
537  if ((nm >= 8.26e-3) && (nm <= 9.24e-3) && (em >= 0.5))
538  irez = 2;
539 
540  // Solar terms
541  ses = ss1 * zns * ss5;
542  sis = ss2 * zns * (sz11 + sz13);
543  sls = -zns * ss3 * (sz1 + sz3 - 14.0 - 6.0 * emsq);
544  sghs = ss4 * zns * (sz31 + sz33 - 6.0);
545  shs = -zns * ss2 * (sz21 + sz23);
546  if ((inclm < 5.2359877e-2) || (inclm > M_PI - 5.2359877e-2))
547  shs = 0.0;
548  if (sinim != 0.0)
549  shs = shs / sinim;
550  sgs = sghs - cosim * shs;
551 
552  // Lunar terms
553  dedt = ses + s1 * znl * s5;
554  didt = sis + s2 * znl * (z11 + z13);
555  dmdt = sls - znl * s3 * (z1 + z3 - 14.0 - 6.0 * emsq);
556  sghl = s4 * znl * (z31 + z33 - 6.0);
557  shll = -znl * s2 * (z21 + z23);
558  if ((inclm < 5.2359877e-2) || (inclm > M_PI - 5.2359877e-2))
559  shll = 0.0;
560  domdt = sgs + sghl;
561  dnodt = shs;
562  if (sinim != 0.0)
563  {
564  domdt = domdt - cosim / sinim * shll;
565  dnodt = dnodt + shll / sinim;
566  }
567 
568  // Calculate deep space resonance effects
569  // Value never used
570  // dndt = 0.0;
571  theta = fmod(gsto + tc * rptim, TWOPI);
572 
573  // Initialize the resonance terms
574  if (irez != 0)
575  {
576  aonv = pow(nm / XKE, X2O3);
577 
578  // Geopotential resonance for 12 hour orbits
579  if (irez == 2)
580  {
581  cosisq = cosim * cosim;
582  // Value never used
583  // emo = em;
584  em = m_eccentricity;
585  emsqo = emsq;
586  emsq = eccsq;
587  eoc = em * emsq;
588  g201 = -0.306 - (em - 0.64) * 0.440;
589 
590  if (em <= 0.65)
591  {
592  g211 = 3.616 - 13.2470 * em + 16.2900 * emsq;
593  g310 = -19.302 + 117.3900 * em - 228.4190 * emsq + 156.5910 * eoc;
594  g322 = -18.9068 + 109.7927 * em - 214.6334 * emsq + 146.5816 * eoc;
595  g410 = -41.122 + 242.6940 * em - 471.0940 * emsq + 313.9530 * eoc;
596  g422 = -146.407 + 841.8800 * em - 1629.014 * emsq + 1083.4350 * eoc;
597  g520 = -532.114 + 3017.977 * em - 5740.032 * emsq + 3708.2760 * eoc;
598  }
599  else
600  {
601  g211 = -72.099 + 331.819 * em - 508.738 * emsq + 266.724 * eoc;
602  g310 = -346.844 + 1582.851 * em - 2415.925 * emsq + 1246.113 * eoc;
603  g322 = -342.585 + 1554.908 * em - 2366.899 * emsq + 1215.972 * eoc;
604  g410 = -1052.797 + 4758.686 * em - 7193.992 * emsq + 3651.957 * eoc;
605  g422 = -3581.690 + 16178.110 * em - 24462.770 * emsq + 12422.520 * eoc;
606  if (em > 0.715)
607  g520 = -5149.66 + 29936.92 * em - 54087.36 * emsq + 31324.56 * eoc;
608  else
609  g520 = 1464.74 - 4664.75 * em + 3763.64 * emsq;
610  }
611  if (em < 0.7)
612  {
613  g533 = -919.22770 + 4988.6100 * em - 9064.7700 * emsq + 5542.21 * eoc;
614  g521 = -822.71072 + 4568.6173 * em - 8491.4146 * emsq + 5337.524 * eoc;
615  g532 = -853.66600 + 4690.2500 * em - 8624.7700 * emsq + 5341.4 * eoc;
616  }
617  else
618  {
619  g533 = -37995.780 + 161616.52 * em - 229838.20 * emsq + 109377.94 * eoc;
620  g521 = -51752.104 + 218913.95 * em - 309468.16 * emsq + 146349.42 * eoc;
621  g532 = -40023.880 + 170470.89 * em - 242699.48 * emsq + 115605.82 * eoc;
622  }
623 
624  sini2 = sinim * sinim;
625  f220 = 0.75 * (1.0 + 2.0 * cosim + cosisq);
626  f221 = 1.5 * sini2;
627  f321 = 1.875 * sinim * (1.0 - 2.0 * cosim - 3.0 * cosisq);
628  f322 = -1.875 * sinim * (1.0 + 2.0 * cosim - 3.0 * cosisq);
629  f441 = 35.0 * sini2 * f220;
630  f442 = 39.3750 * sini2 * sini2;
631  f522 =
632  9.84375 * sinim *
633  (sini2 * (1.0 - 2.0 * cosim - 5.0 * cosisq) + 0.33333333 * (-2.0 + 4.0 * cosim + 6.0 * cosisq));
634  f523 = sinim * (4.92187512 * sini2 * (-2.0 - 4.0 * cosim + 10.0 * cosisq) +
635  6.56250012 * (1.0 + 2.0 * cosim - 3.0 * cosisq));
636  f542 = 29.53125 * sinim * (2.0 - 8.0 * cosim + cosisq * (-12.0 + 8.0 * cosim + 10.0 * cosisq));
637  f543 = 29.53125 * sinim * (-2.0 - 8.0 * cosim + cosisq * (12.0 + 8.0 * cosim - 10.0 * cosisq));
638  xno2 = nm * nm;
639  ainv2 = aonv * aonv;
640  temp1 = 3.0 * xno2 * ainv2;
641  temp = temp1 * root22;
642  d2201 = temp * f220 * g201;
643  d2211 = temp * f221 * g211;
644  temp1 = temp1 * aonv;
645  temp = temp1 * root32;
646  d3210 = temp * f321 * g310;
647  d3222 = temp * f322 * g322;
648  temp1 = temp1 * aonv;
649  temp = 2.0 * temp1 * root44;
650  d4410 = temp * f441 * g410;
651  d4422 = temp * f442 * g422;
652  temp1 = temp1 * aonv;
653  temp = temp1 * root52;
654  d5220 = temp * f522 * g520;
655  d5232 = temp * f523 * g532;
656  temp = 2.0 * temp1 * root54;
657  d5421 = temp * f542 * g521;
658  d5433 = temp * f543 * g533;
659  xlamo = fmod(m_mean_anomaly + m_ra + m_ra - theta - theta, TWOPI);
660  xfact = mdot + dmdt + 2.0 * (nodedot + dnodt - rptim) - m_mean_motion;
661  // Value never used
662  // em = emo;
663  emsq = emsqo;
664  }
665 
666  if (irez == 1)
667  {
668  g200 = 1.0 + emsq * (-2.5 + 0.8125 * emsq);
669  g310 = 1.0 + 2.0 * emsq;
670  g300 = 1.0 + emsq * (-6.0 + 6.60937 * emsq);
671  f220 = 0.75 * (1.0 + cosim) * (1.0 + cosim);
672  f311 = 0.9375 * sinim * sinim * (1.0 + 3.0 * cosim) - 0.75 * (1.0 + cosim);
673  f330 = 1.0 + cosim;
674  f330 = 1.875 * f330 * f330 * f330;
675  del1 = 3.0 * nm * nm * aonv * aonv;
676  del2 = 2.0 * del1 * f220 * g200 * q22;
677  del3 = 3.0 * del1 * f330 * g300 * q33 * aonv;
678  del1 = del1 * f311 * g310 * q31 * aonv;
679  xlamo = fmod(m_mean_anomaly + m_ra + m_arg_perigee - theta, TWOPI);
680  xfact = mdot + xpidot - rptim + dmdt + domdt + dnodt - m_mean_motion;
681  }
682 
683  xli = xlamo;
684  xni = m_mean_motion;
685  atime = 0.0;
686  // Value never used
687  // nm = m_mean_motion + dndt;
688  }
689  }
690 
691  // Set variables if not deep space
692  if (!isimp)
693  {
694  cc1sq = cc1 * cc1;
695  d2 = 4.0 * ao * tsi * cc1sq;
696  temp = d2 * tsi * cc1 / 3.0;
697  d3 = (17.0 * ao + sfour) * temp;
698  d4 = 0.5 * temp * ao * tsi * (221.0 * ao + 31.0 * sfour) * cc1;
699  t3cof = d2 + 2.0 * cc1sq;
700  t4cof = 0.25 * (3.0 * d3 + cc1 * (12.0 * d2 + 10.0 * cc1sq));
701  t5cof = 0.2 * (3.0 * d4 + 12.0 * cc1 * d3 + 6.0 * d2 * d2 + 15.0 * cc1sq * (2.0 * d2 + cc1sq));
702  }
703  }
704 }
705 
707 {
708  KStarsData *data = KStarsData::Instance();
709  return sgp4((data->clock()->utc().djd() - m_tle_jd) * MINPD);
710 }
711 
712 int Satellite::sgp4(double tsince)
713 {
714  KStarsData *data = KStarsData::Instance();
715 
716  int ktr;
717  double am, axnl, aynl, betal, cosim, cnod, cos2u, coseo1 = 0, cosi, cosip, cosisq, cossu, cosu, delm, delomg, em,
718  ecose, el2, eo1, ep, esine, argpm, argpp, argpdf, pl,
719  mrt = 0.0, mvt, rdotl, rl, rvdot, rvdotl, sinim, dndt, sin2u, sineo1 = 0, sini, sinip, sinsu, sinu, snod, su, t2,
720  t3, t4, tem5, temp, temp1, temp2, tempa, tempe, templ, u, ux, uy, uz, vx, vy, vz, inclm, mm, nm, nodem, xinc,
721  xincp, xl, xlm, mp, xmdf, xmx, xmy, nodedf, xnode, nodep, tc, sat_posx, sat_posy, sat_posz, sat_posw, sat_velx,
722  sat_vely, sat_velz, sinlat, obs_posx, obs_posy, obs_posz, obs_posw, /*obs_velx, obs_vely, obs_velz,*/
723  coslat, thetageo, sintheta, costheta, c, sq, achcp, vkmpersec;
724  // double emsq;
725 
726  const double temp4 = 1.5e-12;
727 
728  double jul_utc = data->clock()->utc().djd();
729 
730  vkmpersec = RADIUSEARTHKM * XKE / 60.0;
731 
732  // Update for secular gravity and atmospheric drag
733  xmdf = m_mean_anomaly + mdot * tsince;
734  argpdf = m_arg_perigee + argpdot * tsince;
735  nodedf = m_ra + nodedot * tsince;
736  argpm = argpdf;
737  mm = xmdf;
738  t2 = tsince * tsince;
739  nodem = nodedf + nodecf * t2;
740  tempa = 1.0 - cc1 * tsince;
741  tempe = m_bstar * cc4 * tsince;
742  templ = t2cof * t2;
743 
744  if (!isimp)
745  {
746  delomg = omgcof * tsince;
747  delm = xmcof * (pow((1.0 + eta * cos(xmdf)), 3) - delmo);
748  temp = delomg + delm;
749  mm = xmdf + temp;
750  argpm = argpdf - temp;
751  t3 = t2 * tsince;
752  t4 = t3 * tsince;
753  tempa = tempa - d2 * t2 - d3 * t3 - d4 * t4;
754  tempe = tempe + m_bstar * cc5 * (sin(mm) - sinmao);
755  templ = templ + t3cof * t3 + t4 * (t4cof + tsince * t5cof);
756  }
757 
758  nm = m_mean_motion;
759  em = m_eccentricity;
760  inclm = m_inclination;
761 
762  if (method == 'd')
763  {
764  tc = tsince;
765  // Deep space contributions to mean elements for perturbing third body
766  int iretn;
767  double delt, ft, theta, x2li, x2omi, xl, xldot, xnddt, xndt, xomi;
768 
769  // Define some constants
770  const double fasx2 = 0.13130908;
771  const double fasx4 = 2.8843198;
772  const double fasx6 = 0.37448087;
773  const double g22 = 5.7686396;
774  const double g32 = 0.95240898;
775  const double g44 = 1.8014998;
776  const double g52 = 1.0508330;
777  const double g54 = 4.4108898;
778  const double rptim = 4.37526908801129966e-3; // this equates to 7.29211514668855e-5 rad/sec
779  const double step = 720.0;
780  const double step2 = step * step / 2;
781 
782  // Calculate deep space resonance effects
783  // Value never used
784  // dndt = 0.0;
785  theta = fmod(gsto + tc * rptim, TWOPI);
786  em = em + dedt * tsince;
787 
788  inclm = inclm + didt * tsince;
789  argpm = argpm + domdt * tsince;
790  nodem = nodem + dnodt * tsince;
791  mm = mm + dmdt * tsince;
792 
793  // Update resonances : numerical (euler-maclaurin) integration
794  ft = 0.0;
795  if (irez != 0)
796  {
797  if ((atime == 0.0) || (tsince * atime <= 0.0) || (fabs(tsince) < fabs(atime)))
798  {
799  atime = 0.0;
800  xni = m_mean_motion;
801  xli = xlamo;
802  }
803 
804  if (tsince > 0.0)
805  delt = step;
806  else
807  delt = -step;
808 
809  iretn = 381; // added for do loop
810 
811  while (iretn == 381)
812  {
813  // Near - synchronous resonance terms
814  if (irez != 2)
815  {
816  xndt = del1 * sin(xli - fasx2) + del2 * sin(2.0 * (xli - fasx4)) + del3 * sin(3.0 * (xli - fasx6));
817  xldot = xni + xfact;
818  xnddt = del1 * cos(xli - fasx2) + 2.0 * del2 * cos(2.0 * (xli - fasx4)) +
819  3.0 * del3 * cos(3.0 * (xli - fasx6));
820  xnddt = xnddt * xldot;
821  }
822  else
823  {
824  // Near - half-day resonance terms
825  xomi = m_arg_perigee + argpdot * atime;
826  x2omi = xomi + xomi;
827  x2li = xli + xli;
828  xndt = d2201 * sin(x2omi + xli - g22) + d2211 * sin(xli - g22) + d3210 * sin(xomi + xli - g32) +
829  d3222 * sin(-xomi + xli - g32) + d4410 * sin(x2omi + x2li - g44) + d4422 * sin(x2li - g44) +
830  d5220 * sin(xomi + xli - g52) + d5232 * sin(-xomi + xli - g52) +
831  d5421 * sin(xomi + x2li - g54) + d5433 * sin(-xomi + x2li - g54);
832  xldot = xni + xfact;
833  xnddt = d2201 * cos(x2omi + xli - g22) + d2211 * cos(xli - g22) + d3210 * cos(xomi + xli - g32) +
834  d3222 * cos(-xomi + xli - g32) + d5220 * cos(xomi + xli - g52) +
835  d5232 * cos(-xomi + xli - g52) +
836  2.0 * (d4410 * cos(x2omi + x2li - g44) + d4422 * cos(x2li - g44) +
837  d5421 * cos(xomi + x2li - g54) + d5433 * cos(-xomi + x2li - g54));
838  xnddt = xnddt * xldot;
839  }
840 
841  if (fabs(tsince - atime) >= step)
842  {
843  iretn = 381;
844  }
845  else
846  {
847  ft = tsince - atime;
848  iretn = 0;
849  }
850 
851  if (iretn == 381)
852  {
853  xli = xli + xldot * delt + xndt * step2;
854  xni = xni + xndt * delt + xnddt * step2;
855  atime = atime + delt;
856  }
857  }
858 
859  nm = xni + xndt * ft + xnddt * ft * ft * 0.5;
860  xl = xli + xldot * ft + xndt * ft * ft * 0.5;
861 
862  if (irez != 1)
863  {
864  mm = xl - 2.0 * nodem + 2.0 * theta;
865  dndt = nm - m_mean_motion;
866  }
867  else
868  {
869  mm = xl - nodem - argpm + theta;
870  dndt = nm - m_mean_motion;
871  }
872 
873  nm = m_mean_motion + dndt;
874  }
875  }
876 
877  if (nm <= 0.0)
878  {
879  qDebug() << Q_FUNC_INFO << "Mean motion less than 0.0";
880  return (2);
881  }
882 
883  am = pow((XKE / nm), X2O3) * tempa * tempa;
884  nm = XKE / pow(am, 1.5);
885  em = em - tempe;
886 
887  if ((em >= 1.0) || (em < -0.001))
888  {
889  qDebug() << Q_FUNC_INFO << "Eccentricity >= 1.0 or < -0.001";
890  return (1);
891  }
892 
893  if (em < 1.0e-6)
894  em = 1.0e-6;
895 
896  mm = mm + m_mean_motion * templ;
897  xlm = mm + argpm + nodem;
898  // Value never used
899  // emsq = em * em;
900  // Value never used
901  // temp = 1.0 - emsq;
902 
903  nodem = fmod(nodem, TWOPI);
904  argpm = fmod(argpm, TWOPI);
905  xlm = fmod(xlm, TWOPI);
906  mm = fmod(xlm - argpm - nodem, TWOPI);
907 
908  // Compute extra mean quantities
909  sinim = sin(inclm);
910  cosim = cos(inclm);
911 
912  // Add lunar-solar periodics
913  ep = em;
914  xincp = inclm;
915  argpp = argpm;
916  nodep = nodem;
917  mp = mm;
918  sinip = sinim;
919  cosip = cosim;
920  if (method == 'd')
921  {
922  double alfdp, betdp, cosip, cosop, dalf, dbet, dls, f2, f3, pe, pgh, ph, pinc, pl, sel, ses, sghl, sghs, shll,
923  shs, sil, sinip, sinop, sinzf, sis, sll, sls, xls, xnoh, zf, zm;
924 
925  // Define some constants
926  const double zns = 1.19459e-5;
927  const double zes = 0.01675;
928  const double znl = 1.5835218e-4;
929  const double zel = 0.05490;
930 
931  // Calculate time varying periodics
932  zm = zmos + zns * tsince;
933  zf = zm + 2.0 * zes * sin(zm);
934  sinzf = sin(zf);
935  f2 = 0.5 * sinzf * sinzf - 0.25;
936  f3 = -0.5 * sinzf * cos(zf);
937  ses = se2 * f2 + se3 * f3;
938  sis = si2 * f2 + si3 * f3;
939  sls = sl2 * f2 + sl3 * f3 + sl4 * sinzf;
940  sghs = sgh2 * f2 + sgh3 * f3 + sgh4 * sinzf;
941  shs = sh2 * f2 + sh3 * f3;
942  zm = zmol + znl * tsince;
943 
944  zf = zm + 2.0 * zel * sin(zm);
945  sinzf = sin(zf);
946  f2 = 0.5 * sinzf * sinzf - 0.25;
947  f3 = -0.5 * sinzf * cos(zf);
948  sel = ee2 * f2 + e3 * f3;
949  sil = xi2 * f2 + xi3 * f3;
950  sll = xl2 * f2 + xl3 * f3 + xl4 * sinzf;
951  sghl = xgh2 * f2 + xgh3 * f3 + xgh4 * sinzf;
952  shll = xh2 * f2 + xh3 * f3;
953  pe = ses + sel;
954  pinc = sis + sil;
955  pl = sls + sll;
956  pgh = sghs + sghl;
957  ph = shs + shll;
958 
959  pe = pe - peo;
960  pinc = pinc - pinco;
961  pl = pl - plo;
962  pgh = pgh - pgho;
963  ph = ph - pho;
964  xincp = xincp + pinc;
965  ep = ep + pe;
966  sinip = sin(xincp);
967  cosip = cos(xincp);
968 
969  // Apply periodics directly
970  if (xincp >= 0.2)
971  {
972  ph = ph / sinip;
973  pgh = pgh - cosip * ph;
974  argpp = argpp + pgh;
975  nodep = nodep + ph;
976  mp = mp + pl;
977  }
978  else
979  {
980  // Apply periodics with lyddane modification
981  sinop = sin(nodep);
982  cosop = cos(nodep);
983  alfdp = sinip * sinop;
984  betdp = sinip * cosop;
985  dalf = ph * cosop + pinc * cosip * sinop;
986  dbet = -ph * sinop + pinc * cosip * cosop;
987  alfdp = alfdp + dalf;
988  betdp = betdp + dbet;
989  nodep = fmod(nodep, TWOPI);
990  if (nodep < 0.0)
991  nodep += TWOPI;
992  xls = mp + argpp + cosip * nodep;
993  dls = pl + pgh - pinc * nodep * sinip;
994  xls = xls + dls;
995  xnoh = nodep;
996  nodep = atan2(alfdp, betdp);
997  if ((nodep < 0.0))
998  nodep += TWOPI;
999  if (fabs(xnoh - nodep) > M_PI)
1000  {
1001  if (nodep < xnoh)
1002  nodep += TWOPI;
1003  else
1004  nodep -= TWOPI;
1005  }
1006  mp = mp + pl;
1007  argpp = xls - mp - cosip * nodep;
1008  }
1009 
1010  if (xincp < 0.0)
1011  {
1012  xincp = -xincp;
1013  nodep = nodep + M_PI;
1014  argpp = argpp - M_PI;
1015  }
1016 
1017  if ((ep < 0.0) || (ep > 1.0))
1018  {
1019  qDebug() << Q_FUNC_INFO << "Eccentricity < 0.0 or > 1.0";
1020  return (3);
1021  }
1022  }
1023 
1024  // Long period periodics
1025  if (method == 'd')
1026  {
1027  sinip = sin(xincp);
1028  cosip = cos(xincp);
1029  aycof = -0.5 * J3OJ2 * sinip;
1030  if (fabs(cosip + 1.0) > 1.5e-12)
1031  xlcof = -0.25 * J3OJ2 * sinip * (3.0 + 5.0 * cosip) / (1.0 + cosip);
1032  else
1033  xlcof = -0.25 * J3OJ2 * sinip * (3.0 + 5.0 * cosip) / temp4;
1034  }
1035  axnl = ep * cos(argpp);
1036  temp = 1.0 / (am * (1.0 - ep * ep));
1037  aynl = ep * sin(argpp) + temp * aycof;
1038  xl = mp + argpp + nodep + temp * xlcof * axnl;
1039 
1040  // Solve kepler's equation
1041  u = fmod(xl - nodep, TWOPI);
1042  eo1 = u;
1043  tem5 = 9999.9;
1044  ktr = 1;
1045  while ((fabs(tem5) >= 1.0e-12) && (ktr <= 10))
1046  {
1047  sineo1 = sin(eo1);
1048  coseo1 = cos(eo1);
1049  tem5 = 1.0 - coseo1 * axnl - sineo1 * aynl;
1050  tem5 = (u - aynl * coseo1 + axnl * sineo1 - eo1) / tem5;
1051  if (fabs(tem5) >= 0.95)
1052  tem5 = tem5 > 0.0 ? 0.95 : -0.95;
1053  eo1 = eo1 + tem5;
1054  ktr = ktr + 1;
1055  }
1056 
1057  // Short period preliminary quantities
1058  ecose = axnl * coseo1 + aynl * sineo1;
1059  esine = axnl * sineo1 - aynl * coseo1;
1060  el2 = axnl * axnl + aynl * aynl;
1061  pl = am * (1.0 - el2);
1062 
1063  if (pl < 0.0)
1064  {
1065  qDebug() << Q_FUNC_INFO << "Semi-latus rectum < 0.0";
1066  return (4);
1067  }
1068 
1069  rl = am * (1.0 - ecose);
1070  rdotl = sqrt(am) * esine / rl;
1071  rvdotl = sqrt(pl) / rl;
1072  betal = sqrt(1.0 - el2);
1073  temp = esine / (1.0 + betal);
1074  sinu = am / rl * (sineo1 - aynl - axnl * temp);
1075  cosu = am / rl * (coseo1 - axnl + aynl * temp);
1076  su = atan2(sinu, cosu);
1077  sin2u = (cosu + cosu) * sinu;
1078  cos2u = 1.0 - 2.0 * sinu * sinu;
1079  temp = 1.0 / pl;
1080  temp1 = 0.5 * J2 * temp;
1081  temp2 = temp1 * temp;
1082 
1083  // Update for short period periodics
1084  if (method == 'd')
1085  {
1086  cosisq = cosip * cosip;
1087  con41 = 3.0 * cosisq - 1.0;
1088  x1mth2 = 1.0 - cosisq;
1089  x7thm1 = 7.0 * cosisq - 1.0;
1090  }
1091  mrt = rl * (1.0 - 1.5 * temp2 * betal * con41) + 0.5 * temp1 * x1mth2 * cos2u;
1092  su = su - 0.25 * temp2 * x7thm1 * sin2u;
1093  xnode = nodep + 1.5 * temp2 * cosip * sin2u;
1094  xinc = xincp + 1.5 * temp2 * cosip * sinip * cos2u;
1095  mvt = rdotl - nm * temp1 * x1mth2 * sin2u / XKE;
1096  rvdot = rvdotl + nm * temp1 * (x1mth2 * cos2u + 1.5 * con41) / XKE;
1097 
1098  // Orientation vectors
1099  sinsu = sin(su);
1100  cossu = cos(su);
1101  snod = sin(xnode);
1102  cnod = cos(xnode);
1103  sini = sin(xinc);
1104  cosi = cos(xinc);
1105  xmx = -snod * cosi;
1106  xmy = cnod * cosi;
1107  ux = xmx * sinsu + cnod * cossu;
1108  uy = xmy * sinsu + snod * cossu;
1109  uz = sini * sinsu;
1110  vx = xmx * cossu - cnod * sinsu;
1111  vy = xmy * cossu - snod * sinsu;
1112  vz = sini * cossu;
1113 
1114  // Position and velocity (in km and km/sec)
1115  sat_posx = (mrt * ux) * RADIUSEARTHKM;
1116  sat_posy = (mrt * uy) * RADIUSEARTHKM;
1117  sat_posz = (mrt * uz) * RADIUSEARTHKM;
1118  sat_posw = sqrt(sat_posx * sat_posx + sat_posy * sat_posy + sat_posz * sat_posz);
1119  sat_velx = (mvt * ux + rvdot * vx) * vkmpersec;
1120  sat_vely = (mvt * uy + rvdot * vy) * vkmpersec;
1121  sat_velz = (mvt * uz + rvdot * vz) * vkmpersec;
1122  m_velocity = sqrt(sat_velx * sat_velx + sat_vely * sat_vely + sat_velz * sat_velz);
1123 
1124  // printf("tsince=%.15f\n", tsince);
1125  // printf("sat_posx=%.15f\n", sat_posx);
1126  // printf("sat_posy=%.15f\n", sat_posy);
1127  // printf("sat_posz=%.15f\n", sat_posz);
1128  // printf("sat_velx=%.15f\n", sat_velx);
1129  // printf("sat_vely=%.15f\n", sat_vely);
1130  // printf("sat_velz=%.15f\n", sat_velz);
1131 
1132  if (mrt < 1.0)
1133  {
1134  qDebug() << Q_FUNC_INFO << "Satellite has decayed";
1135  return (6);
1136  }
1137 
1138  // Observer ECI position and velocity
1139  sinlat = sin(data->geo()->lat()->radians());
1140  coslat = cos(data->geo()->lat()->radians());
1141  thetageo = data->geo()->LMST(jul_utc);
1142  sintheta = sin(thetageo);
1143  costheta = cos(thetageo);
1144  c = 1.0 / sqrt(1.0 + F * (F - 2.0) * sinlat * sinlat);
1145  sq = (1.0 - F) * (1.0 - F) * c;
1146  achcp = (RADIUSEARTHKM * c + MEANALT) * coslat;
1147  obs_posx = achcp * costheta;
1148  obs_posy = achcp * sintheta;
1149  obs_posz = (RADIUSEARTHKM * sq + MEANALT) * sinlat;
1150  obs_posw = sqrt(obs_posx * obs_posx + obs_posy * sat_posy + obs_posz * obs_posz);
1151  /*obs_velx = -MFACTOR * obs_posy;
1152  obs_vely = MFACTOR * obs_posx;
1153  obs_velz = 0.;*/
1154 
1155  m_altitude = sat_posw - obs_posw + MEANALT;
1156 
1157  // Az and Dec
1158  double range_posx = sat_posx - obs_posx;
1159  double range_posy = sat_posy - obs_posy;
1160  double range_posz = sat_posz - obs_posz;
1161  m_range = sqrt(range_posx * range_posx + range_posy * range_posy + range_posz * range_posz);
1162  // double range_velx = sat_velx - obs_velx;
1163  // double range_vely = sat_velx - obs_vely;
1164  // double range_velz = sat_velx - obs_velz;
1165 
1166  double top_s = sinlat * costheta * range_posx + sinlat * sintheta * range_posy - coslat * range_posz;
1167  double top_e = -sintheta * range_posx + costheta * range_posy;
1168  double top_z = coslat * costheta * range_posx + coslat * sintheta * range_posy + sinlat * range_posz;
1169 
1170  double azimuth = atan(-top_e / top_s);
1171  if (top_s > 0.)
1172  azimuth += M_PI;
1173  if (azimuth < 0.)
1174  azimuth += TWOPI;
1175  double elevation = arcSin(top_z / m_range);
1176 
1177  // printf("azimuth=%.15f\n\r", azimuth / DEG2RAD);
1178  // printf("elevation=%.15f\n\r", elevation / DEG2RAD);
1179 
1180  setAz(azimuth / DEG2RAD);
1181  setAlt(elevation / DEG2RAD);
1182  HorizontalToEquatorial(data->lst(), data->geo()->lat());
1183 
1184  // is the satellite visible ?
1185  // Find ECI coordinates of the sun
1186  double mjd, year, T, M, L, e, C, O, Lsa, nu, R, eps;
1187 
1188  mjd = jul_utc - 2415020.0;
1189  year = 1900.0 + mjd / 365.25;
1190  T = (mjd + deltaET(year) / (MINPD * 60.0)) / 36525.0;
1191  M = DEG2RAD * (Modulus(358.47583 + Modulus(35999.04975 * T, 360.0) - (0.000150 + 0.0000033 * T) * T * T, 360.0));
1192  L = DEG2RAD * (Modulus(279.69668 + Modulus(36000.76892 * T, 360.0) + 0.0003025 * T * T, 360.0));
1193  e = 0.01675104 - (0.0000418 + 0.000000126 * T) * T;
1194  C = DEG2RAD * ((1.919460 - (0.004789 + 0.000014 * T) * T) * sin(M) + (0.020094 - 0.000100 * T) * sin(2 * M) +
1195  0.000293 * sin(3 * M));
1196  O = DEG2RAD * (Modulus(259.18 - 1934.142 * T, 360.0));
1197  Lsa = Modulus(L + C - DEG2RAD * (0.00569 - 0.00479 * sin(O)), TWOPI);
1198  nu = Modulus(M + C, TWOPI);
1199  R = 1.0000002 * (1.0 - e * e) / (1.0 + e * cos(nu));
1200  eps = DEG2RAD * (23.452294 - (0.0130125 + (0.00000164 - 0.000000503 * T) * T) * T + 0.00256 * cos(O));
1201  R = AU * R;
1202 
1203  double sun_posx = R * cos(Lsa);
1204  double sun_posy = R * sin(Lsa) * cos(eps);
1205  double sun_posz = R * sin(Lsa) * sin(eps);
1206  double sun_posw = R;
1207 
1208  // Calculates satellite's eclipse status and depth
1209  double sd_sun, sd_earth, delta, depth;
1210 
1211  // Determine partial eclipse
1212  sd_earth = arcSin(RADIUSEARTHKM / sat_posw);
1213  double rho_x = sun_posx - sat_posx;
1214  double rho_y = sun_posy - sat_posy;
1215  double rho_z = sun_posz - sat_posz;
1216  double rho_w = sqrt(rho_x * rho_x + rho_y * rho_y + rho_z * rho_z);
1217  sd_sun = arcSin(SR / rho_w);
1218  double earth_x = -1.0 * sat_posx;
1219  double earth_y = -1.0 * sat_posy;
1220  double earth_z = -1.0 * sat_posz;
1221  double earth_w = sat_posw;
1222  delta = PIO2 - arcSin((sun_posx * earth_x + sun_posy * earth_y + sun_posz * earth_z) / (sun_posw * earth_w));
1223  depth = sd_earth - sd_sun - delta;
1224  KSSun *sun = dynamic_cast<KSSun *>(data->skyComposite()->findByName(i18n("Sun")));
1225 
1226  m_is_eclipsed = sd_earth >= sd_sun && depth >= 0;
1227  m_is_visible = !m_is_eclipsed && sun->alt().Degrees() <= -12.0 && elevation >= 0.0;
1228 
1229  return (0);
1230 }
1231 
1233 {
1234  switch (code)
1235  {
1236  case 0:
1237  return i18n("Success");
1238 
1239  case 1:
1240  case 3:
1241  return i18n("Eccentricity >= 1.0 or < -0.001");
1242 
1243  case 2:
1244  return i18n("Mean motion less than 0.0");
1245 
1246  case 4:
1247  return i18n("Semi-latus rectum < 0.0");
1248 
1249  case 6:
1250  return i18n("Satellite has decayed");
1251 
1252  default:
1253  return i18n("Unknown error");
1254  }
1255 }
1256 
1257 double Satellite::arcSin(double arg)
1258 {
1259  if (fabs(arg) >= 1.)
1260  if (arg > 0.)
1261  return PIO2;
1262  else if (arg < 0.)
1263  return -PIO2;
1264  else
1265  return 0.;
1266  else
1267  return (atan(arg / sqrt(1. - arg * arg)));
1268 }
1269 
1270 double Satellite::deltaET(double year)
1271 {
1272  double delta_et;
1273 
1274  delta_et = 26.465 + 0.747622 * (year - 1950) + 1.886913 * sin(TWOPI * (year - 1975) / 33);
1275 
1276  return delta_et;
1277 }
1278 
1279 double Satellite::Modulus(double arg1, double arg2)
1280 {
1281  int i;
1282  double ret_val;
1283 
1284  ret_val = arg1;
1285  i = ret_val / arg2;
1286  ret_val -= i * arg2;
1287 
1288  if (ret_val < 0.0)
1289  ret_val += arg2;
1290 
1291  return ret_val;
1292 }
1293 
1295 {
1296  return m_is_visible;
1297 }
1298 
1300 {
1301  return m_is_selected;
1302 }
1303 
1304 void Satellite::setSelected(bool selected)
1305 {
1306  m_is_selected = selected;
1307 }
1308 
1310 {
1311 #ifndef KSTARS_LITE
1312  pmenu->createSatelliteMenu(this);
1313 #else
1314  Q_UNUSED(pmenu);
1315 #endif
1316 }
1317 
1318 double Satellite::velocity() const
1319 {
1320  return m_velocity;
1321 }
1322 
1323 double Satellite::altitude() const
1324 {
1325  return m_altitude;
1326 }
1327 
1328 double Satellite::range() const
1329 {
1330  return m_range;
1331 }
1332 
1334 {
1335  return m_id;
1336 }
1337 
1339 {
1340  return m_tle;
1341 }
const dms & alt() const
Definition: skypoint.h:281
bool isVisible()
Definition: satellite.cpp:1294
Satellite * clone() const override
Definition: satellite.cpp:96
void setAlt(dms alt)
Sets Alt, the Altitude.
Definition: skypoint.h:194
int toInt(bool *ok, int base) const const
QString id() const
Definition: satellite.cpp:1333
long double djd() const
void initPopupMenu(KSPopupMenu *pmenu) override
Initialize the popup menut.
Definition: satellite.cpp:1309
double altitude() const
Definition: satellite.cpp:1323
QStringRef midRef(int position, int n) const const
CachingDms * lst()
Definition: kstarsdata.h:224
virtual QString name(void) const
Definition: skyobject.h:145
Satellite(const QString &name, const QString &line1, const QString &line2)
Constructor.
Definition: satellite.cpp:48
bool selected()
Definition: satellite.cpp:1299
double velocity() const
Definition: satellite.cpp:1318
Provides necessary information about the Sun.
Definition: kssun.h:23
SkyObject * findByName(const QString &name, bool exact=true) override
Search the children of this SkyMapComposite for a SkyObject whose name matches the argument.
QString i18n(const char *text, const TYPE &arg...)
const CachingDms * lat() const
Definition: geolocation.h:70
void setMag(float m)
Set the object's sorting magnitude.
Definition: skyobject.h:403
double toDouble(bool *ok) const const
GeoLocation * geo()
Definition: kstarsdata.h:230
double LMST(double jd)
double range() const
Definition: satellite.cpp:1328
Q_INVOKABLE SimClock * clock()
Definition: kstarsdata.h:218
void setName(const QString &name)
Set the object's primary name.
Definition: skyobject.h:417
QString tle() const
Definition: satellite.cpp:1338
SkyMapComposite * skyComposite()
Definition: kstarsdata.h:166
void setAz(dms az)
Sets Az, the Azimuth.
Definition: skypoint.h:230
double radians() const
Express the angle in radians.
Definition: dms.h:325
const double & Degrees() const
Definition: dms.h:141
QString sgp4ErrorString(int code)
sgp4ErrorString Get error string associated with sgp4 calculation failure
Definition: satellite.cpp:1232
const QChar at(int position) const const
void setSelected(bool selected)
Select or not the satellite.
Definition: satellite.cpp:1304
void HorizontalToEquatorial(const dms *LST, const dms *lat)
Determine the (RA, Dec) coordinates of the SkyPoint from its (Altitude, Azimuth) coordinates,...
Definition: skypoint.cpp:143
int updatePos()
Update satellite position.
Definition: satellite.cpp:706
QString mid(int position, int n) const const
void setName2(const QString &name2=QString())
Set the object's secondary name.
Definition: skyobject.h:423
const KStarsDateTime & utc() const
Definition: simclock.h:35
void setType(int t)
Set the object's type identifier to the argument.
Definition: skyobject.h:195
void createSatelliteMenu(Satellite *satellite)
Create a popup menu for a satellite.
void setLongName(const QString &longname=QString())
Set the object's long name.
Definition: skyobject.cpp:76
This file is part of the KDE documentation.
Documentation copyright © 1996-2023 The KDE developers.
Generated on Fri Sep 22 2023 03:57:58 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.