KJ135
/
sgp4
mirror of http://172.16.0.12:5000/LAMPSARToolSoftware/sgp4.git
4
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Update after major changes

feature/19
Daniel Warner 2012-10-11 19:52:51 +01:00
parent 69744793e6
commit 2b5e39efe2
37 changed files with 8072 additions and 1826 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.hgignore
View File

@ -1,7 +1,7 @@
# Exact paths to config junk
^INSTALL$
^aclocal.m4$
^configure$
#^configure$
^(config.guess|config.sub|depcomp|install-sh|ltmain.sh|missing)$
^config.(log|status)$
^autom4te.cache$
@ -14,6 +14,7 @@
.(deps|libs)/*$
.cmp$
autom4te.cache/*$
doxygen/*$
Makefile(.in)?$
.(deps|libs)$

2
Makefile.am
View File

@ -1 +1 @@
SUBDIRS = libsgp4 passpredict runtest sattrack
SUBDIRS = libsgp4 sattrack runtest passpredict

25
config.h.in Normal file
View File

@ -0,0 +1,25 @@
/* config.h.in. Generated from configure.ac by autoheader. */
/* Name of package */
#undef PACKAGE
/* Define to the address where bug reports for this package should be sent. */
#undef PACKAGE_BUGREPORT
/* Define to the full name of this package. */
#undef PACKAGE_NAME
/* Define to the full name and version of this package. */
#undef PACKAGE_STRING
/* Define to the one symbol short name of this package. */
#undef PACKAGE_TARNAME
/* Define to the home page for this package. */
#undef PACKAGE_URL
/* Define to the version of this package. */
#undef PACKAGE_VERSION
/* Version number of package */
#undef VERSION

5406
configure vendored Executable file
View File

File diff suppressed because it is too large Load Diff

22
configure.ac
View File

@ -2,8 +2,10 @@ AC_PREREQ([2.61])
AC_INIT([sgp4], 1.0, [info@danrw.com])
AM_INIT_AUTOMAKE
AM_SILENT_RULES([yes])
#AC_CANONICAL_BUILD
AC_CONFIG_SRCDIR([libsgp4/Tle.cpp])
AM_CONFIG_HEADER([config.h])
CFLAGS=" $CFLAGS"
CXXFLAGS=" $CXXFLAGS"
@ -19,11 +21,13 @@ AC_ARG_ENABLE(debug,
enable_debug=no)
if test x$enable_debug = xyes; then
AM_CXXFLAGS="-g -O0 -Wall -Werror"
AM_CXXFLAGS="-g -O0 -Werror -W -Wall -ansi -pedantic -Wpointer-arith -Wwrite-strings -Wno-long-long"
else
AM_CXXFLAGS="-DNDEBUG -O2 -Wall -Werror -fomit-frame-pointer"
AM_CXXFLAGS="-DNDEBUG -O2 -fomit-frame-pointer -Werror -W -Wall -ansi -pedantic -Wpointer-arith -Wwrite-strings -Wno-long-long"
fi
AC_SEARCH_LIBS(clock_gettime, [rt])
AC_SUBST(AM_CXXFLAGS)
AC_CONFIG_FILES([Makefile
@ -33,3 +37,17 @@ AC_CONFIG_FILES([Makefile
sattrack/Makefile])
AC_OUTPUT
echo "---"
echo "Configuration summary for $PACKAGE_NAME version $VERSION"
echo ""
echo " * Installation prefix: $prefix"
echo " * System type: $build_vendor-$build_os"
echo " * Host CPU: $build_cpu"
echo " C Compiler: ${CC}"
echo " C++ Compiler: ${CXX}"
echo " Compiler flags: ${AM_CXXFLAGS}"
echo " Linker flags: ${AM_LDFLAGS}"
echo " Libraries: ${LIBS}"
echo " Debug enabled: $enable_debug"
echo ""
echo "---"

128
libsgp4/CoordGeodetic.h
View File

@ -1,106 +1,129 @@
#ifndef COORDGEODETIC_H_
#define COORDGEODETIC_H_
#include "Globals.h"
#include "Util.h"
#include <string>
#include <sstream>
#include <iomanip>
/**
* Stores a geodetic position
*/
struct CoordGeodetic
{
public:
CoordGeodetic()
: latitude(0.0), longitude(0.0), altitude(0.0)
{
}
/*
* default is in degrees
/**
* Default constructor
*/
CoordGeodetic(double lat, double lon, double alt, bool radians = false)
CoordGeodetic()
: latitude(0.0),
longitude(0.0),
altitude(0.0)
{
if (radians)
{
latitude = lat;
longitude = lon;
}
else
{
latitude = Util::DegreesToRadians(lat);
longitude = Util::DegreesToRadians(lon);
}
altitude = alt;
}
CoordGeodetic(const CoordGeodetic& g)
/**
* Constructor
* @param[in] arg_latitude the latitude in degrees
* @param[in] arg_longitude the longitude in degrees
* @param[in] arg_altitude the altitude in kilometers
*/
CoordGeodetic(
double arg_latitude,
double arg_longitude,
double arg_altitude)
{
latitude = g.latitude;
longitude = g.longitude;
altitude = g.altitude;
latitude = Util::DegreesToRadians(arg_latitude);
longitude = Util::DegreesToRadians(arg_longitude);
altitude = arg_altitude;
}
/**
* Copy constructor
* @param[in] geo object to copy from
*/
CoordGeodetic(const CoordGeodetic& geo)
{
latitude = geo.latitude;
longitude = geo.longitude;
altitude = geo.altitude;
}
/**
* Destructor
*/
virtual ~CoordGeodetic()
{
}
CoordGeodetic& operator=(const CoordGeodetic& g)
/**
* Assignment operator
* @param[in] geo object to copy from
*/
CoordGeodetic& operator=(const CoordGeodetic& geo)
{
if (this != &g)
if (this != &geo)
{
latitude = g.latitude;
longitude = g.longitude;
altitude = g.altitude;
latitude = geo.latitude;
longitude = geo.longitude;
altitude = geo.altitude;
}
return *this;
}
bool operator==(const CoordGeodetic& g) const
/**
* Equality operator
* @param[in] geo the object to compare with
* @returns whether the object is equal
*/
bool operator==(const CoordGeodetic& geo) const
{
return IsEqual(g);
return IsEqual(geo);
}
bool operator!=(const CoordGeodetic& g) const
/**
* Inequality operator
* @param[in] geo the object to compare with
* @returns whether the object is not equal
*/
bool operator!=(const CoordGeodetic& geo) const
{
return !IsEqual(g);
return !IsEqual(geo);
}
/**
* Dump this object to a string
* @returns string
*/
std::string ToString() const
{
std::stringstream ss;
ss << std::right << std::fixed << std::setprecision(2);
ss << std::right << std::fixed << std::setprecision(3);
ss << "Lat: " << std::setw(7) << Util::RadiansToDegrees(latitude);
ss << ", Lon: " << std::setw(7) << Util::RadiansToDegrees(longitude);
ss << ", Alt: " << std::setw(9) << altitude;
return ss.str();
}
/*
* radians (north positive, south negative)
*/
/** latitude in radians (-PI >= latitude < PI) */
double latitude;
/*
* radians (east positive, west negative)
*/
/** latitude in radians (-PI/2 >= latitude <= PI/2) */
double longitude;
/*
* kilometers
*/
/** altitude in kilometers */
double altitude;
protected:
bool IsEqual(const CoordGeodetic& g) const
private:
bool IsEqual(const CoordGeodetic& geo) const
{
if (latitude == g.latitude && longitude == g.longitude &&
altitude == g.altitude)
bool equal = false;
if (latitude == geo.latitude &&
longitude == geo.longitude &&
altitude == geo.altitude)
{
return false;
}
else
{
return true;
equal = false;
}
return equal;
}
};
@ -110,4 +133,3 @@ inline std::ostream& operator<<(std::ostream& strm, const CoordGeodetic& g)
}
#endif

139
libsgp4/CoordTopographic.h
View File

@ -1,99 +1,139 @@
#ifndef COORDTOPOGRAPHIC_H_
#define COORDTOPOGRAPHIC_H_
#include "Globals.h"
#include "Util.h"
#include <iostream>
#include <string>
#include <sstream>
#include <iomanip>
/**
* Stores a topographic position
*/
struct CoordTopographic
{
public:
/**
* Default constructor
*/
CoordTopographic()
: azimuth(0.0), elevation(0.0), range(0.0), range_rate(0.0)
: azimuth(0.0),
elevation(0.0),
range(0.0),
range_rate(0.0)
{
}
CoordTopographic(double az, double el, double rnge, double rnge_rate)
: azimuth(az), elevation(el), range(rnge), range_rate(rnge_rate)
/**
* Constructor
* @param[in] arg_azimuth azimuth in radians
* @param[in] arg_elevation elevation in radians
* @param[in] arg_range range in kilometers
* @param[in] arg_range_rate range rate in kilometers per second
*/
CoordTopographic(
double arg_azimuth,
double arg_elevation,
double arg_range,
double arg_range_rate)
: azimuth(arg_azimuth),
elevation(arg_elevation),
range(arg_range),
range_rate(arg_range_rate)
{
}
CoordTopographic(const CoordTopographic& t)
/**
* Copy constructor
* @param[in] topo object to copy from
*/
CoordTopographic(const CoordTopographic& topo)
{
azimuth = t.azimuth;
elevation = t.elevation;
range = t.range;
range_rate = t.range_rate;
azimuth = topo.azimuth;
elevation = topo.elevation;
range = topo.range;
range_rate = topo.range_rate;
}
/**
* Destructor
*/
virtual ~CoordTopographic()
{
};
}
CoordTopographic& operator=(const CoordTopographic& t)
/**
* Assignment operator
* @param[in] topo object to copy from
*/
CoordTopographic& operator=(const CoordTopographic& topo)
{
if (this != &t) {
azimuth = t.azimuth;
elevation = t.elevation;
range = t.range;
range_rate = t.range_rate;
if (this != &topo)
{
azimuth = topo.azimuth;
elevation = topo.elevation;
range = topo.range;
range_rate = topo.range_rate;
}
return *this;
}
bool operator==(const CoordTopographic& t) const
/**
* Equality operator
* @param[in] topo value to check
* @returns whether the object is equal
*/
bool operator==(const CoordTopographic& topo) const
{
return IsEqual(t);
return IsEqual(topo);
}
bool operator !=(const CoordTopographic& t) const
/**
* Inequality operator
* @param[in] topo the object to compare with
* @returns whether the object is not equal
*/
bool operator !=(const CoordTopographic& topo) const
{
return !IsEqual(t);
return !IsEqual(topo);
}
/**
* Dump this object to a string
* @returns string
*/
std::string ToString() const
{
std::stringstream ss;
ss << std::right << std::fixed << std::setprecision(2);
ss << "Az: " << std::setw(7) << Util::RadiansToDegrees(azimuth);
ss << ", El: " << std::setw(7) << Util::RadiansToDegrees(elevation);
ss << ", Rng: " << std::setw(9) << range;
ss << ", Rng Rt: " << std::setw(6) << range_rate;
ss << std::right << std::fixed << std::setprecision(3);
ss << "Az: " << std::setw(8) << Util::RadiansToDegrees(azimuth);
ss << ", El: " << std::setw(8) << Util::RadiansToDegrees(elevation);
ss << ", Rng: " << std::setw(10) << range;
ss << ", Rng Rt: " << std::setw(7) << range_rate;
return ss.str();
}
/*
* radians
*/
/** azimuth in radians */
double azimuth;
/*
* radians
*/
/** elevations in radians */
double elevation;
/*
* kilometers
*/
/** range in kilometers */
double range;
/*
* kilometers / second
*/
/** range rate in kilometers per second */
double range_rate;
protected:
bool IsEqual(const CoordTopographic& t) const
private:
bool IsEqual(const CoordTopographic& topo) const
{
if (azimuth == t.azimuth && elevation == t.elevation &&
range == t.range && range_rate == t.range_rate)
{
return true;
}
else
{
return false;
}
bool equal = false;
if (azimuth == topo.azimuth &&
elevation == topo.elevation &&
range == topo.range &&
range_rate == topo.range_rate)
{
equal = true;
}
return equal;
}
};
@ -104,4 +144,3 @@ inline std::ostream& operator<<(std::ostream& strm, const CoordTopographic& t)
}
#endif

132
libsgp4/DateTime.cpp Normal file
View File

@ -0,0 +1,132 @@
#include "DateTime.h"
#if 0
bool jd_dmy(int JD, int c_year, int c_month, int c_day)
{
// For the Gregorian calendar:
int a = JD + 32044;
int b = (4 * a + 3) / 146097;
int c = a - (b * 146097) / 4;
// Then, for both calendars:
int d = (4 * c + 3) / 1461;
int e = c - (1461 * d) / 4;
int m = (5 * e + 2) / 153;
int day = e - (153 * m + 2) / 5 + 1;
int month = m + 3 - 12 * (m / 10);
int year = b * 100 + d - 4800 + m / 10;
if (c_year != year || c_month != month || c_day != day)
{
std::cout << year << " " << month << " " << day << std::endl;
return false;
}
else
{
return true;
}
}
int main()
{
for (int year = 1; year <= 9999; year++)
{
for (int month = 1; month <= 12; month++)
{
for (int day = 1; day <= DateTime::DaysInMonth(year, month); day++)
{
int hour = 23;
int minute = 59;
int second = 59;
int microsecond = 999999;
DateTime dt(year, month, day, hour, minute, second, microsecond);
if (dt.Year() != year ||
dt.Month() != month ||
dt.Day() != day ||
dt.Hour() != hour ||
dt.Minute() != minute ||
dt.Second() != second ||
dt.Microsecond() != microsecond)
{
std::cout << "failed" << std::endl;
std::cout << "Y " << dt.Year() << " " << year << std::endl;
std::cout << "M " << dt.Month() << " " << month << std::endl;
std::cout << "D " << dt.Day() << " " << day << std::endl;
std::cout << "H " << dt.Hour() << " " << hour << std::endl;
std::cout << "M " << dt.Minute() << " " << minute << std::endl;
std::cout << "S " << dt.Second() << " " << second << std::endl;
std::cout << "F " << dt.Microsecond() << " " << microsecond << std::endl;
return 0;
}
if (!jd_dmy(dt.Julian() + 0.5, year, month, day))
{
std::cout << "julian" << std::endl;
return 0;
}
}
}
}
for (int hour = 1; hour < 24; hour++)
{
std::cout << hour << std::endl;
for (int minute = 0; minute < 60; minute++)
{
for (int second = 0; second < 60; second++)
{
for (int microsecond = 0; microsecond < 1000000; microsecond += 10000)
{
int year = 1000;
int month = 10;
int day = 23;
DateTime dt(year, month, day, hour, minute, second, microsecond);
if (dt.Year() != year ||
dt.Month() != month ||
dt.Day() != day ||
dt.Hour() != hour ||
dt.Minute() != minute ||
dt.Second() != second ||
dt.Microsecond() != microsecond)
{
std::cout << "failed" << std::endl;
std::cout << "Y " << dt.Year() << " " << year << std::endl;
std::cout << "M " << dt.Month() << " " << month << std::endl;
std::cout << "D " << dt.Day() << " " << day << std::endl;
std::cout << "H " << dt.Hour() << " " << hour << std::endl;
std::cout << "M " << dt.Minute() << " " << minute << std::endl;
std::cout << "S " << dt.Second() << " " << second << std::endl;
std::cout << "F " << dt.Microsecond() << " " << microsecond << std::endl;
return 0;
}
}
}
}
}
jd_dmy(1721425.5, 0, 0, 0);
DateTime d1(1000, 1, 1);
DateTime d2(2000, 1, 1);
DateTime d3(4000, 1, 1);
DateTime d4(6000, 1, 1);
DateTime d5(8000, 1, 1);
std::cout << std::setprecision(20);
std::cout << d1.Julian() << std::endl;
std::cout << d2.Julian() << std::endl;
std::cout << d3.Julian() << std::endl;
std::cout << d4.Julian() << std::endl;
std::cout << d5.Julian() << std::endl;
return 0;
}
#endif

585
libsgp4/DateTime.h Normal file
View File

@ -0,0 +1,585 @@
#ifndef DATETIME_H_
#define DATETIME_H_
#include <iomanip>
#include <iostream>
#include <sstream>
#include "TimeSpan.h"
#include "Util.h"
namespace
{
static int daysInMonth[2][13] = {
// 1 2 3 4 5 6 7 8 9 10 11 12
{0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31},
{0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}
};
static int cumulDaysInMonth[2][13] = {
// 1 2 3 4 5 6 7 8 9 10 11 12
{0, 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334},
{0, 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335}
};
}
class DateTime
{
public:
DateTime()
{
Initialise(1, 1, 1, 0, 0, 0, 0);
}
DateTime(unsigned long long ticks)
: m_encoded(ticks)
{
}
DateTime(int year, double doy)
{
m_encoded = TimeSpan(AbsoluteDays(year, doy) * TicksPerDay).Ticks();
}
DateTime(int year, int month, int day)
{
Initialise(year, month, day, 0, 0, 0, 0);
}
DateTime(int year, int month, int day, int hour, int minute, int second)
{
Initialise(year, month, day, hour, minute, second, 0);
}
DateTime(int year,
int month,
int day,
int hour,
int minute,
int second,
int microsecond)
{
Initialise(year, month, day, hour, minute, second, microsecond);
}
void Initialise(int year,
int month,
int day,
int hour,
int minute,
int second,
int microsecond)
{
if (!IsValidYearMonthDay(year, month, day) ||
hour < 0 || hour > 23 ||
minute < 0 || minute > 59 ||
second < 0 || second > 59 ||
microsecond < 0 || microsecond > 999999)
{
throw 1;
}
m_encoded = TimeSpan(
AbsoluteDays(year, month, day),
hour,
minute,
second,
microsecond).Ticks();
}
static DateTime Now(bool microseconds = false)
{
DateTime dt;
struct timespec ts;
if (clock_gettime(CLOCK_REALTIME, &ts) == 0)
{
if (microseconds)
{
dt = DateTime(UnixEpoch
+ ts.tv_sec * TicksPerSecond
+ ts.tv_nsec / 1000LL * TicksPerMicrosecond);
}
else
{
dt = DateTime(UnixEpoch
+ ts.tv_sec * TicksPerSecond);
}
}
else
{
throw 1;
}
return dt;
}
static bool IsLeapYear(int year)
{
if (!IsValidYear(year))
{
throw 1;
}
return (((year % 4) == 0 && (year % 100) != 0) || (year % 400) == 0);
}
static bool IsValidYear(int year)
{
bool valid = true;
if (year < 1 || year > 9999)
{
valid = false;
}
return valid;
}
static bool IsValidYearMonth(int year, int month)
{
bool valid = true;
if (IsValidYear(year))
{
if (month < 1 || month > 12)
{
valid = false;
}
}
else
{
valid = false;
}
return valid;
}
static bool IsValidYearMonthDay(int year, int month, int day)
{
bool valid = true;
if (IsValidYearMonth(year, month))
{
if (day < 1 || day > DaysInMonth(year, month))
{
valid = false;
}
}
else
{
valid = false;
}
return valid;
}
static int DaysInMonth(int year, int month)
{
if (!IsValidYearMonth(year, month))
{
throw 1;
}
const int* daysInMonthPtr;
if (IsLeapYear(year))
{
daysInMonthPtr = daysInMonth[1];
}
else
{
daysInMonthPtr = daysInMonth[0];
}
return daysInMonthPtr[month];
}
int DayOfYear(int year, int month, int day) const
{
if (!IsValidYearMonthDay(year, month, day))
{
throw 1;
}
int daysThisYear = day;
if (IsLeapYear(year))
{
daysThisYear += cumulDaysInMonth[1][month];
}
else
{
daysThisYear += cumulDaysInMonth[0][month];
}
return daysThisYear;
}
double AbsoluteDays(int year, double doy) const
{
int previousYear = year - 1;
/*
* + days in previous years ignoring leap days
* + Julian leap days before this year
* - minus prior century years
* + plus prior years divisible by 400 days
*/
long long daysSoFar = 365 * previousYear
+ previousYear / 4
- previousYear / 100
+ previousYear / 400;
return daysSoFar + doy - 1.0;
}
int AbsoluteDays(int year, int month, int day) const
{
int previousYear = year - 1;
/*
* days this year (0 - ...)
* + days in previous years ignoring leap days
* + Julian leap days before this year
* - minus prior century years
* + plus prior years divisible by 400 days
*/
int result = DayOfYear(year, month, day) - 1
+ 365 * previousYear
+ previousYear / 4
- previousYear / 100
+ previousYear / 400;
return result;
}
TimeSpan TimeOfDay() const
{
return TimeSpan(Ticks() % TicksPerDay);
}
int DayOfWeek() const
{
/*
* The fixed day 1 (January 1, 1 Gregorian) is Monday.
* 0 Sunday
* 1 Monday
* 2 Tuesday
* 3 Wednesday
* 4 Thursday
* 5 Friday
* 6 Saturday
*/
return (((m_encoded / TicksPerDay) + 1) % 7);
}
bool Equals(const DateTime& dt) const
{
return (m_encoded == dt.m_encoded);
}
int Compare(const DateTime& dt) const
{
int ret = 0;
if (m_encoded < dt.m_encoded)
{
return -1;
}
else if (m_encoded > dt.m_encoded)
{
return 1;
}
return ret;
}
DateTime AddYears(const int years) const
{
return AddMonths(years * 12);
}
DateTime AddMonths(const int months) const
{
int year;
int month;
int day;
FromTicks(year, month, day);
month += months % 12;
year += months / 12;
if (month < 1)
{
month += 12;
--year;
}
else if (month > 12)
{
month -= 12;
++year;
}
int maxday = DaysInMonth(year, month);
day = std::min(day, maxday);
return DateTime(year, month, day).Add(TimeOfDay());
}
DateTime Add(const TimeSpan& t) const
{
return AddTicks(t.Ticks());
}
DateTime AddDays(const double days) const
{
return AddMicroseconds(days * 86400000000.0);
}
DateTime AddHours(const double hours) const
{
return AddMicroseconds(hours * 3600000000.0);
}
DateTime AddMinutes(const double minutes) const
{
return AddMicroseconds(minutes * 60000000.0);
}
DateTime AddSeconds(const double seconds) const
{
return AddMicroseconds(seconds * 1000000.0);
}
DateTime AddMicroseconds(const double microseconds) const
{
long long ticks = microseconds * TicksPerMicrosecond;
return AddTicks(ticks);
}
DateTime AddTicks(long long ticks) const
{
return DateTime(m_encoded + ticks);
}
long long Ticks() const
{
return m_encoded;
}
void FromTicks(int& year, int& month, int& day) const
{
int totalDays = m_encoded / TicksPerDay;
/*
* number of 400 year cycles
*/
int num400 = totalDays / 146097;
totalDays -= num400 * 146097;
/*
* number of 100 year cycles
*/
int num100 = totalDays / 36524;
if (num100 == 4)
{
/*
* last day of the last leap century
*/
num100 = 3;
}
totalDays -= num100 * 36524;
/*
* number of 4 year cycles
*/
int num4 = totalDays / 1461;
totalDays -= num4 * 1461;
/*
* number of years
*/
int num1 = totalDays / 365;
if (num1 == 4)
{
/*
* last day of the last leap olympiad
*/
num1 = 3;
}
totalDays -= num1 * 365;
/*
* find year
*/
year = (num400 * 400) + (num100 * 100) + (num4 * 4) + num1 + 1;
/*
* convert day of year to month/day
*/
const int* daysInMonthPtr;
if (IsLeapYear(year))
{
daysInMonthPtr = daysInMonth[1];
}
else
{
daysInMonthPtr = daysInMonth[0];
}
month = 1;
while (totalDays >= daysInMonthPtr[month] && month <= 12)
{
totalDays -= daysInMonthPtr[month++];
}
day = totalDays + 1;
}
int Year() const
{
int year;
int month;
int day;
FromTicks(year, month, day);
return year;
}
int Month() const
{
int year;
int month;
int day;
FromTicks(year, month, day);
return month;
}
int Day() const
{
int year;
int month;
int day;
FromTicks(year, month, day);
return day;
}
int Hour() const
{
return m_encoded % TicksPerDay / TicksPerHour;
}
int Minute() const
{
return m_encoded % TicksPerHour / TicksPerMinute;
}
int Second() const
{
return m_encoded % TicksPerMinute / TicksPerSecond;
}
int Microsecond() const
{
return m_encoded % TicksPerSecond / TicksPerMicrosecond;
}
double Julian() const
{
TimeSpan ts = TimeSpan(Ticks());
return ts.TotalDays() + 1721425.5;
}
double ToGreenwichSiderealTime() const
{
// t = Julian centuries from 2000 Jan. 1 12h UT1
const double t = (Julian() - 2451545.0) / 36525.0;
// Rotation angle in arcseconds
double theta = 67310.54841
+ (876600.0 * 3600.0 + 8640184.812866) * t
+ 0.093104 * t * t
- 0.0000062 * t * t * t;
// 360.0 / 86400.0 = 1.0 / 240.0
return Util::WrapTwoPI(Util::DegreesToRadians(theta / 240.0));
}
double ToLocalMeanSiderealTime(const double& lon) const
{
return Util::WrapTwoPI(ToGreenwichSiderealTime() + lon);
}
std::string ToString() const
{
std::stringstream ss;
int year;
int month;
int day;
FromTicks(year, month, day);
ss << std::right << std::setfill('0');
ss << std::setw(4) << year << "-";
ss << std::setw(2) << month << "-";
ss << std::setw(2) << day << " ";
ss << std::setw(2) << Hour() << ":";
ss << std::setw(2) << Minute() << ":";
ss << std::setw(2) << Second() << ".";
ss << std::setw(6) << Microsecond() << " UTC";
return ss.str();
}
private:
unsigned long long m_encoded;
};
inline std::ostream& operator<<(std::ostream& strm, const DateTime& dt)
{
return strm << dt.ToString();
}
inline DateTime operator+(const DateTime& dt, TimeSpan ts)
{
long long int res = dt.Ticks() + ts.Ticks();
if (res < 0 || res > MaxValueTicks)
{
throw 1;
}
return DateTime(res);
}
inline DateTime operator-(const DateTime& dt, const TimeSpan& ts)
{
long long int res = dt.Ticks() - ts.Ticks();
if (res < 0 || res > MaxValueTicks)
{
throw 1;
}
return DateTime(res);
}
inline TimeSpan operator-(const DateTime& dt1, const DateTime& dt2)
{
return TimeSpan(dt1.Ticks() - dt2.Ticks());
}
inline bool operator==(const DateTime& dt1, const DateTime& dt2)
{
return dt1.Equals(dt2);
}
inline bool operator>(const DateTime& dt1, const DateTime& dt2)
{
return (dt1.Compare(dt2) > 0);
}
inline bool operator>=(const DateTime& dt1, const DateTime& dt2)
{
return (dt1.Compare(dt2) >= 0);
}
inline bool operator!=(const DateTime& dt1, const DateTime& dt2)
{
return !dt1.Equals(dt2);
}
inline bool operator<(const DateTime& dt1, const DateTime& dt2)
{
return (dt1.Compare(dt2) < 0);
}
inline bool operator<=(const DateTime& dt1, const DateTime& dt2)
{
return (dt1.Compare(dt2) <= 0);
}
#endif

31
libsgp4/DecayedException.h
View File

@ -2,44 +2,65 @@
#define DECAYEDEXCEPTION_H_
#include <exception>
#include <iostream>
#include "Julian.h"
#include "DateTime.h"
#include "Vector.h"
class DecayedException : public std::exception
{
public:
DecayedException(const Julian& dt, const Vector& pos, const Vector& vel)
/**
* Constructor
* @param[in] dt time of the event
* @param[in] pos position of the satellite at dt
* @param[in] vel velocity of the satellite at dt
*/
DecayedException(const DateTime& dt, const Vector& pos, const Vector& vel)
: _dt(dt), _pos(pos), _vel(vel)
{
}
/**
* Destructor
*/
virtual ~DecayedException(void) throw ()
{
}
/**
* @returns the error string
*/
virtual const char* what() const throw ()
{
return "Error: Satellite decayed";
}
Julian GetDate() const
/**
* @returns the date
*/
DateTime GetDateTime() const
{
return _dt;
}
/**
* @returns the position
*/
Vector GetPosition() const
{
return _pos;
}
/**
* @returns the velocity
*/
Vector GetVelocity() const
{
return _vel;
}
private:
Julian _dt;
DateTime _dt;
Vector _pos;
Vector _vel;
};

53
libsgp4/Eci.cpp
View File

@ -1,28 +1,34 @@
#include "Eci.h"
#include "Globals.h"
#include "Util.h"
void Eci::ToEci(const Julian& date, const CoordGeodetic &g)
/**
* Converts a DateTime and Geodetic position to Eci coordinates
* @param[in] dt the date
* @param[in] geo the geodetic position
*/
void Eci::ToEci(const DateTime& dt, const CoordGeodetic &geo)
{
/*
* set date
*/
date_ = date;
m_dt = dt;
static const double mfactor = kTWOPI * (kOMEGA_E / kSECONDS_PER_DAY);
/*
* Calculate Local Mean Sidereal Time for observers longitude
*/
const double theta = date_.ToLocalMeanSiderealTime(g.longitude);
const double theta = m_dt.ToLocalMeanSiderealTime(geo.longitude);
/*
* take into account earth flattening
*/
const double c = 1.0
/ sqrt(1.0 + kF * (kF - 2.0) * pow(sin(g.latitude), 2.0));
/ sqrt(1.0 + kF * (kF - 2.0) * pow(sin(geo.latitude), 2.0));
const double s = pow(1.0 - kF, 2.0) * c;
const double achcp = (kXKMPER * c + g.altitude) * cos(g.latitude);
const double achcp = (kXKMPER * c + geo.altitude) * cos(geo.latitude);
/*
* X position in km
@ -30,10 +36,10 @@ void Eci::ToEci(const Julian& date, const CoordGeodetic &g)
* Z position in km
* W magnitude in km
*/
position_.x = achcp * cos(theta);
position_.y = achcp * sin(theta);
position_.z = (kXKMPER * s + g.altitude) * sin(g.latitude);
position_.w = position_.GetMagnitude();
m_position.x = achcp * cos(theta);
m_position.y = achcp * sin(theta);
m_position.z = (kXKMPER * s + geo.altitude) * sin(geo.latitude);
m_position.w = m_position.GetMagnitude();
/*
* X velocity in km/s
@ -41,27 +47,28 @@ void Eci::ToEci(const Julian& date, const CoordGeodetic &g)
* Z velocity in km/s
* W magnitude in km/s
*/
velocity_.x = -mfactor * position_.y;
velocity_.y = mfactor * position_.x;
velocity_.z = 0.0;
velocity_.w = velocity_.GetMagnitude();
m_velocity.x = -mfactor * m_position.y;
m_velocity.y = mfactor * m_position.x;
m_velocity.z = 0.0;
m_velocity.w = m_velocity.GetMagnitude();
}
/**
* @returns the position in geodetic form
*/
CoordGeodetic Eci::ToGeodetic() const
{
const double theta = Util::AcTan(position_.y, position_.x);
const double theta = Util::AcTan(m_position.y, m_position.x);
// 0 >= lon < 360
// const double lon = Fmod2p(theta - date_.ToGreenwichSiderealTime());
// 180 >= lon < 180
const double lon = Util::WrapNegPosPI(theta - date_.ToGreenwichSiderealTime());
const double lon = Util::WrapNegPosPI(theta
- m_dt.ToGreenwichSiderealTime());
const double r = sqrt((position_.x * position_.x)
+ (position_.y * position_.y));
const double r = sqrt((m_position.x * m_position.x)
+ (m_position.y * m_position.y));
static const double e2 = kF * (2.0 - kF);
double lat = Util::AcTan(position_.z, r);
double lat = Util::AcTan(m_position.z, r);
double phi = 0.0;
double c = 0.0;
int cnt = 0;
@ -71,12 +78,12 @@ CoordGeodetic Eci::ToGeodetic() const
phi = lat;
const double sinphi = sin(phi);
c = 1.0 / sqrt(1.0 - e2 * sinphi * sinphi);
lat = Util::AcTan(position_.z + kXKMPER * c * e2 * sinphi, r);
lat = Util::AcTan(m_position.z + kXKMPER * c * e2 * sinphi, r);
cnt++;
}
while (fabs(lat - phi) >= 1e-10 && cnt < 10);
const double alt = r / cos(lat) - kXKMPER * c;
return CoordGeodetic(lat, lon, alt, true);
return CoordGeodetic(lat, lon, alt);
}

116
libsgp4/Eci.h
View File

@ -3,67 +3,133 @@
#include "CoordGeodetic.h"
#include "Vector.h"
#include "Julian.h"
#include "Globals.h"
#include "DateTime.h"
/**
* A class store to store an Earth-centered inertial position.
* This is only valid for the date specified.
*/
class Eci
{
public:
/*
* in degrees
/**
* @param[in] dt the date to be used for this position
* @param[in] latitude the latitude in degrees
* @param[in] longitude the longitude in degrees
* @param[in] altitude the altitude in kilometers
*/
Eci(const Julian& date, double latitude, double longitude, double altitude)
Eci(const DateTime& dt,
const double latitude,
const double longitude,
const double altitude)
{
ToEci(date, CoordGeodetic(latitude, longitude, altitude));
ToEci(dt, CoordGeodetic(latitude, longitude, altitude));
}
Eci(const Julian& date, const CoordGeodetic& g)
/**
* @param[in] dt the date to be used for this position
* @param[in] geo the position
*/
Eci(const DateTime& dt, const CoordGeodetic& geo)
{
ToEci(date, g);
ToEci(dt, geo);
}
Eci(const Julian &date, const Vector &position)
: date_(date), position_(position)
/**
* @param[in] dt the date to be used for this position
* @param[in] position
*/
Eci(const DateTime &dt, const Vector &position)
: m_dt(dt),
m_position(position)
{
}
Eci(const Julian &date, const Vector &position, const Vector &velocity)
: date_(date), position_(position), velocity_(velocity)
/**
* @param[in] dt the date to be used for this position
* @param[in] position the position
* @param[in] velocity the velocity
*/
Eci(const DateTime &dt, const Vector &position, const Vector &velocity)
: m_dt(dt),
m_position(position),
m_velocity(velocity)
{
}
/**
* Destructor
*/
virtual ~Eci()
{
}
/**
* Equality operator
* @param dt the date to compare
* @returns true if the object matches
*/
bool operator==(const DateTime& dt) const
{
return m_dt == dt;
}
/**
* Inequality operator
* @param dt the date to compare
* @returns true if the object doesn't match
*/
bool operator!=(const DateTime& dt) const
{
return m_dt != dt;
}
/**
* Update this object with a new date and geodetic position
* @param dt new date
* @param geo new geodetic position
*/
void Update(const DateTime& dt, const CoordGeodetic& geo)
{
ToEci(dt, geo);
}
/**
* @returns the position
*/
Vector GetPosition() const
{
return position_;
return m_position;
}
/**
* @returns the velocity
*/
Vector GetVelocity() const
{
return velocity_;
return m_velocity;
}
Julian GetDate() const
/**
* @returns the date
*/
DateTime GetDateTime() const
{
return date_;
return m_dt;
}
/**
* @returns the position in geodetic form
*/
CoordGeodetic ToGeodetic() const;
protected:
void ToEci(const Julian& date, double latitude, double longitude,
double altitude);
void ToEci(const Julian& date, const CoordGeodetic& g);
private:
Julian date_;
Vector position_;
Vector velocity_;
void ToEci(const DateTime& dt, const CoordGeodetic& geo);
DateTime m_dt;
Vector m_position;
Vector m_velocity;
};
#endif

308
libsgp4/Julian.cpp
View File

@ -1,308 +0,0 @@
#include "Globals.h"
#include "Julian.h"
#include "Util.h"
#include <cmath>
#include <ctime>
#include <cassert>
#ifdef WIN32
#include <windows.h>
#else
#include <sys/time.h>
#endif
Julian::Julian()
{
#ifdef WIN32
SYSTEMTIME st;
GetSystemTime(&st);
Initialize(st.wYear,
st.wMonth,
st.wDay,
st.wHour,
st.wMinute,
(double) st.wSecond + (double) st.wMilliseconds / 1000.0);
#else
struct timeval tv;
gettimeofday(&tv, NULL);
struct tm gmt;
gmtime_r(&tv.tv_sec, &gmt);
Initialize(gmt.tm_year + 1900,
gmt.tm_mon + 1,
gmt.tm_mday,
gmt.tm_hour,
gmt.tm_min,
(double) gmt.tm_sec + (double) tv.tv_usec / 1000000.0);
#endif
}
/*
* create julian date given time_t value
*/
Julian::Julian(const time_t t)
{
struct tm ptm;
#if WIN32
if (gmtime_s(&ptm, &t))
{
assert(1);
}
#else
if (gmtime_r(&t, &ptm) == NULL)
{
assert(1);
}
#endif
int year = ptm.tm_year + 1900;
double day = ptm.tm_yday + 1 +
(ptm.tm_hour +
((ptm.tm_min +
(ptm.tm_sec / 60.0)) / 60.0)) / 24.0;
Initialize(year, day);
}
/*
* comparison
*/
bool Julian::operator==(const Julian &date) const
{
return date_ == date.date_ ? true : false;
}
bool Julian::operator!=(const Julian &date) const
{
return !(*this == date);
}
bool Julian::operator>(const Julian &date) const
{
return date_ > date.date_ ? true : false;
}
bool Julian::operator<(const Julian &date) const {
return date_ < date.date_ ? true : false;
}
bool Julian::operator>=(const Julian &date) const
{
return date_ >= date.date_ ? true : false;
}
bool Julian::operator<=(const Julian &date) const
{
return date_ <= date.date_ ? true : false;
}
/*
* assignment
*/
Julian& Julian::operator=(const Julian& b)
{
if (this != &b) {
date_ = b.date_;
}
return (*this);
}
Julian& Julian::operator=(const double b)
{
date_ = b;
return (*this);
}
/*
* arithmetic
*/
Julian Julian::operator +(const Timespan& b) const
{
return Julian(*this) += b;
}
Julian Julian::operator-(const Timespan& b) const
{
return Julian(*this) -= b;
}
Timespan Julian::operator-(const Julian& b) const
{
return Timespan(date_ - b.date_);
}
/*
* compound assignment
*/
Julian & Julian::operator +=(const Timespan& b)
{
date_ += b;
return (*this);
}
Julian & Julian::operator -=(const Timespan& b)
{
date_ -= b;
return (*this);
}
/*
* create julian date from year and day of year
*/
void Julian::Initialize(int year, double day)
{
year--;
int A = (year / 100);
int B = 2 - A + (A / 4);
double new_years = static_cast<int> (365.25 * year) +
static_cast<int> (30.6001 * 14) +
1720994.5 + B;
date_ = new_years + day;
}
/*
* create julian date from individual components
* year: 2004
* mon: 1-12
* day: 1-31
* hour: 0-23
* min: 0-59
* sec: 0-59.99
*/
void Julian::Initialize(int year, int mon, int day,
int hour, int min, double sec)
{
// Calculate N, the day of the year (1..366)
int N;
int F1 = (int) ((275.0 * mon) / 9.0);
int F2 = (int) ((mon + 9.0) / 12.0);
if (IsLeapYear(year))
{
// Leap year
N = F1 - F2 + day - 30;
}
else
{
// Common year
N = F1 - (2 * F2) + day - 30;
}
double dblDay = N + (hour + (min + (sec / 60.0)) / 60.0) / 24.0;
Initialize(year, dblDay);
}
/*
* converts time to time_t
* note: resolution to seconds only
*/
time_t Julian::ToTime() const
{
return static_cast<time_t> ((date_ - 2440587.5) * 86400.0);
}
/*
* Greenwich Mean Sidereal Time
*/
double Julian::ToGreenwichSiderealTime() const
{
#if 0
const double UT = fmod(jul + 0.5, 1.0);
const double TU = (jul - 2451545.0 - UT) / 36525.0;
double GMST = 24110.54841 + TU *
(8640184.812866 + TU * (0.093104 - TU * 6.2e-06));
GMST = fmod(GMST + SEC_PER_DAY * OMEGA_E * UT, SEC_PER_DAY);
if (GMST < 0.0)
{
GMST += SEC_PER_DAY; // "wrap" negative modulo value
}
return (TWOPI * (GMST / SEC_PER_DAY));
#endif
// tut1 = Julian centuries from 2000 Jan. 1 12h UT1
// (since J2000 which is 2451545.0)
// a Julian century is 36525 days
const double tut1 = (date_ - 2451545.0) / 36525.0;
const double tut2 = tut1 * tut1;
const double tut3 = tut2 * tut1;
// Rotation angle in arcseconds
double theta = 67310.54841 + (876600.0 * 3600.0 + 8640184.812866) * tut1
+ 0.093104 * tut2 - 0.0000062 * tut3;
// 360.0 / 86400.0 = 1.0 / 240.0
theta = Util::WrapTwoPI(Util::DegreesToRadians(theta / 240.0));
return theta;
#if 0
static const double C1 = 1.72027916940703639e-2;
static const double THGR70 = 1.7321343856509374;
static const double FK5R = 5.07551419432269442e-15;
/*
* get integer number of days from 0 jan 1970
*/
const double ts70 = date_ - 2433281.5 - 7305.0;
const double ds70 = floor(ts70 + 1.0e-8);
const double tfrac = ts70 - ds70;
/*
* find greenwich location at epoch
*/
const double c1p2p = C1 + kTWOPI;
double gsto = Util::WrapTwoPI(THGR70 + C1 * ds70 + c1p2p * tfrac
+ ts70 * ts70 * FK5R);
return gsto;
#endif
}
/*
* Local Mean Sideral Time
*/
double Julian::ToLocalMeanSiderealTime(const double& lon) const
{
return fmod(ToGreenwichSiderealTime() + lon, kTWOPI);
}
void Julian::ToGregorian(struct DateTimeComponents* datetime) const
{
double jdAdj = GetDate() + 0.5;
int Z = (int) jdAdj;
double F = jdAdj - Z;
int A = 0;
if (Z < 2299161)
{
A = static_cast<int> (Z);
}
else
{
int a = static_cast<int> ((Z - 1867216.25) / 36524.25);
A = static_cast<int> (Z + 1 + a - static_cast<int> (a / 4));
}
int B = A + 1524;
int C = static_cast<int> ((B - 122.1) / 365.25);
int D = static_cast<int> (365.25 * C);
int E = static_cast<int> ((B - D) / 30.6001);
datetime->hours = static_cast<int> (F * 24.0);
F -= datetime->hours / 24.0;
datetime->minutes = static_cast<int> (F * 1440.0);
F -= datetime->minutes / 1440.0;
datetime->seconds = F * 86400.0;
datetime->days = B - D - static_cast<int> (30.6001 * E);
datetime->months = E < 14 ? E - 1 : E - 13;
datetime->years = datetime->months > 2 ? C - 4716 : C - 4715;
}

168
libsgp4/Julian.h
View File

@ -1,168 +0,0 @@
#ifndef JULIAN_H_
#define JULIAN_H_
#include "Globals.h"
#include "Timespan.h"
#include <ctime>
#include <iostream>
#include <string>
#include <sstream>
#include <iomanip>
class Julian
{
public:
Julian();
Julian(const Julian& jul)
{
date_ = jul.date_;
}
Julian(const time_t t);
/*
* create julian date from year and day of year
*/
Julian(int year, double day)
{
Initialize(year, day);
}
/*
* create julian date from individual components
* year: 2004
* mon: 1-12
* day: 1-31
* hour: 0-23
* min: 0-59
* sec: 0-59.99
*/
Julian(int year, int mon, int day, int hour, int min, double sec)
{
Initialize(year, mon, day, hour, min, sec);
}
~Julian()
{
}
// comparison operators
bool operator==(const Julian &date) const;
bool operator!=(const Julian &date) const;
bool operator>(const Julian &date) const;
bool operator<(const Julian &date) const;
bool operator>=(const Julian &date) const;
bool operator<=(const Julian &date) const;
// assignment
Julian& operator =(const Julian& b);
Julian& operator =(const double b);
// arithmetic
Julian operator +(const Timespan& b) const;
Julian operator -(const Timespan& b) const;
Timespan operator -(const Julian& b) const;
// compound assignment
Julian & operator +=(const Timespan& b);
Julian & operator -=(const Timespan& b);
std::string ToString() const
{
std::stringstream ss;
struct DateTimeComponents dt;
ToGregorian(&dt);
ss << std::right << std::fixed;
ss << std::setprecision(6) << std::setfill('0');
ss << std::setw(4) << dt.years << "-";
ss << std::setw(2) << dt.months << "-";
ss << std::setw(2) << dt.days << " ";
ss << std::setw(2) << dt.hours << ":";
ss << std::setw(2) << dt.minutes << ":";
ss << std::setw(9) << dt.seconds << " UTC";
return ss.str();
}
struct DateTimeComponents {
DateTimeComponents()
: years(0), months(0), days(0),
hours(0), minutes(0), seconds(0.0) {
}
int years;
int months;
int days;
int hours;
int minutes;
double seconds;
};
void ToGregorian(struct DateTimeComponents* datetime) const;
time_t ToTime() const;
double ToGreenwichSiderealTime() const;
double ToLocalMeanSiderealTime(const double& lon) const;
double FromJan1_00h_1900() const
{
return date_ - kEPOCH_JAN1_00H_1900;
}
double FromJan1_12h_1900() const
{
return date_ - kEPOCH_JAN1_12H_1900;
}
double FromJan1_12h_2000() const
{
return date_ - kEPOCH_JAN1_12H_2000;
}
double GetDate() const
{
return date_;
}
void AddDay(double day)
{
date_ += day;
}
void AddHour(double hr)
{
date_ += (hr / kHOURS_PER_DAY);
}
void AddMin(double min)
{
date_ += (min / kMINUTES_PER_DAY);
}
void AddSec(double sec)
{
date_ += (sec / kSECONDS_PER_DAY);
}
static bool IsLeapYear(int y)
{
return (y % 4 == 0 && y % 100 != 0) || (y % 400 == 0);
}
protected:
void Initialize(int year, double day);
void Initialize(int year, int mon, int day, int hour, int min, double sec);
/*
* the stored julian date
*/
double date_;
};
inline std::ostream& operator<<(std::ostream& strm, const Julian& j)
{
return strm << j.ToString();
}
#endif

12
libsgp4/Makefile.am
View File

@ -2,24 +2,22 @@ lib_LIBRARIES = libsgp4.a
libsgp4_a_SOURCES = Tle.cpp \
OrbitalElements.cpp \
Observer.cpp \
Timespan.cpp \
Util.cpp \
Eci.cpp \
Julian.cpp \
SGP4.cpp \
SolarPosition.cpp
TimeSpan.cpp \
DateTime.cpp
include_HEADERS = CoordGeodetic.h \
Eci.h \
Julian.h \
OrbitalElements.h \
SatelliteException.h \
Timespan.h \
TleException.h \
Vector.h \
CoordTopographic.h \
Globals.h \
Observer.h \
SGP4.h \
SolarPosition.h \
Tle.h \
Util.h
Util.h \
TimeSpan.h \
DateTime.h

14
libsgp4/Observer.cpp
View File

@ -1,5 +1,7 @@
#include "Observer.h"
#include "CoordTopographic.h"
/*
* calculate lookangle between the observer and the passed in Eci object
*/
@ -9,23 +11,23 @@ CoordTopographic Observer::GetLookAngle(const Eci &eci)
* update the observers Eci to match the time of the Eci passed in
* if necessary
*/
UpdateObserversEci(eci.GetDate());
Update(eci.GetDateTime());
/*
* calculate differences
*/
Vector range_rate = eci.GetVelocity().Subtract(observers_eci_.GetVelocity());
Vector range = eci.GetPosition().Subtract(observers_eci_.GetPosition());
Vector range_rate = eci.GetVelocity() - m_eci.GetVelocity();
Vector range = eci.GetPosition() - m_eci.GetPosition();
range.w = range.GetMagnitude();
/*
* Calculate Local Mean Sidereal Time for observers longitude
*/
double theta = eci.GetDate().ToLocalMeanSiderealTime(geo_.longitude);
double theta = eci.GetDateTime().ToLocalMeanSiderealTime(m_geo.longitude);
double sin_lat = sin(geo_.latitude);
double cos_lat = cos(geo_.latitude);
double sin_lat = sin(m_geo.latitude);
double cos_lat = cos(m_geo.latitude);
double sin_theta = sin(theta);
double cos_theta = cos(theta);

81
libsgp4/Observer.h
View File

@ -2,66 +2,87 @@
#define OBSERVER_H_
#include "CoordGeodetic.h"
#include "CoordTopographic.h"
#include "Julian.h"
#include "Eci.h"
class DateTime;
class CoordTopographic;
class Observer
{
public:
Observer(double latitude, double longitude, double altitude)
: geo_(latitude, longitude, altitude),
observers_eci_(Julian(), geo_)
/**
* Constructor
* @param[in] latitude observers latitude in degrees
* @param[in] longitude observers longitude in degrees
* @param[in] altitude observers altitude in kilometers
*/
Observer(const double latitude,
const double longitude,
const double altitude)
: m_geo(latitude, longitude, altitude),
m_eci(DateTime(), m_geo)
{
}
Observer(const CoordGeodetic &g)
: geo_(g), observers_eci_(Julian(), geo_)
/**
* Constructor
* @param[in] geo the observers position
*/
Observer(const CoordGeodetic &geo)
: m_geo(geo),
m_eci(DateTime(), geo)
{
}
/**
* Destructor
*/
virtual ~Observer()
{
}
void SetLocation(const CoordGeodetic& g)
/**
* Set the observers location
* @param[in] geo the observers position
*/
void SetLocation(const CoordGeodetic& geo)
{
geo_ = g;
observers_eci_ = Eci(observers_eci_.GetDate(), geo_);
m_geo = geo;
m_eci.Update(m_eci.GetDateTime(), m_geo);
}
/**
* Get the observers location
* @returns the observers position
*/
CoordGeodetic GetLocation() const
{
return geo_;
}
Eci GetEciPosition(const Julian &date) const
{
return Eci(date, geo_);
return m_geo;
}
/**
* Get the look angle for the observers position to the object
* @param[in] eci the object to find the look angle to
* @returns the lookup angle
*/
CoordTopographic GetLookAngle(const Eci &eci);
private:
void UpdateObserversEci(const Julian &date)
/**
* @param[in] dt the date to update the observers position for
*/
void Update(const DateTime &dt)
{
/*
* if date has changed, update for new date
*/
if (observers_eci_.GetDate() != date)
if (m_eci != dt)
{
observers_eci_ = Eci(date, geo_);
m_eci.Update(dt, m_geo);
}
}
/*
* the observers position
*/
CoordGeodetic geo_;
/*
* the observers eci for a particular time
*/
Eci observers_eci_;
/** the observers position */
CoordGeodetic m_geo;
/** the observers Eci for a particular time */
Eci m_eci;
};
#endif

2
libsgp4/OrbitalElements.cpp
View File

@ -1,5 +1,7 @@
#include "OrbitalElements.h"
#include "Tle.h"
OrbitalElements::OrbitalElements(const Tle& tle)
{
/*

10
libsgp4/OrbitalElements.h
View File

@ -1,7 +1,10 @@
#ifndef ORBITALELEMENTS_H_
#define ORBITALELEMENTS_H_
#include "Tle.h"
#include "Util.h"
#include "DateTime.h"
class Tle;
class OrbitalElements
{
@ -103,7 +106,7 @@ public:
/*
* EPOCH
*/
Julian Epoch() const
DateTime Epoch() const
{
return epoch_;
}
@ -121,8 +124,7 @@ private:
double recovered_mean_motion_;
double perigee_;
double period_;
Julian epoch_;
DateTime epoch_;
};
#endif

766
libsgp4/SGP4.cpp
View File

File diff suppressed because it is too large Load Diff

113
libsgp4/SGP4.h
View File

@ -22,17 +22,10 @@ public:
void SetTle(const Tle& tle);
Eci FindPosition(double tsince) const;
Eci FindPosition(const Julian& date) const;
Eci FindPosition(const DateTime& date) const;
struct CommonConstants
{
CommonConstants()
: cosio(0.0), sinio(0.0), eta(0.0), t2cof(0.0), a3ovk2(0.0),
x1mth2(0.0), x3thm1(0.0), x7thm1(0.0), aycof(0.0), xlcof(0.0),
xnodcf(0.0), c1(0.0), c4(0.0), omgdot(0.0), xnodot(0.0), xmdot(0.0)
{
}
double cosio;
double sinio;
double eta;
@ -48,17 +41,11 @@ public:
double c4;
double omgdot; // secular rate of omega (radians/sec)
double xnodot; // secular rate of xnode (radians/sec)
double xmdot; // secular rate of xmo (radians/sec)
double xmdot; // secular rate of xmo (radians/sec)
};
struct NearSpaceConstants
{
NearSpaceConstants()
: c5(0.0), omgcof(0.0), xmcof(0.0), delmo(0.0), sinmo(0.0), d2(0.0),
d3(0.0), d4(0.0), t3cof(0.0), t4cof(0.0), t5cof(0.0)
{
}
double c5;
double omgcof;
double xmcof;
@ -74,18 +61,6 @@ public:
struct DeepSpaceConstants
{
DeepSpaceConstants()
: gsto(0.0), zmol(0.0), zmos(0.0), resonance_flag(false),
synchronous_flag(false), sse(0.0), ssi(0.0), ssl(0.0), ssg(0.0),
ssh(0.0), se2(0.0), si2(0.0), sl2(0.0), sgh2(0.0), sh2(0.0), se3(0.0),
si3(0.0), sl3(0.0), sgh3(0.0), sh3(0.0), sl4(0.0), sgh4(0.0), ee2(0.0),
e3(0.0), xi2(0.0), xi3(0.0), xl2(0.0), xl3(0.0), xl4(0.0), xgh2(0.0),
xgh3(0.0), xgh4(0.0), xh2(0.0), xh3(0.0), d2201(0.0), d2211(0.0),
d3210(0.0), d3222(0.0), d4410(0.0), d4422(0.0), d5220(0.0), d5232(0.0),
d5421(0.0), d5433(0.0), del1(0.0), del2(0.0), del3(0.0)
{
}
double gsto;
double zmol;
double zmos;
@ -156,10 +131,6 @@ public:
struct IntegratorValues
{
IntegratorValues() : xndot(0.0), xnddt(0.0), xldot(0.0)
{
}
double xndot;
double xnddt;
double xldot;
@ -167,10 +138,6 @@ public:
struct IntegratorConstants
{
IntegratorConstants() : xfact(0.0), xlamo(0.0)
{
}
/*
* integrator constants
*/
@ -185,10 +152,6 @@ public:
struct IntegratorParams
{
IntegratorParams() : xli(0.0), xni(0.0), atime(0.0)
{
}
/*
* integrator values
*/
@ -203,25 +166,59 @@ public:
private:
void Initialise();
void DeepSpaceInitialise(const double& eosq, const double& sinio, const double& cosio, const double& betao,
const double& theta2, const double& betao2,
const double& xmdot, const double& omgdot, const double& xnodot);
void DeepSpaceCalculateLunarSolarTerms(const double t, double* pe, double* pinc,
double* pl, double* pgh, double* ph) const;
void DeepSpacePeriodics(const double& t, double* em, double* xinc,
double* omgasm, double* xnodes, double* xll) const;
void DeepSpaceSecular(const double& t, double* xll, double* omgasm,
double* xnodes, double* em, double* xinc, double* xn) const;
Eci FindPosition(const Julian& dt, double tsince) const;
Eci FindPositionSDP4(const Julian& dt, double tsince) const;
Eci FindPositionSGP4(const Julian& dt, double tsince) const;
Eci CalculateFinalPositionVelocity(const Julian& dt, const double& e,
const double& a, const double& omega, const double& xl, const double& xnode,
const double& xincl, const double& xlcof, const double& aycof,
const double& x3thm1, const double& x1mth2, const double& x7thm1,
const double& cosio, const double& sinio) const;
void DeepSpaceCalcDotTerms(struct IntegratorValues *values) const;
void DeepSpaceIntegrator(const double delt, const double step2,
Eci FindPositionSDP4(const double tsince) const;
Eci FindPositionSGP4(double tsince) const;
Eci CalculateFinalPositionVelocity(
const double tsince,
const double e,
const double a,
const double omega,
const double xl,
const double xnode,
const double xincl,
const double xlcof,
const double aycof,
const double x3thm1,
const double x1mth2,
const double x7thm1,
const double cosio,
const double sinio) const;
void DeepSpaceInitialise(
const double eosq,
const double sinio,
const double cosio,
const double betao,
const double theta2,
const double betao2,
const double xmdot,
const double omgdot,
const double xnodot);
void DeepSpaceCalculateLunarSolarTerms(
const double tsince,
double& pe,
double& pinc,
double& pl,
double& pgh,
double& ph) const;
void DeepSpacePeriodics(
const double tsince,
double& em,
double& xinc,
double& omgasm,
double& xnodes,
double& xll) const;
void DeepSpaceSecular(
const double tsince,
double& xll,
double& omgasm,
double& xnodes,
double& em,
double& xinc,
double& xn) const;
void DeepSpaceCalcDotTerms(struct IntegratorValues& values) const;
void DeepSpaceIntegrator(
const double delt,
const double step2,
const struct IntegratorValues& values) const;
void Reset();

1
libsgp4/SatelliteException.h
View File

@ -2,7 +2,6 @@
#define SATELLITEEXCEPTION_H_
#include <exception>
#include <iostream>
class SatelliteException : public std::exception
{

2
libsgp4/SolarPosition.cpp
View File

@ -5,7 +5,7 @@
#include <cmath>
Eci SolarPosition::FindPosition(const Julian& j)
Eci SolarPosition::FindPosition(const DateTime& dt)
{
const double mjd = j.FromJan1_12h_1900();
const double year = 1900 + mjd / 365.25;

4
libsgp4/SolarPosition.h
View File

@ -1,7 +1,7 @@
#ifndef SOLARPOSITION_H_
#define SOLARPOSITION_H_
#include "Julian.h"
#include "DateTime.h"
#include "Eci.h"
class SolarPosition
@ -15,7 +15,7 @@ public:
{
}
Eci FindPosition(const Julian& j);
Eci FindPosition(const DateTime& dt);
private:
double Delta_ET(double year) const;

1
libsgp4/TimeSpan.cpp Normal file
View File

@ -0,0 +1 @@
#include "TimeSpan.h"

194
libsgp4/TimeSpan.h Normal file
View File

@ -0,0 +1,194 @@
#ifndef TIMESPAN_H_
#define TIMESPAN_H_
namespace
{
static const long long TicksPerDay = 86400000000LL;
static const long long TicksPerHour = 3600000000LL;
static const long long TicksPerMinute = 60000000LL;
static const long long TicksPerSecond = 1000000LL;
static const long long TicksPerMillisecond = 1000LL;
static const long long TicksPerMicrosecond = 1LL;
static const long long UnixEpoch = 62135596800000000LL;
static const long long MaxValueTicks = 315537897599999999LL;
// 1582-Oct-15
static const long long GregorianStart = 49916304000000000LL;
}
class TimeSpan
{
public:
TimeSpan(long long ticks)
: m_ticks(ticks)
{
}
TimeSpan(int hours, int minutes, int seconds)
{
CalculateTicks(0, hours, minutes, seconds, 0);
}
TimeSpan(int days, int hours, int minutes, int seconds)
{
CalculateTicks(days, hours, minutes, seconds, 0);
}
TimeSpan(int days, int hours, int minutes, int seconds, int microseconds)
{
CalculateTicks(days, hours, minutes, seconds, microseconds);
}
TimeSpan Add(const TimeSpan& ts) const
{
return TimeSpan(m_ticks + ts.m_ticks);
}
TimeSpan Subtract(const TimeSpan& ts) const
{
return TimeSpan(m_ticks - ts.m_ticks);
}
int Compare(const TimeSpan& ts) const
{
int ret = 0;
if (m_ticks < ts.m_ticks)
{
ret = -1;
}
if (m_ticks < ts.m_ticks)
{
ret = 1;
}
return ret;
}
bool Equals(const TimeSpan& ts) const
{
return m_ticks == ts.m_ticks;
}
int Days()
{
return m_ticks / TicksPerDay;
}
int Hours()
{
return (m_ticks % TicksPerDay / TicksPerHour);
}
int Minutes() const
{
return (m_ticks % TicksPerHour / TicksPerMinute);
}
int Seconds() const
{
return (m_ticks % TicksPerMinute / TicksPerSecond);
}
int Milliseconds() const
{
return (m_ticks % TicksPerSecond / TicksPerMillisecond);
}
int Microseconds() const
{
return (m_ticks % TicksPerSecond / TicksPerMicrosecond);
}
long long Ticks() const
{
return m_ticks;
}
double TotalDays() const
{
return m_ticks / static_cast<double>(TicksPerDay);
}
double TotalHours() const
{
return m_ticks / static_cast<double>(TicksPerHour);
}
double TotalMinutes() const
{
return m_ticks / static_cast<double>(TicksPerMinute);
}
double TotalSeconds() const
{
return m_ticks / static_cast<double>(TicksPerSecond);
}
double TotalMilliseconds() const
{
return m_ticks / static_cast<double>(TicksPerMillisecond);
}
double TotalMicroseconds() const
{
return m_ticks / static_cast<double>(TicksPerMicrosecond);
}
private:
long long m_ticks;
void CalculateTicks(int days,
int hours,
int minutes,
int seconds,
int microseconds)
{
m_ticks = days * TicksPerDay +
(hours * 3600LL + minutes * 60LL + seconds) * TicksPerSecond +
microseconds * TicksPerMicrosecond;
}
};
inline TimeSpan operator+(const TimeSpan& ts1, const TimeSpan& ts2)
{
return ts1.Add(ts2);
}
inline TimeSpan operator-(const TimeSpan& ts1, const TimeSpan& ts2)
{
return ts1.Subtract(ts2);
}
inline bool operator==(const TimeSpan& ts1, TimeSpan& ts2)
{
return ts1.Equals(ts2);
}
inline bool operator>(const TimeSpan& ts1, const TimeSpan& ts2)
{
return (ts1.Compare(ts2) > 0);
}
inline bool operator>=(const TimeSpan& ts1, const TimeSpan& ts2)
{
return (ts1.Compare(ts2) >= 0);
}
inline bool operator!=(const TimeSpan& ts1, const TimeSpan& ts2)
{
return !ts1.Equals(ts2);
}
inline bool operator<(const TimeSpan& ts1, const TimeSpan& ts2)
{
return (ts1.Compare(ts2) < 0);
}
inline bool operator<=(const TimeSpan& ts1, const TimeSpan& ts2)
{
return (ts1.Compare(ts2) <= 0);
}
#endif

175
libsgp4/Timespan.cpp
View File

@ -1,175 +0,0 @@
#include "Timespan.h"
#include "Globals.h"
Timespan::Timespan()
: time_span_(0.0) {
}
Timespan::Timespan(const unsigned int days, const unsigned int hours,
const unsigned int minutes, const double seconds) {
SetValue(days, hours, minutes, seconds);
}
Timespan::Timespan(const double b) {
time_span_ = b;
}
Timespan::Timespan(const Timespan& b) {
time_span_ = b.time_span_;
}
Timespan::~Timespan(void) {
}
void Timespan::SetValue(const unsigned int days, const unsigned int hours,
const unsigned int minutes, const double seconds) {
time_span_ = static_cast<double> (days);
AddHours(hours);
AddMinutes(minutes);
AddSeconds(seconds);
}
void Timespan::AddDays(const unsigned int days) {
time_span_ += static_cast<double> (days);
}
void Timespan::AddHours(const unsigned int hours) {
time_span_ += (static_cast<double> (hours) / kHOURS_PER_DAY);
}
void Timespan::AddMinutes(const unsigned int minutes) {
time_span_ += (static_cast<double> (minutes) / kMINUTES_PER_DAY);
}
void Timespan::AddSeconds(const double seconds) {
time_span_ += (seconds / kSECONDS_PER_DAY);
}
double Timespan::GetTotalDays() const {
return time_span_;
}
double Timespan::GetTotalHours() const {
return time_span_ * kHOURS_PER_DAY;
}
double Timespan::GetTotalMinutes() const {
return time_span_ * kMINUTES_PER_DAY;
}
double Timespan::GetTotalSeconds() const {
return time_span_ * kSECONDS_PER_DAY;
}
Timespan& Timespan::operator =(const Timespan& b) {
if (this != &b) {
time_span_ = b.time_span_;
}
return (*this);
}
Timespan Timespan::operator +(const Timespan& b) const {
return Timespan(*this) += b;
}
Timespan Timespan::operator -(const Timespan& b) const {
return Timespan(*this) -= b;
}
Timespan Timespan::operator/(const double b) const {
return Timespan(*this) /= b;
}
Timespan Timespan::operator*(const double b) const {
return Timespan(*this) *= b;
}
Timespan & Timespan::operator+=(const Timespan& b) {
time_span_ += b.time_span_;
return (*this);
}
Timespan & Timespan::operator-=(const Timespan& b) {
time_span_ -= b.time_span_;
return (*this);
}
Timespan & Timespan::operator/=(const double b) {
time_span_ /= b;
return (*this);
}
Timespan & Timespan::operator*=(const double b) {
time_span_ *= b;
return (*this);
}
bool Timespan::operator ==(const Timespan& b) const {
if (time_span_ == b.time_span_)
return true;
else
return false;
}
bool Timespan::operator !=(const Timespan& b) const {
return !(*this == b);
}
bool Timespan::operator>(const Timespan& b) const {
if (time_span_ > b.time_span_)
return true;
else
return false;
}
bool Timespan::operator<(const Timespan& b) const {
if (time_span_ < b.time_span_)
return true;
else
return false;
}
bool Timespan::operator >=(const Timespan& b) const {
if (time_span_ >= b.time_span_)
return true;
else
return false;
}
bool Timespan::operator <=(const Timespan & b) const {
if (time_span_ <= b.time_span_)
return true;
else
return false;
}
double& operator +=(double& a, const Timespan& b) {
a += b.time_span_;
return a;
}
double& operator -=(double& a, const Timespan& b) {
a -= b.time_span_;
return a;
}

57
libsgp4/Timespan.h
View File

@ -1,57 +0,0 @@
#ifndef TIMESPAN_H_
#define TIMESPAN_H_
class Timespan {
public:
Timespan();
Timespan(const unsigned int days, const unsigned int hours,
const unsigned int minutes, const double seconds);
Timespan(const double b);
Timespan(const Timespan& b);
virtual ~Timespan(void);
void SetValue(const unsigned int days, const unsigned int hours,
const unsigned int minutes, const double seconds);
void AddDays(const unsigned int days);
void AddHours(const unsigned int hours);
void AddMinutes(const unsigned int minutes);
void AddSeconds(const double seconds);
double GetTotalDays() const;
double GetTotalHours() const;
double GetTotalMinutes() const;
double GetTotalSeconds() const;
// assignment
Timespan & operator=(const Timespan& b);
// arithmetic
Timespan operator+(const Timespan& b) const;
Timespan operator-(const Timespan& b) const;
Timespan operator/(const double b) const;
Timespan operator*(const double b) const;
// compound arithmetic
Timespan & operator+=(const Timespan& b);
Timespan & operator-=(const Timespan& b);
Timespan & operator/=(const double b);
Timespan & operator*=(const double b);
// comparison
bool operator==(const Timespan& b) const;
bool operator!=(const Timespan& b) const;
bool operator>(const Timespan& b) const;
bool operator<(const Timespan& b) const;
bool operator>=(const Timespan& b) const;
bool operator<=(const Timespan& b) const;
friend double& operator +=(double& a, const Timespan& b);
friend double& operator -=(double& a, const Timespan& b);
private:
/*
* stores value in minutes
*/
double time_span_;
};
#endif

626
libsgp4/Tle.cpp
View File

@ -1,14 +1,9 @@
#include "Tle.h"
#include "Util.h"
#include <locale>
namespace
{
/*
* line 1
*/
static const unsigned int TLE1_COL_NORADNUM = 2;
static const unsigned int TLE1_LEN_NORADNUM = 5;
static const unsigned int TLE1_COL_INTLDESC_A = 9;
@ -21,10 +16,10 @@ namespace
static const unsigned int TLE1_LEN_EPOCH_A = 2;
static const unsigned int TLE1_COL_EPOCH_B = 20;
static const unsigned int TLE1_LEN_EPOCH_B = 12;
static const unsigned int TLE1_COL_MEANMOTIONDT = 33;
static const unsigned int TLE1_LEN_MEANMOTIONDT = 10;
static const unsigned int TLE1_COL_MEANMOTIONDT2 = 44;
static const unsigned int TLE1_LEN_MEANMOTIONDT2 = 8;
static const unsigned int TLE1_COL_MEANMOTIONDT2 = 33;
static const unsigned int TLE1_LEN_MEANMOTIONDT2 = 10;
static const unsigned int TLE1_COL_MEANMOTIONDDT6 = 44;
static const unsigned int TLE1_LEN_MEANMOTIONDDT6 = 8;
static const unsigned int TLE1_COL_BSTAR = 53;
static const unsigned int TLE1_LEN_BSTAR = 8;
static const unsigned int TLE1_COL_EPHEMTYPE = 62;
@ -32,9 +27,6 @@ namespace
static const unsigned int TLE1_COL_ELNUM = 64;
static const unsigned int TLE1_LEN_ELNUM = 4;
/*
* line 2
*/
static const unsigned int TLE2_COL_NORADNUM = 2;
static const unsigned int TLE2_LEN_NORADNUM = 5;
static const unsigned int TLE2_COL_INCLINATION = 8;
@ -51,421 +43,309 @@ namespace
static const unsigned int TLE2_LEN_MEANMOTION = 11;
static const unsigned int TLE2_COL_REVATEPOCH = 63;
static const unsigned int TLE2_LEN_REVATEPOCH = 5;
}
Tle::Tle(const Tle& tle)
{
name_ = tle.name_;
line_one_ = tle.line_one_;
line_two_ = tle.line_two_;
norad_number_ = tle.norad_number_;
international_designator_ = tle.international_designator_;
epoch_ = tle.epoch_;
mean_motion_dot_ = tle.mean_motion_dot_;
mean_motion_dot2_ = tle.mean_motion_dot2_;
bstar_ = tle.bstar_;
inclination_ = tle.inclination_;
right_ascending_node_ = tle.right_ascending_node_;
eccentricity_ = tle.eccentricity_;
argument_perigee_ = tle.argument_perigee_;
mean_anomaly_ = tle.mean_anomaly_;
mean_motion_ = tle.mean_motion_;
orbit_number_ = tle.orbit_number_;
}
/*
* convert a tle raw string into an exponent string
*/
std::string Tle::ExpToDecimal(const std::string& str)
{
static const int START_SIGN = 0;
static const int LENGTH_SIGN = 1;
static const int START_MANTISSA = 1;
static const int LENGTH_MANTISSA = 5;
static const int START_EXP = 6;
static const int LENGTH_EXP = 2;
if ((LENGTH_SIGN + LENGTH_MANTISSA + LENGTH_EXP) != str.length())
{
throw TleException("Invalid string length for exponential conversion.");
}
std::string value = str.substr(START_SIGN, LENGTH_SIGN);
value += ".";
value += str.substr(START_MANTISSA, LENGTH_MANTISSA);
value += "e";
value += str.substr(START_EXP, LENGTH_EXP);
return value;
}
/*
* extract all variables
/**
* Initialise the tle object.
* @exception TleException
*/
void Tle::Initialize()
{
std::string temp;
/*
* trim whitespace
*/
Util::TrimLeft(name_);
Util::TrimRight(name_);
Util::TrimLeft(line_one_);
Util::TrimRight(line_one_);
Util::TrimLeft(line_two_);
Util::TrimRight(line_two_);
/*
* check the two lines are valid
*/
IsValidPair(line_one_, line_two_);
/*
* line 1
*/
temp = ExtractNoradNumber(line_one_, 1);
if (!Util::FromString<unsigned int>(temp, norad_number_))
if (!IsValidLineLength(line_one_))
{
throw TleException("Conversion failed");
throw TleException("Invalid length for line one");
}
/*
* if blank use norad number for name
*/
if (!IsValidLineLength(line_two_))
{
throw TleException("Invalid length for line two");
}
if (line_one_[0] != '1')
{
throw TleException("Invalid line beginning for line one");
}
if (line_two_[0] != '2')
{
throw TleException("Invalid line beginning for line two");
}
unsigned int sat_number_1;
unsigned int sat_number_2;
ExtractInteger(line_one_.substr(TLE1_COL_NORADNUM,
TLE1_LEN_NORADNUM), sat_number_1);
ExtractInteger(line_two_.substr(TLE2_COL_NORADNUM,
TLE2_LEN_NORADNUM), sat_number_2);
if (sat_number_1 != sat_number_2)
{
throw TleException("Satellite numbers do not match");
}
norad_number_ = sat_number_1;
if (name_.empty())
{
name_ = temp;
name_ = line_one_.substr(TLE1_COL_NORADNUM, TLE1_LEN_NORADNUM);
}
international_designator_ = line_one_.substr(TLE1_COL_INTLDESC_A,
int_designator_ = line_one_.substr(TLE1_COL_INTLDESC_A,
TLE1_LEN_INTLDESC_A + TLE1_LEN_INTLDESC_B + TLE1_LEN_INTLDESC_C);
int year;
double day;
if (!Util::FromString<int>(line_one_.substr(TLE1_COL_EPOCH_A,
TLE1_LEN_EPOCH_A), year))
{
throw TleException("Conversion failed");
}
if (!Util::FromString<double>(line_one_.substr(TLE1_COL_EPOCH_B,
TLE1_LEN_EPOCH_B), day))
{
throw TleException("Conversion failed");
}
/*
* generate julian date for epoch
*/
if (year < 57)
year += 2000;
else
year += 1900;
epoch_ = Julian(year, day);
unsigned int year = 0;
double day = 0.0;
if (line_one_[TLE1_COL_MEANMOTIONDT] == '-') {
temp = "-0";
} else
temp = "0";
temp += line_one_.substr(TLE1_COL_MEANMOTIONDT + 1, TLE1_LEN_MEANMOTIONDT);
if (!Util::FromString<double>(temp, mean_motion_dot_))
{
throw TleException("Conversion failed");
}
temp = ExpToDecimal(line_one_.substr(TLE1_COL_MEANMOTIONDT2,
TLE1_LEN_MEANMOTIONDT2));
if (!Util::FromString<double>(temp, mean_motion_dot2_))
{
throw TleException("Conversion failed");
}
temp = ExpToDecimal(line_one_.substr(TLE1_COL_BSTAR,
TLE1_LEN_BSTAR).c_str());
if (!Util::FromString<double>(temp, bstar_))
{
throw TleException("Conversion failed");
}
ExtractInteger(line_one_.substr(TLE1_COL_EPOCH_A,
TLE1_LEN_EPOCH_A), year);
ExtractDouble(line_one_.substr(TLE1_COL_EPOCH_B,
TLE1_LEN_EPOCH_B), 4, day);
ExtractDouble(line_one_.substr(TLE1_COL_MEANMOTIONDT2,
TLE1_LEN_MEANMOTIONDT2), 2, mean_motion_dt2_);
ExtractExponential(line_one_.substr(TLE1_COL_MEANMOTIONDDT6,
TLE1_LEN_MEANMOTIONDDT6), mean_motion_ddt6_);
ExtractExponential(line_one_.substr(TLE1_COL_BSTAR,
TLE1_LEN_BSTAR), bstar_);
/*
* line 2
*/
temp = line_two_.substr(TLE2_COL_INCLINATION, TLE2_LEN_INCLINATION);
Util::TrimLeft(temp);
if (!Util::FromString<double>(temp, inclination_))
{
throw TleException("Conversion failed");
}
temp = line_two_.substr(TLE2_COL_RAASCENDNODE, TLE2_LEN_RAASCENDNODE);
Util::TrimLeft(temp);
if (!Util::FromString<double>(temp, right_ascending_node_))
{
throw TleException("Conversion failed");
}
temp = "0.";
temp += line_two_.substr(TLE2_COL_ECCENTRICITY, TLE2_LEN_ECCENTRICITY);
if (!Util::FromString<double>(temp, eccentricity_))
{
throw TleException("Conversion failed");
}
temp = line_two_.substr(TLE2_COL_ARGPERIGEE, TLE2_LEN_ARGPERIGEE);
Util::TrimLeft(temp);
if (!Util::FromString<double>(temp, argument_perigee_))
{
throw TleException("Conversion failed");
}
temp = line_two_.substr(TLE2_COL_MEANANOMALY, TLE2_LEN_MEANANOMALY);
Util::TrimLeft(temp);
if (!Util::FromString<double>(temp, mean_anomaly_))
{
throw TleException("Conversion failed");
}
temp = line_two_.substr(TLE2_COL_MEANMOTION, TLE2_LEN_MEANMOTION);
Util::TrimLeft(temp);
if (!Util::FromString<double>(temp, mean_motion_))
{
throw TleException("Conversion failed");
}
temp = line_two_.substr(TLE2_COL_REVATEPOCH, TLE2_LEN_REVATEPOCH);
Util::TrimLeft(temp);
if (!Util::FromString<unsigned int>(temp, orbit_number_))
{
throw TleException("Conversion failed");
}
}
/*
* check the two lines have matching norad numbers
* and that the lines themselves are equal
*/
void Tle::IsValidPair(const std::string& line1, const std::string& line2)
{
/*
* validate each line
*/
IsValidLine(line1, 1);
IsValidLine(line2, 2);
/*
* extract norad numbers
*/
std::string norad_1 = ExtractNoradNumber(line1, 1);
std::string norad_2 = ExtractNoradNumber(line2, 2);
/*
* make sure they match
*/
if (norad_1.compare(norad_2) != 0)
{
throw TleException("Norad numbers do not match.");
}
}
/*
* validate a line
*/
void Tle::IsValidLine(const std::string& str, int line_number)
{
/*
* validation patterns
*/
static const std::string line1_pattern = "1 NNNNNC NNNNNXXX NNNNN.NNNNNNNN +.NNNNNNNN +NNNNN-N +NNNNN-N N NNNNN";
static const std::string line2_pattern = "2 NNNNN NNN.NNNN NNN.NNNN NNNNNNN NNN.NNNN NNN.NNNN NN.NNNNNNNNNNNNNN";
/*
* validate line against the pattern
*/
if (1 == line_number)
{
ValidateLine(str, line1_pattern);
}
else if (2 == line_number)
{
ValidateLine(str, line2_pattern);
}
ExtractDouble(line_two_.substr(TLE2_COL_INCLINATION,
TLE2_LEN_INCLINATION), 4, inclination_);
ExtractDouble(line_two_.substr(TLE2_COL_RAASCENDNODE,
TLE2_LEN_RAASCENDNODE), 4, right_ascending_node_);
ExtractDouble(line_two_.substr(TLE2_COL_ECCENTRICITY,
TLE2_LEN_ECCENTRICITY), -1, eccentricity_);
ExtractDouble(line_two_.substr(TLE2_COL_ARGPERIGEE,
TLE2_LEN_ARGPERIGEE), 4, argument_perigee_);
ExtractDouble(line_two_.substr(TLE2_COL_MEANANOMALY,
TLE2_LEN_MEANANOMALY), 4, mean_anomaly_);
ExtractDouble(line_two_.substr(TLE2_COL_MEANMOTION,
TLE2_LEN_MEANMOTION), 3, mean_motion_);
ExtractInteger(line_two_.substr(TLE2_COL_REVATEPOCH,
TLE2_LEN_REVATEPOCH), orbit_number_);
if (year < 57)
year += 2000;
else
{
throw TleException("Invalid line number to check.");
}
/*
* last char in string is modulo 10 checksum
* edited out as checksum isnt consistent
*
* int chk = CheckSum(str);
* if (chk != (str[TLE_LEN_LINE_DATA - 1] - '0'))
* return false;
*/
year += 1900;
epoch_ = DateTime(year, day);
}
/**
* Check
* @param str The string to check
* @returns Whether true of the string has a valid length
*/
bool Tle::IsValidLineLength(const std::string& str)
{
return str.length() == GetLineLength() ? true : false;
}
/*
* validate line given a pattern
/**
* Convert a string containing an integer
* @param[in] str The string to convert
* @param[out] val The result
* @exception TleException on conversion error
*/
void Tle::ValidateLine(const std::string& line, const std::string& pattern)
void Tle::ExtractInteger(const std::string& str, unsigned int& val)
{
/*
* check length of line
*/
if (!IsValidLineLength(line))
std::string temp;
bool found_digit = false;
for (std::string::const_iterator i = str.begin(); i != str.end(); ++i)
{
throw TleException("Invalid line length.");
if (isdigit(*i))
{
found_digit = true;
temp += *i;
}
else if (found_digit)
{
throw TleException("Unexpected non digit");
}
else if (*i != ' ')
{
throw TleException("Invalid character");
}
}
for (size_t i = 0; i < pattern.length(); i++)
if (temp.length() == 0)
{
char ptrn = pattern[i];
char mtch = line[i];
temp += '0';
}
switch (ptrn)
{
case '1':
case '2':
case ' ':
case '.':
{
/*
* should match exactly
*/
if (ptrn != mtch)
{
throw TleException("Invalid character");
}
break;
}
case 'N':
{
/*
* number or ' '
*/
if (!std::isdigit(mtch, std::locale::classic()) && mtch != ' ')
{
throw TleException("Invalid character");
}
break;
}
case '+':
{
/*
* + or - or space or 0 or digit
*/
if (mtch != '+' && mtch != '-'
&& mtch != ' ' && mtch != '0'
&& !std::isdigit(mtch, std::locale::classic()))
{
throw TleException("Invalid character");
}
break;
}
case '-':
{
/*
* + or -
*/
if (mtch != '+' && mtch != '-')
{
throw TleException("Invalid character");
}
break;
}
case 'C':
{
/*
* U or S
*/
if (mtch != 'U' && mtch != 'S')
{
throw TleException("Invalid character");
}
break;
}
case 'X':
{
/*
* alpha or ' '
*/
if (!std::isupper(mtch, std::locale::classic()) &&
!std::isalpha(mtch, std::locale::classic()) &&
mtch != ' ')
{
throw TleException("Invalid character");
}
break;
}
default:
{
throw TleException("Invalid pattern character");
break;
}
}
if (!Util::FromString<unsigned int>(temp, val))
{
throw TleException("Failed to convert value to integer");
}
}
/*
* compute checksum
/**
* Convert a string containing an double
* @param[in] str The string to convert
* @param[in] point_pos The position of the decimal point. (-1 if none)
* @param[out] val The result
* @exception TleException on conversion error
*/
int Tle::CheckSum(const std::string & str)
void Tle::ExtractDouble(const std::string& str, int point_pos, double& val)
{
size_t len = str.size() - 1;
int xsum = 0;
std::string temp;
bool found_digit = false;
for (size_t i = 0; i < len; i++)
for (std::string::const_iterator i = str.begin(); i != str.end(); ++i)
{
char ch = str[i];
/*
* integer part
*/
if (i < str.begin() + point_pos - 1)
{
bool done = false;
if (std::isdigit(ch, std::locale::classic()))
{
xsum += (ch - '0');
if (i == str.begin())
{
if(*i == '-' || *i == '+')
{
/*
* first character could be signed
*/
temp += *i;
done = true;
}
}
if (!done)
{
if (isdigit(*i))
{
found_digit = true;
temp += *i;
}
else if (found_digit)
{
throw TleException("Unexpected non digit");
}
else if (*i != ' ')
{
throw TleException("Invalid character");
}
}
}
else if (ch == '-')
/*
* decimal point
*/
else if (i == str.begin() + point_pos - 1)
{
xsum++;
if (temp.length() == 0)
{
/*
* integer part is blank, so add a '0'
*/
temp += '0';
}
if (*i == '.')
{
/*
* decimal point found
*/
temp += *i;
}
else
{
throw TleException("Failed to find decimal point");
}
}
/*
* fraction part
*/
else
{
if (i == str.begin() && point_pos == -1)
{
/*
* no decimal point expected, add 0. beginning
*/
temp += '0';
temp += '.';
}
/*
* should be a digit
*/
if (isdigit(*i))
{
temp += *i;
}
else
{
throw TleException("Invalid digit");
}
}
}
return (xsum % 10);
if (!Util::FromString<double>(temp, val))
{
throw TleException("Failed to convert value to double");
}
}
std::string Tle::ExtractNoradNumber(const std::string& str, int line_number)
/**
* Convert a string containing an exponential
* @param[in] str The string to convert
* @param[out] val The result
* @exception TleException on conversion error
*/
void Tle::ExtractExponential(const std::string& str, double& val)
{
std::string norad_number;
std::string temp;
/*
* check length
*/
if (!IsValidLineLength(str))
for (std::string::const_iterator i = str.begin(); i != str.end(); ++i)
{
throw TleException("Invalid line length.");
if (i == str.begin())
{
if (*i == '-' || *i == '+' || *i == ' ')
{
if (*i == '-')
{
temp += *i;
}
temp += '0';
temp += '.';
}
else
{
throw TleException("Invalid sign");
}
}
else if (i == str.begin() + str.length() - 2)
{
if (*i == '-' || *i == '+')
{
temp += 'e';
temp += *i;
}
else
{
throw TleException("Invalid exponential sign");
}
}
else
{
if (isdigit(*i))
{
temp += *i;
}
else
{
throw TleException("Invalid digit");
}
}
}
/*
* extract string
*/
if (line_number == 1)
if (!Util::FromString<double>(temp, val))
{
norad_number = str.substr(TLE1_COL_NORADNUM, TLE1_LEN_NORADNUM);
throw TleException("Failed to convert value to double");
}
else if (line_number == 2)
{
norad_number = str.substr(TLE2_COL_NORADNUM, TLE2_LEN_NORADNUM);
}
else
{
throw TleException("Invalid line number to check.");
}
return norad_number;
}

237
libsgp4/Tle.h
View File

@ -1,86 +1,164 @@
#ifndef TLE_H_
#define TLE_H_
#include "Globals.h"
#include "Util.h"
#include "Julian.h"
#include "DateTime.h"
#include "TleException.h"
#include <iostream>
/**
* Tle class
* @details Loads and validates a tle.
*/
class Tle
{
public:
Tle(const std::string& line_one, const std::string& line_two)
: line_one_(line_one), line_two_(line_two)
/**
* @details Initialise given the two lines of a tle
* @param[in] line_one Tle line one
* @param[in] line_two Tle line two
*/
Tle(const std::string& line_one,
const std::string& line_two)
: line_one_(line_one),
line_two_(line_two)
{
Initialize();
}
Tle(const std::string& name, const std::string& line_one,
const std::string& line_two)
: name_(name),line_one_(line_one), line_two_(line_two)
/**
* @details Initialise given the satellite name and the two lines of a tle
* @param[in] name Satellite name
* @param[in] line_one Tle line one
* @param[in] line_two Tle line two
*/
Tle(const std::string& name,
const std::string& line_one,
const std::string& line_two)
: name_(name),
line_one_(line_one),
line_two_(line_two)
{
Initialize();
}
Tle(const Tle& tle);
/**
* Copy constructor
* @param[in] tle Tle object to copy from
*/
Tle(const Tle& tle)
{
name_ = tle.name_;
line_one_ = tle.line_one_;
line_two_ = tle.line_two_;
norad_number_ = tle.norad_number_;
int_designator_ = tle.int_designator_;
epoch_ = tle.epoch_;
mean_motion_dt2_ = tle.mean_motion_dt2_;
mean_motion_ddt6_ = tle.mean_motion_ddt6_;
bstar_ = tle.bstar_;
inclination_ = tle.inclination_;
right_ascending_node_ = tle.right_ascending_node_;
eccentricity_ = tle.eccentricity_;
argument_perigee_ = tle.argument_perigee_;
mean_anomaly_ = tle.mean_anomaly_;
mean_motion_ = tle.mean_motion_;
orbit_number_ = tle.orbit_number_;
}
/**
* Destructor
*/
virtual ~Tle()
{
}
/*
* get raw strings
/**
* Get the satellite name
* @returns the satellite name
*/
std::string GetName() const
{
return name_;
}
/**
* Get the first line of the tle
* @returns the first line of the tle
*/
std::string GetLine1() const
{
return line_one_;
}
/**
* Get the second line of the tle
* @returns the second line of the tle
*/
std::string GetLine2() const
{
return line_two_;
}
/*
* get tle values
/**
* Get the norad number
* @returns the norad number
*/
unsigned int NoradNumber() const
{
return norad_number_;
}
std::string InternationlDesignator() const
/**
* Get the international designator
* @returns the international designator
*/
std::string IntDesignator() const
{
return international_designator_;
return int_designator_;
}
Julian Epoch() const
/**
* Get the tle epoch
* @returns the tle epoch
*/
DateTime Epoch() const
{
return epoch_;
}
double MeanMotionDot() const
/**
* Get the first time derivative of the mean motion divided by two
* @returns the first time derivative of the mean motion divided by two
*/
double MeanMotionDt2() const
{
return mean_motion_dot_;
return mean_motion_dt2_;
}
double MeanMotionDot2() const
/**
* Get the second time derivative of mean motion divided by six
* @returns the second time derivative of mean motion divided by six
*/
double MeanMotionDdt6() const
{
return mean_motion_dot2_;
return mean_motion_ddt6_;
}
/**
* Get the BSTAR drag term
* @returns the BSTAR drag term
*/
double BStar() const
{
return bstar_;
}
/**
* Get the inclination
* @param in_degrees Whether to return the value in degrees or radians
* @returns the inclination
*/
double Inclination(bool in_degrees) const
{
if (in_degrees)
@ -93,6 +171,11 @@ public:
}
}
/**
* Get the right ascension of the ascending node
* @param in_degrees Whether to return the value in degrees or radians
* @returns the right ascension of the ascending node
*/
double RightAscendingNode(const bool in_degrees) const
{
if (in_degrees)
@ -105,11 +188,20 @@ public:
}
}
/**
* Get the eccentricity
* @returns the eccentricity
*/
double Eccentricity() const
{
return eccentricity_;
}
/**
* Get the argument of perigee
* @param in_degrees Whether to return the value in degrees or radians
* @returns the argument of perigee
*/
double ArgumentPerigee(const bool in_degrees) const
{
if (in_degrees)
@ -122,6 +214,11 @@ public:
}
}
/**
* Get the mean anomaly
* @param in_degrees Whether to return the value in degrees or radians
* @returns the mean anomaly
*/
double MeanAnomaly(const bool in_degrees) const
{
if (in_degrees)
@ -134,67 +231,84 @@ public:
}
}
/**
* Get the mean motion
* @returns the mean motion (revolutions per day)
*/
double MeanMotion() const
{
return mean_motion_;
}
/**
* Get the orbit number
* @returns the orbit number
*/
unsigned int OrbitNumber() const
{
return orbit_number_;
}
/*
* helper / validation methods
/**
* Get the expected tle line length
* @returns the tle line length
*/
static unsigned int GetLineLength()
{
return TLE_LEN_LINE_DATA;
}
static void IsValidPair(const std::string& line1, const std::string& line2);
static void IsValidLine(const std::string& str, int line_number);
/**
* Dump this object to a string
* @returns string
*/
std::string ToString() const
{
std::stringstream ss;
ss << std::right << std::fixed;
ss << "Norad Number: " << NoradNumber() << std::endl;
ss << "Int. Designator: " << IntDesignator() << std::endl;
ss << "Epoch: " << Epoch() << std::endl;
ss << "Orbit Number: " << OrbitNumber() << std::endl;
ss << std::setprecision(8);
ss << "Mean Motion Dt2: ";
ss << std::setw(12) << MeanMotionDt2() << std::endl;
ss << "Mean Motion Ddt6: ";
ss << std::setw(12) << MeanMotionDdt6() << std::endl;
ss << "Eccentricity: ";
ss << std::setw(12) << Eccentricity() << std::endl;
ss << "BStar: ";
ss << std::setw(12) << BStar() << std::endl;
ss << "Inclination: ";
ss << std::setw(12) << Inclination(true) << std::endl;
ss << "Right Ascending Node: ";
ss << std::setw(12) << RightAscendingNode(true) << std::endl;
ss << "Argument Perigee: ";
ss << std::setw(12) << ArgumentPerigee(true) << std::endl;
ss << "Mean Anomaly: ";
ss << std::setw(12) << MeanAnomaly(true) << std::endl;
ss << "Mean Motion: ";
ss << std::setw(12) << MeanMotion() << std::endl;
return ss.str();
}
private:
/*
* initialize from raw tle strings
*/
void Initialize();
/*
* format a raw string into an exponent string
*/
static std::string ExpToDecimal(const std::string&);
/*
* validate a line against a pattern
*/
static void ValidateLine(const std::string& line,
const std::string& pattern);
/*
* compute checksum
*/
static int CheckSum(const std::string& str);
static std::string ExtractNoradNumber(const std::string& str,
int line_number);
static bool IsValidLineLength(const std::string& str);
void ExtractInteger(const std::string& str, unsigned int& val);
void ExtractDouble(const std::string& str, int point_pos, double& val);
void ExtractExponential(const std::string& str, double& val);
private:
/*
* raw tle data
*/
std::string name_;
std::string line_one_;
std::string line_two_;
/*
* extracted values all in native units
*/
unsigned int norad_number_;
std::string international_designator_;
Julian epoch_;
double mean_motion_dot_;
double mean_motion_dot2_;
std::string int_designator_;
DateTime epoch_;
double mean_motion_dt2_;
double mean_motion_ddt6_;
double bstar_;
double inclination_;
double right_ascending_node_;
@ -204,13 +318,14 @@ private:
double mean_motion_;
unsigned int orbit_number_;
/*
* line lengths
*/
static const unsigned int TLE_LEN_LINE_DATA = 69;
static const unsigned int TLE_LEN_LINE_NAME = 22;
};
};
inline std::ostream& operator<<(std::ostream& strm, const Tle& t)
{
return strm << t.ToString();
}
#endif

20
libsgp4/TleException.h
View File

@ -2,28 +2,38 @@
#define TLEEXCEPTION_H_
#include <exception>
#include <iostream>
class TleException : public std::exception
{
public:
/**
* Constructor
* @param message Exception message
*/
TleException(const char* message)
: message_(message)
: m_message(message)
{
}
/**
* Destructor
*/
virtual ~TleException(void) throw ()
{
}
/**
* Get the exception message
* @returns the exception message
*/
virtual const char* what() const throw ()
{
return message_.c_str();
return m_message.c_str();
}
private:
std::string message_;
/** the exception message */
std::string m_message;
};
#endif

1
libsgp4/Util.cpp
View File

@ -34,5 +34,4 @@ namespace Util
TrimLeft(s);
TrimRight(s);
}
}

2
libsgp4/Util.h
View File

@ -84,7 +84,7 @@ namespace Util
}
}
}
void TrimLeft(std::string& s);
void TrimRight(std::string& s);
void Trim(std::string& s);

86
libsgp4/Vector.h
View File

@ -10,21 +10,46 @@ struct Vector
{
public:
/**
* Default constructor
*/
Vector()
: x(0.0), y(0.0), z(0.0), w(0.0)
{
}
Vector(double xx, double yy, double zz)
: x(xx), y(yy), z(zz), w(0.0)
/**
* Constructor
* @param arg_x x value
* @param arg_y y value
* @param arg_z z value
*/
Vector(const double arg_x,
const double arg_y,
const double arg_z)
: x(arg_x), y(arg_y), z(arg_z), w(0.0)
{
}
Vector(double xx, double yy, double zz, double ww)
: x(xx), y(yy), z(zz), w(ww)
/**
* Constructor
* @param arg_x x value
* @param arg_y y value
* @param arg_z z value
* @param arg_w w value
*/
Vector(const double arg_x,
const double arg_y,
const double arg_z,
const double arg_w)
: x(arg_x), y(arg_y), z(arg_z), w(arg_w)
{
}
/**
* Copy constructor
* @param v value to copy from
*/
Vector(const Vector& v)
{
x = v.x;
@ -33,10 +58,17 @@ public:
w = v.w;
}
/**
* Destructor
*/
virtual ~Vector()
{
}
/**
* Assignment operator
* @param v value to copy from
*/
Vector& operator=(const Vector& v)
{
if (this != &v)
@ -49,23 +81,31 @@ public:
return *this;
}
double GetMagnitude() const
/**
* Subtract operator
* @param v value to suctract from
*/
Vector operator-(const Vector& v)
{
return sqrt(x * x + y * y + z * z);
}
Vector Subtract(const Vector& v) const
{
/*
* subtract (this) - (v)
* and return result
*/
return Vector(x - v.x,
y - v.y,
z - v.z,
0.0);
}
/**
* Calculates the magnitude of the vector
* @returns magnitude of the vector
*/
double GetMagnitude() const
{
return sqrt(x * x + y * y + z * z);
}
/**
* Calculates the dot product
* @returns dot product
*/
double Dot(const Vector& vec) const
{
return (x * vec.x) +
@ -73,20 +113,28 @@ public:
(z * vec.z);
}
/**
* Converts this vector to a string
* @returns this vector as a string
*/
std::string ToString() const
{
std::stringstream ss;
ss << std::right << std::fixed << std::setprecision(2);
ss << "X: " << std::setw(8) << x;
ss << ", Y: " << std::setw(8) << y;
ss << ", Z: " << std::setw(8) << z;
ss << ", W: " << std::setw(8) << w;
ss << std::right << std::fixed << std::setprecision(3);
ss << "X: " << std::setw(9) << x;
ss << ", Y: " << std::setw(9) << y;
ss << ", Z: " << std::setw(9) << z;
ss << ", W: " << std::setw(9) << w;
return ss.str();
}
/** x value */
double x;
/** y value */
double y;
/** z value */
double z;
/** w value */
double w;
};

321
passpredict/passpredict.cpp
View File

@ -1,17 +1,30 @@
#include <Observer.h>
#include <SGP4.h>
#include <Util.h>
#include <CoordTopographic.h>
#include <CoordGeodetic.h>
#include <cmath>
#include <iostream>
#include <list>
const int kPASS_TIME_STEP = 180; // time step in seconds, when searching for aos/los
double FindMaxElevation(const CoordGeodetic& user_geo, SGP4& sgp4, const Julian& aos, const Julian& los) {
struct PassDetails
{
DateTime aos;
DateTime los;
double max_elevation;
};
double FindMaxElevation(
const CoordGeodetic& user_geo,
SGP4& sgp4,
const DateTime& aos,
const DateTime& los)
{
Observer obs(user_geo);
double max_elevation = -99.9;
double max_elevation = 0.0;
/*
* time step in seconds
*/
@ -25,9 +38,9 @@ double FindMaxElevation(const CoordGeodetic& user_geo, SGP4& sgp4, const Julian&
*/
bool fine_tune = false;
Julian current_time = aos;
while (current_time < los && searching) {
DateTime current_time(aos);
while (current_time < los && searching)
{
/*
* find position
*/
@ -37,16 +50,18 @@ double FindMaxElevation(const CoordGeodetic& user_geo, SGP4& sgp4, const Julian&
/*
* keep updating max elevation
*/
if (topo.elevation > max_elevation) {
if (topo.elevation > max_elevation)
{
max_elevation = topo.elevation;
} else if (!fine_tune) {
}
else if (!fine_tune)
{
/*
* passed max elevation
* max elevation happened in the last 6 minutes
* go back and fine tune max elevation value
*/
current_time.AddSec(-2.0 * time_step);
current_time = current_time.AddSeconds(-2.0 * time_step);
/*
* dont go back before aos
*/
@ -64,7 +79,9 @@ double FindMaxElevation(const CoordGeodetic& user_geo, SGP4& sgp4, const Julian&
*/
max_elevation = -99.9;
} else {
}
else
{
/*
* found max elevation
*/
@ -72,48 +89,72 @@ double FindMaxElevation(const CoordGeodetic& user_geo, SGP4& sgp4, const Julian&
searching = false;
}
if (searching) {
current_time.AddSec(time_step);
if (searching)
{
current_time = current_time.AddSeconds(time_step);
if (current_time > los)
{
current_time = los;
}
}
};
}
return max_elevation;
}
Julian FindCrossingPoint(const CoordGeodetic& user_geo, SGP4& sgp4, const Julian& initial_time1, const Julian& initial_time2, bool finding_aos) {
DateTime FindCrossingPoint(
const CoordGeodetic& user_geo,
SGP4& sgp4,
const DateTime& initial_time1,
const DateTime& initial_time2,
bool finding_aos)
{
Observer obs(user_geo);
bool searching = true;
unsigned int loop_count = 0;
int cnt = 0;
Julian time1 = initial_time1;
Julian time2 = initial_time2;
DateTime time1(initial_time1);
DateTime time2(initial_time2);
Julian middle_time = time1 + (time2 - time1) / 2.0;
while (searching && loop_count < 25) {
double diff = (time2 - time1).TotalSeconds();
if (finding_aos)
{
diff = std::floor(diff);
}
else
{
diff = std::ceil(diff);
}
DateTime middle_time(time1.AddSeconds(diff));
while (searching && cnt < 25)
{
/*
* find position
* calculate satellite position
*/
Eci eci = sgp4.FindPosition(middle_time);
CoordTopographic topo = obs.GetLookAngle(eci);
if (topo.elevation > 0.0) {
if (finding_aos) {
if (topo.elevation > 0.0)
{
if (finding_aos)
{
time2 = middle_time;
} else {
}
else
{
time1 = middle_time;
}
} else {
if (finding_aos) {
}
else
{
if (finding_aos)
{
time1 = middle_time;
} else {
}
else
{
time2 = middle_time;
}
}
@ -121,120 +162,206 @@ Julian FindCrossingPoint(const CoordGeodetic& user_geo, SGP4& sgp4, const Julian
/*
* when two times are within a second, stop
*/
if (((time2.GetDate() - time1.GetDate()) * kSECONDS_PER_DAY) < 0.5) {
if ((time2 - time1).TotalSeconds() < 1.5)
{
searching = false;
}
middle_time = time1 + (time2 - time1) / 2.0;
loop_count++;
};
else
{
diff = (time2 - time1).TotalSeconds();
if (finding_aos)
{
diff = std::floor(diff);
}
else
{
diff = std::ceil(diff);
}
middle_time = time1.AddSeconds(diff);
}
cnt++;
}
return middle_time;
}
void AOSLOS(const CoordGeodetic& user_geo, SGP4& sgp4, const Julian& start_time, const Julian& end_time) {
std::list<struct PassDetails> GeneratePassList(
const CoordGeodetic& user_geo,
SGP4& sgp4,
const DateTime& start_time,
const DateTime& end_time,
const int time_step)
{
std::list<struct PassDetails> pass_list;
Observer obs(user_geo);
Timespan time_step(0, 0, 0, kPASS_TIME_STEP);
DateTime aos_time;
DateTime los_time;
bool first_run = true;
bool end_of_pass = false;
Julian previous_time = start_time;
Julian aos_time;
Julian los_time;
bool found_aos = false;
bool found_los = false;
Julian current_time = start_time;
while (current_time < end_time) {
DateTime previous_time(start_time);
DateTime current_time(start_time);
// find position
while (current_time < end_time)
{
/*
* calculate satellite position
*/
Eci eci = sgp4.FindPosition(current_time);
CoordTopographic topo = obs.GetLookAngle(eci);
if (topo.elevation > 0.0) {
std::cout << std::fixed << current_time << " " << topo.elevation << std::endl;
if (topo.elevation > 0.0)
{
/*
* satellite is above horizon
* and its the first run, so just use start_time
*/
if (first_run) {
aos_time = start_time;
found_aos = true;
found_los = false;
if (start_time == current_time)
{
/*
* not the first iteration and satellite has
* gone above horizon, so find AOS
* satellite was already above the horizon at the start time,
* so use the start time
*/
} else if (!found_aos) {
aos_time = FindCrossingPoint(user_geo, sgp4, previous_time, current_time, true);
found_aos = true;
found_los = false;
aos_time = start_time;
}
} else {
else
{
/*
* find the point at which the satellite crossed the horizon
*/
aos_time = FindCrossingPoint(
user_geo,
sgp4,
previous_time,
current_time,
true);
}
found_aos = true;
found_los = false;
}
else if (found_aos)
{
/*
* if satellite now below horizon and have previously
* found AOS, find LOS
* satellite now below horizon and we have an AOS, so record the
* pass
*/
if (found_aos) {
los_time = FindCrossingPoint(
user_geo,
sgp4,
previous_time,
current_time,
false);
los_time = FindCrossingPoint(user_geo, sgp4, previous_time, current_time, false);
found_aos = false;
found_los = false;
end_of_pass = true;
double max_elevation = FindMaxElevation(user_geo, sgp4, aos_time, los_time);
std::cout << "AOS: " << aos_time << ", LOS: " << los_time << ", Max El: " << Util::RadiansToDegrees(max_elevation) << std::endl;
found_aos = false;
found_los = true;
}
struct PassDetails pd;
pd.aos = aos_time;
pd.los = los_time;
pd.max_elevation = FindMaxElevation(
user_geo,
sgp4,
aos_time,
los_time);
pass_list.push_back(pd);
}
first_run = false;
previous_time = current_time;
// at end of pass, move time along by 20 minutes
if (end_of_pass)
current_time += Timespan(0, 0, 20, 0);
if (found_los)
{
/*
* at the end of the pass move the time along by 20mins
*/
current_time = current_time + TimeSpan(0, 20, 0);
}
else
current_time += time_step;
{
/*
* move the time along by the time step value
*/
current_time = current_time + TimeSpan(0, 0, time_step);
}
// check we dont go past end time
if (current_time > end_time)
{
/*
* dont go past end time
*/
current_time = end_time;
}
end_of_pass = false;
found_los = false;
};
/*
* is satellite still above horizon at end of search period
*/
if (found_aos && !found_los) {
los_time = end_time;
double max_elevation = FindMaxElevation(user_geo, sgp4, aos_time, los_time);
std::cout << "AOS: " << aos_time << ", LOS: " << los_time << ", Max El: " << Util::RadiansToDegrees(max_elevation) << std::endl;
if (found_aos)
{
/*
* satellite still above horizon at end of search period, so use end
* time as los
*/
struct PassDetails pd;
pd.aos = aos_time;
pd.los = end_time;
pd.max_elevation = FindMaxElevation(user_geo, sgp4, aos_time, los_time);
pass_list.push_back(pd);
}
return pass_list;
}
int main() {
CoordGeodetic geo(51.507406923983446, -0.12773752212524414, 0.05);
Tle tle("GIOVE-B ",
"1 32781U 08020A 11158.03814084 .00000088 00000-0 10000-3 0 4620",
"2 32781 55.9142 172.9458 0022365 228.3743 131.4697 1.70953903 19437");
Tle tle("UK-DMC 2 ",
"1 25544U 98067A 12285.65009259 .00017228 00000-0 30018-3 0 4501",
"2 25544 051.6477 262.7396 0017757 155.0745 185.1532 15.50683239796101");
SGP4 sgp4(tle);
Julian start_date;
Julian end_date(start_date);
end_date.AddDay(10.0);
std::cout << tle << std::endl;
/*
* generate 10 day schedule
* generate 1 day schedule
*/
AOSLOS(geo, sgp4, start_date, end_date);
DateTime start_date = DateTime::Now();
DateTime end_date(start_date.AddDays(1.0));
std::list<struct PassDetails> pass_list;
std::list<struct PassDetails>::const_iterator itr;
std::cout << "Start time: " << start_date << std::endl;
std::cout << "End time : " << end_date << std::endl << std::endl;
/*
* generate passes
*/
pass_list = GeneratePassList(geo, sgp4, start_date, end_date, 180);
if (pass_list.begin() == pass_list.end())
{
std::cout << "No passes found" << std::endl;
}
else
{
itr = pass_list.begin();
do
{
std::cout
<< "AOS: " << itr->aos
<< ", LOS: " << itr->los
<< ", Max El: " << Util::RadiansToDegrees(itr->max_elevation)
<< std::endl;
}
while (++itr != pass_list.end());
}
return 0;
}

10
runtest/runtest.cpp
View File

@ -1,8 +1,5 @@
#include <Julian.h>
#include <Tle.h>
#include <SGP4.h>
#include <Globals.h>
#include <Util.h>
#include <Observer.h>
#include <CoordGeodetic.h>
#include <CoordTopographic.h>
@ -187,7 +184,7 @@ void RunTest(const char* infile)
{
if (line.length() >= Tle::GetLineLength())
{
Tle::IsValidLine(line.substr(0, Tle::GetLineLength()), 1);
//Tle::IsValidLine(line.substr(0, Tle::GetLineLength()), 1);
/*
* store line and now read in second line
*/
@ -236,7 +233,7 @@ void RunTest(const char* infile)
{
if (line.length() >= Tle::GetLineLength())
{
Tle::IsValidLine(line.substr(0, Tle::GetLineLength()), 2);
//Tle::IsValidLine(line.substr(0, Tle::GetLineLength()), 2);
Tle tle("Test", line1, line2);
RunTle(tle, start, end, inc);
}
@ -265,6 +262,3 @@ int main()
return 1;
}

45
sattrack/sattrack.cpp
View File

@ -1,3 +1,5 @@
#include <CoordTopographic.h>
#include <CoordGeodetic.h>
#include <Observer.h>
#include <SGP4.h>
@ -5,22 +7,33 @@
int main()
{
Observer obs(51.507406923983446, -0.12773752212524414, 0.05);
Tle tle = Tle("UK-DMC 2 ",
"1 35683U 09041C 11356.17214994 .00000411 00000-0 77254-4 0 6826",
"2 35683 98.0512 251.8127 0001492 79.4611 280.6776 14.69611889128518");
SGP4 sgp4(tle);
Observer obs(51.507406923983446, -0.12773752212524414, 0.05);
Tle tle = Tle("MASAT 1 ",
"1 25544U 98067A 12285.65009259 .00017228 00000-0 30018-3 0 4501",
"2 25544 051.6477 262.7396 0017757 155.0745 185.1532 15.50683239796101");
SGP4 sgp4(tle);
while(1)
{
Julian now;
Eci eci = sgp4.FindPosition(now);
CoordTopographic topo = obs.GetLookAngle(eci);
CoordGeodetic geo = eci.ToGeodetic();
std::cout << now << " ";
std::cout << topo << " ";
std::cout << geo << std::endl;
};
while (true)
{
/*
* current time
*/
DateTime now = DateTime::Now(true);
/*
* calculate satellite position
*/
Eci eci = sgp4.FindPosition(now);
/*
* get look angle for observer to satellite
*/
CoordTopographic topo = obs.GetLookAngle(eci);
/*
* convert satellite position to geodetic coordinates
*/
CoordGeodetic geo = eci.ToGeodetic();
return 0;
std::cout << now << " " << topo << " " << geo << std::endl;
};
return 0;
}