mirror of
https://github.com/gnss-sdr/gnss-sdr
synced 2026-10-04 06:41:49 +00:00
Merge branch 'next' of https://github.com/gnss-sdr/gnss-sdr into fix_FPGA_IP_address_ranges
This commit is contained in:
Marc Majoral
@@ -19,6 +19,7 @@
|
||||
|
||||
#include "gnss_sdr_supl_client.h"
|
||||
#include "GPS_L1_CA.h"
|
||||
#include "MATH_CONSTANTS.h"
|
||||
#include <boost/archive/xml_iarchive.hpp>
|
||||
#include <boost/archive/xml_oarchive.hpp>
|
||||
#include <boost/serialization/map.hpp>
|
||||
@@ -833,7 +834,7 @@ bool Gnss_Sdr_Supl_Client::load_gal_almanac_xml(const std::string& file_name)
|
||||
}
|
||||
catch (std::exception& e)
|
||||
{
|
||||
// Maybe the file is from https://www.gsc-europa.eu/system-status/almanac-data ?
|
||||
// Maybe the file is from https://www.gsc-europa.eu/product-almanacs ?
|
||||
return this->read_gal_almanac_from_gsa(file_name);
|
||||
}
|
||||
LOG(INFO) << "Loaded Galileo almanac map data with " << this->gal_almanac_map.size() << " satellites";
|
||||
@@ -858,6 +859,7 @@ bool Gnss_Sdr_Supl_Client::read_gal_almanac_from_gsa(const std::string& file_nam
|
||||
Galileo_Almanac gal_alm;
|
||||
try
|
||||
{
|
||||
const double sqrtAnominal = 5440.588203494; // square root of Galileo nominal orbit semi-major axis
|
||||
uint32_t prn = static_cast<uint32_t>(std::stoi(almanac.child_value("SVID")));
|
||||
gal_alm.PRN = prn;
|
||||
gal_alm.toa = std::stoi(almanac.child("almanac").child_value("t0a"));
|
||||
@@ -866,7 +868,7 @@ bool Gnss_Sdr_Supl_Client::read_gal_almanac_from_gsa(const std::string& file_nam
|
||||
gal_alm.delta_i = std::stod(almanac.child("almanac").child_value("deltai"));
|
||||
gal_alm.M_0 = std::stod(almanac.child("almanac").child_value("m0"));
|
||||
gal_alm.ecc = std::stod(almanac.child("almanac").child_value("ecc"));
|
||||
gal_alm.sqrtA = std::stod(almanac.child("almanac").child_value("aSqRoot"));
|
||||
gal_alm.sqrtA = std::stod(almanac.child("almanac").child_value("aSqRoot")) + sqrtAnominal;
|
||||
gal_alm.OMEGA_0 = std::stod(almanac.child("almanac").child_value("omega0"));
|
||||
gal_alm.omega = std::stod(almanac.child("almanac").child_value("w"));
|
||||
gal_alm.OMEGAdot = std::stod(almanac.child("almanac").child_value("omegaDot"));
|
||||
|
||||
@@ -6,6 +6,7 @@
|
||||
|
||||
|
||||
set(SYSTEM_PARAMETERS_SOURCES
|
||||
gnss_almanac.cc
|
||||
gnss_ephemeris.cc
|
||||
gnss_satellite.cc
|
||||
gnss_signal.cc
|
||||
|
||||
@@ -36,7 +36,10 @@ public:
|
||||
/*!
|
||||
* Default constructor
|
||||
*/
|
||||
Beidou_Dnav_Almanac() = default;
|
||||
Beidou_Dnav_Almanac()
|
||||
{
|
||||
this->System = 'B';
|
||||
};
|
||||
|
||||
int SV_health{}; //!< SV Health
|
||||
|
||||
|
||||
@@ -36,7 +36,10 @@ public:
|
||||
/*!
|
||||
* Default constructor
|
||||
*/
|
||||
Galileo_Almanac() = default;
|
||||
Galileo_Almanac()
|
||||
{
|
||||
this->System = 'E';
|
||||
};
|
||||
|
||||
int32_t IODa{};
|
||||
int32_t E5b_HS{};
|
||||
|
||||
@@ -0,0 +1,275 @@
|
||||
/*!
|
||||
* \file gnss_almanac.cc
|
||||
* \brief Base class for GNSS almanac storage
|
||||
* \author Carles Fernandez, 2022. cfernandez(at)cttc.es
|
||||
*
|
||||
* -----------------------------------------------------------------------------
|
||||
*
|
||||
* GNSS-SDR is a Global Navigation Satellite System software-defined receiver.
|
||||
* This file is part of GNSS-SDR.
|
||||
*
|
||||
* Copyright (C) 2010-2022 (see AUTHORS file for a list of contributors)
|
||||
* SPDX-License-Identifier: GPL-3.0-or-later
|
||||
*
|
||||
* -----------------------------------------------------------------------------
|
||||
*/
|
||||
|
||||
#include "gnss_almanac.h"
|
||||
#include "MATH_CONSTANTS.h"
|
||||
#include "gnss_frequencies.h"
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <functional>
|
||||
#include <numeric>
|
||||
#include <vector>
|
||||
|
||||
double Gnss_Almanac::check_t(double time) const
|
||||
{
|
||||
const double half_week = 302400.0; // seconds
|
||||
double corrTime = time;
|
||||
if (time > half_week)
|
||||
{
|
||||
corrTime = time - 2.0 * half_week;
|
||||
}
|
||||
else if (time < -half_week)
|
||||
{
|
||||
corrTime = time + 2.0 * half_week;
|
||||
}
|
||||
return corrTime;
|
||||
}
|
||||
|
||||
double Gnss_Almanac::predicted_doppler(double rx_time_s,
|
||||
double lat,
|
||||
double lon,
|
||||
double h,
|
||||
double ve,
|
||||
double vn,
|
||||
double vu,
|
||||
int band) const
|
||||
{
|
||||
const double RE_WGS84 = 6378137.0; //!< earth semimajor axis (WGS84) (m)
|
||||
const double FE_WGS84 = (1.0 / 298.257223563); //!< earth flattening (WGS84)
|
||||
const double lat_rad = lat * D2R;
|
||||
const double lon_rad = lon * D2R;
|
||||
|
||||
const double sinp = sin(lat_rad);
|
||||
const double cosp = cos(lat_rad);
|
||||
const double sinl = sin(lon_rad);
|
||||
const double cosl = cos(lon_rad);
|
||||
|
||||
const double e2 = FE_WGS84 * (2.0 - FE_WGS84);
|
||||
const double v = RE_WGS84 / std::sqrt(1.0 - e2 * sinp * sinp);
|
||||
|
||||
// Position in EFEF
|
||||
const std::vector<double> pos_rx = {(v + h) * cosp * cosl, (v + h) * cosp * sinl, (v * (1.0 - e2) + h) * sinp};
|
||||
|
||||
// Velocity in EFEF
|
||||
const double t = cosp * vu - sinp * vn;
|
||||
const std::vector<double> vel_rx = {cosl * t - sinl * ve, sinl * t + cosl * ve, sinp * vu + cosp * vn};
|
||||
|
||||
std::array<double, 7> sat_pos_vel = {0};
|
||||
satellitePosVelComputation(rx_time_s, sat_pos_vel);
|
||||
const std::vector<double> pos_sat = {sat_pos_vel[0], sat_pos_vel[1], sat_pos_vel[2]};
|
||||
const std::vector<double> vel_sat = {sat_pos_vel[3], sat_pos_vel[4], sat_pos_vel[5]};
|
||||
|
||||
std::vector<double> x_sr = pos_sat;
|
||||
std::transform(x_sr.begin(), x_sr.end(), pos_rx.begin(), x_sr.begin(), std::minus<double>()); // pos_sat - pos_rx
|
||||
|
||||
const double norm_x_sr = std::sqrt(std::inner_product(x_sr.begin(), x_sr.end(), x_sr.begin(), 0.0)); // Euclidean norm
|
||||
|
||||
std::vector<double> v_sr = vel_sat;
|
||||
std::transform(v_sr.begin(), v_sr.end(), vel_rx.begin(), v_sr.begin(), std::minus<double>()); // vel_sat - vel_rx
|
||||
|
||||
const double radial_vel = std::inner_product(v_sr.begin(), v_sr.end(), x_sr.begin(), 0.0) / norm_x_sr;
|
||||
const double predicted_doppler_normalized = -(radial_vel / SPEED_OF_LIGHT_M_S);
|
||||
double predicted_doppler = 0.0;
|
||||
if (this->System == 'E') // Galileo
|
||||
{
|
||||
if (band == 1)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ1;
|
||||
}
|
||||
else if (band == 5)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ5;
|
||||
}
|
||||
else if (band == 6)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ6;
|
||||
}
|
||||
else if (band == 7)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ7;
|
||||
}
|
||||
else if (band == 8)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ8;
|
||||
}
|
||||
else
|
||||
{
|
||||
predicted_doppler = 0.0;
|
||||
}
|
||||
}
|
||||
else if (this->System == 'G') // GPS
|
||||
{
|
||||
if (band == 1)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ1;
|
||||
}
|
||||
else if (band == 2)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ2;
|
||||
}
|
||||
else if (band == 5)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ5;
|
||||
}
|
||||
else
|
||||
{
|
||||
predicted_doppler = 0.0;
|
||||
}
|
||||
}
|
||||
else if (this->System == 'B') // Beidou
|
||||
{
|
||||
if (band == 1)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ1_BDS;
|
||||
}
|
||||
else if (band == 2)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ2_BDS;
|
||||
}
|
||||
else if (band == 3)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ3_BDS;
|
||||
}
|
||||
else
|
||||
{
|
||||
predicted_doppler = 0.0;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
predicted_doppler = 0.0;
|
||||
}
|
||||
return predicted_doppler;
|
||||
}
|
||||
|
||||
|
||||
void Gnss_Almanac::satellitePosVelComputation(double transmitTime, std::array<double, 7>& pos_vel_dtr) const
|
||||
{
|
||||
// Restore semi-major axis
|
||||
const double a = this->sqrtA * this->sqrtA;
|
||||
|
||||
// Computed mean motion
|
||||
double n;
|
||||
if (this->System == 'E')
|
||||
{
|
||||
n = sqrt(GALILEO_GM / (a * a * a));
|
||||
}
|
||||
else if (this->System == 'B')
|
||||
{
|
||||
n = sqrt(BEIDOU_GM / (a * a * a));
|
||||
}
|
||||
else
|
||||
{
|
||||
n = sqrt(GPS_GM / (a * a * a));
|
||||
}
|
||||
|
||||
// Time from ephemeris reference epoch
|
||||
const double tk = check_t(transmitTime - static_cast<double>(this->toa));
|
||||
|
||||
// Mean anomaly
|
||||
const double M = this->M_0 * GNSS_PI + n * tk;
|
||||
|
||||
// Initial guess of eccentric anomaly
|
||||
double E = M;
|
||||
double E_old;
|
||||
double dE;
|
||||
|
||||
// --- Iteratively compute eccentric anomaly -------------------------------
|
||||
for (int32_t ii = 1; ii < 20; ii++)
|
||||
{
|
||||
E_old = E;
|
||||
E = M + this->ecc * sin(E);
|
||||
dE = fmod(E - E_old, 2.0 * GNSS_PI);
|
||||
if (fabs(dE) < 1e-12)
|
||||
{
|
||||
// Necessary precision is reached, exit from the loop
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
const double sek = sin(E);
|
||||
const double cek = cos(E);
|
||||
const double OneMinusecosE = 1.0 - this->ecc * cek;
|
||||
const double sq1e2 = sqrt(1.0 - this->ecc * this->ecc);
|
||||
const double ekdot = n / OneMinusecosE;
|
||||
|
||||
// Compute the true anomaly
|
||||
const double tmp_Y = sq1e2 * sek;
|
||||
const double tmp_X = cek - this->ecc;
|
||||
const double nu = atan2(tmp_Y, tmp_X);
|
||||
|
||||
// Compute angle phi (argument of Latitude)
|
||||
const double phi = nu + this->omega * GNSS_PI;
|
||||
|
||||
const double pkdot = sq1e2 * ekdot / OneMinusecosE;
|
||||
|
||||
// Correct argument of latitude
|
||||
const double suk = sin(phi);
|
||||
const double cuk = cos(phi);
|
||||
|
||||
// Correct radius
|
||||
const double r = a * OneMinusecosE;
|
||||
const double rkdot = a * this->ecc * sek * ekdot;
|
||||
|
||||
// Correct inclination
|
||||
double i;
|
||||
if (this->System == 'E')
|
||||
{
|
||||
i = ((56.0 / 180.0) + this->delta_i) * GNSS_PI;
|
||||
}
|
||||
else
|
||||
{
|
||||
i = (0.3 + this->delta_i) * GNSS_PI;
|
||||
}
|
||||
|
||||
const double sik = sin(i);
|
||||
const double cik = cos(i);
|
||||
|
||||
// Compute the angle between the ascending node and the Greenwich meridian
|
||||
double Omega;
|
||||
double Omega_dot;
|
||||
if (this->System == 'B')
|
||||
{
|
||||
Omega_dot = this->OMEGAdot * GNSS_PI - BEIDOU_OMEGA_EARTH_DOT;
|
||||
Omega = this->OMEGA_0 * GNSS_PI + Omega_dot * tk - BEIDOU_OMEGA_EARTH_DOT * static_cast<double>(this->toa);
|
||||
}
|
||||
else
|
||||
{
|
||||
Omega_dot = this->OMEGAdot * GNSS_PI - GNSS_OMEGA_EARTH_DOT;
|
||||
Omega = this->OMEGA_0 * GNSS_PI + Omega_dot * tk - GNSS_OMEGA_EARTH_DOT * static_cast<double>(this->toa);
|
||||
}
|
||||
|
||||
const double sok = sin(Omega);
|
||||
const double cok = cos(Omega);
|
||||
|
||||
// --- Compute satellite coordinates in Earth-fixed coordinates
|
||||
const double xprime = r * cuk;
|
||||
const double yprime = r * suk;
|
||||
|
||||
pos_vel_dtr[0] = xprime * cok - yprime * cik * sok;
|
||||
pos_vel_dtr[1] = xprime * sok + yprime * cik * cok;
|
||||
pos_vel_dtr[2] = yprime * sik;
|
||||
|
||||
// Satellite's velocity
|
||||
const double xpkdot = rkdot * cuk - yprime * pkdot;
|
||||
const double ypkdot = rkdot * suk + xprime * pkdot;
|
||||
const double tmp = ypkdot * cik;
|
||||
|
||||
pos_vel_dtr[3] = -Omega_dot * pos_vel_dtr[1] + xpkdot * cok - tmp * sok;
|
||||
pos_vel_dtr[4] = Omega_dot * pos_vel_dtr[0] + xpkdot * sok + tmp * cok;
|
||||
pos_vel_dtr[5] = ypkdot * sik;
|
||||
pos_vel_dtr[6] = 0;
|
||||
}
|
||||
@@ -18,6 +18,7 @@
|
||||
#ifndef GNSS_SDR_GNSS_ALMANAC_H
|
||||
#define GNSS_SDR_GNSS_ALMANAC_H
|
||||
|
||||
#include <array>
|
||||
#include <cstdint>
|
||||
|
||||
/** \addtogroup Core
|
||||
@@ -37,6 +38,46 @@ public:
|
||||
*/
|
||||
Gnss_Almanac() = default;
|
||||
|
||||
/*!
|
||||
* \brief Computes prediction of the Doppler shift for a given time and receiver's position and velocity.
|
||||
* \f[
|
||||
* f_{d} = - \mathbf{v} \frac{\mathbf{x}^{T}}{\left| \mathbf{x} \right| } \frac{f_{L}}{c}
|
||||
* \f]
|
||||
* where:
|
||||
* \f[
|
||||
* \mathbf{v} = \mathbf{v}_{sat} - \mathbf{v}_{rx}
|
||||
* \f]
|
||||
* \f[
|
||||
* \mathbf{x} = \mathbf{x}_{sat} - \mathbf{x}_{rx}
|
||||
* \f]
|
||||
* \f[
|
||||
* \left| \mathbf{x} \right| = \sqrt{\mathbf{x}\mathbf{x}^{T}}
|
||||
* \f]
|
||||
*
|
||||
* @param[in] rx_time_s Time of Week in seconds
|
||||
* @param[in] lat Receiver's latitude in degrees
|
||||
* @param[in] lon Receiver's longitude in degrees
|
||||
* @param[in] h Receiver's height in meters
|
||||
* @param[in] ve Receiver's velocity in the East direction [m/s]
|
||||
* @param[in] vn Receiver's velocity in the North direction [m/s]
|
||||
* @param[in] vu Receiver's velocity in the Up direction [m/s]
|
||||
* @param[in] band Signal band for which the Doppler will be computed
|
||||
* (1: L1 C/A, E1B, BI1; 2: L2C, BI2; 3: BI3; 5: L5/E5a; 6: E6B; 7: E5b; 8: E5a+E5b)
|
||||
*/
|
||||
double predicted_doppler(double rx_time_s,
|
||||
double lat,
|
||||
double lon,
|
||||
double h,
|
||||
double ve,
|
||||
double vn,
|
||||
double vu,
|
||||
int band) const;
|
||||
|
||||
/*!
|
||||
* \brief Computes satellite Position and Velocity, in ECEF, for a given time (expressed in seconds of week)
|
||||
*/
|
||||
void satellitePosVelComputation(double transmitTime, std::array<double, 7>& pos_vel_dtr) const;
|
||||
|
||||
uint32_t PRN{}; //!< SV PRN NUMBER
|
||||
double delta_i{}; //!< Inclination Angle at Reference Time (relative to i_0 = 0.30 semi-circles)
|
||||
int32_t toa{}; //!< Almanac data reference time of week [s]
|
||||
@@ -49,6 +90,11 @@ public:
|
||||
double OMEGAdot{}; //!< Rate of Right Ascension [semi-circles/s]
|
||||
double af0{}; //!< Coefficient 0 of code phase offset model [s]
|
||||
double af1{}; //!< Coefficient 1 of code phase offset model [s/s]
|
||||
|
||||
protected:
|
||||
char System{}; //!< Character ID of the GNSS system. 'G': GPS. 'E': Galileo. 'B': BeiDou
|
||||
private:
|
||||
double check_t(double time) const;
|
||||
};
|
||||
|
||||
|
||||
|
||||
@@ -17,10 +17,156 @@
|
||||
|
||||
#include "gnss_ephemeris.h"
|
||||
#include "MATH_CONSTANTS.h"
|
||||
#include "gnss_frequencies.h"
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <functional>
|
||||
#include <numeric>
|
||||
#include <vector>
|
||||
|
||||
|
||||
double Gnss_Ephemeris::sv_clock_drift(double transmitTime)
|
||||
{
|
||||
const double dt = check_t(transmitTime - this->toc);
|
||||
this->dtr = sv_clock_relativistic_term(transmitTime);
|
||||
this->satClkDrift = this->af0 + this->af1 * dt + this->af2 * (dt * dt) + this->dtr;
|
||||
return this->satClkDrift;
|
||||
}
|
||||
|
||||
|
||||
double Gnss_Ephemeris::predicted_doppler(double rx_time_s,
|
||||
double lat,
|
||||
double lon,
|
||||
double h,
|
||||
double ve,
|
||||
double vn,
|
||||
double vu,
|
||||
int band) const
|
||||
{
|
||||
const double RE_WGS84 = 6378137.0; //!< earth semimajor axis (WGS84) (m)
|
||||
const double FE_WGS84 = (1.0 / 298.257223563); //!< earth flattening (WGS84)
|
||||
const double lat_rad = lat * D2R;
|
||||
const double lon_rad = lon * D2R;
|
||||
|
||||
const double sinp = sin(lat_rad);
|
||||
const double cosp = cos(lat_rad);
|
||||
const double sinl = sin(lon_rad);
|
||||
const double cosl = cos(lon_rad);
|
||||
|
||||
const double e2 = FE_WGS84 * (2.0 - FE_WGS84);
|
||||
const double v = RE_WGS84 / std::sqrt(1.0 - e2 * sinp * sinp);
|
||||
|
||||
// Position in EFEF
|
||||
const std::vector<double> pos_rx = {(v + h) * cosp * cosl, (v + h) * cosp * sinl, (v * (1.0 - e2) + h) * sinp};
|
||||
|
||||
// Velocity in EFEF
|
||||
const double t = cosp * vu - sinp * vn;
|
||||
const std::vector<double> vel_rx = {cosl * t - sinl * ve, sinl * t + cosl * ve, sinp * vu + cosp * vn};
|
||||
|
||||
std::array<double, 7> sat_pos_vel = {0};
|
||||
satellitePosVelComputation(rx_time_s, sat_pos_vel);
|
||||
const std::vector<double> pos_sat = {sat_pos_vel[0], sat_pos_vel[1], sat_pos_vel[2]};
|
||||
const std::vector<double> vel_sat = {sat_pos_vel[3], sat_pos_vel[4], sat_pos_vel[5]};
|
||||
|
||||
std::vector<double> x_sr = pos_sat;
|
||||
std::transform(x_sr.begin(), x_sr.end(), pos_rx.begin(), x_sr.begin(), std::minus<double>()); // pos_sat - pos_rx
|
||||
|
||||
const double norm_x_sr = std::sqrt(std::inner_product(x_sr.begin(), x_sr.end(), x_sr.begin(), 0.0)); // Euclidean norm
|
||||
|
||||
std::vector<double> v_sr = vel_sat;
|
||||
std::transform(v_sr.begin(), v_sr.end(), vel_rx.begin(), v_sr.begin(), std::minus<double>()); // vel_sat - vel_rx
|
||||
|
||||
const double radial_vel = std::inner_product(v_sr.begin(), v_sr.end(), x_sr.begin(), 0.0) / norm_x_sr;
|
||||
const double predicted_doppler_normalized = -(radial_vel / SPEED_OF_LIGHT_M_S);
|
||||
double predicted_doppler = 0.0;
|
||||
if (this->System == 'E') // Galileo
|
||||
{
|
||||
if (band == 1)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ1;
|
||||
}
|
||||
else if (band == 5)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ5;
|
||||
}
|
||||
else if (band == 6)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ6;
|
||||
}
|
||||
else if (band == 7)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ7;
|
||||
}
|
||||
else if (band == 8)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ8;
|
||||
}
|
||||
else
|
||||
{
|
||||
predicted_doppler = 0.0;
|
||||
}
|
||||
}
|
||||
else if (this->System == 'G') // GPS
|
||||
{
|
||||
if (band == 1)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ1;
|
||||
}
|
||||
else if (band == 2)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ2;
|
||||
}
|
||||
else if (band == 5)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ5;
|
||||
}
|
||||
else
|
||||
{
|
||||
predicted_doppler = 0.0;
|
||||
}
|
||||
}
|
||||
else if (this->System == 'B') // Beidou
|
||||
{
|
||||
if (band == 1)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ1_BDS;
|
||||
}
|
||||
else if (band == 2)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ2_BDS;
|
||||
}
|
||||
else if (band == 3)
|
||||
{
|
||||
predicted_doppler = predicted_doppler_normalized * FREQ3_BDS;
|
||||
}
|
||||
else
|
||||
{
|
||||
predicted_doppler = 0.0;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
predicted_doppler = 0.0;
|
||||
}
|
||||
return predicted_doppler;
|
||||
}
|
||||
|
||||
|
||||
void Gnss_Ephemeris::satellitePosition(double transmitTime)
|
||||
{
|
||||
std::array<double, 7> pos_vel_dtr = {0};
|
||||
satellitePosVelComputation(transmitTime, pos_vel_dtr);
|
||||
this->satpos_X = pos_vel_dtr[0];
|
||||
this->satpos_Y = pos_vel_dtr[1];
|
||||
this->satpos_Z = pos_vel_dtr[2];
|
||||
this->satvel_X = pos_vel_dtr[3];
|
||||
this->satvel_Y = pos_vel_dtr[4];
|
||||
this->satvel_Z = pos_vel_dtr[5];
|
||||
this->dtr = pos_vel_dtr[6];
|
||||
}
|
||||
|
||||
|
||||
void Gnss_Ephemeris::satellitePosVelComputation(double transmitTime, std::array<double, 7>& pos_vel_dtr) const
|
||||
{
|
||||
// Restore semi-major axis
|
||||
const double a = this->sqrtA * this->sqrtA;
|
||||
@@ -47,10 +193,7 @@ void Gnss_Ephemeris::satellitePosition(double transmitTime)
|
||||
const double n = n0 + this->delta_n;
|
||||
|
||||
// Mean anomaly
|
||||
double M = this->M_0 + n * tk;
|
||||
|
||||
// Reduce mean anomaly to between 0 and 2pi
|
||||
M = fmod((M + 2.0 * GNSS_PI), (2.0 * GNSS_PI));
|
||||
const double M = this->M_0 + n * tk;
|
||||
|
||||
// Initial guess of eccentric anomaly
|
||||
double E = M;
|
||||
@@ -82,10 +225,9 @@ void Gnss_Ephemeris::satellitePosition(double transmitTime)
|
||||
const double nu = atan2(tmp_Y, tmp_X);
|
||||
|
||||
// Compute angle phi (argument of Latitude)
|
||||
double phi = nu + this->omega;
|
||||
const double phi = nu + this->omega;
|
||||
|
||||
// Reduce phi to between 0 and 2*pi rad
|
||||
phi = fmod((phi), (2.0 * GNSS_PI));
|
||||
const double s2pk = sin(2.0 * phi);
|
||||
const double c2pk = cos(2.0 * phi);
|
||||
const double pkdot = sq1e2 * ekdot / OneMinusecosE;
|
||||
@@ -120,53 +262,59 @@ void Gnss_Ephemeris::satellitePosition(double transmitTime)
|
||||
Omega = this->OMEGA_0 + Omega_dot * tk - GNSS_OMEGA_EARTH_DOT * static_cast<double>(this->toe);
|
||||
}
|
||||
|
||||
// Reduce to between 0 and 2*pi rad
|
||||
Omega = fmod((Omega + 2.0 * GNSS_PI), (2.0 * GNSS_PI));
|
||||
const double sok = sin(Omega);
|
||||
const double cok = cos(Omega);
|
||||
|
||||
// --- Compute satellite coordinates in Earth-fixed coordinates
|
||||
const double xprime = r * cuk;
|
||||
const double yprime = r * suk;
|
||||
this->satpos_X = xprime * cok - yprime * cik * sok;
|
||||
this->satpos_Y = xprime * sok + yprime * cik * cok; // ********NOTE: in GALILEO ICD this expression is not correct because it has minus (- sin(u) * r * cos(i) * cos(Omega)) instead of plus
|
||||
this->satpos_Z = yprime * sik;
|
||||
|
||||
pos_vel_dtr[0] = xprime * cok - yprime * cik * sok;
|
||||
pos_vel_dtr[1] = xprime * sok + yprime * cik * cok; // ********NOTE: in GALILEO ICD this expression is not correct because it has minus (- sin(u) * r * cos(i) * cos(Omega)) instead of plus
|
||||
pos_vel_dtr[2] = yprime * sik;
|
||||
|
||||
// Satellite's velocity. Can be useful for Vector Tracking loops
|
||||
const double xpkdot = rkdot * cuk - yprime * ukdot;
|
||||
const double ypkdot = rkdot * suk + xprime * ukdot;
|
||||
const double tmp = ypkdot * cik - this->satpos_Z * ikdot;
|
||||
const double tmp = ypkdot * cik - pos_vel_dtr[2] * ikdot;
|
||||
|
||||
this->satvel_X = -Omega_dot * this->satpos_Y + xpkdot * cok - tmp * sok;
|
||||
this->satvel_Y = Omega_dot * this->satpos_X + xpkdot * sok + tmp * cok;
|
||||
this->satvel_Z = yprime * cik * ikdot + ypkdot * sik;
|
||||
pos_vel_dtr[3] = -Omega_dot * pos_vel_dtr[1] + xpkdot * cok - tmp * sok;
|
||||
pos_vel_dtr[4] = Omega_dot * pos_vel_dtr[0] + xpkdot * sok + tmp * cok;
|
||||
pos_vel_dtr[5] = yprime * cik * ikdot + ypkdot * sik;
|
||||
|
||||
// Time from ephemeris reference clock
|
||||
tk = check_t(transmitTime - this->toc);
|
||||
|
||||
this->dtr = this->af0 + this->af1 * tk + this->af2 * tk * tk;
|
||||
pos_vel_dtr[6] = this->af0 + this->af1 * tk + this->af2 * tk * tk;
|
||||
|
||||
if (this->System == 'E')
|
||||
{
|
||||
this->dtr -= 2.0 * sqrt(GALILEO_GM * a) * this->ecc * sek / (SPEED_OF_LIGHT_M_S * SPEED_OF_LIGHT_M_S);
|
||||
pos_vel_dtr[6] -= 2.0 * sqrt(GALILEO_GM * a) * this->ecc * sek / (SPEED_OF_LIGHT_M_S * SPEED_OF_LIGHT_M_S);
|
||||
}
|
||||
else if (this->System == 'B')
|
||||
{
|
||||
this->dtr -= 2.0 * sqrt(BEIDOU_GM * a) * this->ecc * sek / (SPEED_OF_LIGHT_M_S * SPEED_OF_LIGHT_M_S);
|
||||
pos_vel_dtr[6] -= 2.0 * sqrt(BEIDOU_GM * a) * this->ecc * sek / (SPEED_OF_LIGHT_M_S * SPEED_OF_LIGHT_M_S);
|
||||
}
|
||||
else
|
||||
{
|
||||
this->dtr -= 2.0 * sqrt(GPS_GM * a) * this->ecc * sek / (SPEED_OF_LIGHT_M_S * SPEED_OF_LIGHT_M_S);
|
||||
pos_vel_dtr[6] -= 2.0 * sqrt(GPS_GM * a) * this->ecc * sek / (SPEED_OF_LIGHT_M_S * SPEED_OF_LIGHT_M_S);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
double Gnss_Ephemeris::sv_clock_drift(double transmitTime)
|
||||
double Gnss_Ephemeris::check_t(double time) const
|
||||
{
|
||||
const double dt = check_t(transmitTime - this->toc);
|
||||
this->dtr = sv_clock_relativistic_term(transmitTime);
|
||||
this->satClkDrift = this->af0 + this->af1 * dt + this->af2 * (dt * dt) + this->dtr;
|
||||
return this->satClkDrift;
|
||||
const double half_week = 302400.0; // seconds
|
||||
double corrTime = time;
|
||||
if (time > half_week)
|
||||
{
|
||||
corrTime = time - 2.0 * half_week;
|
||||
}
|
||||
else if (time < -half_week)
|
||||
{
|
||||
corrTime = time + 2.0 * half_week;
|
||||
}
|
||||
return corrTime;
|
||||
}
|
||||
|
||||
|
||||
@@ -184,7 +332,7 @@ double Gnss_Ephemeris::sv_clock_relativistic_term(double transmitTime) const
|
||||
{
|
||||
n0 = sqrt(GALILEO_GM / (a * a * a));
|
||||
}
|
||||
else if (this->System == 'E')
|
||||
else if (this->System == 'B')
|
||||
{
|
||||
n0 = sqrt(BEIDOU_GM / (a * a * a));
|
||||
}
|
||||
@@ -198,9 +346,6 @@ double Gnss_Ephemeris::sv_clock_relativistic_term(double transmitTime) const
|
||||
// Mean anomaly
|
||||
const double M = this->M_0 + n * tk;
|
||||
|
||||
// Reduce mean anomaly to between 0 and 2pi
|
||||
// M = fmod((M + 2.0 * GNSS_PI), (2.0 * GNSS_PI));
|
||||
|
||||
// Initial guess of eccentric anomaly
|
||||
double E = M;
|
||||
double E_old;
|
||||
@@ -235,19 +380,3 @@ double Gnss_Ephemeris::sv_clock_relativistic_term(double transmitTime) const
|
||||
}
|
||||
return dtr_;
|
||||
}
|
||||
|
||||
|
||||
double Gnss_Ephemeris::check_t(double time) const
|
||||
{
|
||||
const double half_week = 302400.0; // seconds
|
||||
double corrTime = time;
|
||||
if (time > half_week)
|
||||
{
|
||||
corrTime = time - 2.0 * half_week;
|
||||
}
|
||||
else if (time < -half_week)
|
||||
{
|
||||
corrTime = time + 2.0 * half_week;
|
||||
}
|
||||
return corrTime;
|
||||
}
|
||||
|
||||
@@ -19,6 +19,7 @@
|
||||
#ifndef GNSS_SDR_GNSS_EPHEMERIS_H
|
||||
#define GNSS_SDR_GNSS_EPHEMERIS_H
|
||||
|
||||
#include <array>
|
||||
#include <cstdint>
|
||||
|
||||
/*!
|
||||
@@ -36,6 +37,34 @@ public:
|
||||
*/
|
||||
double sv_clock_drift(double transmitTime);
|
||||
|
||||
/*!
|
||||
* \brief Computes prediction of the Doppler shift for a given time and receiver's position and velocity.
|
||||
* \f[
|
||||
* f_{d} = - \mathbf{v} \frac{\mathbf{x}^{T}}{\left| \mathbf{x} \right| } \frac{f_{L}}{c}
|
||||
* \f]
|
||||
* where:
|
||||
* \f[
|
||||
* \mathbf{v} = \mathbf{v}_{sat} - \mathbf{v}_{rx}
|
||||
* \f]
|
||||
* \f[
|
||||
* \mathbf{x} = \mathbf{x}_{sat} - \mathbf{x}_{rx}
|
||||
* \f]
|
||||
* \f[
|
||||
* \left| \mathbf{x} \right| = \sqrt{\mathbf{x}\mathbf{x}^{T}}
|
||||
* \f]
|
||||
*
|
||||
* @param[in] rx_time_s Time of Week in seconds
|
||||
* @param[in] lat Receiver's latitude in degrees
|
||||
* @param[in] lon Receiver's longitude in degrees
|
||||
* @param[in] h Receiver's height in meters
|
||||
* @param[in] ve Receiver's velocity in the East direction [m/s]
|
||||
* @param[in] vn Receiver's velocity in the North direction [m/s]
|
||||
* @param[in] vu Receiver's velocity in the Up direction [m/s]
|
||||
* @param[in] band Signal band for which the Doppler will be computed
|
||||
* (1: L1 C/A, E1B, BI1; 2: L2C, BI2; 3: BI3; 5: L5/E5a; 6: E6B; 7: E5b; 8: E5a+E5b)
|
||||
*/
|
||||
double predicted_doppler(double rx_time_s, double lat, double lon, double h, double ve, double vn, double vu, int band) const;
|
||||
|
||||
void satellitePosition(double transmitTime); //!< Computes the ECEF SV coordinates and ECEF velocity
|
||||
|
||||
uint32_t PRN{}; //!< SV ID
|
||||
@@ -83,6 +112,7 @@ protected:
|
||||
char System{}; //!< Character ID of the GNSS system. 'G': GPS. 'E': Galileo. 'B': BeiDou
|
||||
|
||||
private:
|
||||
void satellitePosVelComputation(double transmitTime, std::array<double, 7>& pos_vel_dtr) const;
|
||||
double check_t(double time) const;
|
||||
double sv_clock_relativistic_term(double transmitTime) const;
|
||||
};
|
||||
|
||||
@@ -525,7 +525,7 @@ std::string Gnss_Satellite::what_block(const std::string& system_, uint32_t PRN_
|
||||
}
|
||||
if (system_ == "Galileo")
|
||||
{
|
||||
// Check https://en.wikipedia.org/wiki/List_of_Galileo_satellites and https://www.gsc-europa.eu/system-status/Constellation-Information
|
||||
// Check https://en.wikipedia.org/wiki/List_of_Galileo_satellites and https://www.gsc-europa.eu/system-service-status/constellation-information
|
||||
switch (PRN_)
|
||||
{
|
||||
case 1:
|
||||
|
||||
@@ -38,7 +38,10 @@ public:
|
||||
/*!
|
||||
* Default constructor
|
||||
*/
|
||||
Gps_Almanac() = default;
|
||||
Gps_Almanac()
|
||||
{
|
||||
this->System = 'G';
|
||||
};
|
||||
|
||||
int32_t SV_health{}; //!< SV Health
|
||||
int32_t AS_status{}; //!< Anti-Spoofing Flags and SV Configuration
|
||||
|
||||
@@ -465,18 +465,6 @@ Gps_Ephemeris Gps_Navigation_Message::get_ephemeris() const
|
||||
ephemeris.alert_flag = b_alert_flag;
|
||||
ephemeris.antispoofing_flag = b_antispoofing_flag;
|
||||
|
||||
// These parameters are empty; can be computed later with
|
||||
// ephemeris.sv_clock_drift(double transmitTime);
|
||||
// ephemeris.satellitePosition(double transmitTime);
|
||||
ephemeris.satClkDrift = d_satClkDrift;
|
||||
ephemeris.dtr = d_dtr;
|
||||
ephemeris.satpos_X = d_satpos_X;
|
||||
ephemeris.satpos_Y = d_satpos_Y;
|
||||
ephemeris.satpos_Z = d_satpos_Z;
|
||||
ephemeris.satvel_X = d_satvel_X;
|
||||
ephemeris.satvel_Y = d_satvel_Y;
|
||||
ephemeris.satvel_Z = d_satvel_Z;
|
||||
|
||||
return ephemeris;
|
||||
}
|
||||
|
||||
|
||||
@@ -195,22 +195,6 @@ private:
|
||||
int32_t i_Toa{}; // Almanac reference time [s]
|
||||
int32_t i_WN_A{}; // Modulo 256 of the GPS week number to which the almanac reference time (i_Toa) is referenced
|
||||
|
||||
// clock terms
|
||||
// double d_master_clock{}; // GPS transmission time
|
||||
|
||||
double d_dtr{}; // relativistic clock correction term
|
||||
double d_satClkDrift{};
|
||||
|
||||
// satellite positions
|
||||
double d_satpos_X{}; // Earth-fixed coordinate x of the satellite [m]. Intersection of the IERS Reference Meridian (IRM) and the plane passing through the origin and normal to the Z-axis.
|
||||
double d_satpos_Y{}; // Earth-fixed coordinate y of the satellite [m]. Completes a right-handed, Earth-Centered, Earth-Fixed orthogonal coordinate system.
|
||||
double d_satpos_Z{}; // Earth-fixed coordinate z of the satellite [m]. The direction of the IERS (International Earth Rotation and Reference Systems Service) Reference Pole (IRP).
|
||||
|
||||
// Satellite velocity
|
||||
double d_satvel_X{}; // Earth-fixed velocity coordinate x of the satellite [m]
|
||||
double d_satvel_Y{}; // Earth-fixed velocity coordinate y of the satellite [m]
|
||||
double d_satvel_Z{}; // Earth-fixed velocity coordinate z of the satellite [m]
|
||||
|
||||
// satellite identification info
|
||||
int32_t i_channel_ID{};
|
||||
uint32_t i_satellite_PRN{};
|
||||
|
||||
Reference in New Issue
Block a user