gnav: Adding GLONASS GNAV Ephemeris, Almanac and UTC Model

Generates code for GLONASS GNAV Ephemeris, Almanac and UTC Model with
their respective decoding position indexes for string decoding and
message parsing. Starts developing of satellite position computation
based on its ephemeris and almanac information
This commit is contained in:
Damian Miralles 2017-06-06 07:27:48 -07:00 committed by Damian Miralles
parent 3d5a55276b
commit f8f3574090
7 changed files with 1196 additions and 0 deletions

View File

@ -0,0 +1,176 @@
/*!
* \file GLONASS_L1_CA.h
* \brief Defines system parameters for GLONASS L1 C/A signal and NAV data
* \author Damian Miralles, 2017. dmiralles2009(at)gmail.com
*
* -------------------------------------------------------------------------
*
* Copyright (C) 2010-2015 (see AUTHORS file for a list of contributors)
*
* GNSS-SDR is a software defined Global Navigation
* Satellite Systems receiver
*
* This file is part of GNSS-SDR.
*
* GNSS-SDR is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* GNSS-SDR is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNSS-SDR. If not, see <http://www.gnu.org/licenses/>.
*
* -------------------------------------------------------------------------
*/
#ifndef GNSS_SDR_GLONASS_L1_CA_H_
#define GNSS_SDR_GLONASS_L1_CA_H_
#include <vector>
#include <utility> // std::pair
#include "MATH_CONSTANTS.h"
#include "gnss_frequencies.h"
// Physical constants
const double GLONASS_C_m_s = SPEED_OF_LIGHT; //!< The speed of light, [m/s]
const double GLONASS_C_m_ms = 299792.4580; //!< The speed of light, [m/ms]
const double GLONASS_PI = 3.1415926535898; //!< Pi as defined in IS-GPS-200E
const double GLONASS_TWO_PI = 6.283185307179586; //!< 2Pi as defined in IS-GPS-200E
const double GLONASS_OMEGA_EARTH_DOT = 7.292115e-5; //!< Earth rotation rate, [rad/s]
const double GLONASS_GM = 398600.4418e9; //!< Universal gravitational constant times the mass of the Earth, [km^3/s^2]
const double GLONASS_fM_a = 0.35e9; //!< Gravitational constant of atmosphere [m^3/s^2]
const double GLONASS_SEMI_MAJOR_AXIS = 6378136; //!< Semi-major axis of Earth [m]
const double GLONASS_FLATTENING = 1/29825784; //!< Flattening parameter
const double GLONASS_GRAVITY = 97803284; //!< Equatorial acceleration of gravity [mGal]
const double GLONASS_GRAVITY_CORRECTION = 0.87; //!< Correction to acceleration of gravity at sea-level due to Atmosphere[мGal]
const double GLONASS_J2 = 1082625.75e-9; //!< Second zonal harmonic of the geopotential
const double GLONASS_J4 = -2370.89e-9; //!<Fourth zonal harmonic of the geopotential
const double GLONASS_J6 = 6.08e-9; //!< Sixth zonal harmonic of the geopotential
const double GLONASS_J8 = 1.40e-11; //!< Eighth zonal harmonic of the geopotential
const double GLONASS_U0 = 62636861.4; //!< Normal potential at surface of common terrestrial ellipsoid [m^2/s^2]
const double GLONASS_C20 = -1082.63e-6; //!< Second zonal coefficient of spherical harmonic expansion
const double GLONASS_EARTH_RADIUS = 6378.136; //!< Equatorial radius of Earth [km]
const double GLONASS_EARTH_INCLINATION = 23°26'33'' //!< Mean inclination of ecliptic to equator (23°26'33'').
const double GLONASS_TAU_0 = -334°1946,40;
const double GLONASS_TAU_1 = 4069°0202,52;
const double GLONASS_MOON_Q0 = -63°5343,41;
const double GLONASS_MOON_Q1 = 477198°5056,79;
const double GLONASS_MOON_OMEGA_0 = 259°1059,79;
const double GLONASS_MOON_OMEGA_1 = -1934°0831,23;
const double GLONASS_MOON_GM = 4902.835; //!< Lunar gravitational constant [km^3/s^2]
const double GLONASS_MOON_SEMI_MAJOR_AXIS = 3.84385243e5; //!< Semi-major axis of lunar orbit [km];
const double GLONASS_MOON_ECCENTRICITY = 0.054900489; //!< Eccentricity of lunar orbit
const double GLONASS_MOON_INCLINATION = 5°08'43.4'' //!< Inclination of lunar orbit to ecliptic plane (5°08'43.4'')
const double GLONASS_SUN_OMEGA = 281°1315,0 + 6189, 03Т;
const double GLONASS_SUN_Q0 = 358°2833,04;
const double GLONASS_SUN_Q1 = 129596579,10;
const double GLONASS_SUN_GM = 0.1325263e12; //!< Solar gravitational constant [km^3/s^2]
const double GLONASS_SUN_SEMI_MAJOR_AXIS = 1.49598e8; //!< Semi-major axis of solar orbit [km];
const double GLONASS_SUN_ECCENTRICITY = 0.016719; //!< Eccentricity of solar orbit
// carrier and code frequencies
const double GLONASS_L1_FREQ_HZ = FREQ1_GLO; //!< L1 [Hz]
const double GLONASS_L1_DFREQ_HZ = DFRQ1_GLO; //!< Freq Bias for GLONASS L1 [Hz]
const double GLONASS_L1_CA_CODE_RATE_HZ = 0.511e6; //!< GLONASS L1 C/A code rate [chips/s]
const double GLONASS_L1_CA_CODE_LENGTH_CHIPS = 511.0; //!< GLONASS L1 C/A code length [chips]
const double GLONASS_L1_CA_CODE_PERIOD = 0.001; //!< GLONASS L1 C/A code period [seconds]
const double GLONASS_L1_CA_CHIP_PERIOD = 1.9569e-06; //!< GLONASS L1 C/A chip period [seconds]
const double GLONASS_STARTOFFSET_ms = 68.802; //[ms] Initial sign. travel time (this cannot go here)
// OBSERVABLE HISTORY DEEP FOR INTERPOLATION
const int GLONASS_L1_CA_HISTORY_DEEP = 100;
// NAVIGATION MESSAGE DEMODULATION AND DECODING
#define GLONASS_CA_PREAMBLE {1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0}
const int GLONASS_CA_PREAMBLE_LENGTH_BITS = 30;
const int GLONASS_CA_PREAMBLE_LENGTH_SYMBOLS = 300;
const double GLONASS_CA_PREAMBLE_DURATION_S = 0.3;
const int GLONASS_CA_TELEMETRY_RATE_BITS_SECOND = 50; //!< NAV message bit rate [bits/s]
const int GLONASS_CA_TELEMETRY_SYMBOLS_PER_BIT = 10;
const int GLONASS_CA_TELEMETRY_RATE_SYMBOLS_SECOND = GLONASS_CA_TELEMETRY_RATE_BITS_SECOND*GLONASS_CA_TELEMETRY_SYMBOLS_PER_BIT; //!< NAV message bit rate [symbols/s]
const int GLONASS_GNAV_WORD_LENGTH = 4; //!< \TODO cHECK this, size is not integer bte size
const int GLONASS_GNAV_FRAME_LENGTH = 40; //!< \TODO GPS_WORD_LENGTH x 10 = 40 bytes
const int GLONASS_GNAV_FRAME_BITS = 1725; //!< Number of chips per frame in the GNAV message 15 strings*(85 data bits + 30 time mark bits)[bits]
const int GLONASS_GNAV_FRAME_SECONDS = 30; //!< Subframe duration [seconds]
const int GLONASS_GNAV_FRAME_MS = 30000; //!< Subframe duration [seconds]
const int GLONASS_GNAV_STRING_BITS = 115; //!< Number of bits per string in the GNAV message (85 data bits + 30 time mark bits) [bits]
// GLONASS GNAV NAVIGATION MESSAGE STRUCTURE
// NAVIGATION MESSAGE FIELDS POSITIONS (from IS-GPS-200E Appendix II)
// FRAME 1-4
// COMMON FIELDS
const std::vector<std::pair<int,int>> STRING_ID({{2,4}});
const std::vector<std::pair<int,int>> KX({{78,8}});
//STRING 1
const std::vector<std::pair<int,int>> P1({{8,2}});
const std::vector<std::pair<int,int>> T_K({{10,12}});
const std::vector<std::pair<int,int>> X_N_DOT ({{22,24}});
const std::vector<std::pair<int,int>> X_N_DOT_DOT ({{46,5}});
const std::vector<std::pair<int,int>> X_N({{51,27}});
//STRING 2
const std::vector<std::pair<int,int>> B_N({{6,3}});
const std::vector<std::pair<int,int>> P2({{9,1}});
const std::vector<std::pair<int,int>> T_B({{10,7}});
const std::vector<std::pair<int,int>> Y_N_DOT ({{22,24}});
const std::vector<std::pair<int,int>> Y_N_DOT_DOT ({{46,5}});
const std::vector<std::pair<int,int>> Y_N({{51,27}});
//STRING 3
const std::vector<std::pair<int,int>> P3({{6,1}});
const std::vector<std::pair<int,int>> GAMMA_N({{7,11}});
const std::vector<std::pair<int,int>> P({{19,2}});
const std::vector<std::pair<int,int>> L_N({{21,1}});
const std::vector<std::pair<int,int>> Z_N_DOT ({{22,24}});
const std::vector<std::pair<int,int>> Z_N_DOT_DOT ({{46,5}});
const std::vector<std::pair<int,int>> Z_N({{51,27}});
// STRING 4
const std::vector<std::pair<int,int>> TAU_N({{6,22}});
const std::vector<std::pair<int,int>> DELTA_TAU_N({{28,5}});
const std::vector<std::pair<int,int>> E_N({{33,5}});
const std::vector<std::pair<int,int>> P4 ({{52,1}});
const std::vector<std::pair<int,int>> F_T ({{53,4}});
const std::vector<std::pair<int,int>> N_T({{60,11}});
const std::vector<std::pair<int,int>> N({{71,5}});
const std::vector<std::pair<int,int>> M({{76,2}});
// STRING 5
const std::vector<std::pair<int,int>> N_A({{6,11}});
const std::vector<std::pair<int,int>> TAU_C({{17,32}});
const std::vector<std::pair<int,int>> N_4({{50,5}});
const std::vector<std::pair<int,int>> TAU_GPS({{55,22}});
const std::vector<std::pair<int,int>> L_N({{77,1}});
// STRING 6, 8, 10, 12, 14
const std::vector<std::pair<int,int>> C_N({{6,1}});
const std::vector<std::pair<int,int>> M_N_A({{7,2}});
const std::vector<std::pair<int,int>> n_A({{9,5}});
const std::vector<std::pair<int,int>> TAU_N_A({{14,10}});
const std::vector<std::pair<int,int>> LAMBDA_N_A({{24,21}});
const std::vector<std::pair<int,int>> DELTA_I_N_A({{45,18}});
const std::vector<std::pair<int,int>> EPSILON_N_A({{63,15}});
//STRING 7, 9, 11, 13, 15
const std::vector<std::pair<int,int>> OMEGA_N_A({{6,16}});
const std::vector<std::pair<int,int>> T_LAMBDA_N_A({{22,21}});
const std::vector<std::pair<int,int>> DELTA_T_N_A({{43,22}});
const std::vector<std::pair<int,int>> DELTA_T_DOT_N_A({{65,7}});
const std::vector<std::pair<int,int>> H_N_A({{72,5}});
const std::vector<std::pair<int,int>> L_N({{77,1}});
#endif /* GNSS_SDR_GLONASS_L1_CA_H_ */

View File

@ -0,0 +1,57 @@
/*!
* \file glonass_gnav_almanac.cc
* \brief Interface of a GLONASS GNAV ALMANAC storage
*
* See http://www.gps.gov/technical/icwg/IS-GPS-200E.pdf Appendix II
* \author Damian Miralles , 2017. dmiralles2009(at)gmail.com
*
* -------------------------------------------------------------------------
*
* Copyright (C) 2010-2015 (see AUTHORS file for a list of contributors)
*
* GNSS-SDR is a software defined Global Navigation
* Satellite Systems receiver
*
* This file is part of GNSS-SDR.
*
* GNSS-SDR is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* GNSS-SDR is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNSS-SDR. If not, see <http://www.gnu.org/licenses/>.
*
* -------------------------------------------------------------------------
*/
#include "glonass_gnav_almanac.h"
Glonass_Gnav_Almanac::Glonass_Gnav_Almanac()
{
int i_satellite_freq_channel = 0;
double d_tau_c = 0.0;
double d_tau_gps = 0.0;
double d_N_4 = 0.0;
double d_N_A = 0.0;
double d_n_A = 0.0;
double d_H_n_A = 0.0;
double d_lambda_n_A = 0.0;
double d_t_lambda_n_A = 0.0;
double d_Delta_i_n_A = 0.0;
double d_Delta_T_n_A = 0.0;
double d_Delta_T_n_A_dot = 0.0;
double d_epsilon_n_A = 0.0;
double d_omega_n_A = 0.0;
double d_M_n_A = 0.0;
double d_B1 = 0.0;
double d_B2 = 0.0;
double d_KP = 0.0;
double d_tau_n_A = 0.0;
double d_C_n_A = 0.0;
}

View File

@ -0,0 +1,70 @@
/*!
* \file glonass_gnav_almanac.h
* \brief Interface of a GLONASS GNAV ALMANAC storage
* \author Damian Miralles, 2017. dmiralles2009@gmail.com
*
* -------------------------------------------------------------------------
*
* Copyright (C) 2010-2015 (see AUTHORS file for a list of contributors)
*
* GNSS-SDR is a software defined Global Navigation
* Satellite Systems receiver
*
* This file is part of GNSS-SDR.
*
* GNSS-SDR is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* GNSS-SDR is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNSS-SDR. If not, see <http://www.gnu.org/licenses/>.
*
* -------------------------------------------------------------------------
*/
#ifndef GNSS_SDR_GLONASS_ALMANAC_H_
#define GNSS_SDR_GLONASS_ALMANAC_H_
/*!
* \brief This class is a storage for the GLONASS SV ALMANAC data as described in IS-GPS-200E
* \todo Add proper links for GLONASS ICD
* See http://www.gps.gov/technical/icwg/IS-GPS-200E.pdf Appendix II
*/
class Glonass_Gnav_Almanac
{
public:
int i_satellite_freq_channel; //!< SV PRN NUMBER
double d_tau_c; //!< GLONASS time scale correction to UTC(SU) time. [s]
double d_tau_gps; //!< Correction to GPS time to GLONASS time [day]
double d_N_4; //!< Four year interval number starting from 1996 [4 year interval]
double d_N_A; //!< Calendar day number within the four-year period beginning since the leap year [days]
double d_n_A; //!< Conventional number of satellite within GLONASS space segment [dimensionless]
double d_H_n_A; //!< Carrier frequency number of navigation RF signal transmitted by d_nA satellite [dimensionless]
double d_lambda_n_A; //!< Longitude of the first (within the d_NA day) ascending node of d_nA [semi-circles]
double d_t_lambda_n_A; //!< Time of first ascending node passage [s]
double d_Delta_i_n_A //!< Correction of the mean value of inclination of d_n_A satellite at instant t_lambda_n_A [semi-circles]
double d_Delta_T_n_A; //!< Correction to the mean value of Draconian period of d_n_A satellite at instant t_lambda_n_A[s / orbital period]
double d_Delta_T_n_A_dot; //!< Rate of change of Draconian period of d_n_A satellite at instant t_lambda_n_A [s / orbital period^2]
double d_epsilon_n_A; //!< Eccentricity of d_n_A satellite at instant t_lambda_n_A [dimensionless]
double d_omega_n_A; //!< Argument of preigree of d_n_A satellite at instant t_lambdan_A [semi-circles]
double d_M_n_A; //!< Type of satellite n_A [dimensionless]
double d_B1; //!< Coefficient to determine DeltaUT1 [s]
double d_B2; //!< Coefficient to determine DeltaUT1 [s/msd]
double d_KP; //!< Notification on forthcoming leap second correction of UTC [dimensionless]
double d_tau_n_A; //!< Coarse value of d_n_A satellite time correction to GLONASS time at instant t_lambdan_A[s]
double d_C_n_A; //!< Generalized “unhealthy flag” of n_A satellite at instant of almanac upload [dimensionless]
/*!
* Default constructor
*/
Glonass_Gnav_Almanac();
};
#endif

View File

@ -0,0 +1,600 @@
/*!
* \file gps_ephemeris.cc
* \brief Interface of a GPS EPHEMERIS storage and orbital model functions
*
* See http://www.gps.gov/technical/icwg/IS-GPS-200E.pdf Appendix II
* \author Javier Arribas, 2013. jarribas(at)cttc.es
*
* -------------------------------------------------------------------------
*
* Copyright (C) 2010-2015 (see AUTHORS file for a list of contributors)
*
* GNSS-SDR is a software defined Global Navigation
* Satellite Systems receiver
*
* This file is part of GNSS-SDR.
*
* GNSS-SDR is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* GNSS-SDR is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNSS-SDR. If not, see <http://www.gnu.org/licenses/>.
*
* -------------------------------------------------------------------------
*/
#include "gps_ephemeris.h"
#include <cmath>
#include "GLONASS_L1_CA.h"
#include "gnss_satellite.h"
Gps_Ephemeris::Gps_Ephemeris()
{
i_satellite_freq_channel = 0;
d_m = 0.0; //!< String number within frame [dimensionless]
d_t_k = 0.0; //!< Time referenced to the beginning of the frame within the current day [hours, minutes, seconds]
d_t_b = 0.0; //!< Index of a time interval within current day according to UTC(SU) + 03 hours 00 min. [minutes]
d_M = 0.0; //!< Type of satellite transmitting navigation signal [dimensionless]
d_gamma_n = 0.0; //!< Relative deviation of predicted carrier frequency value of n- satellite from nominal value at the instant tb [dimensionless]
d_tau_n = 0.0; //!< Correction to the nth satellite time (tn) relative to GLONASS time (te),
// satellite positions
d_satpos_X = 0.0; //!< Earth-fixed coordinate x of the satellite in PZ-90.02 coordinate system [km].
d_satpos_Y = 0.0; //!< Earth-fixed coordinate y of the satellite in PZ-90.02 coordinate system [km]
d_satpos_Z = 0.0; //!< Earth-fixed coordinate z of the satellite in PZ-90.02 coordinate system [km]
// Satellite velocity
d_satvel_X = 0.0; //!< Earth-fixed velocity coordinate x of the satellite in PZ-90.02 coordinate system [km/s]
d_satvel_Y = 0.0; //!< Earth-fixed velocity coordinate y of the satellite in PZ-90.02 coordinate system [km/s]
d_satvel_Z = 0.0; //!< Earth-fixed velocity coordinate z of the satellite in PZ-90.02 coordinate system [km/s]
// Satellite acceleration
d_satacc_X = 0.0; //!< Earth-fixed acceleration coordinate x of the satellite in PZ-90.02 coordinate system [km/s^2]
d_satacc_Y = 0.0; //!< Earth-fixed acceleration coordinate y of the satellite in PZ-90.02 coordinate system [km/s^2]
d_satacc_Z = 0.0; //!< Earth-fixed acceleration coordinate z of the satellite in PZ-90.02 coordinate system [km/s^2]
d_B_n = 0.0; //!< Health flag [dimensionless]
d_P = 0.0; //!< Technological parameter of control segment, indication the satellite operation mode in respect of time parameters [dimensionless]
d_N_T = 0.0; //!< Current date, calendar number of day within four-year interval starting from the 1-st of January in a leap year [days]
d_F_T = 0.0; //!< Parameter that provides the predicted satellite user range accuracy at time tb [dimensionless]
d_n = 0.0; //!< Index of the satellite transmitting given navigation signal. It corresponds to a slot number within GLONASS constellation
d_Delta_tau_n = 0.0; //!< Time difference between navigation RF signal transmitted in L2 sub- band and aviation RF signal transmitted in L1 sub-band by nth satellite. [dimensionless]
d_E_n = 0.0; //!< Characterises "age" of a current information [days]
d_P_1 = 0.0; //!< Flag of the immediate data updating.
d_P_2 = 0.0; //!< Flag of oddness ("1") or evenness ("0") of the value of (tb) [dimensionless]
d_P_3 = 0.0; //!< Flag indicating a number of satellites for which almanac is transmitted within given frame: "1" corresponds to 5 satellites and "0" corresponds to 4 satellites [dimensionless]
d_P_4 = 0.0; //!< Flag to show that ephemeris parameters are present. "1" indicates that updated ephemeris or frequency/time parameters have been uploaded by the control segment [dimensionless]
d_l_n = 0.0; //!< Health flag for nth satellite; ln = 0 indicates the n-th satellite is helthy, ln = 1 indicates malfunction of this nth satellite [dimensionless]
// clock terms derived from ephemeris data
d_satClkDrift = 0.0; //!< GLONASS clock error
d_dtr = 0.0;
}
void Glonass_GNAV_Ephemeris::gravitational_perturbations()
{
double T;
double sigma_days;
double tau_prime;
double Omega_moon;
double q_moon;
double q_sun;
double Tau_11; double Tau_12;
double Eta_11; double Eta_12;
double xi_11; double xi_12;
double xi_star;
double eta_star;
double zeta_star;
double E_moon;
double E_sun;
double xi_11;
double xi_12;
double eta_11;
double eta_12;
double zeta_11;
double zeta_12;
double sin_theta_moon;
double cos_theta_moon;
double theta_moon;
double xi_moon_e;
double eta_moon_e;
double zeta_moon_e;
double sin_theta_sun;
double cos_theta_sun;
double theta_sun;
double xi_sun_e;
double eta_sun_e;
double zeta_sun_e;
double r_moon_e;
double mu_bar_moon;
double x_bar_moon;
double y_bar_moon;
double Delta_o_moon;
double mu_bar_sun;
double x_bar_sun;
double y_bar_sun;
double Delta_o_sun;
//<! Directive Cosine computation
/// \TODO Need to define sigma days
T = (27392.375 + sigma_days + (t_e/86400))/(36525);
tau_prime = GLONASS_TAU_0 + GLONASS_TAU_1 * T;
Omega_moon = GLONASS_MOON_OMEGA_0 + GLONASS_MOON_OMEGA_1*T;
q_moon = GLONASS_MOON_Q0 + GLONASS_MOON_Q1*T;
q_sun = GLONASS_SUN_Q0 + GLONASS_SUN_Q1*T;
xi_star = 1 - (cos(Omega_moon)*cos(Omega_moon))*(1 - cos(GLONASS_MOON_INCLINATION));
eta_star = sin(Omega_moon)*(sin(GLONASS_MOON_INCLINATION));
zeta_star = cos(Omega_moon)*(sin(GLONASS_MOON_INCLINATION));
xi_11 = sin(Omega_moon)*cos(Omega_moon)*(1 - cos(GLONASS_MOON_INCLINATION));
xi_12 = 1 - (sin(Omega_moon)*sin(Omega_moon))*(1 - cos(GLONASS_MOON_INCLINATION));
eta_11 = xi_star*cos(GLONASS_EARTH_INCLINATION) - xi_star*sin(GLONASS_EARTH_INCLINATION);
eta_12 = xi_11*cos(GLONASS_EARTH_INCLINATION) + eta_star*sin(GLONASS_EARTH_INCLINATION);
zeta_11 = xi_star*sin(GLONASS_EARTH_INCLINATION) + zeta_star*cos(GLONASS_EARTH_INCLINATION);
zeta_12 = xi_11*sin(GLONASS_EARTH_INCLINATION) + eta_star*cos(GLONASS_EARTH_INCLINATION);
/// \TODO why is this calling sin(E_moon) in the same line where it is being computed
E_moon = q_moon + GLONASS_MOON_ECCENTRICITY*sin(E_moon);
sin_theta_moon = sqrt(1 - GLONASS_MOON_ECCENTRICITY*GLONASS_MOON_ECCENTRICITY)*sin(E_moon)/(1 - GLONASS_MOON_ECCENTRICITY*cos(E_moon));
cos_theta_moon = cos(E_moon - GLONASS_MOON_ECCENTRICITY)/(1 - GLONASS_MOON_ECCENTRICITY*cos(E_moon));
theta_moon = asin(sin_theta_moon);
xi_moon_e = sin(theta_moon + tau_prime)*xi_11 + cos(theta_moon + tau_prime)*xi_12;
eta_moon_e = sin(theta_moon + tau_prime)*eta_11 + cos(theta_moon + tau_prime)*eta_12;
zeta_moon_e = sin(theta_moon + tau_prime)*zeta_11 + cos(theta_moon + tau_prime)*zeta_12;
r_moon_e = GLONASS_MOON_SEMI_MAJOR_AXIS*(1-GLONASS_MOON_ECCENTRICITY*cos(E_moon));
/// \TODO why is this calling sin(E_moon) in the same line where it is being computed
E_sun = q_sun + GLONASS_SUN_ECCENTRICITY*sin(E_sun);
sin_theta_sun = sqrt(1 - GLONASS_SUN_ECCENTRICITY*GLONASS_SUN_ECCENTRICITY)*sin(E_sun)/(1 - GLONASS_SUN_ECCENTRICITY*cos(E_sun));
cos_theta_sun = cos(E_sun - GLONASS_SUN_ECCENTRICITY)/(1 - GLONASS_SUN_ECCENTRICITY*cos(E_sun));
theta_sun = asin(sin_theta_sun);
xi_sun_e = (cos_theta_sun*cos(GLONASS_SUN_OMEGA) - sin_theta_sun*sin(GLONASS_SUN_OMEGA));
eta_sun_e = (sin_theta_sun*cos(GLONASS_SUN_OMEGA) + cos_theta_sun*sin(GLONASS_SUN_OMEGA))*cos(GLONASS_EARTH_INCLINATION);
zeta_sun_e = (sin_theta_sun*cos(GLONASS_SUN_OMEGA) + cos_theta_sun*sin(GLONASS_SUN_OMEGA))*sin(GLONASS_EARTH_INCLINATION);
r_sun_e = GLONASS_SUN_SEMI_MAJOR_AXIS*(1-GLONASS_SUN_ECCENTRICITY*cos(E_sun));
//<! Gravitational computation
mu_bar_moon = GLONASS_MOON_GM/(d_r_moon_e*d_r_moon_e);
x_bar_moon = d_satpos_X/d_r_moon_e;
y_bar_moon = d_satpos_Y/d_r_moon_e;
z_bar_moon = d_satpos_Z/d_r_moon_e;
Delta_o_moon = sqrt((d_xi_moon_e - x_bar_moon)*(d_xi_moon_e - x_bar_moon) + (d_eta_moon_e - y_bar_moon)*(d_eta_moon_e - y_bar_moon) + (d_zeta_moon_e - z_bar_moon)*(d_zeta_moon_e - z_bar_moon));
mu_bar_sun = GLONASS_MOON_GM/(d_r_sun_e*d_r_sun_e);
x_bar_sun = d_satpos_X/d_r_sun_e;
y_bar_sun = d_satpos_Y/d_r_sun_e;
z_bar_sun = d_satpos_Z/d_r_sun_e;
Delta_o_sun = sqrt((d_xi_sun_e - x_bar_sun)*(d_xi_sun_e - x_bar_sun) + (d_eta_sun_e - y_bar_sun)*(d_eta_sun_e - y_bar_sun) + (d_zeta_sun_e - z_bar_sun)*(d_zeta_sun_e - z_bar_sun));
d_Jx_moon = mu_bar_moon*(xi_moon_e - x_bar_moon)/pow(Delta_o_moon,3.0) - xi_moon_e;
d_Jy_moon = mu_bar_moon*(eta_moon_e - y_bar_moon)/pow(Delta_o_moon,3.0) - eta_moon_e;
d_Jz_moon = mu_bar_moon*(zeta_moon_e - z_bar_moon)/pow(Delta_o_moon,3.0) - zeta_moon_e;
d_Jx_sun = mu_bar_sun*(xi_sun_e - x_bar_sun)/pow(Delta_o_sun,3.0) - xi_sun_e;
d_Jy_sun = mu_bar_sun*(eta_sun_e - y_bar_sun)/pow(Delta_o_sun,3.0) - eta_sun_e;
d_Jz_sun = mu_bar_sun*(zeta_sun_e - z_bar_sun)/pow(Delta_o_sun,3.0) - zeta_sun_e;
}
double Glonass_GNAV_Ephemeris::check_t(double time)
{
double corrTime;
double half_day = 43200.0; // seconds
corrTime = time;
if (time > half_day)
{
corrTime = time - 2.0 * half_day;
}
else if (time < -half_day)
{
corrTime = time + 2.0 * half_day;
}
return corrTime;
}
// 20.3.3.3.3.1 User Algorithm for SV Clock Correction.
double Glonass_GNAV_Ephemeris::sv_clock_drift(double transmitTime, double timeCorrUTC)
{
double dt;
dt = check_t(transmitTime - d_t_b);
d_satClkDrift = -(d_tau_n + timeCorrUTC - d_gamma_n * dt);
//Correct satellite group delay and missing relativistic term here
//d_satClkDrift-=d_TGD;
return d_satClkDrift;
}
// compute the relativistic correction term
double Glonass_GNAV_Ephemeris::sv_clock_relativistic_term(double transmitTime)
{
double tk;
double a;
double n;
double n0;
double E;
double E_old;
double dE;
double M;
// Restore semi-major axis
a = d_sqrt_A * d_sqrt_A;
// Time from ephemeris reference epoch
tk = check_t(transmitTime - d_Toe);
// Computed mean motion
n0 = sqrt(GM / (a * a * a));
// Corrected mean motion
n = n0 + d_Delta_n;
// Mean anomaly
M = d_M_0 + n * tk;
// Reduce mean anomaly to between 0 and 2pi
//M = fmod((M + 2.0 * GPS_PI), (2.0 * GPS_PI));
// Initial guess of eccentric anomaly
E = M;
// --- Iteratively compute eccentric anomaly ----------------------------
for (int ii = 1; ii < 20; ii++)
{
E_old = E;
E = M + d_e_eccentricity * sin(E);
dE = fmod(E - E_old, 2.0 * GPS_PI);
if (fabs(dE) < 1e-12)
{
//Necessary precision is reached, exit from the loop
break;
}
}
// Compute relativistic correction term
d_dtr = F * d_e_eccentricity * d_sqrt_A * sin(E);
return d_dtr;
}
double Glonass_GNAV_Ephemeris::simplifiedSatellitePosition(double transmitTime)
{
double dt = 0.0;
double tk = 0.0;
int numberOfIntegrations = 0;
// Find time difference
dt = check_t(transmitTime - d_t_b);
// Calculate clock correction
satClkCorr = -(d_tau_n + d_tau_c - d_gamma_n * dt);
// Find integration time
tk = dt - satClkCorr;
// Check if to integrate forward or backward
if tk < 0
tau = -60;
else
tau = 60;
end
// x,y,z coordinates to meters
x0 = d_satpos_X * 1e3;
y0 = d_satpos_Y * 1e3;
z0 = d_satpos_Z * 1e3;
// x,y,z velocities to meters/s
Vx0 = d_satvel_X * 1e3;
Vy0 = d_satvel_Y * 1e3;
Vz0 = d_satvel_Z * 1e3;
// x,y,z accelerations to meters/sec^2
Ax0 = d_satacc_X * 1e3;
Ay0 = d_satacc_Y * 1e3;
Az0 = d_satacc_Z * 1e3;
for(numberOfIntegrations = tau; numberOfIntegrations < tk+tau; numberOfIntegrations += tau)
{
// Check if last integration step. If last integration step, make one more step that has the remaining time length...
if(fabs(numberOfIntegrations) > fabs(tk))
{
// if there is more time left to integrate...
if mod(tk,tau) ~= 0
{
tau = mod(time,tau);
}
else // otherwise make a zero step.
{
tau = 0;
}
}
// Runge-Kutta Integration Stage K1
r0 = sqrt(x0*x0 + y0*y0 + z0*z0);
r0_2 = r0*r0;
r0_3 = r0*r0*r0;
r0_5 = r0*r0*r0*r0*r0;
dVx_1 = - GLONASS_GM / r0_3 * x0 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a * GLONASS_a) / r0_5 * x0 * (1 - 5 * (z0*z0) / r0_2) + (GLONASS_OMEGA_EARTH_DOT*GLONASS_OMEGA_EARTH_DOT) * x0 + 2 * GLONASS_OMEGA_EARTH_DOT * Vy0 + Ax0;
dVy_1 = - GLONASS_GM / r0_3 * y0 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a * GLONASS_a) / r0_5 * y0 * (1 - 5 * (z0*z0) / r0_2) + (GLONASS_OMEGA_EARTH_DOT*GLONASS_OMEGA_EARTH_DOT) * y0 - 2 * GLONASS_OMEGA_EARTH_DOT * Vx0 + Ay0;
dVz_1 = - GLONASS_GM / r0_3 * z0 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a * GLONASS_a) / r0_5 * z0 * (3 - 5 * (z0*z0) / r0_2) + Az0;
dx_1 = Vx0;
dy_1 = Vy0;
dz_1 = Vz0;
// Runge-Kutta Integration Stage K2
Vx_1 = Vx0 + 1/2 * tau * dVx_1;
Vy_1 = Vy0 + 1/2 * tau * dVy_1;
Vz_1 = Vz0 + 1/2 * tau * dVz_1;
x_1 = x0 + 1/2 * tau * dx_1;
y_1 = y0 + 1/2 * tau * dy_1;
z_1 = z0 + 1/2 * tau * dz_1;
r1 = sqrt( x_1*x_1 + y_1*y_1 + z_1*z_1 );
r1_2 = r1*r1;
r1_3 = r1*r1*r1;
r1_5 = r1*r1*r1*r1*r1;
dVx_2 = - GLONASS_GM / r1_3 * x_1 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a*GLONASS_a) / r1_5 * x_1 * (1 - 5 * (z_1*z_1) / r1_2) + (GLONASS_OMEGA_EARTH_DOT * GLONASS_OMEGA_EARTH_DOT) * x_1 + 2 * GLONASS_OMEGA_EARTH_DOT * Vy_1 + Ax0;
dVy_2 = - GLONASS_GM / r1_3 * y_1 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a*GLONASS_a) / r1_5 * y_1 * (1 - 5 * (z_1*z_1) / r1_2) + (GLONASS_OMEGA_EARTH_DOT * GLONASS_OMEGA_EARTH_DOT) * y_1 - 2 * GLONASS_OMEGA_EARTH_DOT * Vx_1 + Ay0;
dVz_2 = - GLONASS_GM / r1_3 * z_1 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a*GLONASS_a) / r1_5 * z_1 * (3 - 5 * (z_1*z_1) / r1_2) + Az0;
dx_2 = Vx_1;
dy_2 = Vy_1;
dz_2 = Vz_1;
// Runge-Kutta Integration Stage K2
Vx_2 = Vx0 + 1/2 * tau * dVx_2;
Vy_2 = Vy0 + 1/2 * tau * dVy_2;
Vz_2 = Vz0 + 1/2 * tau * dVz_2;
x_2 = x0 + 1/2 * tau * dx_2;
y_2 = y0 + 1/2 * tau * dy_2;
z_2 = z0 + 1/2 * tau * dz_2;
r2 = sqrt( x_2*x_2 + y_2*y_2 + z_2*z_2 );
r2_2 = r2*r2;
r2_3 = r2*r2*r2;
r2_5 = r2*r2*r2*r2*r2;
dVx_3 = - GLONASS_GM / r2_3 * x_2 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a*GLONASS_a) / r2_5 * x_2 * (1 - 5 * (z_2*z_2) / r2_2) + (GLONASS_OMEGA_EARTH_DOT * GLONASS_OMEGA_EARTH_DOT) * x_2 + 2 * GLONASS_OMEGA_EARTH_DOT * Vy_2 + Ax0;
dVy_3 = - GLONASS_GM / r2_3 * y_2 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a*GLONASS_a) / r2_5 * y_2 * (1 - 5 * (z_2*z_2) / r2_2) + (GLONASS_OMEGA_EARTH_DOT * GLONASS_OMEGA_EARTH_DOT) * y_2 - 2 * GLONASS_OMEGA_EARTH_DOT * Vx_2 + Ay0;
dVz_3 = - GLONASS_GM / r2_3 * z_2 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a*GLONASS_a) / r2_5 * z_2 * (3 - 5 * (z_2*z_2) / r2_2) + Az0;
dx_3 = Vx_2;
dy_3 = Vy_2;
dz_3 = Vz_2;
// Runge-Kutta Integration Stage K3
Vx_3 = Vx0 + tau * dVx_3;
Vy_3 = Vy0 + tau * dVy_3;
Vz_3 = Vz0 + tau * dVz_3;
x_3 = x0 + tau * dx_3;
y_3 = y0 + tau * dy_3;
z_3 = z0 + tau * dz_3;
r3 = sqrt( x_3*x_3 + y_3*y_3 + z_3*z_3 );
r3_2 = r3*r3;
r3_3 = r3*r3*r3;
r3_5 = r3*r3*r3*r3*r3;
dVx_4 = - GLONASS_GM / r3_3 * x_3 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a*GLONASS_a) / r3_5 * x_3 * (1 - 5 * (z_3*z_3) / r3_2) + (GLONASS_OMEGA_EARTH_DOT * GLONASS_OMEGA_EARTH_DOT) * x_3 + 2 * GLONASS_OMEGA_EARTH_DOT * Vy_3 + Ax0;
dVy_4 = - GLONASS_GM / r3_3 * y_3 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a*GLONASS_a) / r3_5 * y_3 * (1 - 5 * (z_3*z_3) / r3_2) + (GLONASS_OMEGA_EARTH_DOT * GLONASS_OMEGA_EARTH_DOT) * y_3 - 2 * GLONASS_OMEGA_EARTH_DOT * Vx_3 + Ay0;
dVz_4 = - GLONASS_GM / r3_3 * z_3 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a*GLONASS_a) / r3_5 * z_3 * (3 - 5 * (z_3*z_3) / r3_2) + Az0;
dx_4 = Vx_3;
dy_4 = Vy_3;
dz_4 = Vz_3;
// Final results showcased here
Vx0 = Vx0 + 1/6 * tau * ( dVx_1 + 2 * dVx_2 + 2 * dVx_3 + dVx_4 );
Vy0 = Vy0 + 1/6 * tau * ( dVy_1 + 2 * dVy_2 + 2 * dVy_3 + dVy_4 );
Vz0 = Vz0 + 1/6 * tau * ( dVz_1 + 2 * dVz_2 + 2 * dVz_3 + dVz_4 );
x0 = x0 + 1/6 * tau * ( dx_1 + 2 * dx_2 + 2 * dx_3 + dx_4 );
y0 = y0 + 1/6 * tau * ( dy_1 + 2 * dy_2 + 2 * dy_3 + dy_4 );
z0 = z0 + 1/6 * tau * ( dz_1 + 2 * dz_2 + 2 * dz_3 + dz_4 );
}
// Reasign position, velocities and accelerations for next integration
d_satpos_X = x0;
d_satpos_Y = y0;
d_satpos_Z = z0;
d_satvel_X = Vx0;
d_satvel_Y = Vy0;
d_satvel_Z = Vz0;
d_satacc_X = d_satacc_X; // No change in accelerations reported over interval
d_satacc_Y = d_satacc_Y; // No change in accelerations reported over interval
d_satacc_Z = d_satacc_Z; // No change in accelerations reported over interval
// Time from ephemeris reference clock
tk = check_t(transmitTime - d_Toc);
double dtr_s = d_A_f0 + d_A_f1 * tk + d_A_f2 * tk * tk;
/* relativity correction */
dtr_s -= 2.0 * sqrt(GM * GLONASS_a) * d_e_eccentricity * sin(E) / (GPS_C_m_s * GPS_C_m_s);
return dtr_s;
}
double Glonass_GNAV_Ephemeris::satellitePosition(double transmitTime)
{
double dt = 0.0;
double tk = 0.0;
int numberOfIntegrations = 0;
// Find time difference
dt = check_t(transmitTime - d_t_b);
// Calculate clock correction
satClkCorr = -(d_tau_n + d_tau_c - d_gamma_n * dt);
// Find integration time
tk = dt - satClkCorr;
// Check if to integrate forward or backward
if tk < 0
tau = -60;
else
tau = 60;
end
// Coordinates transformation to an inertial reference frame
// x,y,z coordinates to meters
x0 = d_satpos_X * 1e3;
y0 = d_satpos_Y * 1e3;
z0 = d_satpos_Z * 1e3;
// x,y,z velocities to meters/s
Vx0 = d_satvel_X * 1e3;
Vy0 = d_satvel_Y * 1e3;
Vz0 = d_satvel_Z * 1e3;
// x,y,z accelerations to meters/sec^2
Ax0 = d_satacc_X * 1e3;
Ay0 = d_satacc_Y * 1e3;
Az0 = d_satacc_Z * 1e3;
for(numberOfIntegrations = tau; numberOfIntegrations < tk+tau; numberOfIntegrations += tau)
{
// Check if last integration step. If last integration step, make one more step that has the remaining time length...
if(fabs(numberOfIntegrations) > fabs(tk))
{
// if there is more time left to integrate...
if mod(tk,tau) ~= 0
{
tau = mod(time,tau);
}
else // otherwise make a zero step.
{
tau = 0;
}
}
// Runge-Kutta Integration Stage K1
r0 = sqrt(x0*x0 + y0*y0 + z0*z0);
r0_2 = r0*r0;
r0_3 = r0*r0*r0;
r0_5 = r0*r0*r0*r0*r0;
dx_1 = Vx0;
dy_1 = Vy0;
dz_1 = Vz0;
dVx_1 = - GLONASS_GM / r0_3 * x0 + 3/2 * GLONASS_C20 * GLONASS_GM * (GLONASS_EARTH_RADIUS * GLONASS_EARTH_RADIUS) / r0_5 * x0 * (1 - 5 * 1 / (z0*z0)) + Jx_moon + Jx_sun;
dVy_1 = - GLONASS_GM / r0_3 * y0 + 3/2 * GLONASS_C20 * GLONASS_GM * (GLONASS_EARTH_RADIUS * GLONASS_EARTH_RADIUS) / r0_5 * y0 * (1 - 5 * 1 / (z0*z0)) + Jy_moon + Jy_sun;
dVz_1 = - GLONASS_GM / r0_3 * z0 + 3/2 * GLONASS_C20 * GLONASS_GM * (GLONASS_EARTH_RADIUS * GLONASS_EARTH_RADIUS) / r0_5 * z0 * (3 - 5 * 1 / (z0*z0)) + Jz_moon + Jz_sun;
dx_1 = Vx0;
dy_1 = Vy0;
dz_1 = Vz0;
// Runge-Kutta Integration Stage K2
Vx_1 = Vx0 + 1/2 * tau * dVx_1;
Vy_1 = Vy0 + 1/2 * tau * dVy_1;
Vz_1 = Vz0 + 1/2 * tau * dVz_1;
x_1 = x0 + 1/2 * tau * dx_1;
y_1 = y0 + 1/2 * tau * dy_1;
z_1 = z0 + 1/2 * tau * dz_1;
r1 = sqrt( x_1*x_1 + y_1*y_1 + z_1*z_1 );
r1_2 = r1*r1;
r1_3 = r1*r1*r1;
r1_5 = r1*r1*r1*r1*r1;
dVx_2 = - GLONASS_GM / r1_3 * x_1 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a*GLONASS_a) / r1_5 * x_1 * (1 - 5 * (z_1*z_1) / r1_2) + (GLONASS_OMEGA_EARTH_DOT * GLONASS_OMEGA_EARTH_DOT) * x_1 + 2 * GLONASS_OMEGA_EARTH_DOT * Vy_1 + Ax0;
dVy_2 = - GLONASS_GM / r1_3 * y_1 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a*GLONASS_a) / r1_5 * y_1 * (1 - 5 * (z_1*z_1) / r1_2) + (GLONASS_OMEGA_EARTH_DOT * GLONASS_OMEGA_EARTH_DOT) * y_1 - 2 * GLONASS_OMEGA_EARTH_DOT * Vx_1 + Ay0;
dVz_2 = - GLONASS_GM / r1_3 * z_1 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a*GLONASS_a) / r1_5 * z_1 * (3 - 5 * (z_1*z_1) / r1_2) + Az0;
dx_2 = Vx_1;
dy_2 = Vy_1;
dz_2 = Vz_1;
// Runge-Kutta Integration Stage K2
Vx_2 = Vx0 + 1/2 * tau * dVx_2;
Vy_2 = Vy0 + 1/2 * tau * dVy_2;
Vz_2 = Vz0 + 1/2 * tau * dVz_2;
x_2 = x0 + 1/2 * tau * dx_2;
y_2 = y0 + 1/2 * tau * dy_2;
z_2 = z0 + 1/2 * tau * dz_2;
r2 = sqrt( x_2*x_2 + y_2*y_2 + z_2*z_2 );
r2_2 = r2*r2;
r2_3 = r2*r2*r2;
r2_5 = r2*r2*r2*r2*r2;
dVx_3 = - GLONASS_GM / r2_3 * x_2 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a*GLONASS_a) / r2_5 * x_2 * (1 - 5 * (z_2*z_2) / r2_2) + (GLONASS_OMEGA_EARTH_DOT * GLONASS_OMEGA_EARTH_DOT) * x_2 + 2 * GLONASS_OMEGA_EARTH_DOT * Vy_2 + Ax0;
dVy_3 = - GLONASS_GM / r2_3 * y_2 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a*GLONASS_a) / r2_5 * y_2 * (1 - 5 * (z_2*z_2) / r2_2) + (GLONASS_OMEGA_EARTH_DOT * GLONASS_OMEGA_EARTH_DOT) * y_2 - 2 * GLONASS_OMEGA_EARTH_DOT * Vx_2 + Ay0;
dVz_3 = - GLONASS_GM / r2_3 * z_2 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a*GLONASS_a) / r2_5 * z_2 * (3 - 5 * (z_2*z_2) / r2_2) + Az0;
dx_3 = Vx_2;
dy_3 = Vy_2;
dz_3 = Vz_2;
// Runge-Kutta Integration Stage K3
Vx_3 = Vx0 + tau * dVx_3;
Vy_3 = Vy0 + tau * dVy_3;
Vz_3 = Vz0 + tau * dVz_3;
x_3 = x0 + tau * dx_3;
y_3 = y0 + tau * dy_3;
z_3 = z0 + tau * dz_3;
r3 = sqrt( x_3*x_3 + y_3*y_3 + z_3*z_3 );
r3_2 = r3*r3;
r3_3 = r3*r3*r3;
r3_5 = r3*r3*r3*r3*r3;
dVx_4 = - GLONASS_GM / r3_3 * x_3 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a*GLONASS_a) / r3_5 * x_3 * (1 - 5 * (z_3*z_3) / r3_2) + (GLONASS_OMEGA_EARTH_DOT * GLONASS_OMEGA_EARTH_DOT) * x_3 + 2 * GLONASS_OMEGA_EARTH_DOT * Vy_3 + Ax0;
dVy_4 = - GLONASS_GM / r3_3 * y_3 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a*GLONASS_a) / r3_5 * y_3 * (1 - 5 * (z_3*z_3) / r3_2) + (GLONASS_OMEGA_EARTH_DOT * GLONASS_OMEGA_EARTH_DOT) * y_3 - 2 * GLONASS_OMEGA_EARTH_DOT * Vx_3 + Ay0;
dVz_4 = - GLONASS_GM / r3_3 * z_3 - 3/2 * GLONASS_J2 * GLONASS_GM * (GLONASS_a*GLONASS_a) / r3_5 * z_3 * (3 - 5 * (z_3*z_3) / r3_2) + Az0;
dx_4 = Vx_3;
dy_4 = Vy_3;
dz_4 = Vz_3;
// Final results showcased here
Vx0 = Vx0 + 1/6 * tau * ( dVx_1 + 2 * dVx_2 + 2 * dVx_3 + dVx_4 );
Vy0 = Vy0 + 1/6 * tau * ( dVy_1 + 2 * dVy_2 + 2 * dVy_3 + dVy_4 );
Vz0 = Vz0 + 1/6 * tau * ( dVz_1 + 2 * dVz_2 + 2 * dVz_3 + dVz_4 );
x0 = x0 + 1/6 * tau * ( dx_1 + 2 * dx_2 + 2 * dx_3 + dx_4 );
y0 = y0 + 1/6 * tau * ( dy_1 + 2 * dy_2 + 2 * dy_3 + dy_4 );
z0 = z0 + 1/6 * tau * ( dz_1 + 2 * dz_2 + 2 * dz_3 + dz_4 );
}
// Time from ephemeris reference clock
tk = check_t(transmitTime - d_Toc);
double dtr_s = d_A_f0 + d_A_f1 * tk + d_A_f2 * tk * tk;
/* relativity correction */
dtr_s -= 2.0 * sqrt(GM * GLONASS_a) * d_e_eccentricity * sin(E) / (GPS_C_m_s * GPS_C_m_s);
return dtr_s;
}

View File

@ -0,0 +1,170 @@
/*!
* \file glonass_gnav_ephemeris.h
* \brief Interface of a GLONASS EPHEMERIS storage
* \author Damian Miralles, 2017. dmiralles2009(at)gmail.com
*
* -------------------------------------------------------------------------
*
* Copyright (C) 2010-2015 (see AUTHORS file for a list of contributors)
*
* GNSS-SDR is a software defined Global Navigation
* Satellite Systems receiver
*
* This file is part of GNSS-SDR.
*
* GNSS-SDR is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* GNSS-SDR is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNSS-SDR. If not, see <http://www.gnu.org/licenses/>.
*
* -------------------------------------------------------------------------
*/
#ifndef GNSS_SDR_GLONASS_GNAV_EPHEMERIS_H_
#define GNSS_SDR_GLONASS_GNAV_EPHEMERIS_H_
#include <map>
#include <string>
#include "boost/assign.hpp"
#include <boost/serialization/nvp.hpp>
/*!
* \brief This class is a storage and orbital model functions for the GLONASS SV ephemeris data as described in GLONASS ICD Edition 5.1
*
* See http://www.gps.gov/technical/icwg/IS-GPS-200E.pdf Appendix II
*/
class Glonass_GNAV_Ephemeris
{
private:
/*
* Accounts for the beginning or end of week crossover
*
* See paragraph 20.3.3.3.3.1 (IS-GPS-200E)
* \param[in] - time in seconds
* \param[out] - corrected time, in seconds
*/
double check_t(double time);
double gravitational_perturbations(const double mu);
double d_Jx_moon; //!< Moon gravitational perturbation
double d_Jx_sun; //!< Sun gravitational perturbation
public:
unsigned int i_satellite_freq_channel; // SV Frequency Channel Number
double d_m; //!< String number within frame [dimensionless]
double d_t_k; //!< Time referenced to the beginning of the frame within the current day [hours, minutes, seconds]
double d_t_b; //!< Index of a time interval within current day according to UTC(SU) + 03 hours 00 min. [minutes]
double d_M; //!< Type of satellite transmitting navigation signal [dimensionless]
double d_gamma_n; //!< Relative deviation of predicted carrier frequency value of n- satellite from nominal value at the instant tb [dimensionless]
double d_tau_n; //!< Correction to the nth satellite time (tn) relative to GLONASS time (te),
// satellite positions
double d_satpos_X; //!< Earth-fixed coordinate x of the satellite in PZ-90.02 coordinate system [km].
double d_satpos_Y; //!< Earth-fixed coordinate y of the satellite in PZ-90.02 coordinate system [km]
double d_satpos_Z; //!< Earth-fixed coordinate z of the satellite in PZ-90.02 coordinate system [km]
// Satellite velocity
double d_satvel_X; //!< Earth-fixed velocity coordinate x of the satellite in PZ-90.02 coordinate system [km/s]
double d_satvel_Y; //!< Earth-fixed velocity coordinate y of the satellite in PZ-90.02 coordinate system [km/s]
double d_satvel_Z; //!< Earth-fixed velocity coordinate z of the satellite in PZ-90.02 coordinate system [km/s]
// Satellite acceleration
double d_satacc_X; //!< Earth-fixed acceleration coordinate x of the satellite in PZ-90.02 coordinate system [km/s^2]
double d_satacc_Y; //!< Earth-fixed acceleration coordinate y of the satellite in PZ-90.02 coordinate system [km/s^2]
double d_satacc_Z; //!< Earth-fixed acceleration coordinate z of the satellite in PZ-90.02 coordinate system [km/s^2]
double d_B_n; //!< Health flag [dimensionless]
double d_P; //!< Technological parameter of control segment, indication the satellite operation mode in respect of time parameters [dimensionless]
double d_N_T; //!< Current date, calendar number of day within four-year interval starting from the 1-st of January in a leap year [days]
double d_F_T; //!< Parameter that provides the predicted satellite user range accuracy at time tb [dimensionless]
double d_n; //!< Index of the satellite transmitting given navigation signal. It corresponds to a slot number within GLONASS constellation
double d_Delta_tau_n; //!< Time difference between navigation RF signal transmitted in L2 sub- band and aviation RF signal transmitted in L1 sub-band by nth satellite. [dimensionless]
double d_E_n; //!< Characterises "age" of a current information [days]
double d_P_1; //!< Flag of the immediate data updating.
double d_P_2; //!< Flag of oddness ("1") or evenness ("0") of the value of (tb) [dimensionless]
double d_P_3; //!< Flag indicating a number of satellites for which almanac is transmitted within given frame: "1" corresponds to 5 satellites and "0" corresponds to 4 satellites [dimensionless]
double d_P_4; //!< Flag to show that ephemeris parameters are present. "1" indicates that updated ephemeris or frequency/time parameters have been uploaded by the control segment [dimensionless]
double d_l_n; //!< Health flag for nth satellite; ln = 0 indicates the n-th satellite is helthy, ln = 1 indicates malfunction of this nth satellite [dimensionless]
// clock terms derived from ephemeris data
double d_satClkDrift; //!< GPS clock error
double d_dtr; //!< relativistic clock correction term
template<class Archive>
/*!
* \brief Serialize is a boost standard method to be called by the boost XML serialization. Here is used to save the ephemeris data on disk file.
*/
void serialize(Archive& archive, const unsigned int version)
{
using boost::serialization::make_nvp;
if(version){};
archive & make_nvp("i_satellite_freq_channel", i_satellite_freq_channel); //!< SV PRN frequency channel number
archive & make_nvp("d_m", d_m); //!< String number within frame [dimensionless]
archive & make_nvp("d_t_k", d_t_k); //!< Time referenced to the beginning of the frame within the current day [hours, minutes, seconds]
archive & make_nvp("d_t_b", d_t_b); //!< Index of a time interval within current day according to UTC(SU) + 03 hours 00 min. [minutes]
archive & make_nvp("d_M", d_M); //!< Type of satellite transmitting navigation signal [dimensionless]
archive & make_nvp("d_gamma_n", d_gamma_n); //!< Relative deviation of predicted carrier frequency value of n- satellite from nominal value at the instant tb [dimensionless]
archive & make_nvp("d_tau_n", d_tau_n); //!< Correction to the nth satellite time (tn) relative to GLONASS time (te)
archive & make_nvp("d_B_n", d_B_n); //!< Health flag [dimensionless]
archive & make_nvp("d_P", d_P); //!< Technological parameter of control segment, indication the satellite operation mode in respect of time parameters [dimensionless]
archive & make_nvp("d_N_T", d_N_T); //!< Current date, calendar number of day within four-year interval starting from the 1-st of January in a leap year [days]
archive & make_nvp("d_F_T", d_F_T); //!< Parameter that provides the predicted satellite user range accuracy at time tb [dimensionless]
archive & make_nvp("d_n", d_n); //!< Index of the satellite transmitting given navigation signal. It corresponds to a slot number within GLONASS constellation
archive & make_nvp("d_Delta_tau_n", d_Delta_tau_n);//!< Time difference between navigation RF signal transmitted in L2 sub- band and aviation RF signal transmitted in L1 sub-band by nth satellite. [dimensionless]
archive & make_nvp("d_E_n", d_E_n); //!< Characterises "age" of a current information [days]
archive & make_nvp("d_P_1", d_P_1); //!< Flag of the immediate data updating.
archive & make_nvp("d_P_2", d_P_2); //!< Flag of oddness ("1") or evenness ("0") of the value of (tb) [dimensionless]
archive & make_nvp("d_P_3", d_P_3); //!< Flag indicating a number of satellites for which almanac is transmitted within given frame: "1" corresponds to 5 satellites and "0" corresponds to 4 satellites [dimensionless]
archive & make_nvp("d_P_4", d_P_4); //!< Flag to show that ephemeris parameters are present. "1" indicates that updated ephemeris or frequency/time parameters have been uploaded by the control segment [dimensionless]
archive & make_nvp("d_l_n", d_l_n); //!< Health flag for nth satellite; ln = 0 indicates the n-th satellite is helthy, ln = 1 indicates malfunction of this nth satellite [dimensionless]
}
/*!
* \brief Compute the ECEF SV coordinates and ECEF velocity
* Implementation of Algorithm A.3.1.2 in GLONASS ICD v5.1
* and compute the clock bias term including relativistic effect (return value)
* \param transmitTime Time of ephemeris transmission
* \return clock bias of satellite
*/
double simplifiedSatellitePosition(double transmitTime);
/*!
* \brief Compute the ECEF SV coordinates and ECEF velocity
* Implementation of Algorithm A.3.1.1 in GLONASS ICD v5.1
* and compute the clock bias term including relativistic effect (return value)
* \param transmitTime Time of ephemeris transmission
* \return clock bias of satellite
*/
double satellitePosition(double transmitTime);
/*!
* \brief Sets (\a d_satClkDrift)and returns the clock drift in seconds according to the User Algorithm for SV Clock Correction
* (IS-GPS-200E, 20.3.3.3.3.1)
*/
double sv_clock_drift(double transmitTime);
/*!
* \brief Sets (\a d_dtr) and returns the clock relativistic correction term in seconds according to the User Algorithm for SV Clock Correction
* (IS-GPS-200E, 20.3.3.3.3.1)
*/
double sv_clock_relativistic_term(double transmitTime);
/*!
* Default constructor
*/
Glonass_GNAV_Ephemeris();
};
#endif

View File

@ -0,0 +1,48 @@
/*
* \file glonass_gnav_utc_model.h
* \brief Interface of a GLONASS GNAV UTC MODEL storage
* \author Damian Miralles, 2017. dmiralles2009(at).gmail.com
*
* -------------------------------------------------------------------------
*
* Copyright (C) 2010-2015 (see AUTHORS file for a list of contributors)
*
* GNSS-SDR is a software defined Global Navigation
* Satellite Systems receiver
*
* This file is part of GNSS-SDR.
*
* GNSS-SDR is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* GNSS-SDR is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNSS-SDR. If not, see <http://www.gnu.org/licenses/>.
*
* -------------------------------------------------------------------------
*/
#include "glonass_gnav_utc_model.h"
#include <cmath>
Glonass_Gnav_Utc_Model::Glonass_Gnav_Utc_Model()
{
valid = false;
d_tau_c = 0.0;
}
double Glonass_Gnav_Utc_Model::utc_time(double glonass_time_corrected)
{
double t_utc;
// GLONASS Time is relative to UTC Moscow, so we simply add its time difference
t_utc = glonass_time_corrected + 3*3600 + d_tau_c;
return t_utc;
}

View File

@ -0,0 +1,75 @@
/*!
* \file glonass_gnav_utc_model.h
* \brief Interface of a GLONASS GNAV UTC MODEL storage
* \author Damian Miralles, 2017. dmiralles2009(at)gmail.com
*
* -------------------------------------------------------------------------
*
* Copyright (C) 2010-2015 (see AUTHORS file for a list of contributors)
*
* GNSS-SDR is a software defined Global Navigation
* Satellite Systems receiver
*
* This file is part of GNSS-SDR.
*
* GNSS-SDR is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* GNSS-SDR is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNSS-SDR. If not, see <http://www.gnu.org/licenses/>.
*
* -------------------------------------------------------------------------
*/
#ifndef GNSS_SDR_GLONASS_GNAV_UTC_MODEL_H_
#define GNSS_SDR_GLONASS_GNAV_UTC_MODEL_H_
#include "boost/assign.hpp"
#include <boost/serialization/nvp.hpp>
#include <iostream>
/*!
* \brief This class is a storage for the GLONASS GNAV UTC MODEL data as described in GLONASS ICD (Edition 5.1)
* See http://russianspacesystems.ru/wp-content/uploads/2016/08/ICD_GLONASS_eng_v5.1.pdf
*/
class Glonass_Gnav_Utc_Model
{
public:
bool valid;
// UTC
double d_tau_c; //!< GLONASS time scale correction to UTC(SU) time. [s]
/*!
* Default constructor
*/
Glonass_Gnav_Utc_Model();
/*!
* \brief Computes the Coordinated Universal Time (UTC) and
* returns it in [s] (GLONASS ICD (Edition 5.1) Section 3.3.3 GLONASS Time)
*/
double utc_time(double glonass_time_corrected);
template<class Archive>
/*
* \brief Serialize is a boost standard method to be called by the boost XML serialization. Here is used to save the ephemeris data on disk file.
*/
void serialize(Archive& archive, const unsigned int version)
{
using boost::serialization::make_nvp;
if(version){};
archive & make_nvp("valid",valid);
archive & make_nvp("d_tau_c",d_tau_c);
}
};
#endif