Here we go through the code that reads RINEX format Navigation and Observation files from a CORS location and calculates the User position for all the pseudorange measurements in the file.
This code uses the same routines to calculate the user position as the GPS Hardware receiver the main difference is that it uses pseudorange and ADR data from a CORS file and runs on a Linux PC and not directly in a GPS receiver.
The program needs to be supplied the name and location of an Observation file it assumes there is an associated Navigation file in the same directory with the same name but a different file extension, ending in 'n' instead of 'o'.
The program begins in the standard manner by defining include files that contain definitions and functions used by the body of the program, followed by the definition of memory structures used to hold the GPS variables required by the program.
/********************************************************************************
* C routines for calculating positions of satellites
* at a given time of day.Starting with satellite frame data.
*
*******************************************************************************/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
#include "gnss_types.h"
#include "gnss_gps.h"
#include "fileio.h"
#include "position.h"
#include "time.h"
#include "List.h"
#include "constants.h"
GPS_structEphemeris ephs[N_SATELLITES];
GPS_structAlmanac almanacs[N_SATELLITES];
GNSS_structMeasurement rx_meas[N_SATELLITES];
GNSS_structRxPVT rx_pos;
GNSS_structKlobuchar iono_data;
GPS_structUTC utc_conv;
time_t gps_start_time;
int gps_week,gps_time=0;
const double L1_LAMBDA = (double)SPEED_OF_LIGHT/(double)(154*10.23E6);
/******************************************************************************
*
* Main Routine to calculate positions of satellites at the current time
* with Ephemeris data read from a Rinex Navigation file downloaded from
* the internet or generated by a GPS Rx
* and Measurement data read from a Rinex Observation file downloaded from
* the internet or generated by a GPS Rx.
*
*******************************************************************************
*/
int main (int argc, char **argv)
{
struct tm *gmt;
time_t gmt_time_secs;
int satno,prn,count,chan,i;
int sats[N_CHANNELS];
RinexOH_ rnxoh;
RinexOD_ rnxod;
GPS_structEphemeris eph,*eph_ptr,*next_eph;
GNSS_structSatPVT sat_pos;
GNSS_structPVT *spvt;
List_ *eph_list[N_SATELLITES];
FILE *fdata;
char data_file[256];
int gps_wday,no_sats;
double user_pos[4];
double user_xyz[4];
double theta;
double phi;
double horz_error,avg_horz_error,avg_vert_error;
double pos_error,avg_pos_error,no_obs;
double max_pos_error=0.0;
double min_pos_error=10000.0;
The program begins by checking that a filename has been supplied for the RINEX files to be read and then proceeeds to initialise the time and date structures needed to hold the time for the calculations, then it set the position back to zero.
if(argc != 2)
{
fprintf(stderr,"Correct usage is eph_user_pos rinex09o\n");
exit(0);
}
setenv("TZ","GMT",1);
tzset();
time(&gmt_time_secs);
gmt = gmtime(&gmt_time_secs);
init_gps_start_time();
// printf("GPS start time is %.26s",asctime(gmtime(&gps_start_time)));
for(i=0;i<4;i++)
{
user_pos[i]=0.0;
user_xyz[i]=0.0;
}
strcpy(data_file,argv[1]);
if((fdata=fopen(data_file,"rb"))==NULL)
{
fprintf(stderr,"Cannot open file %s \n",data_file);
exit(1);
}
count = rinexO_hdr(fdata,&rnxoh);
fclose(fdata);
gps_time = difftime(rnxoh.time_1st_obs,gps_start_time);
gps_week = gps_time/604800;
gps_time %= 604800;
gmt_time_secs = rnxoh.time_1st_obs+43200;
gmt = gmtime(&gmt_time_secs);
gps_wday = gmt->tm_wday;
gmt_time_secs = rnxoh.time_1st_obs;
printf("GPS time is %.26s",asctime(gmtime(&gmt_time_secs)));
data_file[strlen(data_file)-1] = 'n';
printf("Navigation file is %s \n",data_file);
if((fdata=fopen(data_file,"rb"))==NULL)
{
fprintf(stderr,"Cannot open file %s \n",data_file);
exit(1);
}
for(satno=0;satno<N_SATELLITES;satno++)
eph_list[satno] = NULL;
count = rinexN_hdr(fdata, &iono_data,&utc_conv);
printf("UTC GPS week is %d\n",count);
To begin with Navigation and Observation file headers are read and the GPS start time of the Observations derived.
while(feof(fdata) == 0)
{
prn = rinexN_dtl(fdata, &eph);
if(prn > 0 && prn <= 32)
{
// printf("Sat %d Ephemeris is %.24s toe is %lf\n",
// prn,asctime(gmtime(&eph.toc)),eph.toe);
eph_ptr = calloc(1,sizeof(GPS_structEphemeris));
*eph_ptr = eph;
satno = prn - 1;
eph_list[satno] = AddToList(eph_list[satno],eph_ptr);
}
}
fclose(fdata);
for(satno=0;satno<N_SATELLITES;satno++)
{
eph_ptr = getNext(eph_list[satno]);
if(eph_ptr != NULL) ephs[satno] = *eph_ptr;
}
Now all the Ephemeris parameters in the Navigation file are read and loaded into memory. Then the structures that will hold the Ephemeris parameters used in the actual user position calculation are initialised.
strcpy(data_file,argv[1]);
if((fdata=fopen(data_file,"rb"))==NULL)
{
fprintf(stderr,"Cannot open file %s \n",data_file);
exit(2);
}
count = rinexO_hdr(fdata, &rnxoh);
printf("GPS week is %d\n",count);
rx_pos.pvt.clk = rnxoh.upos[0];
rx_pos.pvt.x = rnxoh.upos[1];
rx_pos.pvt.y = rnxoh.upos[2];
rx_pos.pvt.z = rnxoh.upos[3];
rx_pos.isKalValid = 0;
rx_pos.isValid = 1;
for(chan=0;chan<N_CHANNELS;chan++)
rx_meas[chan].flags.isAdrValid = 0;
while(feof(fdata) == 0)
{
no_sats = rinexO_dtl(fdata, &rnxod, &rnxoh);
printf("\nNo of Sats observed is %d Observation time %.26s TOW %d\n",
no_sats,asctime(gmtime(&rnxod.too)),rnxod.gps_time);
//
// Single Frequency un-smoothed pseudorange
//
if(no_sats >= 4)
{
for(chan=0,no_sats=0;chan<rnxod.no_sats;chan++)
{
prn = rnxod.sv[chan];
satno = prn - 1;
eph_ptr = &ephs[chan];
if(eph_ptr == NULL || prn != eph_ptr->prn)
eph_ptr = viewCurr(eph_list[satno]);
next_eph = viewNext(eph_list[satno]);
while(fabs(difftime(eph_ptr->toe,rnxod.gps_time)) >
fabs(difftime(next_eph->toe,rnxod.gps_time)))
{
eph_ptr = getNext(eph_list[satno]);
next_eph = viewNext(eph_list[satno]);
if(next_eph != NULL)
{
printf("Sat %d Ephemeris toe is %d Next toe is %d too is %d\n",
prn,eph_ptr->toe,next_eph->toe,rnxod.gps_time);
}
else
{
printf("Reached end of Ephemeris List %d\n",prn);
}
}
if(eph_ptr != NULL)
{
ephs[chan] = *eph_ptr;
eph_sat_pos(chan,rnxod.gps_time+rnxod.clk_bias);
// printf("Chan %d Sat %2.2d position is x=%lf,y=%lf,z=%lf\n",
// chan,prn,rx_meas[chan].sat_pos.pvt.x,rx_meas[chan].sat_pos.pvt.y,
// rx_meas[chan].sat_pos.pvt.z);fflush(stdout);
theta = sat_el_angle(&rx_meas[chan].sat_pos);
phi = sat_az_angle(&rx_meas[chan].sat_pos);
// printf("Satellite %d elevation is %lf azimuth is %lf\n",
// ephs[chan].prn,theta*180.0/PI,phi*180.0/PI);
if(rnxod.lli[chan][rnxoh.C1]&0x3>0
||rnxod.ss[chan][rnxoh.C1]<5||ephs[chan].health>0)
{
// printf("Sat %d lli %d ss %d health %d\n",prn,rnxod.lli[chan],rnxod.ss[chan],ephs[chan].health);
}
// else
// {
rx_meas[chan].id = prn;
rx_meas[chan].psr = rnxod.C1[chan];
rx_meas[chan].psr_smoothed = rnxod.C1[chan];
printf("Sat %d pseudorange is %lf\n",prn,rnxod.C1[chan]);
rx_meas[chan].tow = (double)rnxod.gps_time + rnxod.clk_bias;
rx_meas[chan].corrections.prcIono = 0.0;
rx_meas[chan].corrections.prcTropoDry = 0.0;
if(rnxod.C1[chan] > 20000000.0) sats[no_sats++] = chan;
// }
rx_meas[chan].channel = chan;
// printf("Sat %d Ephemeris is %.26s",prn,asctime(gmtime(&sats[satno].eph->toc)));
}
}
//
// Looped over Observation channel next is position calculation
//
rnxod.no_sats = no_sats;
// printf("Ephemeris is %.26s",asctime(gmtime(&sats[satno].eph.toc)));
if(no_sats >= 4)
{
printf("User position is %lf,%lf,%lf,%lf\n",
rnxoh.upos[0],rnxoh.upos[1],rnxoh.upos[2],rnxoh.upos[3]);
rx_pos.isValid = 0;
rx_pos.pvt.clk = 0.0;
rx_pos.pvt.x = 0.0;
rx_pos.pvt.y = 0.0;
rx_pos.pvt.z = 0.0;
pos_calc_lsqf(no_sats,sats);
user_xyz[0] = rx_pos.pvt.clk;
user_xyz[1] = rx_pos.pvt.x;
user_xyz[2] = rx_pos.pvt.y;
user_xyz[3] = rx_pos.pvt.z;
printf("Calc position is %lf,%lf,%lf,%lf \n",
rx_pos.pvt.clk,rx_pos.pvt.x,rx_pos.pvt.y,rx_pos.pvt.z);
pos_error = sqrt(vdist(user_xyz,rnxoh.upos));
printf("Position error is %lf no Sats used is %d\n",
pos_error,no_sats);fflush(stdout);
}
//
// Calculate error for current position calculation
//
if(no_sats >= 4)
{
for(i=0;i<4;i++)user_pos[i] += user_xyz[i];
avg_pos_error += pos_error;
no_obs ++;
// hvdelta(rnxoh.upos,user_xyz);
// horz_error = user_xyz[1]*user_xyz[1]+user_xyz[2]*user_xyz[2];
// horz_error = sqrt(horz_error);
// avg_horz_error += horz_error;
// avg_vert_error += user_xyz[3];
if(min_pos_error &;gt; pos_error)min_pos_error = pos_error;
if(max_pos_error &;lt; pos_error)max_pos_error = pos_error;
} // Calc position if enough valid Sats
} // Don't do anything if less than 4 sats read from file this time.
} // End of File read loop
for(i=0;i<4;i++) user_pos[i] /= no_obs;
printf("\nAvg error's are posn %lf horizontal %lf vertical %lf\n",
avg_pos_error/no_obs,avg_horz_error/no_obs,avg_vert_error/no_obs);
printf("Averaged calc posn is %lf,%lf,%lf,%lf\n",
user_pos[0],user_pos[1],user_pos[2],user_pos[3]);
printf("Actual user posn is %lf,%lf,%lf,%lf\n",
rnxoh.upos[0],rnxoh.upos[1],rnxoh.upos[2],rnxoh.upos[3]);
// ecef_to_llh(rnxoh.upos);
// printf("User position is %lf, %lf, %lf\n",
// rnxoh.upos[1]*180.0/PI,rnxoh.upos[2]*180.0/PI,rnxoh.upos[3]);
// wgs_to_nad(user_pos);
// printf("Averaged NAD calc posn is %lf %lf %lf\n",
// user_pos[1],user_pos[2],user_pos[3]);
printf("Min pos error is %lf Max pos error is %lf\n",
min_pos_error,max_pos_error);
ecef_to_llh(user_pos);
printf("Averaged calc posn is %lf %lf %lf\n",
user_pos[1]*180.0/PI,user_pos[2]*180.0/PI,user_pos[3]);
fclose(fdata);
exit(0);
}
Back to Top of Page