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
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
48Satellite::Satellite(const QString &name, const QString &line1, const QString &line2)
49{
50 //m_name = name;
51 m_number = line1.mid(2, 5).toInt();
52 m_class = line1.at(7);
53 m_id = line1.mid(9, 8);
54 m_epoch = line1.mid(18, 14).toDouble();
55 m_first_deriv = line1.mid(33, 10).toDouble() / (XPDOTP * MINPD);
56 m_second_deriv =
57 line1.mid(44, 6).toDouble() * (1.0e-5 / pow(10.0, line1.mid(51, 1).toDouble())) / (XPDOTP * MINPD * MINPD);
58 m_bstar = line1.mid(53, 6).toDouble() * 1.0e-5 / pow(10.0, line1.mid(60, 1).toDouble());
59 m_ephem_type = line1.mid(62, 1).toInt();
60 m_elem_number = line1.mid(64, 4).toInt();
61 m_inclination = line2.mid(8, 8).toDouble() * DEG2RAD;
62 m_ra = line2.mid(17, 8).toDouble() * DEG2RAD;
63 m_eccentricity = line2.mid(26, 7).toDouble() * 1.0e-7;
64 m_arg_perigee = line2.mid(34, 8).toDouble() * DEG2RAD;
65 m_mean_anomaly = line2.mid(43, 8).toDouble() * DEG2RAD;
66 m_mean_motion = line2.mid(52, 11).toDouble() * TWOPI / MINPD;
67 m_nb_revolution = line2.mid(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
96Satellite *Satellite::clone() const
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
102void 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
706int Satellite::updatePos()
707{
708 KStarsData *data = KStarsData::Instance();
709 return sgp4((data->clock()->utc().djd() - m_tle_jd) * MINPD);
710}
711
712int 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
1232QString Satellite::sgp4ErrorString(int code)
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
1257double 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
1270double 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
1279double 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
1294bool Satellite::isVisible()
1295{
1296 return m_is_visible;
1297}
1298
1299bool Satellite::selected()
1300{
1301 return m_is_selected;
1302}
1303
1304void Satellite::setSelected(bool selected)
1305{
1306 m_is_selected = selected;
1307}
1308
1309void Satellite::initPopupMenu(KSPopupMenu *pmenu)
1310{
1311#ifndef KSTARS_LITE
1312 pmenu->createSatelliteMenu(this);
1313#else
1314 Q_UNUSED(pmenu);
1315#endif
1316}
1317
1318double Satellite::velocity() const
1319{
1320 return m_velocity;
1321}
1322
1323double Satellite::altitude() const
1324{
1325 return m_altitude;
1326}
1327
1328double Satellite::range() const
1329{
1330 return m_range;
1331}
1332
1333QString Satellite::id() const
1334{
1335 return m_id;
1336}
1337
1338QString Satellite::tle() const
1339{
1340 return m_tle;
1341}
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
KSSun
Child class of KSPlanetBase; encapsulates information about the Sun.
Definition kssun.h:24
KStarsData
KStarsData is the backbone of KStars.
Definition kstarsdata.h:72
KStarsData::lst
CachingDms * lst()
Definition kstarsdata.h:224
KStarsData::clock
Q_INVOKABLE SimClock * clock()
Definition kstarsdata.h:218
KStarsData::geo
GeoLocation * geo()
Definition kstarsdata.h:230
KStarsData::skyComposite
SkyMapComposite * skyComposite()
Definition kstarsdata.h:166
KStarsDateTime::djd
long double djd() const
Definition kstarsdatetime.h:163
Satellite
Represents an artificial satellites.
Definition satellite.h:23
Satellite::updatePos
int updatePos()
Update satellite position.
Definition satellite.cpp:706
Satellite::setSelected
void setSelected(bool selected)
Select or not the satellite.
Definition satellite.cpp:1304
Satellite::range
double range() const
Definition satellite.cpp:1328
Satellite::id
QString id() const
Definition satellite.cpp:1333
Satellite::altitude
double altitude() const
Definition satellite.cpp:1323
Satellite::sgp4ErrorString
QString sgp4ErrorString(int code)
sgp4ErrorString Get error string associated with sgp4 calculation failure
Definition satellite.cpp:1232
Satellite::velocity
double velocity() const
Definition satellite.cpp:1318
Satellite::initPopupMenu
void initPopupMenu(KSPopupMenu *pmenu) override
Initialize the popup menut.
Definition satellite.cpp:1309
Satellite::tle
QString tle() const
Definition satellite.cpp:1338
Satellite::isVisible
bool isVisible()
Definition satellite.cpp:1294
Satellite::Satellite
Satellite(const QString &name, const QString &line1, const QString &line2)
Constructor.
Definition satellite.cpp:48
Satellite::selected
bool selected()
Definition satellite.cpp:1299
Satellite::clone
Satellite * clone() const override
Definition satellite.cpp:96
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:423
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:403
SkyObject::name
virtual QString name(void) const
Definition skyobject.h:145
SkyObject::setName
void setName(const QString &name)
Set the object's primary name.
Definition skyobject.h:417
SkyObject::setType
void setType(int t)
Set the object's type identifier to the argument.
Definition skyobject.h:195
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-2024 The KDE developers.
Generated on Mon Nov 4 2024 16:38:44 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.