VE7UWU/sdrangel
1
0
Fork 0
mirror of https://github.com/f4exb/sdrangel.git synced 2024-11-13 20:01:46 -05:00
Code Issues Projects Releases Wiki Activity
cc0e0c0f28
sdrangel/plugins/feature/satellitetracker/satellitetrackersgp4.cpp

536 lines
20 KiB
C++
Raw Blame History

///////////////////////////////////////////////////////////////////////////////////
// Copyright (C) 2021-2023 Jon Beniston, M7RCE <jon@beniston.com> //
// //
// This program is free software; you can redistribute it and/or modify //
// it under the terms of the GNU General Public License as published by //
// the Free Software Foundation as version 3 of the License, or //
// (at your option) any later version. //
// //
// This program is distributed in the hope that it will be useful, //
// but WITHOUT ANY WARRANTY; without even the implied warranty of //
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
// GNU General Public License V3 for more details. //
// //
// You should have received a copy of the GNU General Public License //
// along with this program. If not, see <http://www.gnu.org/licenses/>. //
///////////////////////////////////////////////////////////////////////////////////
#include <cmath>
#include <CoordTopocentric.h>
#include <CoordGeodetic.h>
#include <Observer.h>
#include <SGP4.h>
#include <QTimeZone>
#include "util/units.h"
#include "satellitetrackersgp4.h"
using namespace libsgp4;
// Convert QGP4 DateTime to Qt QDataTime
static QDateTime dateTimeToQDateTime(DateTime dt)
{
QDateTime qdt(QDate(dt.Year(), dt.Month(), dt.Day()), QTime(dt.Hour(), dt.Minute(), dt.Second(), (int)(dt.Microsecond()/1000.0)), Qt::UTC);
return qdt;
}
// Convert Qt QDataTime to QGP4 DateTime
static DateTime qDateTimeToDateTime(QDateTime qdt)
{
QDateTime utc = qdt.toUTC();
QDate date = utc.date();
QTime time = utc.time();
DateTime dt;
dt.Initialise(date.year(), date.month(), date.day(), time.hour(), time.minute(), time.second(), time.msec() * 1000);
return dt;
}
// Get ground track
// Throws SatelliteException, DecayedException and TleException
void getGroundTrack(QDateTime dateTime,
const QString& tle0, const QString& tle1, const QString& tle2,
int steps, bool forward,
QList<QGeoCoordinate *>& coordinates,
QList<QDateTime *>& coordinateDateTimes)
{
Tle tle = Tle(tle0.toStdString(), tle1.toStdString(), tle2.toStdString());
SGP4 sgp4(tle);
OrbitalElements ele(tle);
double periodMins;
double timeStep;
// For 3D map, we want to quantize to minutes, so we replace previous
// position data, rather than insert additional positions alongside the old
// which can result is the camera view jumping around
dateTime = QDateTime(dateTime.date(), QTime(dateTime.time().hour(), dateTime.time().minute()), dateTime.timeZone());
// Note 2D map doesn't support paths wrapping around Earth several times
// So we just have a slight overlap here, with the future track being longer
DateTime currentTime = qDateTimeToDateTime(dateTime);
DateTime endTime;
if (forward)
{
periodMins = ele.Period() * 0.9;
endTime = currentTime.AddMinutes(periodMins);
timeStep = periodMins / (steps * 0.9);
}
else
{
periodMins = ele.Period() * 0.4;
endTime = currentTime.AddMinutes(-periodMins);
timeStep = -periodMins / (steps * 0.4);
}
// Quantize time step to 30 seconds
timeStep *= 2.0;
if (timeStep > 0.0) {
timeStep = std::max(timeStep, 1.0);
} else if (timeStep < 0.0) {
timeStep = std::min(timeStep, -1.0);
}
timeStep = round(timeStep);
timeStep /= 2.0;
while ((forward && (currentTime < endTime)) || (!forward && (currentTime > endTime)))
{
// Calculate satellite position
Eci eci = sgp4.FindPosition(currentTime);
// Convert satellite position to geodetic coordinates (lat and long)
CoordGeodetic geo = eci.ToGeodetic();
QGeoCoordinate *coord = new QGeoCoordinate(Units::radiansToDegrees(geo.latitude),
Units::radiansToDegrees(geo.longitude),
geo.altitude * 1000.0);
coordinates.append(coord);
QDateTime *coordDateTime = new QDateTime(dateTimeToQDateTime(currentTime));
coordinateDateTimes.append(coordDateTime);
// 2D map is stretched at poles, so use finer steps
if (std::abs(Units::radiansToDegrees(geo.latitude)) >= 70)
currentTime = currentTime.AddMinutes(timeStep/4);
else
currentTime = currentTime.AddMinutes(timeStep);
}
}
// Find azimuth and elevation points during a pass
void getPassAzEl(QLineSeries* azimuth, QLineSeries* elevation, QLineSeries* polar,
const QString& tle0, const QString& tle1, const QString& tle2,
double latitude, double longitude, double altitude,
const QDateTime& aos, const QDateTime& los)
{
try
{
Tle tle = Tle(tle0.toStdString(), tle1.toStdString(), tle2.toStdString());
SGP4 sgp4(tle);
Observer obs(latitude, longitude, altitude);
DateTime aosTime = qDateTimeToDateTime(aos);
DateTime losTime = qDateTimeToDateTime(los);
DateTime currentTime(aosTime);
int steps = 150; // Needs to be high enough, so rotator intersect with satellite position
double timeStep = (losTime - aosTime).TotalSeconds() / steps;
if (timeStep <= 0.0)
{
qDebug() << "getPassAzEl: AOS is the same as or after LOS";
return;
}
while (currentTime <= losTime)
{
// Calculate satellite position
Eci eci = sgp4.FindPosition(currentTime);
// Calculate angle to satellite from antenna
CoordTopocentric topo = obs.GetLookAngle(eci);
// Save azimuth and elevation in series
QDateTime qdt = dateTimeToQDateTime(currentTime);
if (azimuth != nullptr)
azimuth->append(qdt.toMSecsSinceEpoch(), Units::radiansToDegrees(topo.azimuth));
if (elevation != nullptr)
elevation->append(qdt.toMSecsSinceEpoch(), Units::radiansToDegrees(topo.elevation));
if (polar != nullptr)
polar->append(Units::radiansToDegrees(topo.azimuth), 90.0-Units::radiansToDegrees(topo.elevation));
currentTime = currentTime.AddSeconds(timeStep);
}
}
catch (SatelliteException& se)
{
qDebug() << se.what();
}
catch (DecayedException& de)
{
qDebug() << de.what();
}
catch (TleException& tlee)
{
qDebug() << tlee.what();
}
}
// Get whether a pass passes through 0 degreees
bool getPassesThrough0Deg(const QString& tle0, const QString& tle1, const QString& tle2,
double latitude, double longitude, double altitude,
QDateTime& aos, QDateTime& los)
{
try
{
Tle tle = Tle(tle0.toStdString(), tle1.toStdString(), tle2.toStdString());
SGP4 sgp4(tle);
Observer obs(latitude, longitude, altitude);
DateTime aosTime = qDateTimeToDateTime(aos);
DateTime losTime = qDateTimeToDateTime(los);
DateTime currentTime(aosTime);
int steps = 20;
double timeStep = (losTime - aosTime).TotalSeconds() / steps;
double prevAz;
for (int i = 0; i < steps; i++)
{
// Calculate satellite position
Eci eci = sgp4.FindPosition(currentTime);
// Calculate angle to satellite from antenna
CoordTopocentric topo = obs.GetLookAngle(eci);
double az = Units::radiansToDegrees(topo.azimuth);
if (i == 0)
prevAz = az;
// Does it cross 0 degrees?
if (((prevAz > 270.0) && (az < 90.0)) || ((prevAz < 90.0) && (az >= 270.0)))
return true;
prevAz = az;
currentTime = currentTime.AddSeconds(timeStep);
}
}
catch (SatelliteException& se)
{
qDebug() << se.what();
}
catch (DecayedException& de)
{
qDebug() << de.what();
}
catch (TleException& tlee)
{
qDebug() << tlee.what();
}
return false;
}
// Find maximum elevation in a pass
static double findMaxElevation(Observer& obs1, SGP4& sgp4, const DateTime& aos, const DateTime& los)
{
Observer obs(obs1.GetLocation());
bool running;
double timeStep = (los - aos).TotalSeconds() / 9.0;
DateTime currentTime(aos);
DateTime time1(aos);
DateTime time2(los);
double maxElevation;
do
{
running = true;
maxElevation = -INFINITY;
while (running && (currentTime < time2))
{
Eci eci = sgp4.FindPosition(currentTime);
CoordTopocentric topo = obs.GetLookAngle(eci);
if (topo.elevation > maxElevation)
{
maxElevation = topo.elevation;
currentTime = currentTime.AddSeconds(timeStep);
if (currentTime > time2)
currentTime = time2;
}
else
running = false;
}
time1 = currentTime.AddSeconds(-2.0 * timeStep);
time2 = currentTime;
currentTime = time1;
timeStep = (time2 - time1).TotalSeconds() / 9.0;
}
while (timeStep > 1.0);
return Units::radiansToDegrees(maxElevation);
}
// Find the time at which the satellite crossed the minimum elevation required for AOS or LOS
static DateTime findCrossingPoint(Observer& obs, SGP4& sgp4, const DateTime& initialTime1, const DateTime& initialTime2, double minElevation, bool findingAOS)
{
bool running;
int cnt;
DateTime time1(initialTime1);
DateTime time2(initialTime2);
DateTime middleTime;
running = true;
cnt = 0;
while (running && (cnt++ < 16))
{
middleTime = time1.AddSeconds((time2 - time1).TotalSeconds() / 2.0);
Eci eci = sgp4.FindPosition(middleTime);
CoordTopocentric topo = obs.GetLookAngle(eci);
if (topo.elevation > minElevation)
{
if (findingAOS)
time2 = middleTime;
else
time1 = middleTime;
}
else
{
if (findingAOS)
time1 = middleTime;
else
time2 = middleTime;
}
if ((time2 - time1).TotalSeconds() < 1.0)
{
running = false;
int us = middleTime.Microsecond();
middleTime = middleTime.AddMicroseconds(-us);
middleTime = middleTime.AddSeconds(findingAOS ? 1 : -1);
}
}
running = true;
cnt = 0;
while (running && (cnt++ < 6))
{
Eci eci = sgp4.FindPosition(middleTime);
CoordTopocentric topo = obs.GetLookAngle(eci);
if (topo.elevation > minElevation)
middleTime = middleTime.AddSeconds(findingAOS ? -1 : 1);
else
running = false;
}
return middleTime;
}
// Find when AOS occured, by stepping backwards
static DateTime findAOSBackwards(Observer& obs, SGP4& sgp4, DateTime& startTime,
int predictionPeriod, double minElevation, bool& aosUnknown)
{
DateTime previousTime(startTime);
DateTime currentTime(startTime);
DateTime endTime(startTime.AddDays(-predictionPeriod));
while (currentTime >= endTime)
{
Eci eci = sgp4.FindPosition(currentTime);
CoordTopocentric topo = obs.GetLookAngle(eci);
if (topo.elevation < minElevation)
{
aosUnknown = false;
return findCrossingPoint(obs, sgp4, currentTime, previousTime, minElevation, true);
}
previousTime = currentTime;
currentTime = currentTime - TimeSpan(0, 0, 180);
}
aosUnknown = true;
return currentTime;
}
bool inPassWindow(DateTime dateTime, QTime passStartTime, QTime passEndTime, bool utc)
{
// Don't compare seconds as not currently settable in GUI
QDateTime qdt = dateTimeToQDateTime(dateTime);
if (!utc)
qdt = qdt.toLocalTime();
QTime qt(qdt.time().hour(), qdt.time().minute());
passStartTime = QTime(passStartTime.hour(), passStartTime.minute());
passEndTime = QTime(passEndTime.hour(), passEndTime.minute());
// If passEndTime is before passStartTime, then we allow overnight passes
if (passEndTime > passStartTime)
{
return (qt >= passStartTime) && (qt <= passEndTime);
}
else
{
return (qt <= passEndTime) || (qt >= passStartTime);
}
}
// Create a list of satellite passes, between the given start and end times, that exceed the specified minimum elevation
// We return an uninitalised QDateTime if AOS or LOS is outside of predictionPeriod
static QList<SatellitePass> createPassList(Observer& obs, SGP4& sgp4, DateTime& startTime,
int predictionPeriod, double minAOSElevation, double minPassElevationDeg,
QTime passStartTime, QTime passEndTime, bool utc,
int noOfPasses)
{
QList<SatellitePass> passes;
bool aos = false;
bool aosUnknown = true;
double aosAz;
double losAz;
DateTime previousTime(startTime);
DateTime currentTime(startTime);
DateTime endTime(startTime.AddDays(predictionPeriod));
DateTime aosTime;
DateTime losTime;
while (currentTime < endTime)
{
bool endOfPass = false;
Eci eci = sgp4.FindPosition(currentTime);
CoordTopocentric topo = obs.GetLookAngle(eci);
if (!aos && (topo.elevation > minAOSElevation))
{
if (startTime == currentTime)
{
// AOS is before startTime
aosTime = findAOSBackwards(obs, sgp4, startTime, predictionPeriod, minAOSElevation, aosUnknown);
}
else
{
aosTime = findCrossingPoint(obs, sgp4, previousTime, currentTime, minAOSElevation, true);
aosUnknown = false;
}
aos = true;
eci = sgp4.FindPosition(aosTime);
topo = obs.GetLookAngle(eci);
aosAz = Units::radiansToDegrees(topo.azimuth);
}
else if (aos && (topo.elevation < minAOSElevation))
{
aos = false;
endOfPass = true;
losTime = findCrossingPoint(obs, sgp4, previousTime, currentTime, minAOSElevation, false);
eci = sgp4.FindPosition(losTime);
topo = obs.GetLookAngle(eci);
losAz = Units::radiansToDegrees(topo.azimuth);
double maxElevationDeg = findMaxElevation(obs, sgp4, aosTime, losTime);
if ((maxElevationDeg >= minPassElevationDeg)
&& inPassWindow(aosTime, passStartTime, passEndTime, utc)
&& inPassWindow(losTime, passStartTime, passEndTime, utc))
{
SatellitePass pass;
pass.m_aos = aosUnknown ? QDateTime() : dateTimeToQDateTime(aosTime);
pass.m_los = dateTimeToQDateTime(losTime);
pass.m_maxElevation = maxElevationDeg;
pass.m_aosAzimuth = aosAz;
pass.m_losAzimuth = losAz;
pass.m_northToSouth = std::min(360.0-aosAz, aosAz-0.0) < std::min(360.0-losAz, losAz-0.0);
passes.append(pass);
noOfPasses--;
if (noOfPasses <= 0)
return passes;
}
}
previousTime = currentTime;
if (endOfPass)
currentTime = currentTime + TimeSpan(0, 30, 0); // 30 minutes - no orbit likely to be that fast
else
currentTime = currentTime + TimeSpan(0, 0, 180);
if (currentTime > endTime)
currentTime = endTime;
}
if (aos)
{
// Pass still in progress at end time
Eci eci = sgp4.FindPosition(currentTime);
CoordTopocentric topo = obs.GetLookAngle(eci);
losAz = Units::radiansToDegrees(topo.azimuth);
double maxElevationDeg = findMaxElevation(obs, sgp4, aosTime, losTime);
if ((maxElevationDeg >= minPassElevationDeg)
&& inPassWindow(aosTime, passStartTime, passEndTime, utc)
&& inPassWindow(losTime, passStartTime, passEndTime, utc))
{
SatellitePass pass;
pass.m_aos = aosUnknown ? QDateTime() : dateTimeToQDateTime(aosTime);
pass.m_los = QDateTime();
pass.m_aosAzimuth = aosAz;
pass.m_losAzimuth = losAz;
pass.m_maxElevation = maxElevationDeg;
pass.m_northToSouth = std::min(360.0-aosAz, aosAz-0.0) < std::min(360.0-losAz, losAz-0.0);
passes.append(pass);
}
}
return passes;
}
void getSatelliteState(QDateTime dateTime,
const QString& tle0, const QString& tle1, const QString& tle2,
double latitude, double longitude, double altitude,
int predictionPeriod, int minAOSElevationDeg, int minPassElevationDeg,
QTime passStartTime, QTime passFinishTime, bool utc,
int noOfPasses, int groundTrackSteps, SatelliteState *satState)
{
try
{
Tle tle = Tle(tle0.toStdString(), tle1.toStdString(), tle2.toStdString());
SGP4 sgp4(tle);
Observer obs(latitude, longitude, altitude);
DateTime dt = qDateTimeToDateTime(dateTime);
// Calculate satellite position
Eci eci = sgp4.FindPosition(dt);
// Calculate angle to satellite from antenna
CoordTopocentric topo = obs.GetLookAngle(eci);
// Convert satellite position to geodetic coordinates (lat and long)
CoordGeodetic geo = eci.ToGeodetic();
satState->m_latitude = Units::radiansToDegrees(geo.latitude);
satState->m_longitude = Units::radiansToDegrees(geo.longitude);
satState->m_altitude = geo.altitude;
satState->m_azimuth = Units::radiansToDegrees(topo.azimuth);
satState->m_elevation = Units::radiansToDegrees(topo.elevation);
satState->m_range = topo.range;
satState->m_rangeRate = topo.range_rate;
OrbitalElements ele(tle);
satState->m_speed = eci.Velocity().Magnitude();
satState->m_period = ele.Period();
if (noOfPasses > 0)
{
satState->m_passes.clear();
satState->m_passes = createPassList(obs, sgp4, dt, predictionPeriod,
Units::degreesToRadians((double)minAOSElevationDeg),
minPassElevationDeg,
passStartTime, passFinishTime, utc,
noOfPasses);
}
qDeleteAll(satState->m_groundTrack);
satState->m_groundTrack.clear();
qDeleteAll(satState->m_groundTrackDateTime);
satState->m_groundTrackDateTime.clear();
qDeleteAll(satState->m_predictedGroundTrack);
satState->m_predictedGroundTrack.clear();
qDeleteAll(satState->m_predictedGroundTrackDateTime);
satState->m_predictedGroundTrackDateTime.clear();
getGroundTrack(dateTime, tle0, tle1, tle2, groundTrackSteps, false, satState->m_groundTrack, satState->m_groundTrackDateTime);
getGroundTrack(dateTime, tle0, tle1, tle2, groundTrackSteps, true, satState->m_predictedGroundTrack, satState->m_predictedGroundTrackDateTime);
}
catch (SatelliteException& se)
{
qDebug() << "getSatelliteState:SatelliteException " << satState->m_name << ": " << se.what();
}
catch (DecayedException& de)
{
qDebug() << "getSatelliteState:DecayedException " << satState->m_name << ": " << de.what();
}
catch (TleException& tlee)
{
qDebug() << "getSatelliteState:TleException " << satState->m_name << ": " << tlee.what() << "\n" << tle0 << "\n" << tle1 << "\n" << tle2;
}
}
Reference in New Issue View Git Blame Copy Permalink