Kstars

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

KDE's Doxygen guidelines are available online.