ogr2rnx.c

The OpenGPSRec receiver software runs on the real time operating system RTAI-Linux. It compiles with

/********************************************************************
** @source
**
** @author Copyright (C) 2000,2001 Antonio Tabernero
** @version 1.48
** @last modified Ago 28 2001 Antonio Tabernero. 
** @last modified Mar 22 2002 Antonio Tabernero. 
** @last modified Apr  3 2002 Antonio Tabernero. 
** @@
**
** This program is free software; you can redistribute it and/or
** modify it under the terms of the GNU Library General Public
** License as published by the Free Software Foundation; either
** version 2 of the License, or (at your option) any later version.
**
** This program 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
** Library General Public License for more details.
**
** You should have received a copy of the GNU Library General Public
** License along with this library; if not, write to the
** Free Software Foundation, Inc., 59 Temple Place - Suite 330,
** Boston, MA  02111-1307, USA.
********************************************************************/

/*
 * Changes:
 *
 *  18.07.03 - dump at 1 Hz 
 *
 */


#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#include <time.h>
#include <sys/types.h>

#define VERSION 1.48

/* ******************************************************************
1.3  Can use records 0x0e and 0x11 to get aprox time and position 
     instead of requiring record 0x33 (which seems to be missing 
     in some models/firmwares)

     Output Doppler data (using the flag +doppler) when applied to 
     binary files logged with the 1.2 (or newer) version of async.

1.4  Added generation of RINEX ephemeris files (-nav option)
     Added decoding of the navigation message
     Added option to get input from stdin
     Some bugs corrected

1.45 -monitor, -V, -sf  options
1.48 version published. Minor cosmetic changes from 1.45
1.49 run source thru indent -bli0 --no-tabs
     generate_nav: check parity of all WORDs before writing nav file
     patch file header with ionospheric and UTC data (option -f only)
******************************************************************* */


#define AS_BYTE   0
#define AS_CHAR   1
#define AS_UINT   2
#define AS_ULONG  3
#define AS_FLOAT  4
#define AS_DOUBLE 5

#define NEVER_USED 0
#define NO_DUMP    1
#define DUMP       2
#define DUMP_PHASE_ONLY 3

#define c 2.99792458e8          // WGS84 Speed of light
#define L1 1575420000           // L1 frequency
#define lambda (c/L1)           // L1 wavelength
#define WGS84_PI 3.1415926535898



typedef unsigned char BOOLEAN;
typedef unsigned char BYTE;
typedef unsigned long ULONG;
typedef unsigned short int UINT;
typedef short int INT;

BOOLEAN VC_format = 1;
BOOLEAN STDIN = 0;

BOOLEAN ETREX;
BOOLEAN RESET_CLOCK;
BYTE ONLY_STATS, SELECTED_SV, ONE_SAT, DIF_RECORDS;
int SELECTED_SF, SELECTED_PAGE;
BYTE VERBOSE, VERBOSE_NAV;
BYTE NAV_GENERATION,
  MONITOR_NAV, PARSE_RECORDS, RINEX_GENERATION, VERIFY_TIME_TAGS, NO_SNR;
BYTE OPT1, RELAX;

double USER_XYZ[3];
double IONO[8], UTC_A0, UTC_A1,
       UTC_DT_LS, UTC_DT_LSF;
ULONG  UTC_TOT, UTC_WN_T, UTC_WN_LSF, UTC_DAY_NO;

int USER_DATE[3];
BOOLEAN GIVEN_XYZ, GIVEN_DATE;


int L1_FACTOR, GET_THIS;
int SELECTED_RECORDS[16];
BYTE RINEX_FILE;
BYTE OBS_MASK, N_OBS;
UINT INTERVAL;
long START, LAST, ELAPSED;

char DATAFILE[60];
char COMMAND_LINE[256];
char location[6];
char marker[32];


BYTE LSB[32] = { 7, 6, 5, 4, 3, 2, 1, 0,
  15, 14, 13, 12, 11, 10, 9, 8,
  23, 22, 21, 20, 19, 18, 17, 16,
  31, 30, 29, 28, 27, 26, 25, 24
};

ULONG NAV_WORD;

// data bits: bits 30 -> 7 within ULONG NAV_WORD
#define d1  ( (NAV_WORD>>29) & 1 )
#define d2  ( (NAV_WORD>>28) & 1 )

#define d3  ( (NAV_WORD>>27) & 1 )
#define d4  ( (NAV_WORD>>26) & 1 )
#define d5  ( (NAV_WORD>>25) & 1 )
#define d6  ( (NAV_WORD>>24) & 1 )
#define d7  ( (NAV_WORD>>23) & 1 )
#define d8  ( (NAV_WORD>>22) & 1 )
#define d9  ( (NAV_WORD>>21) & 1 )
#define d10  ( (NAV_WORD>>20) & 1 )

#define d11  ( (NAV_WORD>>19) & 1 )
#define d12  ( (NAV_WORD>>18) & 1 )
#define d13  ( (NAV_WORD>>17) & 1 )
#define d14  ( (NAV_WORD>>16) & 1 )
#define d15  ( (NAV_WORD>>15) & 1 )
#define d16  ( (NAV_WORD>>14) & 1 )
#define d17  ( (NAV_WORD>>13) & 1 )
#define d18  ( (NAV_WORD>>12) & 1 )
#define d19  ( (NAV_WORD>>11) & 1 )
#define d20  ( (NAV_WORD>>10) & 1 )
#define d21  ( (NAV_WORD>>9) & 1 )
#define d22  ( (NAV_WORD>>8) & 1 )
#define d23  ( (NAV_WORD>>7) & 1 )
#define d24  ( (NAV_WORD>>6) & 1 )

// Parity bits: bits 6 -> 1 within ULONG NAV_WORD
#define D25  ( (NAV_WORD>>5) & 1 )
#define D26  ( (NAV_WORD>>4) & 1 )
#define D27  ( (NAV_WORD>>3) & 1 )
#define D28  ( (NAV_WORD>>2) & 1 )
#define D29  ( (NAV_WORD>>1) & 1 )
#define D30  ( (NAV_WORD>>0) & 1 )

// Previous word parity bits 29-30: bits 32 and 31 of ULONG NAV_WORD
#define P29  ( (NAV_WORD>>31) & 1 )
#define P30  ( (NAV_WORD>>30) & 1 )


BYTE OK[3] = { 1, 1, 1 };
BYTE RESET[30] =
  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0,
  0, 0, 0, 0, 0
};

BYTE frame[32][30];
BOOLEAN check_frame[32][30];

int pages[64] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
  16, 17, 18, 19, 20, 21, 22, 23, 24, 2, 3, 4, 5, 7, 8, 9,
  10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
  -1, -1, -1, 25, 13, 14, 15, 17, 18, 1, 19, 20, 22, 23, 12, 25
};

// //// STRUCTS of RECORDS with (more or less) known fields

typedef struct
{
  double llh[3];
  double xyz[3];
} type_rec0x11;

typedef struct
{
  int week;
  int garmin_wdays;
  double tow;
} type_rec0x0e;

typedef struct
{
  float delta_pr;
  float f1;
  double pr;
  float f2;
  BYTE sv;
} type_rec0x16;

typedef struct
{
  BYTE sv;
  BYTE elev;
  UINT frac;
  UINT db;
  BYTE tracked;
  BYTE flag;
} type_rec0x1a;

typedef struct
{
  float altitud;
  float epe[3];
  UINT fix;
  double tow;
  double pos[2];
  float vel[3];
  float h_ellip;
  UINT leap;
  ULONG wdays;
} type_rec0x33;


typedef struct
{
  ULONG c50;
  BYTE uk[4];
  BYTE sv;
} type_rec0x36;


typedef struct
{
  ULONG c_phase;
  long int tracked;             // BYTE tracked; UINT cc2; BYTE flag; 
  UINT delta_f;
  ULONG int_phase;
  double pr;
  ULONG c511;
  UINT db;
  double tow;
  BYTE sv;
} type_rec0x38;


typedef struct
{
  BYTE used;
  UINT db;
  double prange;
  double phase;
  float doppler;
  BYTE elev;
  long int tracked;
  ULONG c511;
  float last36;
} rinex_obs;


// Global variables used to study rates of change, etc.

type_rec0x38 rec38, last38;
float mean, sigma, mean_q;
float phase_rate;
float mean2 = 0;
float sigma2 = 0;
int n;
int n2 = 0;


void llh2xyz (double *, double *);


type_rec0x11
process_0x11 (BYTE * record)    // Position
{
  type_rec0x11 rec;

  memcpy (&rec.llh[0], record, 8);
  memcpy (&rec.llh[1], record + 8, 8);


  rec.llh[0] *= (180 / WGS84_PI);
  rec.llh[1] *= (180 / WGS84_PI);
  rec.llh[2] = 0.0;
  llh2xyz (rec.llh, rec.xyz);

  if (VERBOSE)
  {
    printf ("Record 0x11 :: "),
      printf ("Lat %.4f  Long %.4f\n", rec.llh[0], rec.llh[1]);
    printf ("XYZ(h=0.0) = [ %.0f %.0f %.0f]\n", rec.xyz[0], rec.xyz[1],
            rec.xyz[2]);
  }

  return rec;
}

double
julday (year, month, day, hour, min, sec)
     UINT year, month, day, hour, min;
     double sec;
{
  double jd;

  if (month <= 2)
  {
    year -= 1;
    month += 12;
  }
  jd = floor (365.25 * year) + floor (30.6001 * (month + 1)) + day;
  jd += (hour + (min + sec / 60.0) / 60.0) / 24 + 1720981.5;
  return jd;
}

double
gps_time (double julday, int *week)
{
  double tow, days_gps;
  int wn;
  // int weekday;

  // weekday=floor(julday+0.5)%7 + 1 ; // 0 Sunday, 1 Monday, ...

  // Correct at least for the 13 leap seconds as of writting this 
  julday += (13.0 / 86400.0);

  days_gps = (julday - 2444244.5);
  wn = floor (days_gps / 7);

  tow = (days_gps - wn * 7) * 86400;

  *week = wn;

  return tow;
}

type_rec0x0e
process_0x0e (BYTE * record)
{
  type_rec0x0e rec;
  UINT month, day, year, hour, min, sec;
  double jd;

  /*
   * typedef struct {BYTE month;BYTE day;UINT year;UINT hour;BYTE
   * min;BYTE sec;} UTC; UTC utc; utc.month=record[0];
   * utc.day=record[1]; utc.year=*((short int*)(record+2));
   * utc.hour=*((short int*)(record+4)); utc.min=record[6];
   * utc.sec=record[7]; 
   */

  month = (UINT) record[0];
  day = (UINT) record[1];
  year = *((short int *) (record + 2));
  hour = *((short int *) (record + 4));
  min = (UINT) record[6];
  sec = (UINT) record[7];

  jd = julday (year, month, day, hour, min, (double) sec);
  rec.tow = gps_time (jd, &rec.week);
  rec.garmin_wdays = (rec.week - 521) * 7;


  if (VERBOSE)
  {
    printf ("%02u/%02u/%4d ", day, month, year);
    printf ("%02d:%02d:%02d -> JD %f -> ", hour, min, sec, jd);
    printf ("GPS Week %4d ToW %6.0f sec  Garmin Weekdays %d\n",
            rec.week, rec.tow, rec.garmin_wdays);
  }

  return rec;
}

void
process_0x12_gps38 (BYTE * record)
{
  printf ("PRN %2d ", record[3]);
  printf ("TOW %14.7f ", *((double *) (record + 18)));
  printf ("Pseudorange %14.3f\n", *((double *) (record + 26)));
  return;
}




void
process_0x14 (BYTE * record)
{

  printf ("Fix %d %d\n", *((UINT *) record), *((UINT *) (record + 2)));
  printf ("TOW4 %.7f  TOW12 %.7f\n", *((double *) (record + 4)),
          *((double *) (record + 12)));

  printf ("Posicion %.4f %.4f\n",
          *((double *) (record + 28)) * 180 / WGS84_PI,
          *((double *) (record + 36)) * 180 / WGS84_PI);

  printf ("Vel  %f %f %f\n", *((float *) (record + 48)),
          *((float *) (record + 52)), *((float *) (record + 56)));
}



type_rec0x36
process_0x36 (BYTE * record)
{
  type_rec0x36 rec;
  int i;

  rec.c50 = *((ULONG *) (record));      // 50 Hz counter
  for (i = 0; i < 4; i++)
    rec.uk[i] = record[4 + i];  // 30 bit nav word

  rec.sv = record[8];           // PRN-1

  /*
   * memcpy(&rec.c50,record,4); memcpy(&rec.uk,record+4,4);
   * memcpy(&rec.sv,record+8,1); 
   */

  if (ONE_SAT && rec.sv != SELECTED_SV)
    return rec;

  if (VERBOSE)
  {
    printf
      ("0x36 ----------------------------------------------------------------\n");
    printf ("PRN %02d: C (50 Hz) %8lu = TOW %8.1f ", rec.sv + 1,
            rec.c50, rec.c50 / 50.0);
    for (i = 0; i < 4; i++)
      printf ("%02x ", rec.uk[i]);
    printf ("\n");
  }

  return rec;
}


type_rec0x33
process_0x33 (BYTE * record)
{
  type_rec0x33 rec;

  memcpy (&rec.altitud, record, 4);
  record += 4;
  memcpy (&rec.epe, record, 12);
  record += 12;
  memcpy (&rec.fix, record, 2);
  record += 2;
  memcpy (&rec.tow, record, 8);
  record += 8;
  memcpy (&rec.pos, record, 16);
  record += 16;
  memcpy (&rec.vel, record, 12);
  record += 12;
  memcpy (&rec.h_ellip, record, 4);
  record += 4;
  memcpy (&rec.leap, record, 2);
  record += 2;
  memcpy (&rec.wdays, record, 4);
  record += 4;


  if (VERBOSE)
  {
    printf
      ("0x33 ----------------------------------------------------------------\n");
    printf ("Pos:  (N,E,H )  (%f %f %f)  H_ellip=%.0f  FIX %dD\n",
            rec.pos[0] * 180 / WGS84_PI, rec.pos[1] * 180 / WGS84_PI,
            rec.altitud, rec.h_ellip, rec.fix);
    printf ("Vel:  (N,E,Up)  (%f %f %f)  ", rec.vel[0], rec.vel[1],
            rec.vel[2]);
    printf ("Epe  (dX,dY,dH) (%.1f %.1f %.1f)\n", rec.epe[0],
            rec.epe[1], rec.epe[2]);
    printf ("Wdays %ld  TOW %13.6f  Leap %d\n", rec.wdays, rec.tow, rec.leap);
  }

  return rec;
}


void
process_0x37 (BYTE * record)
{
  double pr, tow;
  ULONG c511;
  BYTE svid, unknown[4];
  UINT cc, countdown, c1, c2;
  float TR;

  memcpy (&c1, record, 2);
  memcpy (unknown, record + 2, 2);
  memcpy (&c2, record + 4, 2);
  memcpy (&cc, record + 6, 2);
  memcpy (&countdown, record + 8, 2);
  TR = ((float) 1024.0 - countdown) / 16.0;
  memcpy (&pr, record + 10, 8);
  memcpy (&c511, record + 18, 4);
  memcpy (&tow, record + 22, 8);
  memcpy (unknown + 2, record + 30, 2);
  svid = record[32];


  if (ONE_SAT && svid != SELECTED_SV)
    return;

  printf
    ("0x37 ----------------------------------------------------------------\n");
  printf ("PRN %02d: PR %10.1f TOW %10.3f (c511 %10lu) CC %5d\n",
          svid + 1, pr, tow, c511, cc);
  printf ("       Countdown %5d (%2.0f'') C1 %04x C2 %04x\n", countdown,
          TR, c1, c2);
  printf ("       UNKNOWN: ");
  printf ("%02x %02x | %02x %02x\n", unknown[0], unknown[1], unknown[2],
          unknown[3]);
}

void
get_0x1a_info (BYTE * record, type_rec0x1a chan[])
{
  BYTE *ptr;
  int k;

  for (k = 0, ptr = record; k < 12; k++, ptr += 8)
  {
    if (ETREX)
    {
      chan[k].sv = ptr[5];
      chan[k].elev = ptr[4];
      chan[k].frac = *((UINT *) (ptr + 0));
      chan[k].db = *((UINT *) (ptr + 2));
      chan[k].tracked = ptr[6];
      chan[k].flag = ptr[7];
    }
    else
    {
      chan[k].sv = ptr[0];
      chan[k].elev = ptr[1];
      chan[k].frac = *((UINT *) (ptr + 2));
      chan[k].db = *((UINT *) (ptr + 4));
      chan[k].tracked = ptr[6];
      chan[k].flag = ptr[7];
    }
  }

  return;
}

void
process_0x1a (BYTE * record)
{
  int k;
  BYTE *rec;
  BYTE svid, elev, trk, flag;
  unsigned int sq;
  float frac_phase;


  if (!ONE_SAT)
  {
    printf
      ("0x1a ---------------------------------------------------------------\n");
    printf ("Chann -> PRN | Elev | Q Sgnl |  Phase  | Track | Flag |\n");
  }

  for (k = 0, rec = record; k < 12; k++, rec += 8)
  {
    if (ETREX)
    {
      svid = rec[5];
      elev = rec[4];
      sq = *((UINT *) (rec + 2));
      frac_phase = *((UINT *) (rec + 0)) / 2048.0;
      trk = rec[6];
      flag = rec[7];
    }
    else
    {
      svid = rec[0];
      elev = rec[1];
      sq = *((UINT *) (rec + 4));
      frac_phase = *((UINT *) (rec + 2)) / 2048.0;