1 /***************************************************************************
2 * Name: gSatWrapper.hpp
3 *
4 * Description: Wrapper over gSatTEME class.
5 * This class allow use Satellite orbit calculation module (gSAt) in
6 * Stellarium 'native' mode using Stellarium objects.
7 *
8 ***************************************************************************/
9 /***************************************************************************
10 * Copyright (C) 2006 by J.L. Canales *
11 * jlcanales.gasco@gmail.com *
12 * *
13 * This program is free software; you can redistribute it and/or modify *
14 * it under the terms of the GNU General Public License as published by *
15 * the Free Software Foundation; either version 2 of the License, or *
16 * (at your option) any later version. *
17 * *
18 * This program is distributed in the hope that it will be useful, *
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
21 * GNU General Public License for more details. *
22 * *
23 * You should have received a copy of the GNU General Public License *
24 * along with this program; if not, write to the *
25 * Free Software Foundation, Inc., *
26 * 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA. *
27 ***************************************************************************/
28
29 #include "gsatellite/stdsat.h"
30 #include "gsatellite/mathUtils.hpp"
31
32 #include "gSatWrapper.hpp"
33 #include "StelApp.hpp"
34 #include "StelCore.hpp"
35 #include "StelUtils.hpp"
36
37 #include "SolarSystem.hpp"
38 #include "StelModuleMgr.hpp"
39
40 #include <QDebug>
41 #include <QByteArray>
42
gSatWrapper(QString designation,QString tle1,QString tle2)43 gSatWrapper::gSatWrapper(QString designation, QString tle1,QString tle2)
44 {
45 // The TLE library actually modifies the TLE strings, which is annoying (because
46 // when we get updates, we want to check if there has been a change by using ==
47 // with the original. Thus we make a copy to send to the TLE library.
48 QByteArray t1(tle1.toLatin1().data()), t2(tle2.toLatin1().data());
49
50 // Also, the TLE library expects no more than 130 characters length input. We
51 // shouldn't have sane input with a TLE longer than about 80, but just in case
52 // we have a mal-formed input, we will truncate here to be safe
53 t1.truncate(130);
54 t2.truncate(130);
55
56 pSatellite = new gSatTEME(designation.toLatin1().data(), t1.data(), t2.data());
57 setEpoch(StelApp::getInstance().getCore()->getJD());
58 }
59
~gSatWrapper()60 gSatWrapper::~gSatWrapper()
61 {
62 if (pSatellite != Q_NULLPTR)
63 delete pSatellite;
64 }
65
getTEMEPos() const66 Vec3d gSatWrapper::getTEMEPos() const
67 {
68 Vec3d returnedVector;
69 if (pSatellite != Q_NULLPTR)
70 returnedVector = pSatellite->getPos();
71 else
72 qWarning() << "gSatWrapper::getTEMEPos Method called without pSatellite initialized";
73
74 return returnedVector;
75 }
76
getTEMEVel() const77 Vec3d gSatWrapper::getTEMEVel() const
78 {
79 Vec3d returnedVector;
80 if (pSatellite != Q_NULLPTR)
81 returnedVector = pSatellite->getVel();
82 else
83 qWarning() << "gSatWrapper::getTEMEVel Method called without pSatellite initialized";
84
85 return returnedVector;
86 }
87
getSubPoint() const88 Vec3d gSatWrapper::getSubPoint() const
89 {
90 Vec3d returnedVector;
91 if (pSatellite != Q_NULLPTR)
92 returnedVector = pSatellite->getSubPoint();
93 else
94 qWarning() << "gSatWrapper::getTEMEVel Method called without pSatellite initialized";
95
96 return returnedVector;
97 }
98
setEpoch(double ai_julianDaysEpoch)99 void gSatWrapper::setEpoch(double ai_julianDaysEpoch)
100 {
101 epoch = ai_julianDaysEpoch;
102 if (pSatellite)
103 pSatellite->setEpoch(epoch);
104 }
105
calcObserverECIPosition(Vec3d & ao_position,Vec3d & ao_velocity)106 void gSatWrapper::calcObserverECIPosition(Vec3d& ao_position, Vec3d& ao_velocity)
107 {
108 if (epoch != lastCalcObserverECIPosition)
109 {
110 StelLocation loc = StelApp::getInstance().getCore()->getCurrentLocation();
111
112 double radLatitude = loc.latitude * KDEG2RAD;
113 double theta = epoch.toThetaLMST(loc.longitude * KDEG2RAD);
114 double r;
115 double c,sq;
116
117 /* Reference: Explanatory supplement to the Astronomical Almanac 1992, page 209-210. */
118 /* Elipsoid earth model*/
119 /* c = Nlat/a */
120 c = 1/std::sqrt(1 + __f*(__f - 2)*Sqr(sin(radLatitude)));
121 sq = Sqr(1 - __f)*c;
122
123 r = (KEARTHRADIUS*c + (loc.altitude/1000))*cos(radLatitude);
124 ao_position[0] = r * cos(theta);/*kilometers*/
125 ao_position[1] = r * sin(theta);
126 ao_position[2] = (KEARTHRADIUS*sq + (loc.altitude/1000))*sin(radLatitude);
127 ao_velocity[0] = -KMFACTOR*ao_position[1];/*kilometers/second*/
128 ao_velocity[1] = KMFACTOR*ao_position[0];
129 ao_velocity[2] = 0;
130
131 lastCalcObserverECIPosition=epoch;
132 }
133 }
134
getAltAz() const135 Vec3d gSatWrapper::getAltAz() const
136 {
137 StelLocation loc = StelApp::getInstance().getCore()->getCurrentLocation();
138 Vec3d topoSatPos;
139
140 const double radLatitude = loc.latitude * KDEG2RAD;
141 const double theta = epoch.toThetaLMST(loc.longitude * KDEG2RAD);
142 const double sinRadLatitude = sin(radLatitude);
143 const double cosRadLatitude = cos(radLatitude);
144 const double sinTheta = sin(theta);
145 const double cosTheta = cos(theta);
146
147 // This now only updates if required.
148 calcObserverECIPosition(observerECIPos, observerECIVel);
149
150 Vec3d satECIPos = getTEMEPos();
151 Vec3d slantRange = satECIPos - observerECIPos;
152
153 //top_s
154 topoSatPos[0] = (sinRadLatitude * cosTheta*slantRange[0]
155 + sinRadLatitude* sinTheta*slantRange[1]
156 - cosRadLatitude* slantRange[2]);
157 //top_e
158 topoSatPos[1] = ((-1.0)* sinTheta*slantRange[0]
159 + cosTheta*slantRange[1]);
160
161 //top_z
162 topoSatPos[2] = (cosRadLatitude * cosTheta*slantRange[0]
163 + cosRadLatitude * sinTheta*slantRange[1]
164 + sinRadLatitude *slantRange[2]);
165
166 return topoSatPos;
167 }
168
getSlantRange(double & ao_slantRange,double & ao_slantRangeRate) const169 void gSatWrapper::getSlantRange(double &ao_slantRange, double &ao_slantRangeRate) const
170 {
171 calcObserverECIPosition(observerECIPos, observerECIVel);
172
173 Vec3d satECIPos = getTEMEPos();
174 Vec3d satECIVel = getTEMEVel();
175 Vec3d slantRange = satECIPos - observerECIPos;
176 Vec3d slantRangeVelocity = satECIVel - observerECIVel;
177
178 ao_slantRange = slantRange.length();
179 ao_slantRangeRate = slantRange.dot(slantRangeVelocity)/ao_slantRange;
180 }
181
182 // Does the actual computation only if necessary (0.1s apart) and caches the result in a static variable.
updateSunECIPos()183 void gSatWrapper::updateSunECIPos()
184 {
185 // All positions in ECI system are positions referenced in a StelCore::EquinoxEq system centered in the earth centre
186 calcObserverECIPosition(observerECIPos, observerECIVel);
187
188 static const SolarSystem *solsystem = (SolarSystem*)StelApp::getInstance().getModuleMgr().getModule("SolarSystem");
189 Vec3d sunEquinoxEqPos = solsystem->getSun()->getEquinoxEquatorialPos(StelApp::getInstance().getCore());
190
191 //sunEquinoxEqPos is measured in AU. we need measure it in Km
192 sunECIPos.set(sunEquinoxEqPos[0]*AU, sunEquinoxEqPos[1]*AU, sunEquinoxEqPos[2]*AU);
193 sunECIPos = sunECIPos + observerECIPos; //Change ref system centre
194 }
195
getSunECIPos()196 Vec3d gSatWrapper::getSunECIPos()
197 {
198 if (epoch != lastSunECIepoch)
199 {
200 updateSunECIPos();
201 lastSunECIepoch=epoch;
202 }
203 return sunECIPos;
204 }
205
206 // Operation getVisibilityPredict
207 // @brief This operation predicts the satellite visibility conditions.
getVisibilityPredict() const208 gSatWrapper::Visibility gSatWrapper::getVisibilityPredict() const
209 {
210 Vec3d satAltAzPos = getAltAz();
211
212 if (satAltAzPos[2] > 0)
213 {
214 Vec3d satECIPos = getTEMEPos();
215 static const SolarSystem *solsystem = (SolarSystem*)StelApp::getInstance().getModuleMgr().getModule("SolarSystem");
216 Vec3d sunAltAzPos = solsystem->getSun()->getAltAzPosGeometric(StelApp::getInstance().getCore());
217 Vec3d sunECIPos = getSunECIPos();
218
219 if (sunAltAzPos[2] > 0.0)
220 {
221 return RADAR_SUN;
222 }
223 else
224 {
225 double sunSatAngle = sunECIPos.angle(satECIPos);
226 double Dist = satECIPos.length()*cos(sunSatAngle - (M_PI/2));
227
228 if (Dist > KEARTHRADIUS)
229 {
230 return VISIBLE;
231 }
232 else
233 {
234 return RADAR_NIGHT;
235 }
236 }
237 }
238 else
239 return NOT_VISIBLE;
240 }
241
getPhaseAngle() const242 double gSatWrapper::getPhaseAngle() const
243 {
244 Vec3d sunECIPos = getSunECIPos();
245 return sunECIPos.angle(getTEMEPos());
246 }
247
getOrbitalPeriod() const248 double gSatWrapper::getOrbitalPeriod() const
249 {
250 return pSatellite->getPeriod();
251 }
252
getOrbitalInclination() const253 double gSatWrapper::getOrbitalInclination() const
254 {
255 return pSatellite->getInclination();
256 }
257
getPerigeeApogeeAltitudes() const258 Vec2d gSatWrapper::getPerigeeApogeeAltitudes() const
259 {
260 return pSatellite->getPerigeeApogee();
261 }
262
263 gTime gSatWrapper::epoch;
264 gTime gSatWrapper::lastSunECIepoch=0.0; // store last time of computation to avoid all-1 computations.
265 gTime gSatWrapper::lastCalcObserverECIPosition;
266
267 Vec3d gSatWrapper::sunECIPos; // enough to have this once.
268 Vec3d gSatWrapper::observerECIPos;
269 Vec3d gSatWrapper::observerECIVel;
270