File indexing completed on 2024-04-21 03:44:42
0001 /* 0002 SPDX-FileCopyrightText: 2009 Jerome SONRIER <jsid@emor3j.fr.eu.org> 0003 0004 SPDX-License-Identifier: GPL-2.0-or-later 0005 */ 0006 0007 #include "satellite.h" 0008 0009 #include "ksplanetbase.h" 0010 #ifndef KSTARS_LITE 0011 #include "kspopupmenu.h" 0012 #endif 0013 #include "kstarsdata.h" 0014 #include "kssun.h" 0015 #include "Options.h" 0016 #include "skymapcomposite.h" 0017 0018 #include <QDebug> 0019 0020 #include <cmath> 0021 #include <typeinfo> 0022 0023 // Define some constants 0024 // WGS-72 constants 0025 #define RADIUSEARTHKM 6378.135 // Earth radius (km) 0026 #define XKE 0.07436691613317 // 60.0 / sqrt(RADIUSEARTHKM^3/MU) 0027 #define J2 0.001082616 // The second gravitational zonal harmonic of the Earth 0028 #define J4 -0.00000165597 // The fourth gravitational zonal harmonic of the Earth 0029 #define J3OJ2 -2.34506972e-3 // J3 / J2 0030 0031 // Mathematical constants 0032 #define TWOPI 6.2831853071795864769 // 2*PI 0033 #define PIO2 1.5707963267948966192 // PI/2 0034 #define X2O3 .66666666666666666667 // 2/3 0035 #define DEG2RAD 1.745329251994330e-2 // Deg -> Rad 0036 0037 // Other constants 0038 #define MINPD 1440 // Minutes per day 0039 #define MEANALT 0.84 // Mean altitude (km) 0040 #define SR 6.96000e5 // Solar radius - km (IAU 76) 0041 #define AU 1.49597870691e8 // Astronomical unit - km (IAU 76) 0042 #define XPDOTP 229.1831180523293 // 1440.0 / (2.0 * pi) 0043 #define SS 1.0122292801892716288 // Parameter for the SGP4 density function 0044 #define QZMS2T 1.8802791590152706439e-9 // (( 120.0 - 78.0) / RADIUSEARTHKM )^4 0045 #define F 3.35281066474748e-3 // Flattening factor 0046 #define MFACTOR 7.292115e-5 0047 0048 Satellite::Satellite(const QString &name, const QString &line1, const QString &line2) 0049 { 0050 //m_name = name; 0051 m_number = line1.midRef(2, 5).toInt(); 0052 m_class = line1.at(7); 0053 m_id = line1.mid(9, 8); 0054 m_epoch = line1.midRef(18, 14).toDouble(); 0055 m_first_deriv = line1.midRef(33, 10).toDouble() / (XPDOTP * MINPD); 0056 m_second_deriv = 0057 line1.midRef(44, 6).toDouble() * (1.0e-5 / pow(10.0, line1.midRef(51, 1).toDouble())) / (XPDOTP * MINPD * MINPD); 0058 m_bstar = line1.midRef(53, 6).toDouble() * 1.0e-5 / pow(10.0, line1.midRef(60, 1).toDouble()); 0059 m_ephem_type = line1.midRef(62, 1).toInt(); 0060 m_elem_number = line1.midRef(64, 4).toInt(); 0061 m_inclination = line2.midRef(8, 8).toDouble() * DEG2RAD; 0062 m_ra = line2.midRef(17, 8).toDouble() * DEG2RAD; 0063 m_eccentricity = line2.midRef(26, 7).toDouble() * 1.0e-7; 0064 m_arg_perigee = line2.midRef(34, 8).toDouble() * DEG2RAD; 0065 m_mean_anomaly = line2.midRef(43, 8).toDouble() * DEG2RAD; 0066 m_mean_motion = line2.midRef(52, 11).toDouble() * TWOPI / MINPD; 0067 m_nb_revolution = line2.midRef(63, 5).toInt(); 0068 m_tle = name + "\n" + line1 + "\n" + line2; 0069 0070 setName(name); 0071 setName2(name); 0072 setLongName(name + " (" + m_id + ')'); 0073 setType(SkyObject::SATELLITE); 0074 setMag(0.0); 0075 0076 m_is_selected = Options::selectedSatellites().contains(name); 0077 0078 // Convert TLE epoch to Julian date 0079 double day = modf(m_epoch * 1.e-3, &m_epoch_year) * 1.e3; 0080 if (m_epoch_year < 57.) 0081 m_epoch_year += 2000.; 0082 else 0083 m_epoch_year += 1900.; 0084 double year = m_epoch_year - 1.; 0085 long i = year / 100; 0086 long A = i; 0087 i = A / 4; 0088 long B = 2 - A + i; 0089 i = 365.25 * year; 0090 i += 30.6001 * 14; 0091 m_tle_jd = i + 1720994.5 + B + day; 0092 0093 init(); 0094 } 0095 0096 Satellite *Satellite::clone() const 0097 { 0098 Q_ASSERT(typeid(this) == typeid(static_cast<const Satellite *>(this))); // Ensure we are not slicing a derived class 0099 return new Satellite(*this); 0100 } 0101 0102 void Satellite::init() 0103 { 0104 double ao, cosio, sinio, cosio2, omeosq, posq, rp, rteosq, eccsq, con42, cnodm, snodm, cosim, sinim, cosomm, sinomm, 0105 cc1sq, cc2, cc3, coef, coef1, cosio4, day, em, emsq, eeta, etasq, gam, inclm, nm, perige, pinvsq, psisq, 0106 qzms24, rtemsq, s1, s2, s3, s4, s5, s6, s7, sfour, ss1(0), ss2(0), ss3(0), ss4(0), ss5(0), ss6(0), ss7(0), 0107 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, 0108 temp, temp1, temp2, temp3, tsi, xpidot, xhdot1, z1, z2, z3, z11, z12, z13, z21, z22, z23, z31, z32, z33, ak, d1, 0109 del, adel, po, ds70, ts70, tfrac, c1, thgr70, fk5r, c1p2p; 0110 // double dndt; 0111 0112 // Init near earth variables 0113 isimp = false; 0114 aycof = 0.; 0115 con41 = 0.; 0116 cc1 = 0.; 0117 cc4 = 0.; 0118 cc5 = 0.; 0119 d2 = 0.; 0120 d3 = 0.; 0121 d4 = 0.; 0122 delmo = 0.; 0123 eta = 0.; 0124 argpdot = 0.; 0125 omgcof = 0.; 0126 sinmao = 0.; 0127 t = 0.; 0128 t2cof = 0.; 0129 t3cof = 0.; 0130 t4cof = 0.; 0131 t5cof = 0.; 0132 x1mth2 = 0.; 0133 x7thm1 = 0.; 0134 mdot = 0.; 0135 nodedot = 0.; 0136 xlcof = 0.; 0137 xmcof = 0.; 0138 nodecf = 0.; 0139 0140 // Init deep space variables 0141 irez = 0; 0142 d2201 = 0.; 0143 d2211 = 0.; 0144 d3210 = 0.; 0145 d3222 = 0.; 0146 d4410 = 0.; 0147 d4422 = 0.; 0148 d5220 = 0.; 0149 d5232 = 0.; 0150 d5421 = 0.; 0151 d5433 = 0.; 0152 dedt = 0.; 0153 del1 = 0.; 0154 del2 = 0.; 0155 del3 = 0.; 0156 didt = 0.; 0157 dmdt = 0.; 0158 dnodt = 0.; 0159 domdt = 0.; 0160 e3 = 0.; 0161 ee2 = 0.; 0162 peo = 0.; 0163 pgho = 0.; 0164 pho = 0.; 0165 pinco = 0.; 0166 plo = 0.; 0167 se2 = 0.; 0168 se3 = 0.; 0169 sgh2 = 0.; 0170 sgh3 = 0.; 0171 sgh4 = 0.; 0172 sh2 = 0.; 0173 sh3 = 0.; 0174 si2 = 0.; 0175 si3 = 0.; 0176 sl2 = 0.; 0177 sl3 = 0.; 0178 sl4 = 0.; 0179 gsto = 0.; 0180 xfact = 0.; 0181 xgh2 = 0.; 0182 xgh3 = 0.; 0183 xgh4 = 0.; 0184 xh2 = 0.; 0185 xh3 = 0.; 0186 xi2 = 0.; 0187 xi3 = 0.; 0188 xl2 = 0.; 0189 xl3 = 0.; 0190 xl4 = 0.; 0191 xlamo = 0.; 0192 zmol = 0.; 0193 zmos = 0.; 0194 atime = 0.; 0195 xli = 0.; 0196 xni = 0.; 0197 0198 method = 'n'; 0199 0200 m_is_visible = false; 0201 0202 // Divisor for divide by zero check on inclination 0203 const double temp4 = 1.5e-12; 0204 0205 /*----- Initializes variables for sgp4 -----*/ 0206 // Calculate auxiliary epoch quantities 0207 eccsq = m_eccentricity * m_eccentricity; 0208 omeosq = 1.0 - eccsq; 0209 rteosq = sqrt(omeosq); 0210 cosio = cos(m_inclination); 0211 cosio2 = cosio * cosio; 0212 0213 // Un-kozai the mean motion 0214 ak = pow(XKE / m_mean_motion, X2O3); 0215 d1 = 0.75 * J2 * (3.0 * cosio2 - 1.0) / (rteosq * omeosq); 0216 del = d1 / (ak * ak); 0217 adel = ak * (1.0 - del * del - del * (1.0 / 3.0 + 134.0 * del * del / 81.0)); 0218 del = d1 / (adel * adel); 0219 m_mean_motion = m_mean_motion / (1.0 + del); 0220 0221 ao = pow(XKE / m_mean_motion, X2O3); 0222 sinio = sin(m_inclination); 0223 po = ao * omeosq; 0224 con42 = 1.0 - 5.0 * cosio2; 0225 con41 = -con42 - (2.0 * cosio2); 0226 posq = po * po; 0227 rp = ao * (1.0 - m_eccentricity); 0228 method = 'n'; 0229 0230 // Find sidereal time 0231 ts70 = m_tle_jd - 2433281.5 - 7305.0; 0232 ds70 = floor(ts70 + 1.0e-8); 0233 tfrac = ts70 - ds70; 0234 // find greenwich location at epoch 0235 c1 = 1.72027916940703639e-2; 0236 thgr70 = 1.7321343856509374; 0237 fk5r = 5.07551419432269442e-15; 0238 c1p2p = c1 + TWOPI; 0239 gsto = fmod(thgr70 + c1 * ds70 + c1p2p * tfrac + ts70 * ts70 * fk5r, TWOPI); 0240 if (gsto < 0.0) 0241 gsto = gsto + TWOPI; 0242 0243 if ((omeosq >= 0.0) || (m_mean_motion >= 0.0)) 0244 { 0245 if (rp < (220.0 / RADIUSEARTHKM + 1.0)) 0246 isimp = true; 0247 sfour = SS; 0248 qzms24 = QZMS2T; 0249 perige = (rp - 1.0) * RADIUSEARTHKM; 0250 0251 // For perigees below 156 km, s and qoms2t are altered 0252 if (perige < 156.0) 0253 { 0254 sfour = perige - 78.0; 0255 if (perige < 98.0) 0256 sfour = 20.0; 0257 qzms24 = pow(((120.0 - sfour) / RADIUSEARTHKM), 4.0); 0258 sfour = sfour / RADIUSEARTHKM + 1.0; 0259 } 0260 pinvsq = 1.0 / posq; 0261 0262 tsi = 1.0 / (ao - sfour); 0263 eta = ao * m_eccentricity * tsi; 0264 etasq = eta * eta; 0265 eeta = m_eccentricity * eta; 0266 psisq = fabs(1.0 - etasq); 0267 coef = qzms24 * pow(tsi, 4.0); 0268 coef1 = coef / pow(psisq, 3.5); 0269 cc2 = coef1 * m_mean_motion * 0270 (ao * (1.0 + 1.5 * etasq + eeta * (4.0 + etasq)) + 0271 0.375 * J2 * tsi / psisq * con41 * (8.0 + 3.0 * etasq * (8.0 + etasq))); 0272 cc1 = m_bstar * cc2; 0273 cc3 = 0.0; 0274 if (m_eccentricity > 1.0e-4) 0275 cc3 = -2.0 * coef * tsi * J3OJ2 * m_mean_motion * sinio / m_eccentricity; 0276 x1mth2 = 1.0 - cosio2; 0277 cc4 = 2.0 * m_mean_motion * coef1 * ao * omeosq * 0278 (eta * (2.0 + 0.5 * etasq) + m_eccentricity * (0.5 + 2.0 * etasq) - 0279 J2 * tsi / (ao * psisq) * 0280 (-3.0 * con41 * (1.0 - 2.0 * eeta + etasq * (1.5 - 0.5 * eeta)) + 0281 0.75 * x1mth2 * (2.0 * etasq - eeta * (1.0 + etasq)) * cos(2.0 * m_arg_perigee))); 0282 cc5 = 2.0 * coef1 * ao * omeosq * (1.0 + 2.75 * (etasq + eeta) + eeta * etasq); 0283 0284 cosio4 = cosio2 * cosio2; 0285 temp1 = 1.5 * J2 * pinvsq * m_mean_motion; 0286 temp2 = 0.5 * temp1 * J2 * pinvsq; 0287 temp3 = -0.46875 * J4 * pinvsq * pinvsq * m_mean_motion; 0288 mdot = m_mean_motion + 0.5 * temp1 * rteosq * con41 + 0289 0.0625 * temp2 * rteosq * (13.0 - 78.0 * cosio2 + 137.0 * cosio4); 0290 argpdot = -0.5 * temp1 * con42 + 0.0625 * temp2 * (7.0 - 114.0 * cosio2 + 395.0 * cosio4) + 0291 temp3 * (3.0 - 36.0 * cosio2 + 49.0 * cosio4); 0292 xhdot1 = -temp1 * cosio; 0293 nodedot = xhdot1 + (0.5 * temp2 * (4.0 - 19.0 * cosio2) + 2.0 * temp3 * (3.0 - 7.0 * cosio2)) * cosio; 0294 xpidot = argpdot + nodedot; 0295 omgcof = m_bstar * cc3 * cos(m_arg_perigee); 0296 xmcof = 0.0; 0297 if (m_eccentricity > 1.0e-4) 0298 xmcof = -X2O3 * coef * m_bstar / eeta; 0299 nodecf = 3.5 * omeosq * xhdot1 * cc1; 0300 t2cof = 1.5 * cc1; 0301 // Do not divide by zero 0302 if (fabs(1.0 + cosio) > 1.5e-12) 0303 xlcof = -0.25 * J3OJ2 * sinio * (3.0 + 5.0 * cosio) / (1.0 + cosio); 0304 else 0305 xlcof = -0.25 * J3OJ2 * sinio * (3.0 + 5.0 * cosio) / temp4; 0306 aycof = -0.5 * J3OJ2 * sinio; 0307 delmo = pow((1.0 + eta * cos(m_mean_anomaly)), 3); 0308 sinmao = sin(m_mean_anomaly); 0309 x7thm1 = 7.0 * cosio2 - 1.0; 0310 0311 // Deep space initialization 0312 if ((TWOPI / m_mean_motion) >= 225.0) 0313 { 0314 method = 'd'; 0315 isimp = true; 0316 tc = 0.0; 0317 inclm = m_inclination; 0318 0319 // Init deep space common variables 0320 // Define some constants 0321 const double zes = 0.01675; 0322 const double zel = 0.05490; 0323 const double c1ss = 2.9864797e-6; 0324 const double c1l = 4.7968065e-7; 0325 const double zsinis = 0.39785416; 0326 const double zcosis = 0.91744867; 0327 const double zcosgs = 0.1945905; 0328 const double zsings = -0.98088458; 0329 0330 int lsflg; 0331 double a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, betasq, cc, ctem, stem, x1, x2, x3, x4, x5, x6, x7, x8, 0332 xnodce, xnoi, zcosg, zcosgl, zcosh, zcoshl, zcosi, zcosil, zsing, zsingl, zsinh, zsinhl, zsini, zsinil, 0333 zx, zy; 0334 0335 nm = m_mean_motion; 0336 em = m_eccentricity; 0337 snodm = sin(m_ra); 0338 cnodm = cos(m_ra); 0339 sinomm = sin(m_arg_perigee); 0340 cosomm = cos(m_arg_perigee); 0341 sinim = sin(m_inclination); 0342 cosim = cos(m_inclination); 0343 emsq = em * em; 0344 betasq = 1.0 - emsq; 0345 rtemsq = sqrt(betasq); 0346 0347 // Initialize lunar solar terms 0348 peo = 0.0; 0349 pinco = 0.0; 0350 plo = 0.0; 0351 pgho = 0.0; 0352 pho = 0.0; 0353 day = m_tle_jd - 2433281.5 + 18261.5 + tc / 1440.0; 0354 xnodce = fmod(4.5236020 - 9.2422029e-4 * day, TWOPI); 0355 stem = sin(xnodce); 0356 ctem = cos(xnodce); 0357 zcosil = 0.91375164 - 0.03568096 * ctem; 0358 zsinil = sqrt(1.0 - zcosil * zcosil); 0359 zsinhl = 0.089683511 * stem / zsinil; 0360 zcoshl = sqrt(1.0 - zsinhl * zsinhl); 0361 gam = 5.8351514 + 0.0019443680 * day; 0362 zx = 0.39785416 * stem / zsinil; 0363 zy = zcoshl * ctem + 0.91744867 * zsinhl * stem; 0364 zx = atan2(zx, zy); 0365 zx = gam + zx - xnodce; 0366 zcosgl = cos(zx); 0367 zsingl = sin(zx); 0368 0369 // Solar terms 0370 zcosg = zcosgs; 0371 zsing = zsings; 0372 zcosi = zcosis; 0373 zsini = zsinis; 0374 zcosh = cnodm; 0375 zsinh = snodm; 0376 cc = c1ss; 0377 xnoi = 1.0 / nm; 0378 0379 for (lsflg = 1; lsflg <= 2; ++lsflg) 0380 { 0381 a1 = zcosg * zcosh + zsing * zcosi * zsinh; 0382 a3 = -zsing * zcosh + zcosg * zcosi * zsinh; 0383 a7 = -zcosg * zsinh + zsing * zcosi * zcosh; 0384 a8 = zsing * zsini; 0385 a9 = zsing * zsinh + zcosg * zcosi * zcosh; 0386 a10 = zcosg * zsini; 0387 a2 = cosim * a7 + sinim * a8; 0388 a4 = cosim * a9 + sinim * a10; 0389 a5 = -sinim * a7 + cosim * a8; 0390 a6 = -sinim * a9 + cosim * a10; 0391 0392 x1 = a1 * cosomm + a2 * sinomm; 0393 x2 = a3 * cosomm + a4 * sinomm; 0394 x3 = -a1 * sinomm + a2 * cosomm; 0395 x4 = -a3 * sinomm + a4 * cosomm; 0396 x5 = a5 * sinomm; 0397 x6 = a6 * sinomm; 0398 x7 = a5 * cosomm; 0399 x8 = a6 * cosomm; 0400 0401 z31 = 12.0 * x1 * x1 - 3.0 * x3 * x3; 0402 z32 = 24.0 * x1 * x2 - 6.0 * x3 * x4; 0403 z33 = 12.0 * x2 * x2 - 3.0 * x4 * x4; 0404 z1 = 3.0 * (a1 * a1 + a2 * a2) + z31 * emsq; 0405 z2 = 6.0 * (a1 * a3 + a2 * a4) + z32 * emsq; 0406 z3 = 3.0 * (a3 * a3 + a4 * a4) + z33 * emsq; 0407 z11 = -6.0 * a1 * a5 + emsq * (-24.0 * x1 * x7 - 6.0 * x3 * x5); 0408 z12 = -6.0 * (a1 * a6 + a3 * a5) + emsq * (-24.0 * (x2 * x7 + x1 * x8) - 6.0 * (x3 * x6 + x4 * x5)); 0409 z13 = -6.0 * a3 * a6 + emsq * (-24.0 * x2 * x8 - 6.0 * x4 * x6); 0410 z21 = 6.0 * a2 * a5 + emsq * (24.0 * x1 * x5 - 6.0 * x3 * x7); 0411 z22 = 6.0 * (a4 * a5 + a2 * a6) + emsq * (24.0 * (x2 * x5 + x1 * x6) - 6.0 * (x4 * x7 + x3 * x8)); 0412 z23 = 6.0 * a4 * a6 + emsq * (24.0 * x2 * x6 - 6.0 * x4 * x8); 0413 z1 = z1 + z1 + betasq * z31; 0414 z2 = z2 + z2 + betasq * z32; 0415 z3 = z3 + z3 + betasq * z33; 0416 s3 = cc * xnoi; 0417 s2 = -0.5 * s3 / rtemsq; 0418 s4 = s3 * rtemsq; 0419 s1 = -15.0 * em * s4; 0420 s5 = x1 * x3 + x2 * x4; 0421 s6 = x2 * x3 + x1 * x4; 0422 s7 = x2 * x4 - x1 * x3; 0423 0424 // Lunar terms 0425 if (lsflg == 1) 0426 { 0427 ss1 = s1; 0428 ss2 = s2; 0429 ss3 = s3; 0430 ss4 = s4; 0431 ss5 = s5; 0432 ss6 = s6; 0433 ss7 = s7; 0434 sz1 = z1; 0435 sz2 = z2; 0436 sz3 = z3; 0437 sz11 = z11; 0438 sz12 = z12; 0439 sz13 = z13; 0440 sz21 = z21; 0441 sz22 = z22; 0442 sz23 = z23; 0443 sz31 = z31; 0444 sz32 = z32; 0445 sz33 = z33; 0446 zcosg = zcosgl; 0447 zsing = zsingl; 0448 zcosi = zcosil; 0449 zsini = zsinil; 0450 zcosh = zcoshl * cnodm + zsinhl * snodm; 0451 zsinh = snodm * zcoshl - cnodm * zsinhl; 0452 cc = c1l; 0453 } 0454 } 0455 0456 zmol = fmod(4.7199672 + 0.22997150 * day - gam, TWOPI); 0457 zmos = fmod(6.2565837 + 0.017201977 * day, TWOPI); 0458 0459 // Solar terms 0460 se2 = 2.0 * ss1 * ss6; 0461 se3 = 2.0 * ss1 * ss7; 0462 si2 = 2.0 * ss2 * sz12; 0463 si3 = 2.0 * ss2 * (sz13 - sz11); 0464 sl2 = -2.0 * ss3 * sz2; 0465 sl3 = -2.0 * ss3 * (sz3 - sz1); 0466 sl4 = -2.0 * ss3 * (-21.0 - 9.0 * emsq) * zes; 0467 sgh2 = 2.0 * ss4 * sz32; 0468 sgh3 = 2.0 * ss4 * (sz33 - sz31); 0469 sgh4 = -18.0 * ss4 * zes; 0470 sh2 = -2.0 * ss2 * sz22; 0471 sh3 = -2.0 * ss2 * (sz23 - sz21); 0472 0473 // Lunar terms 0474 ee2 = 2.0 * s1 * s6; 0475 e3 = 2.0 * s1 * s7; 0476 xi2 = 2.0 * s2 * z12; 0477 xi3 = 2.0 * s2 * (z13 - z11); 0478 xl2 = -2.0 * s3 * z2; 0479 xl3 = -2.0 * s3 * (z3 - z1); 0480 xl4 = -2.0 * s3 * (-21.0 - 9.0 * emsq) * zel; 0481 xgh2 = 2.0 * s4 * z32; 0482 xgh3 = 2.0 * s4 * (z33 - z31); 0483 xgh4 = -18.0 * s4 * zel; 0484 xh2 = -2.0 * s2 * z22; 0485 xh3 = -2.0 * s2 * (z23 - z21); 0486 0487 // Apply deep space long period periodic contributions to the mean elements 0488 // double f2, f3, sinzf, zf, zm; 0489 double ses, sghl, sghs, shll, shs, sis, sls; 0490 0491 // Define some constants 0492 const double zns = 1.19459e-5; 0493 const double znl = 1.5835218e-4; 0494 0495 // Calculate time varying periodics 0496 // These results are never used, should we remove them? 0497 // zm = zmos; 0498 // zf = zm + 2.0 * zes * sin(zm); 0499 // sinzf = sin(zf); 0500 // f2 = 0.5 * sinzf * sinzf - 0.25; 0501 // f3 = -0.5 * sinzf * cos(zf); 0502 // ses = se2 * f2 + se3 * f3; 0503 // sis = si2 * f2 + si3 * f3; 0504 // sls = sl2 * f2 + sl3 * f3 + sl4 * sinzf; 0505 // sghs = sgh2 * f2 + sgh3 * f3 + sgh4 * sinzf; 0506 // shs = sh2 * f2 + sh3 * f3; 0507 // zm = zmol; 0508 0509 // zf = zm + 2.0 * zel * sin(zm); 0510 // sinzf = sin(zf); 0511 // f2 = 0.5 * sinzf * sinzf - 0.25; 0512 // f3 = -0.5 * sinzf * cos(zf); 0513 // sghl = xgh2 * f2 + xgh3 * f3 + xgh4 * sinzf; 0514 // shll = xh2 * f2 + xh3 * f3; 0515 0516 // Deep space contributions to mean motion dot due to geopotential resonance with half day and one day orbits 0517 double ainv2, aonv = 0.0, cosisq, eoc, f220, f221, f311, f321, f322, f330, f441, f442, f522, f523, f542, 0518 f543, g200, g201, g211, g300, g310, g322, g410, g422, g520, g521, g532, g533, sgs, sini2, 0519 temp, temp1, theta, xno2, emsqo; 0520 // double emo; 0521 0522 // Define some constant 0523 const double q22 = 1.7891679e-6; 0524 const double q31 = 2.1460748e-6; 0525 const double q33 = 2.2123015e-7; 0526 const double root22 = 1.7891679e-6; 0527 const double root44 = 7.3636953e-9; 0528 const double root54 = 2.1765803e-9; 0529 const double rptim = 4.37526908801129966e-3; // this equates to 7.29211514668855e-5 rad/sec 0530 const double root32 = 3.7393792e-7; 0531 const double root52 = 1.1428639e-7; 0532 0533 // Deep space initialization 0534 irez = 0; 0535 if ((nm < 0.0052359877) && (nm > 0.0034906585)) 0536 irez = 1; 0537 if ((nm >= 8.26e-3) && (nm <= 9.24e-3) && (em >= 0.5)) 0538 irez = 2; 0539 0540 // Solar terms 0541 ses = ss1 * zns * ss5; 0542 sis = ss2 * zns * (sz11 + sz13); 0543 sls = -zns * ss3 * (sz1 + sz3 - 14.0 - 6.0 * emsq); 0544 sghs = ss4 * zns * (sz31 + sz33 - 6.0); 0545 shs = -zns * ss2 * (sz21 + sz23); 0546 if ((inclm < 5.2359877e-2) || (inclm > M_PI - 5.2359877e-2)) 0547 shs = 0.0; 0548 if (sinim != 0.0) 0549 shs = shs / sinim; 0550 sgs = sghs - cosim * shs; 0551 0552 // Lunar terms 0553 dedt = ses + s1 * znl * s5; 0554 didt = sis + s2 * znl * (z11 + z13); 0555 dmdt = sls - znl * s3 * (z1 + z3 - 14.0 - 6.0 * emsq); 0556 sghl = s4 * znl * (z31 + z33 - 6.0); 0557 shll = -znl * s2 * (z21 + z23); 0558 if ((inclm < 5.2359877e-2) || (inclm > M_PI - 5.2359877e-2)) 0559 shll = 0.0; 0560 domdt = sgs + sghl; 0561 dnodt = shs; 0562 if (sinim != 0.0) 0563 { 0564 domdt = domdt - cosim / sinim * shll; 0565 dnodt = dnodt + shll / sinim; 0566 } 0567 0568 // Calculate deep space resonance effects 0569 // Value never used 0570 // dndt = 0.0; 0571 theta = fmod(gsto + tc * rptim, TWOPI); 0572 0573 // Initialize the resonance terms 0574 if (irez != 0) 0575 { 0576 aonv = pow(nm / XKE, X2O3); 0577 0578 // Geopotential resonance for 12 hour orbits 0579 if (irez == 2) 0580 { 0581 cosisq = cosim * cosim; 0582 // Value never used 0583 // emo = em; 0584 em = m_eccentricity; 0585 emsqo = emsq; 0586 emsq = eccsq; 0587 eoc = em * emsq; 0588 g201 = -0.306 - (em - 0.64) * 0.440; 0589 0590 if (em <= 0.65) 0591 { 0592 g211 = 3.616 - 13.2470 * em + 16.2900 * emsq; 0593 g310 = -19.302 + 117.3900 * em - 228.4190 * emsq + 156.5910 * eoc; 0594 g322 = -18.9068 + 109.7927 * em - 214.6334 * emsq + 146.5816 * eoc; 0595 g410 = -41.122 + 242.6940 * em - 471.0940 * emsq + 313.9530 * eoc; 0596 g422 = -146.407 + 841.8800 * em - 1629.014 * emsq + 1083.4350 * eoc; 0597 g520 = -532.114 + 3017.977 * em - 5740.032 * emsq + 3708.2760 * eoc; 0598 } 0599 else 0600 { 0601 g211 = -72.099 + 331.819 * em - 508.738 * emsq + 266.724 * eoc; 0602 g310 = -346.844 + 1582.851 * em - 2415.925 * emsq + 1246.113 * eoc; 0603 g322 = -342.585 + 1554.908 * em - 2366.899 * emsq + 1215.972 * eoc; 0604 g410 = -1052.797 + 4758.686 * em - 7193.992 * emsq + 3651.957 * eoc; 0605 g422 = -3581.690 + 16178.110 * em - 24462.770 * emsq + 12422.520 * eoc; 0606 if (em > 0.715) 0607 g520 = -5149.66 + 29936.92 * em - 54087.36 * emsq + 31324.56 * eoc; 0608 else 0609 g520 = 1464.74 - 4664.75 * em + 3763.64 * emsq; 0610 } 0611 if (em < 0.7) 0612 { 0613 g533 = -919.22770 + 4988.6100 * em - 9064.7700 * emsq + 5542.21 * eoc; 0614 g521 = -822.71072 + 4568.6173 * em - 8491.4146 * emsq + 5337.524 * eoc; 0615 g532 = -853.66600 + 4690.2500 * em - 8624.7700 * emsq + 5341.4 * eoc; 0616 } 0617 else 0618 { 0619 g533 = -37995.780 + 161616.52 * em - 229838.20 * emsq + 109377.94 * eoc; 0620 g521 = -51752.104 + 218913.95 * em - 309468.16 * emsq + 146349.42 * eoc; 0621 g532 = -40023.880 + 170470.89 * em - 242699.48 * emsq + 115605.82 * eoc; 0622 } 0623 0624 sini2 = sinim * sinim; 0625 f220 = 0.75 * (1.0 + 2.0 * cosim + cosisq); 0626 f221 = 1.5 * sini2; 0627 f321 = 1.875 * sinim * (1.0 - 2.0 * cosim - 3.0 * cosisq); 0628 f322 = -1.875 * sinim * (1.0 + 2.0 * cosim - 3.0 * cosisq); 0629 f441 = 35.0 * sini2 * f220; 0630 f442 = 39.3750 * sini2 * sini2; 0631 f522 = 0632 9.84375 * sinim * 0633 (sini2 * (1.0 - 2.0 * cosim - 5.0 * cosisq) + 0.33333333 * (-2.0 + 4.0 * cosim + 6.0 * cosisq)); 0634 f523 = sinim * (4.92187512 * sini2 * (-2.0 - 4.0 * cosim + 10.0 * cosisq) + 0635 6.56250012 * (1.0 + 2.0 * cosim - 3.0 * cosisq)); 0636 f542 = 29.53125 * sinim * (2.0 - 8.0 * cosim + cosisq * (-12.0 + 8.0 * cosim + 10.0 * cosisq)); 0637 f543 = 29.53125 * sinim * (-2.0 - 8.0 * cosim + cosisq * (12.0 + 8.0 * cosim - 10.0 * cosisq)); 0638 xno2 = nm * nm; 0639 ainv2 = aonv * aonv; 0640 temp1 = 3.0 * xno2 * ainv2; 0641 temp = temp1 * root22; 0642 d2201 = temp * f220 * g201; 0643 d2211 = temp * f221 * g211; 0644 temp1 = temp1 * aonv; 0645 temp = temp1 * root32; 0646 d3210 = temp * f321 * g310; 0647 d3222 = temp * f322 * g322; 0648 temp1 = temp1 * aonv; 0649 temp = 2.0 * temp1 * root44; 0650 d4410 = temp * f441 * g410; 0651 d4422 = temp * f442 * g422; 0652 temp1 = temp1 * aonv; 0653 temp = temp1 * root52; 0654 d5220 = temp * f522 * g520; 0655 d5232 = temp * f523 * g532; 0656 temp = 2.0 * temp1 * root54; 0657 d5421 = temp * f542 * g521; 0658 d5433 = temp * f543 * g533; 0659 xlamo = fmod(m_mean_anomaly + m_ra + m_ra - theta - theta, TWOPI); 0660 xfact = mdot + dmdt + 2.0 * (nodedot + dnodt - rptim) - m_mean_motion; 0661 // Value never used 0662 // em = emo; 0663 emsq = emsqo; 0664 } 0665 0666 if (irez == 1) 0667 { 0668 g200 = 1.0 + emsq * (-2.5 + 0.8125 * emsq); 0669 g310 = 1.0 + 2.0 * emsq; 0670 g300 = 1.0 + emsq * (-6.0 + 6.60937 * emsq); 0671 f220 = 0.75 * (1.0 + cosim) * (1.0 + cosim); 0672 f311 = 0.9375 * sinim * sinim * (1.0 + 3.0 * cosim) - 0.75 * (1.0 + cosim); 0673 f330 = 1.0 + cosim; 0674 f330 = 1.875 * f330 * f330 * f330; 0675 del1 = 3.0 * nm * nm * aonv * aonv; 0676 del2 = 2.0 * del1 * f220 * g200 * q22; 0677 del3 = 3.0 * del1 * f330 * g300 * q33 * aonv; 0678 del1 = del1 * f311 * g310 * q31 * aonv; 0679 xlamo = fmod(m_mean_anomaly + m_ra + m_arg_perigee - theta, TWOPI); 0680 xfact = mdot + xpidot - rptim + dmdt + domdt + dnodt - m_mean_motion; 0681 } 0682 0683 xli = xlamo; 0684 xni = m_mean_motion; 0685 atime = 0.0; 0686 // Value never used 0687 // nm = m_mean_motion + dndt; 0688 } 0689 } 0690 0691 // Set variables if not deep space 0692 if (!isimp) 0693 { 0694 cc1sq = cc1 * cc1; 0695 d2 = 4.0 * ao * tsi * cc1sq; 0696 temp = d2 * tsi * cc1 / 3.0; 0697 d3 = (17.0 * ao + sfour) * temp; 0698 d4 = 0.5 * temp * ao * tsi * (221.0 * ao + 31.0 * sfour) * cc1; 0699 t3cof = d2 + 2.0 * cc1sq; 0700 t4cof = 0.25 * (3.0 * d3 + cc1 * (12.0 * d2 + 10.0 * cc1sq)); 0701 t5cof = 0.2 * (3.0 * d4 + 12.0 * cc1 * d3 + 6.0 * d2 * d2 + 15.0 * cc1sq * (2.0 * d2 + cc1sq)); 0702 } 0703 } 0704 } 0705 0706 int Satellite::updatePos() 0707 { 0708 KStarsData *data = KStarsData::Instance(); 0709 return sgp4((data->clock()->utc().djd() - m_tle_jd) * MINPD); 0710 } 0711 0712 int Satellite::sgp4(double tsince) 0713 { 0714 KStarsData *data = KStarsData::Instance(); 0715 0716 int ktr; 0717 double am, axnl, aynl, betal, cosim, cnod, cos2u, coseo1 = 0, cosi, cosip, cosisq, cossu, cosu, delm, delomg, em, 0718 ecose, el2, eo1, ep, esine, argpm, argpp, argpdf, pl, 0719 mrt = 0.0, mvt, rdotl, rl, rvdot, rvdotl, sinim, dndt, sin2u, sineo1 = 0, sini, sinip, sinsu, sinu, snod, su, t2, 0720 t3, t4, tem5, temp, temp1, temp2, tempa, tempe, templ, u, ux, uy, uz, vx, vy, vz, inclm, mm, nm, nodem, xinc, 0721 xincp, xl, xlm, mp, xmdf, xmx, xmy, nodedf, xnode, nodep, tc, sat_posx, sat_posy, sat_posz, sat_posw, sat_velx, 0722 sat_vely, sat_velz, sinlat, obs_posx, obs_posy, obs_posz, obs_posw, /*obs_velx, obs_vely, obs_velz,*/ 0723 coslat, thetageo, sintheta, costheta, c, sq, achcp, vkmpersec; 0724 // double emsq; 0725 0726 const double temp4 = 1.5e-12; 0727 0728 double jul_utc = data->clock()->utc().djd(); 0729 0730 vkmpersec = RADIUSEARTHKM * XKE / 60.0; 0731 0732 // Update for secular gravity and atmospheric drag 0733 xmdf = m_mean_anomaly + mdot * tsince; 0734 argpdf = m_arg_perigee + argpdot * tsince; 0735 nodedf = m_ra + nodedot * tsince; 0736 argpm = argpdf; 0737 mm = xmdf; 0738 t2 = tsince * tsince; 0739 nodem = nodedf + nodecf * t2; 0740 tempa = 1.0 - cc1 * tsince; 0741 tempe = m_bstar * cc4 * tsince; 0742 templ = t2cof * t2; 0743 0744 if (!isimp) 0745 { 0746 delomg = omgcof * tsince; 0747 delm = xmcof * (pow((1.0 + eta * cos(xmdf)), 3) - delmo); 0748 temp = delomg + delm; 0749 mm = xmdf + temp; 0750 argpm = argpdf - temp; 0751 t3 = t2 * tsince; 0752 t4 = t3 * tsince; 0753 tempa = tempa - d2 * t2 - d3 * t3 - d4 * t4; 0754 tempe = tempe + m_bstar * cc5 * (sin(mm) - sinmao); 0755 templ = templ + t3cof * t3 + t4 * (t4cof + tsince * t5cof); 0756 } 0757 0758 nm = m_mean_motion; 0759 em = m_eccentricity; 0760 inclm = m_inclination; 0761 0762 if (method == 'd') 0763 { 0764 tc = tsince; 0765 // Deep space contributions to mean elements for perturbing third body 0766 int iretn; 0767 double delt, ft, theta, x2li, x2omi, xl, xldot, xnddt, xndt, xomi; 0768 0769 // Define some constants 0770 const double fasx2 = 0.13130908; 0771 const double fasx4 = 2.8843198; 0772 const double fasx6 = 0.37448087; 0773 const double g22 = 5.7686396; 0774 const double g32 = 0.95240898; 0775 const double g44 = 1.8014998; 0776 const double g52 = 1.0508330; 0777 const double g54 = 4.4108898; 0778 const double rptim = 4.37526908801129966e-3; // this equates to 7.29211514668855e-5 rad/sec 0779 const double step = 720.0; 0780 const double step2 = step * step / 2; 0781 0782 // Calculate deep space resonance effects 0783 // Value never used 0784 // dndt = 0.0; 0785 theta = fmod(gsto + tc * rptim, TWOPI); 0786 em = em + dedt * tsince; 0787 0788 inclm = inclm + didt * tsince; 0789 argpm = argpm + domdt * tsince; 0790 nodem = nodem + dnodt * tsince; 0791 mm = mm + dmdt * tsince; 0792 0793 // Update resonances : numerical (euler-maclaurin) integration 0794 ft = 0.0; 0795 if (irez != 0) 0796 { 0797 if ((atime == 0.0) || (tsince * atime <= 0.0) || (fabs(tsince) < fabs(atime))) 0798 { 0799 atime = 0.0; 0800 xni = m_mean_motion; 0801 xli = xlamo; 0802 } 0803 0804 if (tsince > 0.0) 0805 delt = step; 0806 else 0807 delt = -step; 0808 0809 iretn = 381; // added for do loop 0810 0811 while (iretn == 381) 0812 { 0813 // Near - synchronous resonance terms 0814 if (irez != 2) 0815 { 0816 xndt = del1 * sin(xli - fasx2) + del2 * sin(2.0 * (xli - fasx4)) + del3 * sin(3.0 * (xli - fasx6)); 0817 xldot = xni + xfact; 0818 xnddt = del1 * cos(xli - fasx2) + 2.0 * del2 * cos(2.0 * (xli - fasx4)) + 0819 3.0 * del3 * cos(3.0 * (xli - fasx6)); 0820 xnddt = xnddt * xldot; 0821 } 0822 else 0823 { 0824 // Near - half-day resonance terms 0825 xomi = m_arg_perigee + argpdot * atime; 0826 x2omi = xomi + xomi; 0827 x2li = xli + xli; 0828 xndt = d2201 * sin(x2omi + xli - g22) + d2211 * sin(xli - g22) + d3210 * sin(xomi + xli - g32) + 0829 d3222 * sin(-xomi + xli - g32) + d4410 * sin(x2omi + x2li - g44) + d4422 * sin(x2li - g44) + 0830 d5220 * sin(xomi + xli - g52) + d5232 * sin(-xomi + xli - g52) + 0831 d5421 * sin(xomi + x2li - g54) + d5433 * sin(-xomi + x2li - g54); 0832 xldot = xni + xfact; 0833 xnddt = d2201 * cos(x2omi + xli - g22) + d2211 * cos(xli - g22) + d3210 * cos(xomi + xli - g32) + 0834 d3222 * cos(-xomi + xli - g32) + d5220 * cos(xomi + xli - g52) + 0835 d5232 * cos(-xomi + xli - g52) + 0836 2.0 * (d4410 * cos(x2omi + x2li - g44) + d4422 * cos(x2li - g44) + 0837 d5421 * cos(xomi + x2li - g54) + d5433 * cos(-xomi + x2li - g54)); 0838 xnddt = xnddt * xldot; 0839 } 0840 0841 if (fabs(tsince - atime) >= step) 0842 { 0843 iretn = 381; 0844 } 0845 else 0846 { 0847 ft = tsince - atime; 0848 iretn = 0; 0849 } 0850 0851 if (iretn == 381) 0852 { 0853 xli = xli + xldot * delt + xndt * step2; 0854 xni = xni + xndt * delt + xnddt * step2; 0855 atime = atime + delt; 0856 } 0857 } 0858 0859 nm = xni + xndt * ft + xnddt * ft * ft * 0.5; 0860 xl = xli + xldot * ft + xndt * ft * ft * 0.5; 0861 0862 if (irez != 1) 0863 { 0864 mm = xl - 2.0 * nodem + 2.0 * theta; 0865 dndt = nm - m_mean_motion; 0866 } 0867 else 0868 { 0869 mm = xl - nodem - argpm + theta; 0870 dndt = nm - m_mean_motion; 0871 } 0872 0873 nm = m_mean_motion + dndt; 0874 } 0875 } 0876 0877 if (nm <= 0.0) 0878 { 0879 qDebug() << Q_FUNC_INFO << "Mean motion less than 0.0"; 0880 return (2); 0881 } 0882 0883 am = pow((XKE / nm), X2O3) * tempa * tempa; 0884 nm = XKE / pow(am, 1.5); 0885 em = em - tempe; 0886 0887 if ((em >= 1.0) || (em < -0.001)) 0888 { 0889 qDebug() << Q_FUNC_INFO << "Eccentricity >= 1.0 or < -0.001"; 0890 return (1); 0891 } 0892 0893 if (em < 1.0e-6) 0894 em = 1.0e-6; 0895 0896 mm = mm + m_mean_motion * templ; 0897 xlm = mm + argpm + nodem; 0898 // Value never used 0899 // emsq = em * em; 0900 // Value never used 0901 // temp = 1.0 - emsq; 0902 0903 nodem = fmod(nodem, TWOPI); 0904 argpm = fmod(argpm, TWOPI); 0905 xlm = fmod(xlm, TWOPI); 0906 mm = fmod(xlm - argpm - nodem, TWOPI); 0907 0908 // Compute extra mean quantities 0909 sinim = sin(inclm); 0910 cosim = cos(inclm); 0911 0912 // Add lunar-solar periodics 0913 ep = em; 0914 xincp = inclm; 0915 argpp = argpm; 0916 nodep = nodem; 0917 mp = mm; 0918 sinip = sinim; 0919 cosip = cosim; 0920 if (method == 'd') 0921 { 0922 double alfdp, betdp, cosip, cosop, dalf, dbet, dls, f2, f3, pe, pgh, ph, pinc, pl, sel, ses, sghl, sghs, shll, 0923 shs, sil, sinip, sinop, sinzf, sis, sll, sls, xls, xnoh, zf, zm; 0924 0925 // Define some constants 0926 const double zns = 1.19459e-5; 0927 const double zes = 0.01675; 0928 const double znl = 1.5835218e-4; 0929 const double zel = 0.05490; 0930 0931 // Calculate time varying periodics 0932 zm = zmos + zns * tsince; 0933 zf = zm + 2.0 * zes * sin(zm); 0934 sinzf = sin(zf); 0935 f2 = 0.5 * sinzf * sinzf - 0.25; 0936 f3 = -0.5 * sinzf * cos(zf); 0937 ses = se2 * f2 + se3 * f3; 0938 sis = si2 * f2 + si3 * f3; 0939 sls = sl2 * f2 + sl3 * f3 + sl4 * sinzf; 0940 sghs = sgh2 * f2 + sgh3 * f3 + sgh4 * sinzf; 0941 shs = sh2 * f2 + sh3 * f3; 0942 zm = zmol + znl * tsince; 0943 0944 zf = zm + 2.0 * zel * sin(zm); 0945 sinzf = sin(zf); 0946 f2 = 0.5 * sinzf * sinzf - 0.25; 0947 f3 = -0.5 * sinzf * cos(zf); 0948 sel = ee2 * f2 + e3 * f3; 0949 sil = xi2 * f2 + xi3 * f3; 0950 sll = xl2 * f2 + xl3 * f3 + xl4 * sinzf; 0951 sghl = xgh2 * f2 + xgh3 * f3 + xgh4 * sinzf; 0952 shll = xh2 * f2 + xh3 * f3; 0953 pe = ses + sel; 0954 pinc = sis + sil; 0955 pl = sls + sll; 0956 pgh = sghs + sghl; 0957 ph = shs + shll; 0958 0959 pe = pe - peo; 0960 pinc = pinc - pinco; 0961 pl = pl - plo; 0962 pgh = pgh - pgho; 0963 ph = ph - pho; 0964 xincp = xincp + pinc; 0965 ep = ep + pe; 0966 sinip = sin(xincp); 0967 cosip = cos(xincp); 0968 0969 // Apply periodics directly 0970 if (xincp >= 0.2) 0971 { 0972 ph = ph / sinip; 0973 pgh = pgh - cosip * ph; 0974 argpp = argpp + pgh; 0975 nodep = nodep + ph; 0976 mp = mp + pl; 0977 } 0978 else 0979 { 0980 // Apply periodics with lyddane modification 0981 sinop = sin(nodep); 0982 cosop = cos(nodep); 0983 alfdp = sinip * sinop; 0984 betdp = sinip * cosop; 0985 dalf = ph * cosop + pinc * cosip * sinop; 0986 dbet = -ph * sinop + pinc * cosip * cosop; 0987 alfdp = alfdp + dalf; 0988 betdp = betdp + dbet; 0989 nodep = fmod(nodep, TWOPI); 0990 if (nodep < 0.0) 0991 nodep += TWOPI; 0992 xls = mp + argpp + cosip * nodep; 0993 dls = pl + pgh - pinc * nodep * sinip; 0994 xls = xls + dls; 0995 xnoh = nodep; 0996 nodep = atan2(alfdp, betdp); 0997 if ((nodep < 0.0)) 0998 nodep += TWOPI; 0999 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 1232 QString 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 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 1294 bool Satellite::isVisible() 1295 { 1296 return m_is_visible; 1297 } 1298 1299 bool Satellite::selected() 1300 { 1301 return m_is_selected; 1302 } 1303 1304 void Satellite::setSelected(bool selected) 1305 { 1306 m_is_selected = selected; 1307 } 1308 1309 void Satellite::initPopupMenu(KSPopupMenu *pmenu) 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 1333 QString Satellite::id() const 1334 { 1335 return m_id; 1336 } 1337 1338 QString Satellite::tle() const 1339 { 1340 return m_tle; 1341 }