/*!
* \file gps_navigation_message.cc
* \brief Implementation of a GPS NAV Data message decoder as described in IS-GPS-200E
*
* See http://www.gps.gov/technical/icwg/IS-GPS-200E.pdf Appendix II
* \author Javier Arribas, 2011. jarribas(at)cttc.es
*
* -------------------------------------------------------------------------
*
* Copyright (C) 2010-2012 (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 .
*
* -------------------------------------------------------------------------
*/
#include "gps_navigation_message.h"
#include
#include "boost/date_time/posix_time/posix_time.hpp"
#define num_of_slices(x) sizeof(x)/sizeof(bits_slice)
void Gps_Navigation_Message::reset()
{
b_update_tow_flag=false;
b_valid_ephemeris_set_flag=false;
d_TOW=0;
d_TOW_SF1 = 0;
d_TOW_SF2 = 0;
d_TOW_SF3 = 0;
d_TOW_SF4 = 0;
d_TOW_SF5 = 0;
d_IODE_SF2 = 0;
d_IODE_SF3 = 0;
d_Crs = 0;
d_Delta_n = 0;
d_M_0 = 0;
d_Cuc = 0;
d_e_eccentricity = 0;
d_Cus = 0;
d_sqrt_A = 0;
d_Toe = 0;
d_Toc = 0;
d_Cic = 0;
d_OMEGA0 = 0;
d_Cis = 0;
d_i_0 = 0;
d_Crc = 0;
d_OMEGA = 0;
d_OMEGA_DOT = 0;
d_IDOT = 0;
i_code_on_L2 = 0;
i_GPS_week = 0;
b_L2_P_data_flag = false;
i_SV_accuracy = 0;
i_SV_health = 0;
d_TGD = 0;
d_IODC = -1;
i_AODO = 0;
b_fit_interval_flag = false;
d_spare1 = 0;
d_spare2 = 0;
d_A_f0 = 0;
d_A_f1 = 0;
d_A_f2 = 0;
//clock terms
//d_master_clock=0;
d_dtr = 0;
d_satClkCorr = 0;
// satellite positions
d_satpos_X = 0;
d_satpos_Y = 0;
d_satpos_Z = 0;
// info
i_channel_ID = 0;
i_satellite_PRN = 0;
// time synchro
d_subframe_timestamp_ms = 0;
// flags
b_alert_flag = false;
b_integrity_status_flag = false;
b_antispoofing_flag = false;
// Ionosphere and UTC
d_alpha0 = 0;
d_alpha1 = 0;
d_alpha2 = 0;
d_alpha3 = 0;
d_beta0 = 0;
d_beta1 = 0;
d_beta2 = 0;
d_beta3 = 0;
d_A1 = 0;
d_A0 = 0;
d_t_OT = 0;
i_WN_T = 0;
d_DeltaT_LS = 0;
i_WN_LSF = 0;
i_DN = 0;
d_DeltaT_LSF= 0;
//Almanac
d_Toa = 0;
i_WN_A = 0;
for (int i=1; i < 32; i++ )
{
almanacHealth[i] = 0;
}
// Satellite velocity
d_satvel_X = 0;
d_satvel_Y = 0;
d_satvel_Z = 0;
//Plane A (info from http://www.navcen.uscg.gov/?Do=constellationStatus)
satelliteBlock[9] = "IIA";
satelliteBlock[31] = "IIR-M";
satelliteBlock[8] = "IIA";
satelliteBlock[7] = "IIR-M";
satelliteBlock[27] = "IIA";
//Plane B
satelliteBlock[16] = "IIR";
satelliteBlock[25] = "IIF";
satelliteBlock[28] = "IIR";
satelliteBlock[12] = "IIR-M";
satelliteBlock[30] = "IIA";
//Plane C
satelliteBlock[29] = "IIR-M";
satelliteBlock[3] = "IIA";
satelliteBlock[19] = "IIR";
satelliteBlock[17] = "IIR-M";
satelliteBlock[6] = "IIA";
//Plane D
satelliteBlock[2] = "IIR";
satelliteBlock[1] = "IIF";
satelliteBlock[21] = "IIR";
satelliteBlock[4] = "IIA";
satelliteBlock[11] = "IIR";
satelliteBlock[24] = "IIA"; // Decommissioned from active service on 04 Nov 2011
//Plane E
satelliteBlock[20] = "IIR";
satelliteBlock[22] = "IIR";
satelliteBlock[5] = "IIR-M";
satelliteBlock[18] = "IIR";
satelliteBlock[32] = "IIA";
satelliteBlock[10] = "IIA";
//Plane F
satelliteBlock[14] = "IIR";
satelliteBlock[15] = "IIR-M";
satelliteBlock[13] = "IIR";
satelliteBlock[23] = "IIR";
satelliteBlock[26] = "IIA";
}
Gps_Navigation_Message::Gps_Navigation_Message()
{
reset();
}
void Gps_Navigation_Message::print_gps_word_bytes(unsigned int GPS_word)
{
std::cout << " Word =";
std::cout << std::bitset<32>(GPS_word);
std::cout << std::endl;
}
bool Gps_Navigation_Message::read_navigation_bool(std::bitset bits, const bits_slice *slices)
{
bool value;
if (bits[GPS_SUBFRAME_BITS - slices[0].position] == 1)
{
value = true;
}
else
{
value = false;
}
return value;
}
unsigned long int Gps_Navigation_Message::read_navigation_unsigned(std::bitset bits, const bits_slice *slices, int num_of_slices)
{
unsigned long int value = 0;
for (int i=0; i bits, const bits_slice *slices, int num_of_slices)
{
signed long int value = 0;
// Discriminate between 64 bits and 32 bits compiler
int long_int_size_bytes = sizeof(signed long int);
if (long_int_size_bytes == 8) // if a long int takes 8 bytes, we are in a 64 bits system
{
// read the MSB and perform the sign extension
if (bits[GPS_SUBFRAME_BITS - slices[0].position] == 1)
{
value^=0xFFFFFFFFFFFFFFFF; //64 bits variable
}
else
{
value&=0;
}
for (int i=0; i half_week)
{
corrTime = time - 2*half_week;
}
else if (time < -half_week)
{
corrTime = time + 2*half_week;
}
return corrTime;
}
// 20.3.3.3.3.1 User Algorithm for SV Clock Correction.
double Gps_Navigation_Message::sv_clock_correction(double transmitTime)
{
double dt;
dt = check_t(transmitTime - d_Toc);
d_satClkCorr = (d_A_f2 * dt + d_A_f1) * dt + d_A_f0 + d_dtr;
double correctedTime = transmitTime - d_satClkCorr;
return correctedTime;
}
void Gps_Navigation_Message::satellitePosition(double transmitTime)
{
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
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*GPS_PI), (2*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*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);
// 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
phi = fmod((phi), (2*GPS_PI));
// Correct argument of latitude
u = phi + d_Cuc * cos(2*phi) + d_Cus * sin(2*phi);
// Correct radius
r = a * (1 - d_e_eccentricity*cos(E)) + d_Crc * cos(2*phi) + d_Crs * sin(2*phi);
// Correct inclination
i = d_i_0 + d_IDOT * tk + d_Cic * cos(2*phi) + d_Cis * sin(2*phi);
// Compute the angle between the ascending node and the Greenwich meridian
Omega = d_OMEGA0 + (d_OMEGA_DOT - OMEGA_EARTH_DOT)*tk - OMEGA_EARTH_DOT * d_Toe;
// Reduce to between 0 and 2*pi rad
Omega = fmod((Omega + 2*GPS_PI), (2*GPS_PI));
// --- 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
double Omega_dot = d_OMEGA_DOT - OMEGA_EARTH_DOT;
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);
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);
}
int Gps_Navigation_Message::subframe_decoder(char *subframe)
{
int subframe_ID = 0;
int SV_data_ID = 0;
int SV_page = 0;
//double tmp_TOW;
unsigned int gps_word;
// UNPACK BYTES TO BITS AND REMOVE THE CRC REDUNDANCE
std::bitset subframe_bits;
std::bitset word_bits;
for (int i=0; i<10; i++)
{
memcpy(&gps_word, &subframe[i*4], sizeof(char)*4);
word_bits = std::bitset<(GPS_WORD_BITS+2)>(gps_word);
for (int j=0; j= 0) // is not in the past
{
//Detect if the effectivity time and user's time is within six hours = 6 * 60 *60 = 21600 s
int secondOfLeapSecondEvent = i_DN * 24 * 60 * 60;
if (weeksToLeapSecondEvent > 0)
{
t_utc_daytime = fmod(gpstime_corrected - Delta_t_UTC, 86400);
}
else //we are in the same week than the leap second event
{
if (abs(gpstime_corrected - secondOfLeapSecondEvent) > 21600)
{
/* 20.3.3.5.2.4a
* Whenever the effectivity time indicated by the WN_LSF and the DN values
* is not in the past (relative to the user's present time), and the user's
* present time does not fall in the time span which starts at six hours prior
* to the effectivity time and ends at six hours after the effectivity time,
* the UTC/GPS-time relationship is given by
*/
t_utc_daytime = fmod(gpstime_corrected - Delta_t_UTC, 86400);
}
else
{
/* 20.3.3.5.2.4b
* Whenever the user's current time falls within the time span of six hours
* prior to the effectivity time to six hours after the effectivity time,
* proper accommodation of the leap second event with a possible week number
* transition is provided by the following expression for UTC:
*/
int W = fmod(gpstime_corrected - Delta_t_UTC - 43200, 86400) + 43200;
t_utc_daytime = fmod(W, 86400 + d_DeltaT_LSF - d_DeltaT_LS);
//implement something to handle a leap second event!
}
if ( (gpstime_corrected - secondOfLeapSecondEvent) > 21600)
{
Delta_t_UTC = d_DeltaT_LSF + d_A0 + d_A1 * (gpstime_corrected - d_t_OT + 604800*(double)(i_GPS_week - i_WN_T));
t_utc_daytime = fmod(gpstime_corrected - Delta_t_UTC, 86400);
}
}
}
else // the effectivity time is in the past
{
/* 20.3.3.5.2.4c
* Whenever the effectivity time of the leap second event, as indicated by the
* WNLSF and DN values, is in the "past" (relative to the user's current time),
* and the user�s current time does not fall in the time span as given above
* in 20.3.3.5.2.4b,*/
Delta_t_UTC = d_DeltaT_LSF + d_A0 + d_A1 * (gpstime_corrected - d_t_OT + 604800 * (double)(i_GPS_week - i_WN_T));
t_utc_daytime = fmod(gpstime_corrected - Delta_t_UTC, 86400);
}
double secondsOfWeekBeforeToday = 43200 * floor(gpstime_corrected / 43200);
t_utc = secondsOfWeekBeforeToday + t_utc_daytime;
return t_utc;
}
bool Gps_Navigation_Message::satellite_validation()
{
bool flag_data_valid = false;
b_valid_ephemeris_set_flag = false;
// First Step:
// check Issue Of Ephemeris Data (IODE IODC..) to find a possible interrupted reception
// and check if the data have been filled (!=0)
if (d_TOW_SF1 != 0 and d_TOW_SF2 != 0 and d_TOW_SF3 != 0)
{
if (d_IODE_SF2 == d_IODE_SF3 and d_IODC == d_IODE_SF2 and d_IODC!=-1)
{
flag_data_valid = true;
b_valid_ephemeris_set_flag=true;
}
}
return flag_data_valid;
}