2013-03-11 18:29:33 +00:00
|
|
|
/*!
|
|
|
|
* \file gps_ephemeris.cc
|
2013-03-14 12:52:32 +00:00
|
|
|
* \brief Interface of a GPS EPHEMERIS storage and orbital model functions
|
2013-03-11 18:29:33 +00:00
|
|
|
*
|
2020-01-25 12:07:03 +00:00
|
|
|
* See https://www.gps.gov/technical/icwg/IS-GPS-200K.pdf Appendix II
|
2013-03-11 18:29:33 +00:00
|
|
|
* \author Javier Arribas, 2013. jarribas(at)cttc.es
|
|
|
|
*
|
|
|
|
* -------------------------------------------------------------------------
|
|
|
|
*
|
2019-07-26 10:38:20 +00:00
|
|
|
* Copyright (C) 2010-2019 (see AUTHORS file for a list of contributors)
|
2013-03-11 18:29:33 +00:00
|
|
|
*
|
|
|
|
* GNSS-SDR is a software defined Global Navigation
|
|
|
|
* Satellite Systems receiver
|
|
|
|
*
|
|
|
|
* This file is part of GNSS-SDR.
|
|
|
|
*
|
2020-02-08 00:20:02 +00:00
|
|
|
* SPDX-License-Identifier: GPL-3.0-or-later
|
2013-03-11 18:29:33 +00:00
|
|
|
*
|
|
|
|
* -------------------------------------------------------------------------
|
|
|
|
*/
|
|
|
|
|
|
|
|
#include "gps_ephemeris.h"
|
2018-02-26 02:15:53 +00:00
|
|
|
#include "GPS_L1_CA.h"
|
2018-12-02 13:32:22 +00:00
|
|
|
#include "gnss_satellite.h"
|
2018-02-26 02:15:53 +00:00
|
|
|
#include <cmath>
|
2013-03-11 18:29:33 +00:00
|
|
|
|
2020-06-10 08:15:59 +00:00
|
|
|
|
2013-03-11 18:29:33 +00:00
|
|
|
Gps_Ephemeris::Gps_Ephemeris()
|
|
|
|
{
|
2016-05-02 10:06:23 +00:00
|
|
|
auto gnss_sat = Gnss_Satellite();
|
2018-03-03 01:03:39 +00:00
|
|
|
std::string _system("GPS");
|
2018-08-13 23:13:07 +00:00
|
|
|
for (uint32_t i = 1; i < 33; i++)
|
2016-05-02 10:06:23 +00:00
|
|
|
{
|
|
|
|
satelliteBlock[i] = gnss_sat.what_block(_system, i);
|
|
|
|
}
|
2013-03-14 12:52:32 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
double Gps_Ephemeris::check_t(double time)
|
|
|
|
{
|
2013-07-04 13:47:40 +00:00
|
|
|
double corrTime;
|
2018-03-03 01:03:39 +00:00
|
|
|
double half_week = 302400.0; // seconds
|
2013-07-04 13:47:40 +00:00
|
|
|
corrTime = time;
|
|
|
|
if (time > half_week)
|
|
|
|
{
|
2015-05-13 14:40:46 +00:00
|
|
|
corrTime = time - 2.0 * half_week;
|
2013-07-04 13:47:40 +00:00
|
|
|
}
|
|
|
|
else if (time < -half_week)
|
|
|
|
{
|
2015-05-13 14:40:46 +00:00
|
|
|
corrTime = time + 2.0 * half_week;
|
2013-07-04 13:47:40 +00:00
|
|
|
}
|
|
|
|
return corrTime;
|
2013-03-14 12:52:32 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
// 20.3.3.3.3.1 User Algorithm for SV Clock Correction.
|
2013-06-12 15:19:32 +00:00
|
|
|
double Gps_Ephemeris::sv_clock_drift(double transmitTime)
|
2013-03-14 12:52:32 +00:00
|
|
|
{
|
2018-03-03 01:03:39 +00:00
|
|
|
// double dt;
|
|
|
|
// dt = check_t(transmitTime - d_Toc);
|
|
|
|
//
|
2018-08-13 23:13:07 +00:00
|
|
|
// for (int32_t i = 0; i < 2; i++)
|
2018-03-03 01:03:39 +00:00
|
|
|
// {
|
|
|
|
// dt -= d_A_f0 + d_A_f1 * dt + d_A_f2 * (dt * dt);
|
|
|
|
// }
|
|
|
|
// d_satClkDrift = d_A_f0 + d_A_f1 * dt + d_A_f2 * (dt * dt);
|
2017-04-12 15:04:51 +00:00
|
|
|
|
2013-07-04 13:47:40 +00:00
|
|
|
double dt;
|
|
|
|
dt = check_t(transmitTime - d_Toc);
|
2017-04-12 15:04:51 +00:00
|
|
|
d_satClkDrift = d_A_f0 + d_A_f1 * dt + d_A_f2 * (dt * dt) + sv_clock_relativistic_term(transmitTime);
|
2019-08-18 20:16:13 +00:00
|
|
|
// Correct satellite group delay
|
2018-03-03 01:03:39 +00:00
|
|
|
d_satClkDrift -= d_TGD;
|
2017-01-30 18:03:18 +00:00
|
|
|
|
2013-07-04 13:47:40 +00:00
|
|
|
return d_satClkDrift;
|
2013-03-11 18:29:33 +00:00
|
|
|
}
|
2013-07-04 13:47:40 +00:00
|
|
|
|
2013-11-23 13:05:38 +00:00
|
|
|
|
2013-06-12 15:19:32 +00:00
|
|
|
// compute the relativistic correction term
|
|
|
|
double Gps_Ephemeris::sv_clock_relativistic_term(double transmitTime)
|
|
|
|
{
|
2013-07-04 13:47:40 +00:00
|
|
|
double tk;
|
|
|
|
double a;
|
|
|
|
double n;
|
|
|
|
double n0;
|
|
|
|
double E;
|
|
|
|
double E_old;
|
|
|
|
double dE;
|
|
|
|
double M;
|
|
|
|
|
|
|
|
// Restore semi-major axis
|
2015-05-13 14:40:46 +00:00
|
|
|
a = d_sqrt_A * d_sqrt_A;
|
2013-07-04 13:47:40 +00:00
|
|
|
|
|
|
|
// Time from ephemeris reference epoch
|
|
|
|
tk = check_t(transmitTime - d_Toe);
|
|
|
|
|
|
|
|
// Computed mean motion
|
2020-07-05 18:20:02 +00:00
|
|
|
n0 = sqrt(GPS_GM / (a * a * a));
|
2013-07-04 13:47:40 +00:00
|
|
|
// Corrected mean motion
|
|
|
|
n = n0 + d_Delta_n;
|
|
|
|
// Mean anomaly
|
|
|
|
M = d_M_0 + n * tk;
|
|
|
|
|
|
|
|
// Reduce mean anomaly to between 0 and 2pi
|
2020-07-05 18:20:02 +00:00
|
|
|
// M = fmod((M + 2.0 * GNSS_PI), (2.0 * GNSS_PI));
|
2013-07-04 13:47:40 +00:00
|
|
|
|
|
|
|
// Initial guess of eccentric anomaly
|
|
|
|
E = M;
|
|
|
|
|
|
|
|
// --- Iteratively compute eccentric anomaly ----------------------------
|
2018-08-13 23:13:07 +00:00
|
|
|
for (int32_t ii = 1; ii < 20; ii++)
|
2013-07-04 13:47:40 +00:00
|
|
|
{
|
2018-03-03 01:03:39 +00:00
|
|
|
E_old = E;
|
|
|
|
E = M + d_e_eccentricity * sin(E);
|
2020-07-05 18:20:02 +00:00
|
|
|
dE = fmod(E - E_old, 2.0 * GNSS_PI);
|
2013-07-04 13:47:40 +00:00
|
|
|
if (fabs(dE) < 1e-12)
|
|
|
|
{
|
2019-08-18 20:16:13 +00:00
|
|
|
// Necessary precision is reached, exit from the loop
|
2013-07-04 13:47:40 +00:00
|
|
|
break;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
// Compute relativistic correction term
|
2020-07-05 18:20:02 +00:00
|
|
|
d_dtr = GPS_F * d_e_eccentricity * d_sqrt_A * sin(E);
|
2013-07-04 13:47:40 +00:00
|
|
|
return d_dtr;
|
2013-06-12 15:19:32 +00:00
|
|
|
}
|
2013-03-14 12:52:32 +00:00
|
|
|
|
2013-07-04 13:47:40 +00:00
|
|
|
|
2017-01-30 18:03:18 +00:00
|
|
|
double Gps_Ephemeris::satellitePosition(double transmitTime)
|
2013-03-14 12:52:32 +00:00
|
|
|
{
|
2013-07-04 13:47:40 +00:00
|
|
|
double tk;
|
|
|
|
double a;
|
|
|
|
double n;
|
|
|
|
double n0;
|
|
|
|
double M;
|
|
|
|
double E;
|
|
|
|
double E_old;
|
|
|
|
double dE;
|
|
|
|
double nu;
|
|
|
|
double phi;
|
|
|
|
double u;
|
|
|
|
double r;
|
|
|
|
double i;
|
|
|
|
double Omega;
|
|
|
|
|
|
|
|
// Find satellite's position ----------------------------------------------
|
|
|
|
|
|
|
|
// Restore semi-major axis
|
2017-01-30 19:26:50 +00:00
|
|
|
a = d_sqrt_A * d_sqrt_A;
|
2013-07-04 13:47:40 +00:00
|
|
|
|
|
|
|
// Time from ephemeris reference epoch
|
|
|
|
tk = check_t(transmitTime - d_Toe);
|
|
|
|
|
|
|
|
// Computed mean motion
|
2020-07-05 18:20:02 +00:00
|
|
|
n0 = sqrt(GPS_GM / (a * a * a));
|
2013-07-04 13:47:40 +00:00
|
|
|
|
|
|
|
// Corrected mean motion
|
|
|
|
n = n0 + d_Delta_n;
|
|
|
|
|
|
|
|
// Mean anomaly
|
|
|
|
M = d_M_0 + n * tk;
|
|
|
|
|
|
|
|
// Reduce mean anomaly to between 0 and 2pi
|
2020-07-05 18:20:02 +00:00
|
|
|
// M = fmod((M + 2.0 * GNSS_PI), (2.0 * GNSS_PI));
|
2013-07-04 13:47:40 +00:00
|
|
|
|
|
|
|
// Initial guess of eccentric anomaly
|
|
|
|
E = M;
|
|
|
|
|
|
|
|
// --- Iteratively compute eccentric anomaly ----------------------------
|
2018-08-13 23:13:07 +00:00
|
|
|
for (int32_t ii = 1; ii < 20; ii++)
|
2013-07-04 13:47:40 +00:00
|
|
|
{
|
2018-03-03 01:03:39 +00:00
|
|
|
E_old = E;
|
|
|
|
E = M + d_e_eccentricity * sin(E);
|
2020-07-05 18:20:02 +00:00
|
|
|
dE = fmod(E - E_old, 2.0 * GNSS_PI);
|
2013-07-04 13:47:40 +00:00
|
|
|
if (fabs(dE) < 1e-12)
|
|
|
|
{
|
2019-08-18 20:16:13 +00:00
|
|
|
// Necessary precision is reached, exit from the loop
|
2013-07-04 13:47:40 +00:00
|
|
|
break;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
// Compute the true anomaly
|
|
|
|
double tmp_Y = sqrt(1.0 - d_e_eccentricity * d_e_eccentricity) * sin(E);
|
|
|
|
double tmp_X = cos(E) - d_e_eccentricity;
|
|
|
|
nu = atan2(tmp_Y, tmp_X);
|
|
|
|
|
|
|
|
// Compute angle phi (argument of Latitude)
|
|
|
|
phi = nu + d_OMEGA;
|
|
|
|
|
|
|
|
// Reduce phi to between 0 and 2*pi rad
|
2020-07-05 18:20:02 +00:00
|
|
|
// phi = fmod((phi), (2.0 * GNSS_PI));
|
2013-07-04 13:47:40 +00:00
|
|
|
|
|
|
|
// Correct argument of latitude
|
2018-03-03 01:03:39 +00:00
|
|
|
u = phi + d_Cuc * cos(2.0 * phi) + d_Cus * sin(2.0 * phi);
|
2013-07-04 13:47:40 +00:00
|
|
|
|
|
|
|
// Correct radius
|
2018-03-03 01:03:39 +00:00
|
|
|
r = a * (1.0 - d_e_eccentricity * cos(E)) + d_Crc * cos(2.0 * phi) + d_Crs * sin(2.0 * phi);
|
2013-07-04 13:47:40 +00:00
|
|
|
|
|
|
|
// Correct inclination
|
2017-01-30 19:26:50 +00:00
|
|
|
i = d_i_0 + d_IDOT * tk + d_Cic * cos(2.0 * phi) + d_Cis * sin(2.0 * phi);
|
2013-07-04 13:47:40 +00:00
|
|
|
|
|
|
|
// Compute the angle between the ascending node and the Greenwich meridian
|
2020-07-05 18:20:02 +00:00
|
|
|
Omega = d_OMEGA0 + (d_OMEGA_DOT - GNSS_OMEGA_EARTH_DOT) * tk - GNSS_OMEGA_EARTH_DOT * d_Toe;
|
2013-07-04 13:47:40 +00:00
|
|
|
|
|
|
|
// Reduce to between 0 and 2*pi rad
|
2020-07-05 18:20:02 +00:00
|
|
|
// Omega = fmod((Omega + 2.0 * GNSS_PI), (2.0 * GNSS_PI));
|
2013-07-04 13:47:40 +00:00
|
|
|
|
|
|
|
// --- Compute satellite coordinates in Earth-fixed coordinates
|
|
|
|
d_satpos_X = cos(u) * r * cos(Omega) - sin(u) * r * cos(i) * sin(Omega);
|
|
|
|
d_satpos_Y = cos(u) * r * sin(Omega) + sin(u) * r * cos(i) * cos(Omega);
|
|
|
|
d_satpos_Z = sin(u) * r * sin(i);
|
|
|
|
|
|
|
|
// Satellite's velocity. Can be useful for Vector Tracking loops
|
2020-07-05 18:20:02 +00:00
|
|
|
double Omega_dot = d_OMEGA_DOT - GNSS_OMEGA_EARTH_DOT;
|
2018-03-03 01:03:39 +00:00
|
|
|
d_satvel_X = -Omega_dot * (cos(u) * r + sin(u) * r * cos(i)) + d_satpos_X * cos(Omega) - d_satpos_Y * cos(i) * sin(Omega);
|
2013-07-04 13:47:40 +00:00
|
|
|
d_satvel_Y = Omega_dot * (cos(u) * r * cos(Omega) - sin(u) * r * cos(i) * sin(Omega)) + d_satpos_X * sin(Omega) + d_satpos_Y * cos(i) * cos(Omega);
|
|
|
|
d_satvel_Z = d_satpos_Y * sin(i);
|
2017-01-30 18:03:18 +00:00
|
|
|
|
|
|
|
// Time from ephemeris reference clock
|
|
|
|
tk = check_t(transmitTime - d_Toc);
|
|
|
|
|
2017-01-30 19:26:50 +00:00
|
|
|
double dtr_s = d_A_f0 + d_A_f1 * tk + d_A_f2 * tk * tk;
|
2017-01-30 18:03:18 +00:00
|
|
|
|
|
|
|
/* relativity correction */
|
2020-07-05 18:20:02 +00:00
|
|
|
dtr_s -= 2.0 * sqrt(GPS_GM * a) * d_e_eccentricity * sin(E) / (SPEED_OF_LIGHT_M_S * SPEED_OF_LIGHT_M_S);
|
2017-01-30 18:03:18 +00:00
|
|
|
|
|
|
|
return dtr_s;
|
2013-03-14 12:52:32 +00:00
|
|
|
}
|