/*
BJfmtl_rnx.c
*/

/*
This code is largely extracted from Larry Romans's code BJfmt.c

It should conform to vanilla ANSI C.

It is compiled on HP-UX with cc -Aa -o BJfmtl_rnx BJfmtl_rnx.c -lm

I intend to continue developing and maintaining it in the foreseeable future, and would very much appreciate hearing of any problems
 or suggestions.

- Byron

Dr. Byron A. Iijima         Ionospheric and Atmospheric Remote Sensing Group
Jet Propulsion Laboratory M/S 138-308                 Phone: +1 818 354 8503
4800 Oak Grove Drive                                    Fax: +1 818 393 5115
Pasadena, CA 91109  USA                   Email: Byron.A.Iijima@jpl.nasa.gov

*/

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

#define MAXLINE 512
#define dnint(A) ( floor((A)+.5) )
#define J2000 946728000
#define loop(I,N) for(I=0;I<N;I++)
#define ARG_IS(A) (!strcmp(arg, A))

void print_rnx(void);
void print_header(void);
int snr2n(short);
void cmdline(void);
void argprob(char *);
void do_str(char *, int);
void do_int(int *);
void do_uint(unsigned int *);
void do_short(short *);
void do_real(double *);
void do_vec(double *);
void do_help(void);

int obsdcmp(int *,int *);

typedef struct obsd_str {
unsigned int t;
int prn;
short snr[3];
double ph[3];
double pr[3];
int flag; /* 1=print, 2=discard, 3=save */
} obsd_str;

obsd_str obsd[20];
int nbuf,ibuf;
unsigned int toutlast = 0;

int NARGS, iarg;
char **ARGV, *arg, *pname;

unsigned char q = 0;

/* for header */

char satname[9];
char pgm[21], runby[21], date[21],
  marker[61], observer[21], agency[41],
  recnum[21], rectype[21], recvers[21],
  antnum[21], anttype[21];
double pos[3], ant[3];

/* other control parameters */

int delconflict=1;
unsigned int use_samprate=10;
double p1cutfac=-1.;
short p1cut=-1;
short p2cut=-1;
int noedit, do_snrs;

struct tm *pT;

/**********************************************************************************************/

main(int argc, char *argv[]) {

unsigned int antid,obstype,samprate,ch0,ch1,ch2;
unsigned int tlast=0;
int dlength;
char line[MAXLINE];
obsd_str obsd1;

struct tm *pTT;
time_t tp;
tp=time(NULL);
pTT = gmtime(&tp);
sprintf(date, "%04d-%02d-%02d %02d:%02d:%02.0f",
      pTT->tm_year + 1900, pTT->tm_mon + 1, pTT->tm_mday,
      pTT->tm_hour, pTT->tm_min, (double)pTT->tm_sec);

NARGS = argc;
ARGV = argv;
cmdline();

nbuf=20;
ibuf=0;

while(fgets(line,MAXLINE,stdin)) {
  if(q) printf("Hello,%s",line);
/*
40095820 OBSD prn22 03 83 0a 57 CA: 00 297 8978434.728137 23355782.575862  P2: 02 55 6996152.443221 23355791.037362  P1: 01 57 8978434.491994 23355781.606868 
*/
  if(sscanf(line,"%u OBSD prn%d %x %x %x %d CA: %x %hd %lf %lf P2: %x %hd %lf %lf P1: %x %hd %lf %lf",
   &obsd1.t,&obsd1.prn,&antid,&obstype,&samprate,&dlength,
   &ch0,&obsd1.snr[0],&obsd1.ph[0],&obsd1.pr[0],
   &ch1,&obsd1.snr[1],&obsd1.ph[1],&obsd1.pr[1],
   &ch2,&obsd1.snr[2],&obsd1.ph[2],&obsd1.pr[2]
   )!=18) continue;
  if(samprate != use_samprate) continue;
  if(p2cut>=0 && obsd1.snr[1]<p2cut) continue;
  if(p1cut>=0 && obsd1.snr[2]<p1cut) continue;
  if(p1cutfac>=0. && obsd1.snr[2]<(short)(p1cutfac*obsd1.snr[0]*obsd1.snr[0]/1000.)) continue;
  if(obsd1.t<tlast) fprintf(stderr,"t<tlast %u < %u\n",obsd1.t,tlast);
  tlast=obsd1.t;
  memcpy((void *)(obsd+ibuf),(void *)&obsd1,sizeof(obsd_str));
  if(q) printf("obsdt=%u\n",obsd[ibuf].t);
  ibuf++;
  if(ibuf>=nbuf-1) print_rnx();
}

while(ibuf>0) print_rnx();

}

/*************************************************************/

/* print most recent data */
/* collapse remaining records */

void print_rnx(void) {

unsigned int tmin;
int indx[20];
int i,nout,nin,ndiscrep;
static int first=1;
time_t ttmp;

/* Sort records by time and prn */

for(i=0;i<ibuf;i++) indx[i]=i;
qsort((void *) indx, (size_t) ibuf, sizeof(int), (int (*)(void*,void*))obsdcmp);
/*
fprintf(stderr,"Hello %d\n",ibuf);
for(i=0;i<ibuf;i++) {
  fprintf(stderr,"%d %d\n",obsd[indx[i]].t,obsd[indx[i]].prn);
}
*/

/* Find minimum time and all records for that time */

tmin=obsd[indx[0]].t;
nout=0;
for(i=0;i<ibuf;i++) {
  if(obsd[i].t==tmin) { nout++; obsd[i].flag=1; }
  else obsd[i].flag=3;
}
/*
fprintf(stderr,"tmin %d nout %d\n",tmin,nout);
*/

/* Look for discrepant records */

ndiscrep=0;
for(i=0;i<nout-1;i++) {
  if(obsd[indx[i]].prn==obsd[indx[i+1]].prn) { 
    fprintf(stderr,"Conflicting records at time %d, prn %d\n",(int)tmin,obsd[indx[i]].prn);
    if(delconflict) {
      if(obsd[indx[i]].flag==1) { obsd[indx[i]].flag=2; ndiscrep++; }
      if(obsd[indx[i+1]].flag==1) { obsd[indx[i+1]].flag=2; ndiscrep++; }
    }
  }
}

/* Print out tmin records */

ttmp=(time_t) (J2000 + tmin);
pT=gmtime(&ttmp);

if(first) { print_header(); first=0; }

if(tmin<=toutlast) {
  fprintf(stderr,"tmin<toutlast %d < %d\n",tmin,toutlast);
}
else {
  printf(" %02d %2d %2d %2d %2d %10.7f %2d %2d",
   pT->tm_year%100,pT->tm_mon+1,pT->tm_mday,pT->tm_hour,pT->tm_min,(double)pT->tm_sec,0,nout-ndiscrep);
  for(i=0;i<nout;i++) if(obsd[indx[i]].flag==1) printf(" %2d",obsd[indx[i]].prn); 
  printf("\n");

  for(i=0;i<nout;i++) {
    if(obsd[indx[i]].flag==2) continue;
    printf("%14.3f %1d",obsd[indx[i]].ph[0],snr2n(obsd[indx[i]].snr[0]));
    printf("%14.3f %1d",obsd[indx[i]].ph[1],snr2n(obsd[indx[i]].snr[1]));
    printf("%14.3f  ",obsd[indx[i]].pr[2]);
    printf("%14.3f  ",obsd[indx[i]].pr[1]);
    printf("%14.3f  ",obsd[indx[i]].pr[0]);
    printf("\n");
  }
  toutlast=tmin;
}

/* Collapse remaining records */

nin=0;
for(i=0;i<ibuf;i++) {
  if(obsd[i].flag==3) {
    memcpy(obsd+nin,obsd+i,sizeof(obsd_str));
    nin++;
  }
}
ibuf=nin;

/*
fprintf(stderr,"nin %d\n",nin);
*/
return;
}

int obsdcmp(int *i, int *j) {
  if(obsd[*i].t<obsd[*j].t) return -1;
  if(obsd[*i].t>obsd[*j].t) return 1;
  if(obsd[*i].prn<obsd[*j].prn) return -1;
  if(obsd[*i].prn>obsd[*j].prn) return 1;
  return 0;
}

/*************************************************************/

int snr2n(short s)
{
  if (s > 10) {
    if (s > 100) return( (s > 500) ? 9 : 8 );
    else return( (s > 50) ? 7 : 6 );
  }
  else {
    if (s > 3) return( (s > 5) ? 5 : 4 );
    else return( (s > 1) ? 3 : 2 );
  }
}


/*************************************************************/

void print_header(void) {

  int i;

  strncpy(pgm, "BJfmt", 20);
  strncpy(runby, getenv("USER"), 20);
  strncpy(marker, "MARKER", 60);
  strncpy(observer, "OBSERVER", 20);
  strncpy(agency, "AGENCY", 40);
  strncpy(recnum, "RECNUM", 20);
  strncpy(rectype, "RECTYPE", 20);
  strncpy(recvers, "RECVERS", 20);
  strncpy(antnum, "ANTNUM", 20);
  strncpy(anttype, "ANTTYPE", 20);
  strncpy(satname, "CHAMP", 8);
  for(i=0;i<3;i++) pos[i] = ant[i] = 0.0;

    printf("     2              OBSERVATION DATA    GPS                 RINEX VERSION / TYPE\n");
    printf("%-20s%-20s%-20sPGM / RUN BY / DATE\n", pgm, runby, date);
    printf("%-20s%-20s%-20sREC # / TYPE / VERS\n", recnum, rectype, recvers);
    printf("%-20s%-20s                    ANT # / TYPE\n", antnum, anttype);
    printf("%14.4f%14.4f%14.4f                  APPROX POSITION XYZ\n", pos[0], pos[1], pos[2]);
    printf("%14.4f%14.4f%14.4f                  ANTENNA: DELTA H/E/N\n", ant[0], ant[1], ant[2]);
    printf("     1     1     0                                          WAVELENGTH FACT L1/2\n");

/*
    if (do_snrs) printf("     7    L1    L2    C1    P1    P2    S1    S2            # / TYPES OF OBSERV\n");
    else printf("     5    L1    L2    C1    P1    P2                        # / TYPES OF OBSERV\n");
*/

    printf("     5    L1    L2    P1    P2    C1                        # / TYPES OF OBSERV\n");
    printf("%6d                                                      INTERVAL\n", (int)use_samprate);
    printf("  %04d    %02d    %02d    %02d    %02d   %09.6f      GPS         TIME OF FIRST OBS\n",
      pT->tm_year + 1900, pT->tm_mon + 1, pT->tm_mday,
      pT->tm_hour, pT->tm_min, (double)pT->tm_sec);
    printf("SNR is mapped to signal strength [0,1,2-9]                  COMMENT\n");
    printf(" SNR: >500 >100  >50  >10   >5   >3  >1  >0  bad  n/a       COMMENT\n");
    printf(" sig:    9    8    7    6    5    4   3   2    1    0       COMMENT\n");
    printf("                                                            END OF HEADER\n");

}

void cmdline(void)
{
  int i;
/*
  tt = (time_t) (time(&tt));
  pT = localtime(&tt);
  yr = (pT->tm_year % 100) + 1900;
  if (yr < 1980) yr += 100;
*/
  
  strncpy(pgm, "BJfmt", 20);
  strncpy(runby, getenv("USER"), 20);
/*
  sprintf(date, "%04d-%02d-%02d %02d:%02d:%02.0f",
      yr, pT->tm_mon + 1, pT->tm_mday,
      pT->tm_hour, pT->tm_min, (double)pT->tm_sec);
*/
  strncpy(marker, "MARKER", 60);
  strncpy(observer, "OBSERVER", 20);
  strncpy(agency, "AGENCY", 40);
  strncpy(recnum, "RECNUM", 20);
  strncpy(rectype, "RECTYPE", 20);
  strncpy(recvers, "RECVERS", 20);
  strncpy(antnum, "ANTNUM", 20);
  strncpy(anttype, "ANTTYPE", 20);
  strncpy(satname, "CHAMP", 8);
  loop(i,3) pos[i] = ant[i] = 0.0;

  pname = ARGV[0];
  iarg = 1;
  while (iarg < NARGS) {
    arg = ARGV[iarg++];
    if (ARG_IS("-h") || ARG_IS("-H")) do_help();
    else if (ARG_IS("-noedit")) noedit = 1;
    else if (ARG_IS("-delconflict")) do_int(&delconflict);
    else if (ARG_IS("-p1cut")) do_short(&p1cut);
    else if (ARG_IS("-p2cut")) do_short(&p2cut);
    else if (ARG_IS("-p1cutfac")) do_real(&p1cutfac);
    else if (ARG_IS("-samprate")) do_uint(&use_samprate);
    else if (ARG_IS("-snrs")) do_snrs = 1;
    else if (ARG_IS("-run_by")) do_str(runby, 20);
    else if (ARG_IS("-marker")) do_str(marker, 60);
    else if (ARG_IS("-observer")) do_str(observer, 20);
    else if (ARG_IS("-agency")) do_str(agency, 40);
    else if (ARG_IS("-recnum")) do_str(recnum, 20);
    else if (ARG_IS("-rectype")) do_str(rectype, 20);
    else if (ARG_IS("-recvers")) do_str(recvers, 20);
    else if (ARG_IS("-antnum")) do_str(antnum, 20);
    else if (ARG_IS("-anttype")) do_str(anttype, 20);
    else {
      fprintf(stderr, "\n%s: Fatal error: unknown flag \"%s\"\n\n", pname, arg);
      fprintf(stderr, "For help type:  %s -h\n\n", pname);
      exit(1);
    }
  }
}

void argprob(char *s)
{
  fprintf(stderr, "\n%s: Fatal problem processing command-line flag \"%s\"\n\n",
     pname, s);
  fprintf(stderr, "For help type:  %s -h\n\n", pname);
  exit(1);
}

void do_str(char *s, int n)
{
  char *argval;
  if (iarg == NARGS) argprob(arg);
  argval = ARGV[iarg++];
  strncpy(s, argval, n+1);
  s[n] = '\0';
  if (strlen(argval) > n) {
    fprintf(stderr,
"\n%s (Warning): Length of argument \"%s\" for flag \"%s\" exceeds allowance %d; truncated to \"%s\"\n",
      pname, argval, arg, n, s);
  }
}

void do_int(int *v)
{
  if (iarg == NARGS) argprob(arg);
  if (sscanf(ARGV[iarg], "%d", v) != 1) argprob(arg);
  iarg++;
}

void do_uint(unsigned int *v)
{
  if (iarg == NARGS) argprob(arg);
  if (sscanf(ARGV[iarg], "%u", v) != 1) argprob(arg);
  iarg++;
}

void do_short(short *v)
{
  if (iarg == NARGS) argprob(arg);
  if (sscanf(ARGV[iarg], "%hd", v) != 1) argprob(arg);
  iarg++;
}


void do_real(double *v)
{
  if (iarg == NARGS) argprob(arg);
  if (sscanf(ARGV[iarg], "%lf", v) != 1) argprob(arg);
  iarg++;
}

void do_vec(double *v)
{
  int i;
  char *argval;
  loop (i,3) {
    if (iarg + i == NARGS) argprob(arg);
    argval = ARGV[iarg + i];
    if (sscanf(argval, "%lf", v + i) != 1) argprob(arg);
  }
  iarg += 3;
}

/************************************************* help *************************************************/


void do_help(void)
{

fprintf(stderr,"BJfmtl_rnx reads output from BJfmtl and writes RINEX 2.1 format output\n");
fprintf(stderr,"\n");
fprintf(stderr,"An example of usage is\n");
fprintf(stderr,"\n");
fprintf(stderr,"cat file | BJfmtl -obs -weekbdry | BJfmtl_rnx -p1cutfac 0.5 -p2cut 5 > site0010.01o\n");
fprintf(stderr,"\n");
fprintf(stderr,"where \"file\" is binary BlackJack data with TRSR wrapper, and \"site0010.01o\" is a rinex format file\n");
fprintf(stderr,"\n");
fprintf(stderr,"---------------------------------------------------------\n");
fprintf(stderr,"\n");
fprintf(stderr,"Some command line flags\n");
fprintf(stderr,"\n");
fprintf(stderr,"-h                 : help\n");
fprintf(stderr,"-p1cut <float x>   : delete rinex record if P1 snr < x (default -1)\n");
fprintf(stderr,"-p2cut <float x>   : delete rinex record if P2 snr < x (default -1)\n");
fprintf(stderr,"-p1cutfac <float x>: delete rinex record if P1 snr < x*(CA snr)*(CA snr)/1000. (default -1)\n");
fprintf(stderr,"-delconflict <0 or 1>: delete conflicting rinex records (when two records for same prn have same time tag). 1 indicates delete. (default 1)\n");
fprintf(stderr,"-samprate <int x>  : use observable data with sample rate x only (default 10)\n");
fprintf(stderr,"\n");
fprintf(stderr,"SNRs are one second SNRV.\n");
fprintf(stderr,"\n");
fprintf(stderr,"The recommended values of p1cutfac and p2cut are 0.5 and 5 respectively.\n");
fprintf(stderr,"\n");
fprintf(stderr,"---------------------------------------------------------\n");
fprintf(stderr,"\n");
fprintf(stderr,"Command line flags for specifying rinex header:\n");
fprintf(stderr,"\n");
fprintf(stderr,"-run_by USER\n");
fprintf(stderr,"-marker MARKER\n");
fprintf(stderr,"-observer OBSERVER\n");
fprintf(stderr,"-agency AGENCY\n");
fprintf(stderr,"-recnum RECNUM\n");
fprintf(stderr,"-rectype RECTYPE\n");
fprintf(stderr,"-recvers RECVERS\n");
fprintf(stderr,"-antnum ANTNUM\n");
fprintf(stderr,"-anttype ANTTYPE\n");
fprintf(stderr,"\n");
  exit(0);

}


/************************************************* end *************************************************/