1 // SPDX-License-Identifier: LGPL-2.1-or-later
2 //
3 // SPDX-FileCopyrightText: 2011 Guillaume Martres <smarter@ubuntu.com>
4 //
5
6 #include "SatellitesTLEItem.h"
7
8 #include "MarbleClock.h"
9 #include "MarbleDebug.h"
10 #include "MarbleGlobal.h"
11 #include "GeoPainter.h"
12 #include "GeoDataCoordinates.h"
13 #include "GeoDataPlacemark.h"
14 #include "GeoDataStyle.h"
15 #include "GeoDataTrack.h"
16
17 #include <sgp4ext.h>
18
19 #include <QFile>
20 #include <QDateTime>
21 #include <QAction>
22 #include <QColor>
23
24 #include <cmath>
25
26 namespace Marble {
27
SatellitesTLEItem(const QString & name,elsetrec satrec,const MarbleClock * clock)28 SatellitesTLEItem::SatellitesTLEItem( const QString &name,
29 elsetrec satrec,
30 const MarbleClock *clock )
31 : TrackerPluginItem( name ),
32 m_satrec( satrec ),
33 m_track( new GeoDataTrack() ),
34 m_clock( clock )
35 {
36 double tumin, mu, xke, j2, j3, j4, j3oj2;
37 double radiusearthkm;
38 getgravconst( wgs84, tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 );
39 m_earthSemiMajorAxis = radiusearthkm;
40
41 setDescription();
42
43 placemark()->setVisualCategory(GeoDataPlacemark::Satellite);
44 placemark()->setZoomLevel( 0 );
45 placemark()->setGeometry( m_track );
46
47 update();
48 }
49
setDescription()50 void SatellitesTLEItem::setDescription()
51 {
52 QFile templateFile(QStringLiteral(":/marble/satellites/satellite.html"));
53 if (!templateFile.open(QIODevice::ReadOnly)) {
54 placemark()->setDescription(QObject::tr("No info available."));
55 return;
56 }
57 QString html = templateFile.readAll();
58
59 html.replace("%name%", name());
60 html.replace("%noradId%", QString::number(m_satrec.satnum));
61 html.replace("%perigee%", QString::number(perigee(), 'f', 2));
62 html.replace("%apogee%", QString::number(apogee(), 'f', 2));
63 html.replace("%inclination%", QString::number(inclination(), 'f', 2));
64 html.replace("%period%", QString::number(period(), 'f', 2));
65 html.replace("%semiMajorAxis%", QString::number(semiMajorAxis(), 'f', 2));
66
67 placemark()->setDescription( html );
68 }
69
update()70 void SatellitesTLEItem::update()
71 {
72 if( !isEnabled() ) {
73 return;
74 }
75
76 QDateTime startTime = m_clock->dateTime();
77 QDateTime endTime = startTime;
78 if( isTrackVisible() ) {
79 startTime = startTime.addSecs( -2 * 60 );
80 endTime = startTime.addSecs( period() );
81 }
82
83 m_track->removeBefore( startTime );
84 m_track->removeAfter( endTime );
85
86 addPointAt( m_clock->dateTime() );
87
88 // time interval between each point in the track, in seconds
89 double step = period() / 100.0;
90
91 for ( double i = startTime.toTime_t(); i < endTime.toTime_t(); i += step ) {
92 // No need to add points in this interval
93 if ( i >= m_track->firstWhen().toTime_t() ) {
94 i = m_track->lastWhen().toTime_t() + step;
95 }
96
97 addPointAt( QDateTime::fromTime_t( i ) );
98 }
99 }
100
addPointAt(const QDateTime & dateTime)101 void SatellitesTLEItem::addPointAt( const QDateTime &dateTime )
102 {
103 // in minutes
104 double timeSinceEpoch = (double)( dateTime.toTime_t() -
105 timeAtEpoch().toTime_t() ) / 60.0;
106
107 double r[3], v[3];
108 sgp4( wgs84, m_satrec, timeSinceEpoch, r, v );
109
110 GeoDataCoordinates coordinates = fromTEME(
111 r[0], r[1], r[2], gmst( timeSinceEpoch ) );
112 if ( m_satrec.error != 0 ) {
113 return;
114 }
115
116 m_track->addPoint( dateTime, coordinates);
117 }
118
timeAtEpoch() const119 QDateTime SatellitesTLEItem::timeAtEpoch() const
120 {
121 int year = m_satrec.epochyr + ( m_satrec.epochyr < 57 ? 2000 : 1900 );
122
123 int month, day, hours, minutes;
124 double seconds;
125 days2mdhms( year, m_satrec.epochdays, month, day, hours , minutes, seconds );
126
127 int ms = fmod(seconds * 1000.0, 1000.0);
128
129 return QDateTime( QDate( year, month, day ),
130 QTime( hours, minutes, (int)seconds, ms ),
131 Qt::UTC );
132 }
133
period() const134 double SatellitesTLEItem::period() const
135 {
136 // no := mean motion (rad / min)
137 return 60 * (2 * M_PI / m_satrec.no);
138 }
139
apogee() const140 double SatellitesTLEItem::apogee() const
141 {
142 return m_satrec.alta * m_earthSemiMajorAxis;
143 }
144
perigee() const145 double SatellitesTLEItem::perigee() const
146 {
147 return m_satrec.altp * m_earthSemiMajorAxis;
148 }
149
semiMajorAxis() const150 double SatellitesTLEItem::semiMajorAxis() const
151 {
152
153 return m_satrec.a * m_earthSemiMajorAxis;
154 }
155
inclination() const156 double SatellitesTLEItem::inclination() const
157 {
158 return m_satrec.inclo / M_PI * 180;
159 }
160
fromTEME(double x,double y,double z,double gmst) const161 GeoDataCoordinates SatellitesTLEItem::fromTEME( double x,
162 double y,
163 double z,
164 double gmst ) const
165 {
166 double lon = atan2( y, x );
167 // Rotate the angle by gmst (the origin goes from the vernal equinox
168 // point to the Greenwich Meridian)
169 lon = GeoDataCoordinates::normalizeLon( fmod(lon - gmst, 2 * M_PI) );
170
171 double lat = atan2( z, sqrt( x*x + y*y ) );
172
173 //TODO: determine if this is worth the extra precision
174 // Algorithm from https://celestrak.com/columns/v02n03/
175 //TODO: demonstrate it.
176 double a = m_earthSemiMajorAxis;
177 double planetRadius = sqrt( x*x + y*y );
178 double latp = lat;
179 double C;
180 for ( int i = 0; i < 3; i++ ) {
181 C = 1 / sqrt( 1 - square( m_satrec.ecco * sin( latp ) ) );
182 lat = atan2( z + a * C * square( m_satrec.ecco ) * sin( latp ), planetRadius );
183 }
184
185 double alt = planetRadius / cos( lat ) - a * C;
186
187 lat = GeoDataCoordinates::normalizeLat( lat );
188
189 return GeoDataCoordinates( lon, lat, alt * 1000 );
190 }
191
gmst(double minutesP) const192 double SatellitesTLEItem::gmst( double minutesP ) const
193 {
194 // Earth rotation rate in rad/min, from sgp4io.cpp
195 double rptim = 4.37526908801129966e-3;
196 return fmod( m_satrec.gsto + rptim * minutesP, 2 * M_PI );
197 }
198
square(double x)199 double SatellitesTLEItem::square( double x )
200 {
201 return x * x;
202 }
203
204 } // namespace Marble
205