📄 ppt.cpp
字号:
// PPT.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
#include "fstream.h"
#include "math.h"
#include "ErrorCorrection.h"
#include "GlobleDefine.h"
#include "Phase.h"
#include "RinexNav.h"
#include "RinexObs.h"
#include "SP3.h"
#include "TimeCorTran.h"
#include "iomanip.h"
#include "Matrixcsl.h"
#include "PreClk.h"
#include "Lambda.h"
//#include "MatrixLYQ.h"
#include "FileProcess.h"
#include "PPT.h"
#define ConFile "ConFile.ini"
int main(int argc, char* argv[])
{
printf(" ******************************************************************\n");
printf(" * *\n");
printf(" * Welcome using TriP-T GPS Data Processor 1.0 (2006) *\n");
printf(" * *\n");
printf(" * Developed by Dr.Xiaohong Zhang *\n");
printf(" * *\n");
printf(" * School of Geodesy and Geomatics,Wuhan University,China *\n");
printf(" * *\n");
printf(" * *\n");
printf(" * All Copyright resered! *\n");
printf(" * *\n");
printf(" ******************************************************************\n\n");
CONPARAMETAR Par;
char fname_inputobs [81];
char fname_inputnav [81];
char fname_innavsp3 [81];
char fname_inClk [81];
char fname_outputpos[81];
char fname_outputclk[81];
if(!ReadIni (ConFile,fname_inputobs, fname_inputnav,fname_innavsp3,fname_inClk,fname_outputpos,fname_outputclk,Par))
return 0;
if(Par.ModelPPP)
{
if(!PPPMatchFile(fname_inputobs,fname_innavsp3,fname_inClk,fname_outputpos,fname_outputclk,Par))
return 0;
}
else
{
if(!SPPMatchFile(fname_inputobs,fname_inputnav,fname_inClk,fname_outputpos,fname_outputclk,Par))
return 0;
}
cerr<<" *********************The Program Ends and Succeed!*************************\n";
return 1;
}
int SPP(char * fname_inputobs,char * fname_inputnav,char * fname_inClk,
char * fname_outputpos,char * fname_outputclk,CONPARAMETAR Par)
{
ifstream inObsFile(fname_inputobs,ios::in);
ifstream inRinexNavFile(fname_inputnav,ios::in);
ifstream inClkFile(fname_inClk,ios::in);
ofstream outPosFile(fname_outputpos,ios::out);
ofstream outClkFile(fname_outputclk,ios::out);
/********Reading the Obs file and store the observation*****************/
printf(" Process the OBS file %s \n and Nav file %s...\n",fname_inputobs,fname_inputnav);
/********Reading the Obs file and store the observation*****************/
printf(" Reading the Obs file and store the observation...\n");
OBSREC *obsrec;
if((obsrec=new OBSREC[3000])==NULL)
cerr<<"Error allocating memory!"<<endl;
OBSHEAD head;
int obstype[MAXTYPE];
int TYPENUM=0;
READEPOOBSHEAD(inObsFile,outPosFile,head,obstype,TYPENUM);
int TOTALOBSNUM=0;
while(1)
{ if(inObsFile.peek()==EOF)
break;
READEPOOBSREC(inObsFile,outPosFile,obsrec[TOTALOBSNUM].epohead,
obsrec[TOTALOBSNUM].eporec,TYPENUM);
TOTALOBSNUM++;
}
printf(" The Total number of observations are: %d \n",TOTALOBSNUM);
/***************************************************************************/
//transform the format of the record to detect the cycle slip
printf(" Transform the format of the OBS record ...\n");
SATOBSREC *SatObsRec;
SatObsRec = new SATOBSREC[31];
for (int i=0;i<31;i++)
SatObsRec[i].PObsEpoRec= new OBSRECPERSAT*[TOTALOBSNUM];
int *iEpoNumofSat =new int[MAXNAVSAT];
SatObsRecIni(SatObsRec,obsrec,iEpoNumofSat,TOTALOBSNUM);
printf(" Detect the cycle slip with MW and Geometry-free method...\n");
if(obsrec[0].eporec[0].lOBS[obstype[P1]]&&obsrec[0].eporec[0].lOBS[obstype[P2]])
CycleSlipDetectionMW(outClkFile,SatObsRec,iEpoNumofSat,obstype);
else cout<<" There are not P1 or P2 code for WM cycle slip detection!";
CycleSlipDetectionIon(outClkFile,SatObsRec,iEpoNumofSat,obstype);
printf(" Smooth the pseudorange with the carrier phase...\n");
PhaSmoPseudo(SatObsRec,obstype,iEpoNumofSat);
////////////////////////////////////////////////////////////////////////
//read the Precise clock file
printf(" Read the Precise clock file and store the record...\n");
PRECLK PreClk;
bool isClkInterval30S=false;
if (strstr(fname_inClk,"cod")||strstr(fname_inClk,"COD"))
{
isClkInterval30S=true;
printf(" The Precise clock file is from COD with the Interval 30s.\n");
}
if(isClkInterval30S)
{ PreClk.pRecClk =new SIGCLK[MAXCLKRECORDS];
PreClk.pSavClk =new SATCLK[MAXCLKRECORDS*10];
PreClk.pDifToReference=new SIGCLK[MAXCLKRECORDS*10];
READCODCLKFILE(inClkFile,head.Markername,PreClk);
}
else
{
PreClk.pRecClk =new SIGCLK[MAXCLKRECORDS];
PreClk.pSavClk =new SATCLK[MAXCLKRECORDS];
PreClk.pDifToReference=new SIGCLK[MAXCLKRECORDS];
READCLKFILE(inClkFile,head.Markername,PreClk);
}
if(!head.dapproxmiatexyz.lX||!head.dapproxmiatexyz.lY||!head.dapproxmiatexyz.lZ)
{
head.dapproxmiatexyz.lX=PreClk.RecPos[0];
head.dapproxmiatexyz.lY=PreClk.RecPos[1];
head.dapproxmiatexyz.lZ=PreClk.RecPos[2];
}
//RECANTCORRECTION(head.dAPPROXX,head.dAPPROXY,head.dAPPROXZ,)
/********************Navigation file reading and storing ********************/
printf(" Read the Navigation file and store the Nav record...\n");
NAV *nav;
if((nav=new NAV[500])==NULL)
cerr<<" Error allocating memory!"<<endl;
int TOTALNAVNUM;
NAVHEAD header;
READNAV(inRinexNavFile,outPosFile,header,nav, TOTALNAVNUM);
printf(" The total nav number is %d.\n",TOTALNAVNUM);
SATPOS *satpos=new SATPOS[TOTALOBSNUM*MAXNAVSAT];
long double *dSatClk=new long double[TOTALOBSNUM*MAXNAVSAT];
/*****************************************************************************/
long double *dRawRecClk=new long double[TOTALOBSNUM] ;
long double *dCorrecedRecclk=new long double[TOTALOBSNUM] ;
//the above code is storing the clock data with the clock jump
long double **pPosXYZDif=TwoArrayAlloc(3,TOTALOBSNUM);
// the core part of processing the obs and Eph epoch by epoch
printf(" Process The Observables by Epoch\n to get the final receiver coordinates and clock ...\n");
outPosFile.flags(ios::fixed|ios::showpoint);
outPosFile<<" The initial position in WGS84 cordinate system are : \n "<<head.dapproxmiatexyz.lX
<<" "<<head.dapproxmiatexyz.lY<<" "<<head.dapproxmiatexyz.lZ<<endl;
GPSTIME *gpst=new GPSTIME[TOTALOBSNUM];
for(i=0;i<TOTALOBSNUM;i++)
{ if(i>=536)
int csl=i;
int num=obsrec[i].epohead.iSATNUM;
long double **lL0,**lA,**lP,**lDELTA,**lQx;
int demention=1;
if(!Par.ModelHolding)
demention=4;
lA=TwoArrayAlloc(num,demention);
lP=TwoArrayAlloc(num,num);
lL0=TwoArrayAlloc(num,1);
lDELTA=TwoArrayAlloc(demention,1);
lQx =TwoArrayAlloc(demention,demention);
// CO-efficient Matrix is used for caculate the TDop of each epoch
// lQX= inv(A*AT)
RECPOS *recpos=new RECPOS[2];
long double *ldRecClk=new long double[2];
long double lGDop;
long double lPDop;
long double lTDop;
// deliberately transfer the value of recpos[0] from approximate to 0 to test this part of code
// the initial coordinate of the receiver will have enormous infulence on the final position! NT
recpos[0].xyz.lX=head.dapproxmiatexyz.lX,recpos[0].xyz.lY=head.dapproxmiatexyz.lY,
recpos[0].xyz.lZ=head.dapproxmiatexyz.lZ;
ldRecClk[0]=0;
RECXYZ_BLH(recpos[0]);
UTCCOVGPS(&obsrec[i].epohead.TIME,&gpst[i]);
UTCTIME *Temtime=&obsrec[i].epohead.TIME;
outPosFile<<setw(4)<<i+1<<" "<<2000+Temtime->iYEAR<<"-"<<Temtime->iMONTH<<"-"<<setw(2)<<Temtime->iDAY;
outPosFile<<" "<<setw(2)<<Temtime->iHOUR<<":"<<setw(2)<<Temtime->iMINUTE<<":" <<setw(2)<<(int)Temtime->dSECOND<<" ";
printf(" Epoch %d %d-%d-%d %d:%d:%d.\n",i+1,2000+Temtime->iYEAR,Temtime->iMONTH,Temtime->iDAY,Temtime->iHOUR,Temtime->iMINUTE,(int)Temtime->dSECOND);
while(1)
{ int j;
for(j=0;j<num;j++)
{ if(j==4)
int csl=j;
int k=obsrec[i].eporec[j].iPRN-1;
int inear =FindNearToe(gpst[i].dWEEKSECONDS,obsrec[i].eporec[j].iPRN,TOTALNAVNUM,nav);
GETFINALSVPOSFRONAV(gpst[i].dWEEKSECONDS,obsrec[i].eporec[j].iPRN,inear,nav,
satpos[i*MAXNAVSAT+k],recpos[0].xyz.lX,recpos[0].xyz.lY,recpos[0].xyz.lZ);
long double *dSatClkCofVec=new long double[3];
dSatClkCofVec[0]=nav[inear].dDELTASV;
dSatClkCofVec[1]=nav[inear].dDELTASV1;
dSatClkCofVec[2]=nav[inear].dDELTASV2;
GPSTIME eocgpst;
UTCCOVGPS(&nav[inear].UTC,&eocgpst);
double t=satpos[i*MAXNAVSAT+k].time.dWEEKSECONDS-eocgpst.dWEEKSECONDS;
dSatClk[i*MAXNAVSAT+k]=dSatClkCofVec[0]+dSatClkCofVec[1]*t+dSatClkCofVec[2]*t*t;
delete[]dSatClkCofVec;
//////////////////////////////////////////////////////////////////////////////
//compute the weight of the different satellite
lP[j][j]=1.0;
SATXYZ_AE(satpos[i*MAXNAVSAT+k],recpos[0]);
if(satpos[i*MAXNAVSAT+k].dE<30)
lP[j][j]=sin(satpos[i*MAXNAVSAT+k].dE/180*PI);
if(!obsrec[i].eporec[j].lOBS[obstype[P1]]||!obsrec[i]
.eporec[j].lOBS[obstype[P2]])
lP[j][j]=0;
//just for phase smooth psudorange
//compute the various error correction
int pst,tst,rhst;
pst =PS*pow(1-(2.26e-5)*(recpos->blh.lH-0.0),5.225);
tst =TS-(recpos->blh.lH-0.0)*0.065;
rhst =RHS*exp(-6.396e-4*(recpos->blh.lH-0.0));
long double ldiono=-1.54573*(obsrec[i].eporec[j]
.lOBS[obstype[P1]]-obsrec[i].eporec[j].lOBS[obstype[P2]]);
long double ldtrop=SIMPLEHOPFIELD(Par.TropModel,satpos[i*MAXNAVSAT+k].dE,
PS,TS,RHS,recpos->blh);
/*
long double ldtrop=SIMPLEHOPFIELD(Par.TropModel,satpos[i*MAXNAVSAT+k].dE,
pst,tst,rhst,recpos->blh);*/
long double ldRelCor=0;
ldRelCor=RELATIVECORRECTION(satpos[i*MAXNAVSAT+k]);
//relative correction
long double ldSatCor=dSatClk[i*MAXNAVSAT+k]*VC;
long double ldPseduRow=obsrec[i].eporec[j].lOBS[P1]
+ldSatCor-ldRecClk[0]*VC-ldiono-ldtrop-ldRelCor;
long double ldDistance=SSDistance(recpos[0].xyz.lX,recpos[0].xyz.lY
,recpos[0].xyz.lZ,satpos[i*MAXNAVSAT+k].xyz.lX
,satpos[i*MAXNAVSAT+k].xyz.lY,satpos[i*MAXNAVSAT+k].xyz.lZ);
long double lL[MAXCHA],lM[MAXCHA],lN[MAXCHA];
lL[j]=(satpos[i*MAXNAVSAT+k].xyz.lX-recpos[0].xyz.lX)/ldDistance;
lM[j]=(satpos[i*MAXNAVSAT+k].xyz.lY-recpos[0].xyz.lY)/ldDistance;
lN[j]=(satpos[i*MAXNAVSAT+k].xyz.lZ-recpos[0].xyz.lZ)/ldDistance;
lA[j][0]=-lL[j];
lA[j][1]=-lM[j];
lA[j][2]=-lN[j];
lA[j][3]=1;
lL0[j][0]=ldPseduRow-ldDistance;
}
LEASTSQUAREAJUST(lDELTA,lA,num,4,1,lL0,lP,lQx);
recpos[0].xyz.lX+=lDELTA[0][0];
recpos[0].xyz.lY+=lDELTA[1][0];
recpos[0].xyz.lZ+=lDELTA[2][0];
ldRecClk[0]+=lDELTA[3][0]/VC;
// recpos[1]=recpos[0];
// if((fabs(lDELTA[0][0])+fabs(lDELTA[1][0])+fabs(lDELTA[2][0]))<0.1)
if((fabs(lDELTA[3][0])<0.1))
break;
}
dRawRecClk[i] =ldRecClk[0];
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -