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

Updated Eci

feature/19
Daniel Warner 2011-12-13 22:37:13 +00:00
parent 0d17be39fa
commit 4aab9bf601
9 changed files with 104 additions and 79 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
CoordGeodetic.h
View File

@ -16,11 +16,21 @@ public:
}
/*
* radians
* default is in degrees
*/
CoordGeodetic(double lat, double lon, double alt)
: latitude(lat), longitude(lon), altitude(alt)
CoordGeodetic(double lat, double lon, double alt, bool radians = false)
{
if (radians)
{
latitude = lat;
longitude = lon;
}
else
{
latitude = DegreesToRadians(lat);
longitude = DegreesToRadians(lon);
}
altitude = alt;
}
CoordGeodetic(const CoordGeodetic& g)

49
Eci.cpp
View File

@ -1,26 +1,26 @@
#include "Eci.h"
#include "Globals.h"
Eci::Eci(const Julian &date, const CoordGeodetic &geo)
: date_(date) {
void Eci::ToEci(const Julian& date, const CoordGeodetic &g)
{
/*
* set date
*/
date_ = date;
static const double mfactor = kTWOPI * (kOMEGA_E / kSECONDS_PER_DAY);
const double latitude = geo.latitude;
const double longitude = geo.longitude;
const double altitude = geo.altitude;
/*
* Calculate Local Mean Sidereal Time for observers longitude
*/
const double theta = date_.ToLocalMeanSiderealTime(longitude);
const double theta = date_.ToLocalMeanSiderealTime(g.longitude);
/*
* take into account earth flattening
*/
const double c = 1.0 / sqrt(1.0 + kF * (kF - 2.0) * pow(sin(latitude), 2.0));
const double c = 1.0 /
sqrt(1.0 + kF * (kF - 2.0) * pow(sin(g.latitude), 2.0));
const double s = pow(1.0 - kF, 2.0) * c;
const double achcp = (kXKMPER * c + altitude) * cos(latitude);
const double achcp = (kXKMPER * c + g.altitude) * cos(g.latitude);
/*
* X position in km
@ -30,7 +30,7 @@ Eci::Eci(const Julian &date, const CoordGeodetic &geo)
*/
position_.x = achcp * cos(theta);
position_.y = achcp * sin(theta);
position_.z = (kXKMPER * s + altitude) * sin(latitude);
position_.z = (kXKMPER * s + g.altitude) * sin(g.latitude);
position_.w = position_.GetMagnitude();
/*
@ -45,21 +45,8 @@ Eci::Eci(const Julian &date, const CoordGeodetic &geo)
velocity_.w = velocity_.GetMagnitude();
}
Eci::Eci(const Julian &date, const Vector &position)
: date_(date), position_(position) {
}
Eci::Eci(const Julian &date, const Vector &position, const Vector &velocity)
: date_(date), position_(position), velocity_(velocity) {
}
Eci::~Eci(void) {
}
CoordGeodetic Eci::ToGeodetic() const {
CoordGeodetic Eci::ToGeodetic() const
{
const double theta = AcTan(position_.y, position_.x);
// 0 >= lon < 360
@ -67,7 +54,9 @@ CoordGeodetic Eci::ToGeodetic() const {
// 180 >= lon < 180
const double lon = fmod(theta - date_.ToGreenwichSiderealTime(), kPI);
const double r = sqrt((position_.x * position_.x) + (position_.y * position_.y));
const double r = sqrt((position_.x * position_.x) +
(position_.y * position_.y));
static const double e2 = kF * (2.0 - kF);
double lat = AcTan(position_.z, r);
@ -75,13 +64,15 @@ CoordGeodetic Eci::ToGeodetic() const {
double c = 0.0;
int cnt = 0;
do {
do
{
phi = lat;
const double sinphi = sin(phi);
c = 1.0 / sqrt(1.0 - e2 * sinphi * sinphi);
lat = AcTan(position_.z + kXKMPER * c * e2 * sinphi, r);
cnt++;
} while (fabs(lat - phi) >= 1e-10 && cnt < 10);
}
while (fabs(lat - phi) >= 1e-10 && cnt < 10);
const double alt = r / cos(lat) - kXKMPER * c;

50
Eci.h
View File

@ -4,31 +4,61 @@
#include "CoordGeodetic.h"
#include "Vector.h"
#include "Julian.h"
#include "Globals.h"
class Eci {
class Eci
{
public:
Eci() {
};
Eci(const Julian &date, const CoordGeodetic &geo);
Eci(const Julian &date, const Vector &position);
Eci(const Julian &date, const Vector &position, const Vector &velocity);
virtual ~Eci(void);
/*
* in degrees
*/
Eci(const Julian& date, double latitude, double longitude, double altitude)
{
ToEci(date, CoordGeodetic(latitude, longitude, altitude));
}
Vector GetPosition() const {
Eci(const Julian& date, const CoordGeodetic& g)
{
ToEci(date, g);
}
Eci(const Julian &date, const Vector &position)
: date_(date), position_(position)
{
}
Eci(const Julian &date, const Vector &position, const Vector &velocity)
: date_(date), position_(position), velocity_(velocity)
{
}
virtual ~Eci()
{
}
Vector GetPosition() const
{
return position_;
}
Vector GetVelocity() const {
Vector GetVelocity() const
{
return velocity_;
}
Julian GetDate() const {
Julian GetDate() const
{
return date_;
}
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_;

16
Observer.cpp
View File

@ -5,19 +5,15 @@
/*
* in degrees!
*/
Observer::Observer(const double latitude, const double longitude, const double altitude) {
geo_.latitude = DegreesToRadians(latitude);
geo_.longitude = DegreesToRadians(longitude);
geo_.altitude = altitude;
observers_eci_ = Eci(Julian(), geo_);
Observer::Observer(const double latitude, const double longitude, const double altitude)
: geo_(DegreesToRadians(latitude), DegreesToRadians(longitude), altitude),
observers_eci_(Julian(), geo_)
{
}
Observer::Observer(const CoordGeodetic &geo)
: geo_(geo) {
observers_eci_ = Eci(Julian(), geo_);
: geo_(geo), observers_eci_(Julian(), geo_)
{
}
Observer::~Observer(void) {

9
PassPredict.cpp
View File

@ -9,7 +9,6 @@ const int kPASS_TIME_STEP = 180; // time step in seconds, when searching for aos
double FindMaxElevation(const CoordGeodetic& user_geo, SGP4& sgp4, const Julian& aos, const Julian& los) {
Observer obs(user_geo);
Eci eci;
double max_elevation = -99.9;
/*
@ -31,7 +30,7 @@ double FindMaxElevation(const CoordGeodetic& user_geo, SGP4& sgp4, const Julian&
/*
* find position
*/
sgp4.FindPosition(&eci, current_time);
Eci eci = sgp4.FindPosition(current_time);
CoordTopographic topo = obs.GetLookAngle(eci);
/*
@ -85,7 +84,6 @@ double FindMaxElevation(const CoordGeodetic& user_geo, SGP4& sgp4, const Julian&
Julian FindCrossingPoint(const CoordGeodetic& user_geo, SGP4& sgp4, const Julian& initial_time1, const Julian& initial_time2, bool finding_aos) {
Observer obs(user_geo);
Eci eci;
bool searching = true;
unsigned int loop_count = 0;
@ -100,7 +98,7 @@ Julian FindCrossingPoint(const CoordGeodetic& user_geo, SGP4& sgp4, const Julian
/*
* find position
*/
sgp4.FindPosition(&eci, middle_time);
Eci eci = sgp4.FindPosition(middle_time);
CoordTopographic topo = obs.GetLookAngle(eci);
if (topo.elevation > 0.0) {
@ -136,7 +134,6 @@ Julian FindCrossingPoint(const CoordGeodetic& user_geo, SGP4& sgp4, const Julian
void AOSLOS(const CoordGeodetic& user_geo, SGP4& sgp4, const Julian& start_time, const Julian& end_time) {
Observer obs(user_geo);
Eci eci;
Timespan time_step(0, 0, 0, kPASS_TIME_STEP);
@ -153,7 +150,7 @@ void AOSLOS(const CoordGeodetic& user_geo, SGP4& sgp4, const Julian& start_time,
while (current_time < end_time) {
// find position
sgp4.FindPosition(&eci, current_time);
Eci eci = sgp4.FindPosition(current_time);
CoordTopographic topo = obs.GetLookAngle(eci);
if (topo.elevation > 0.0) {

3
RunTest.cpp
View File

@ -23,7 +23,6 @@ void RunTle(Tle tle, double start, double end, double inc) {
while (running) {
try {
double val;
Eci eci;
if (first_run && current != 0.0) {
/*
* make sure first run is always as zero
@ -35,7 +34,7 @@ void RunTle(Tle tle, double start, double end, double inc) {
*/
val = current;
}
model.FindPosition(&eci, val);
Eci eci = model.FindPosition(val);
Vector position = eci.GetPosition();
Vector velocity = eci.GetVelocity();

24
SGP4.cpp
View File

@ -178,22 +178,22 @@ void SGP4::Initialize() {
first_run_ = false;
}
void SGP4::FindPosition(Eci* eci, double tsince) const {
Eci SGP4::FindPosition(double tsince) const {
if (use_deep_space_)
FindPositionSDP4(eci, tsince);
return FindPositionSDP4(tsince);
else
FindPositionSGP4(eci, tsince);
return FindPositionSGP4(tsince);
}
void SGP4::FindPosition(Eci* eci, const Julian& date) const {
Eci SGP4::FindPosition(const Julian& date) const {
Timespan diff = date - elements_.Epoch();
FindPosition(eci, diff.GetTotalMinutes());
return FindPosition(diff.GetTotalMinutes());
}
void SGP4::FindPositionSDP4(Eci* eci, double tsince) const {
Eci SGP4::FindPositionSDP4(double tsince) const {
/*
* the final values
@ -286,7 +286,7 @@ void SGP4::FindPositionSDP4(Eci* eci, double tsince) const {
/*
* using calculated values, find position and velocity
*/
CalculateFinalPositionVelocity(eci, tsince, e,
return CalculateFinalPositionVelocity(tsince, e,
a, omega, xl, xnode,
xincl, perturbed_xlcof, perturbed_aycof,
perturbed_x3thm1, perturbed_x1mth2, perturbed_x7thm1,
@ -294,7 +294,7 @@ void SGP4::FindPositionSDP4(Eci* eci, double tsince) const {
}
void SGP4::FindPositionSGP4(Eci* eci, double tsince) const {
Eci SGP4::FindPositionSGP4(double tsince) const {
/*
* the final values
@ -363,7 +363,7 @@ void SGP4::FindPositionSGP4(Eci* eci, double tsince) const {
* using calculated values, find position and velocity
* we can pass in constants from Initialize() as these dont change
*/
CalculateFinalPositionVelocity(eci, tsince, e,
return CalculateFinalPositionVelocity(tsince, e,
a, omega, xl, xnode,
xincl, common_consts_.xlcof, common_consts_.aycof,
common_consts_.x3thm1, common_consts_.x1mth2, common_consts_.x7thm1,
@ -371,7 +371,7 @@ void SGP4::FindPositionSGP4(Eci* eci, double tsince) const {
}
void SGP4::CalculateFinalPositionVelocity(Eci* eci, const double& tsince, const double& e,
Eci SGP4::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,
@ -522,7 +522,9 @@ void SGP4::CalculateFinalPositionVelocity(Eci* eci, const double& tsince, const
Julian julian(elements_.Epoch());
julian.AddMin(tsince);
(*eci) = Eci(julian, position, velocity);
Eci eci(julian, position, velocity);
return eci;
}
static inline double EvaluateCubicPolynomial(const double x, const double constant,

13
SGP4.h
View File

@ -6,14 +6,15 @@
#include "Eci.h"
#include "SatelliteException.h"
class SGP4 {
class SGP4
{
public:
SGP4(const Tle& tle);
virtual ~SGP4(void);
void SetTle(const Tle& tle);
void FindPosition(Eci* eci, double tsince) const;
void FindPosition(Eci* eci, const Julian& date) const;
Eci FindPosition(double tsince) const;
Eci FindPosition(const Julian& date) const;
struct CommonConstants {
@ -196,9 +197,9 @@ private:
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;
void FindPositionSDP4(Eci* eci, double tsince) const;
void FindPositionSGP4(Eci* eci, double tsince) const;
void CalculateFinalPositionVelocity(Eci* eci, const double& tsince, const double& e,
Eci FindPositionSDP4(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,

3
SatTrack.cpp
View File

@ -10,11 +10,10 @@ int main() {
"1 25544U 98067A 11146.36888985 .00025753 00000-0 16912-3 0 4201",
"2 25544 51.6504 272.6534 0003891 329.5510 71.2188 15.75539412717473");
SGP4 sgp4(tle);
Eci eci;
while(1) {
Julian now;
sgp4.FindPosition(&eci, now);
Eci eci = sgp4.FindPosition(now);
CoordTopographic topo = obs.GetLookAngle(eci);
CoordGeodetic geo = eci.ToGeodetic();
std::cout << now << " ";