/*
readrinex3.c : extract data from a rinex file (version 3.00)

 ex : [Observables, epochs, svG, svE, svR, svS, apcoords]=readrinex3(fic);
 
 This routine extracts quickly the data from a rinex file (version 3.00);  
 GPS (G), Galileo (E),GLONASS (R) and SBAS (S) datas are extracted; all the observations are gathered into
 one structure. 
 The epoch flag is extracted.
 A description of the rinex 3.00 format can be found on  : ftp://igscb.jpl.nasa.gov/pub/data/format/rinex300.pdf
 
 Input arguments
    fic   : name of the file to be read

 Output arguments:
    Observables  : structure containing the different observables (C1C,C5X,...)
    epochs       : matrix (nb_epochs * 7); form of the column : <Year;Month;Day;HH;MM;SS;epoch flag>
    svG          : vector (column); the ith element gives the PRN of the sat in the ith column of the GPS Observables
    svE          : vector (column); the ith element gives the PRN of the sat in the ith column of the GALILEO Observables
    svR          : vector (column); the ith element gives the PRN of the sat in the ith column of the GLONASS Observables
    svS          : vector (column); the ith element gives the PRN of the sat in the ith column of the SBAS Observables
    apcoords     : vector (column); gives the approx position XYZ
	
 ----------------------------------------------------------------------
 File.....: readrinex3.c
 Date.....: 18-MAY-2011
 Author...: Emmanuel SCHIELIN
            (inspired from readrinex.c, Eric CALAIS)
 ----------------------------------------------------------------------
*/
 

#include <stdio.h>
#include <string.h>
#include <math.h>
#include "mex.h"


FILE *fIn;

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	char *FileName;
    char **Observables;
	char GetString[200];
    int NumObsG, NumObsE,NumObsR, NumObsS;
	int i,k,FileNameLength,NumObs,NumSats,NumEpochs;
    int ns_col;
    int NumSatsE, NumSatsG,NumSatsR, NumSatsS;
	int ColumnIndexG[99]={0},SatColumnG[99]={0},SatCountG[99]={0};
    int ColumnIndexE[99]={0},SatColumnE[99]={0},SatCountE[99]={0};
    int ColumnIndexR[99]={0},SatColumnR[99]={0},SatCountR[99]={0};
    int ColumnIndexS[99]={0},SatColumnS[99]={0},SatCountS[99]={0};
    //number(99) can be modified according to the constellations
	double XYZ[3], *pSatellites;
	double **pObs, *pEpoch;
	mxArray **mxObs;
	
	float version;

	if (nrhs!=1)
		mexErrMsgTxt("Usage: [Observables,epochs,svG,svE,svR,svS,apcoords]=readrinex(filename).");
    
    if (nlhs!=7)
		mexErrMsgTxt("Usage: [Observables,epochs,svG,svE,svR,svS,apcoords]=readrinex(filename).");
    
	if (mxIsChar(prhs[0])==0)
		mexErrMsgTxt("First argument must be a rinex file name.");

	FileNameLength = (mxGetM(prhs[0])*mxGetN(prhs[0]))+1;
	FileName =(char *)mxCalloc(FileNameLength,sizeof(char));

	mxGetString(prhs[0],FileName,FileNameLength);

	/* Try to open file. */

		fIn = fopen(FileName, "r");
		if (fIn==NULL)
			mexErrMsgTxt("Could not open file.");

	/* Check for valid RINEX file.*/

		fgets(GetString, 82, fIn);
		sscanf(GetString,"%f",&version);
		printf("\n Reading : '%s', \n RINEX version %4.2f\n",FileName,version);
		if(strstr(GetString,"RINEX VERSION / TYPE") == 0){
			mexErrMsgTxt("Invalid RINEX header.");
		}
        if(strstr(GetString,"3.00") == 0){
			printf("\n%s",GetString);
            printf("Beware, this script can only read Rinex 3.00 version. \n");
			mexErrMsgTxt("Beware, this script can only read Rinex 3.00 version. \n");
		}

	/* Determine number/type of observables, number/type of satellites. */
        NumSatsG = 0;   NumObsG = 0;
        NumSatsE = 0;   NumObsE = 0;
        NumSatsR = 0;   NumObsR = 0;
        NumSatsS = 0;   NumObsS = 0;
		while (fgets(GetString, 82, fIn))
		{
			if(GetString[68]=='/') /* Look for the '/'(=47 in ascii) in '# / TYPES OF OBSERV' */
            /*Rq = in C, table indice begins at 0, so GetString[62] gives us the symbol of the 63rd column... */
			{
                if(GetString[0]=='G')  /* Look for the 'G'=GPS sat in 'TYPES OF OBSERV' */
                    NumObsG=atoi(&GetString[4]);

                if(GetString[0]=='E')  /* Look for the 'E'=GALILEO sat in 'TYPES OF OBSERV' */
                    NumObsE=atoi(&GetString[4]);
                
                if(GetString[0]=='R')  /* Look for the 'R'=GLONASS sat in 'TYPES OF OBSERV' */
                    NumObsR=atoi(&GetString[4]);

                if(GetString[0]=='S')  /* Look for the 'S'=SBAS sat in 'TYPES OF OBSERV' */
                    NumObsS=atoi(&GetString[4]);
			}
            
			if(GetString[62]=='P') /* Look for the 'P' in 'APPROX POSITION XYZ' */
			{
				XYZ[0]=atof(&GetString[0]);
				XYZ[1]=atof(&GetString[15]);
				XYZ[2]=atof(&GetString[29]);
				printf("Approx. position (XYZ, m)= %.4f %.4f %.4f\n",XYZ[0],XYZ[1],XYZ[2]);
			}
            if(GetString[60]=='#') /*Look for the '#' in '# OF SATELLITES' */    
                sscanf(GetString,"%d",&NumSats);
            
            if(GetString[3]=='G' && GetString[60]==80) {  /*Look for the 'G' before the PRN of the satellite */
                SatCountG[NumSatsG]=atoi(&GetString[4]);    // SatCountG(sat num column in obs vect) = PRN sat
                ColumnIndexG[SatCountG[NumSatsG]]=NumSatsG; // ColumnIndexG(PRN sat) = column number of this PRN in the Obs vector
                NumSatsG++;      /* We count the number of GPS sats*/ 
            }
            if(GetString[3]=='E' && GetString[60]==80) {
                SatCountE[NumSatsE]=atoi(&GetString[4]);
                ColumnIndexE[SatCountE[NumSatsE]]=NumSatsE; 
                NumSatsE++;   
            }
            if(GetString[3]=='R' && GetString[60]==80) {  /*Look for the 'G' before the PRN of the satellite */
                SatCountR[NumSatsR]=atoi(&GetString[4]);    
                ColumnIndexR[SatCountR[NumSatsR]]=NumSatsR; 
                NumSatsR++;  
            }
            if(GetString[3]=='S' && GetString[60]==80) { /*Look for the 'S' before the PRN of the satellite */
                SatCountS[NumSatsS]=atoi(&GetString[4]);
                ColumnIndexS[SatCountS[NumSatsS]]=NumSatsS; 
                NumSatsS++; 
            }
                        
			if(GetString[62]=='D') /* Look for the 1st 'D' in 'END OF HEADER' */
            {   // Reservation of the observables = sum of all (GPS,GALILEO,GLONASS,SBAS) obs
                NumObs = NumObsE + NumObsG + NumObsR + NumObsS;               
 
                printf("Found %d GPS, %d GALILEO, %d GLONASS and %d SBAS observables\n",NumObsG,NumObsE,NumObsR,NumObsS);
                printf("Found %d observables\n",NumObs);
				break;
            }
		}//END while
 
        /* New loop to fill the Observables. */
        Observables = (char **)mxCalloc(NumObs,sizeof(char *));
        fseek(fIn,0,SEEK_SET);
		while (fgets(GetString, 82, fIn))
        {
               if(GetString[0]=='G' && GetString[74]=='T') { //GPS, Look for the 'T' in 'OBS TYPES' to be sure it's really the ligne of an observation ...
                    for(i=0;i<NumObsG;i++) {
                        Observables[i] = (char *)mxCalloc(3,sizeof(char));
                        sscanf(&GetString[7+i*4],"%3s",Observables[i]);
                    }
                }
               if(GetString[0]=='E' && GetString[74]=='T') { //GALILEO
                    for(i=NumObsG;i<NumObsE+NumObsG;i++) {
                        Observables[i] = (char *)mxCalloc(3,sizeof(char));
                        sscanf(&GetString[7+(i-NumObsG)*4],"%3s",Observables[i]);
                    }
                }
              if(GetString[0]=='R' && GetString[74]=='T') { //GLONASS
                    for(i=NumObsE+NumObsG;i<NumObsE+NumObsG+NumObsR;i++) {
                        Observables[i] = (char *)mxCalloc(3,sizeof(char));
                        sscanf(&GetString[7+(i-NumObsE-NumObsG)*4],"%3s",Observables[i]);
                    }
                }
              if(GetString[0]=='S' && GetString[74]=='T') { //SBAS
                    for(i=NumObsE+NumObsG+NumObsR;i<NumObsE+NumObsG+NumObsR+NumObsS;i++) {
                        Observables[i] = (char *)mxCalloc(3,sizeof(char));
                        sscanf(&GetString[7+(i-NumObsE-NumObsG-NumObsR)*4],"%3s",Observables[i]);
                    }
                }
			if(GetString[62]=='D') /* Look for the 1st 'D' in 'END OF HEADER' */ 
				break;
        }
	
	/* Count up number of epochs. */

		NumEpochs=0;
		while(fgets(GetString, 82, fIn))
		{
            if(GetString[0]=='>') /*Look for the 1st '>' behind every epoch */
                NumEpochs++;
		}
        printf("Found %d epochs \n",NumEpochs);
		printf("Found %d satellites (%d GPS, %d GALILEO, %d GLONASS, %d SBAS) \n",NumSats,NumSatsG,NumSatsE,NumSatsR,NumSatsS);

        
	/* Allocate memory to hold outputs */

		mxObs = (mxArray **)mxCalloc(NumObs,sizeof(mxArray *));
		pObs = (double **)mxCalloc(NumObs,sizeof(double *));
		plhs[0]=mxCreateStructMatrix(1,1,NumObs,Observables);   /* Obs : output1 */
		plhs[1]=mxCreateDoubleMatrix(NumEpochs,7,mxREAL);       /* epochs : output2 */
		plhs[2]=mxCreateDoubleMatrix(NumSatsG,1,mxREAL);        /* sv_G : output3 */
        plhs[3]=mxCreateDoubleMatrix(NumSatsE,1,mxREAL);        /* sv_E : output4 */
		plhs[4]=mxCreateDoubleMatrix(NumSatsR,1,mxREAL);        /* sv_R : output5 */
        plhs[5]=mxCreateDoubleMatrix(NumSatsS,1,mxREAL);        /* sv_S : output6 */
		plhs[6]=mxCreateDoubleMatrix(3,1,mxREAL);               /* apcoords : output7 */

		pSatellites=mxGetPr(plhs[2]);
		for (i=0;i<NumSatsG;i++)
			pSatellites[i]=SatCountG[i]; /*pSatellites = pointer to the first element of plhs[2]*/
        
        pSatellites=mxGetPr(plhs[3]);
		for (i=0;i<NumSatsE;i++)
			pSatellites[i]=SatCountE[i];

        pSatellites=mxGetPr(plhs[4]);
		for (i=0;i<NumSatsR;i++)
			pSatellites[i]=SatCountR[i];
        
        pSatellites=mxGetPr(plhs[5]);
		for (i=0;i<NumSatsS;i++)
			pSatellites[i]=SatCountS[i];
       
		memcpy(mxGetPr(plhs[6]),XYZ,3*sizeof(double));
		pEpoch = mxGetPr(plhs[1]);

		for(i=0;i<NumObsG;i++)  
		{ /*Allocation of the first part of Obs corresponding to GPS sats*/
			mxObs[i] = mxCreateDoubleMatrix(NumEpochs,NumSatsG,mxREAL);
			pObs[i] = mxGetPr(mxObs[i]);
			memset(pObs[i],255,NumEpochs*NumSatsG*sizeof(double));
		} 
        for(i=NumObsG;i<NumObsE+NumObsG;i++)  
		{ /*Allocation of the first part of Obs corresponding to GALILEO sats*/
			mxObs[i] = mxCreateDoubleMatrix(NumEpochs,NumSatsE,mxREAL);
			pObs[i] = mxGetPr(mxObs[i]);
			memset(pObs[i],255,NumEpochs*NumSatsE*sizeof(double));
		}
        for(i=NumObsE+NumObsG;i<NumObsE+NumObsG+NumObsR;i++)  
		{ /*Allocation of the first part of Obs corresponding to GLONASS sats*/
			mxObs[i] = mxCreateDoubleMatrix(NumEpochs,NumSatsR,mxREAL);
			pObs[i] = mxGetPr(mxObs[i]);
			memset(pObs[i],255,NumEpochs*NumSatsR*sizeof(double));
		}
        for(i=NumObsE+NumObsG+NumObsR;i<NumObsE+NumObsG+NumObsR+NumObsS;i++)  
		{ /*Allocation of the first part of Obs corresponding to SBAS sats*/
			mxObs[i] = mxCreateDoubleMatrix(NumEpochs,NumSatsS,mxREAL);
			pObs[i] = mxGetPr(mxObs[i]);
			memset(pObs[i],255,NumEpochs*NumSatsS*sizeof(double));
		}
	/* Get observables. */

		k=-1;
		fseek(fIn,0,SEEK_SET);
		while (fgets(GetString, 82, fIn))
			if(GetString[62]==68) /* Look for the 1st 'D' in 'END OF HEADER' */ 
				break;
        while (fgets(GetString, 200, fIn))
        {
            if(GetString[0]=='>')  /*Look for the 1st '>' behind every epoch */
            {
                k++;
                pEpoch[k]=atoi(&GetString[2]);              // year
                pEpoch[k+NumEpochs]=atoi(&GetString[7]);    // month
                pEpoch[k+2*NumEpochs]=atoi(&GetString[10]); // day
                pEpoch[k+3*NumEpochs]=atoi(&GetString[13]); // hh
                pEpoch[k+4*NumEpochs]=atoi(&GetString[16]); // mm
                pEpoch[k+5*NumEpochs]=atof(&GetString[19]); // ss (float!)
                pEpoch[k+6*NumEpochs]=atoi(&GetString[31]); // Epoch flag
            }                  
            else {      
                if(GetString[0]=='G') // for GPS
                {
                    ns_col = (strlen(GetString) - 4)/16;		//number of Observation for this sat and this epoch
                    for(i=0;i<ns_col;i++) 
                        pObs[i][k+ColumnIndexG[atoi(&GetString[1])]*NumEpochs]=atof(&GetString[4+i*16]);
                }
                if(GetString[0]=='E') // for GALILEO 
                {
                    ns_col = (strlen(GetString) - 4)/16;		//number of Observation for this sat and this epoch
                    for(i=0;i<ns_col;i++)
                        pObs[i+NumObsG][k+ColumnIndexE[atoi(&GetString[1])]*NumEpochs]=atof(&GetString[4+i*16]);
                }  
                if(GetString[0]=='R') // for GLONASS
                {
                    ns_col = (strlen(GetString) - 4)/16;		//number of Observation for this sat and this epoch
                    for(i=0;i<ns_col;i++) 
                        pObs[i+NumObsG+NumObsE][k+ColumnIndexR[atoi(&GetString[1])]*NumEpochs]=atof(&GetString[4+i*16]);
                }
                if(GetString[0]=='S') // for SBAS
                {
                    ns_col = (strlen(GetString) - 4)/16;		//number of Observation for this sat and this epoch
                    for(i=0;i<ns_col;i++)
                        pObs[i+NumObsG+NumObsE+NumObsR][k+ColumnIndexS[atoi(&GetString[1])]*NumEpochs]=atof(&GetString[4+i*16]);
                }
            } // END else
        }// END while
 
	fclose(fIn);

	/* Output observables to structure. */
	
	for(i=0;i<NumObs;i++)
        mxSetFieldByNumber(plhs[0],0,i,mxObs[i]);

}// END function