📄 ppt.cpp
字号:
pPosXYZDif[0][i]=recpos[0].xyz.lX-head.dapproxmiatexyz.lX;
pPosXYZDif[1][i]=recpos[0].xyz.lY-head.dapproxmiatexyz.lY;
pPosXYZDif[2][i]=recpos[0].xyz.lZ-head.dapproxmiatexyz.lZ;
lTDop =pow(lQx[demention-1][demention-1],0.5);
lPDop =pow((lQx[0][0]+lQx[1][1]+lQx[2][2]),0.5);
lGDop =pow((lQx[0][0]+lQx[1][1]+lQx[2][2]+lQx[3][3]),0.5);
//store the epoch receiver clock error to the dRawRecClk
outPosFile<<" "<<setprecision(3)<<setw(7)
<<recpos[0].xyz.lX-head.dapproxmiatexyz.lX<<" "<<setw(7)<<recpos[0].xyz.lY-head.dapproxmiatexyz.lY
<<" "<<setw(7)<<recpos[0].xyz.lZ-head.dapproxmiatexyz.lZ;
outPosFile<<" "<<setprecision(4)<<setw(15)<<recpos[0].xyz.lX<<" "<<setw(15)<<recpos[0].xyz.lY<<" "<<setw(15)<<recpos[0].xyz.lZ
<<" "<<setw(12)<<setprecision(10)<<ldRecClk[0]<<" ";
outPosFile<<" "<<setw(15)<<setprecision(11)<<ldRecClk[0]<<" "<<lTDop<<" "<<lPDop<<" "<<lGDop<<endl;
TwoArrayFree(lDELTA);
TwoArrayFree(lA) ;
TwoArrayFree(lL0);
TwoArrayFree(lP) ;
TwoArrayFree(lQx);
delete[]recpos ;
delete[]ldRecClk ;
}
/*
outPosFile<<" X coordinate： ";
GETPOSRMS(outPosFile,pPosXYZDif[0],TOTALOBSNUM);
outPosFile<<" Y coordinate： ";
GETPOSRMS(outPosFile,pPosXYZDif[1],TOTALOBSNUM);
outPosFile<<" Z coordinate： ";
GETPOSRMS(outPosFile,pPosXYZDif[2],TOTALOBSNUM);*/
// process the clock error record
/* if(!CORRECTCLK(TOTALOBSNUM,dRawRecClk,dCorrecedRecclk))
{ printf(" The clock has dected with clock jump by TRIP-T \n We correct the jump with our software.\n ");
outClkFile<<" After the clock jump being detected and corrected :"<<endl;
GETCLKRMS(outClkFile,dCorrecedRecclk,TOTALOBSNUM);
}
else
GETCLKRMS(outClkFile,dRawRecClk,TOTALOBSNUM);*/
if(PreClk.Makername)
{ int EpoPerGroup=10; //clock devided to groups at the interval of 30s
int ClkNum=TOTALOBSNUM/EpoPerGroup;
long double * lGroupClk=new long double[ClkNum];
//Clk30SecTo5Min(outClkFile,dRawRecClk,lGroupClk,ClkNum,EpoPerGroup);
for(int i=0;i<ClkNum;i++)
lGroupClk[i]=dRawRecClk[EpoPerGroup*i];
GETWCLKRMS(outClkFile,lGroupClk,PreClk.pRecClk,PreClk.pDifToReference,ClkNum);
// outClkFile<<" The true value of clock from the IGS CLK Precise clock files are :"<<endl;
// for (i=0;i<MAXCLKRECORDS;i++)
// outClkFile<<setw(4)<<i<<" "<<setw(12)<<PreClk.pRecClk[i].lClk[0]<<endl;
delete[]lGroupClk;
}
/************************************************************************/
delete[]iEpoNumofSat;
SatObsRecDesTroy(SatObsRec);
delete[]dRawRecClk;
delete[]dCorrecedRecclk;
delete[]satpos;
delete[]dSatClk;
delete[]gpst;
// delete[]SatObsRec; //have deleted before in SatObsRecDesTroy
delete[]obsrec;
TwoArrayFree(pPosXYZDif);
inObsFile.close();
outPosFile.close();
inRinexNavFile.close();
outPosFile.close();
outClkFile.close();
return 1;
}
int PPP(char * fname_inputobs,char * fname_innavsp3,char * fname_inClk,
char * fname_outputpos,char * fname_outputclk,CONPARAMETAR Par)
{
ifstream inObsFile(fname_inputobs,ios::in);
ifstream inSp3NavFile(fname_innavsp3,ios::in);
ifstream inClkFile(fname_inClk,ios::in);
ofstream outPosFile(fname_outputpos,ios::out);
ofstream outClkFile(fname_outputclk,ios::out);
printf(" Process the OBS file %s \n and SP3 file %s\n and clk file %s .....\n",fname_inputobs,fname_innavsp3,fname_inClk);
/********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++;
if(TOTALOBSNUM==2826)
int csl=8;
}
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,)
/********************sp3 Precise Eph reading and storing ********************/
printf(" Process the SP3/Eph file and store the satellite record...\n");
NAVSP3HEAD NavSp3Head;
READNAVSP3HEAD(inSp3NavFile,outPosFile,NavSp3Head);
int iEpoSatNum = NavSp3Head.iRealSatNumber;
int iEpoNum =NavSp3Head.iEpoNumber;
NAVSP3EPO *NavSp3Epo;
NavSp3Epo=new NAVSP3EPO[iEpoNum];
READNAVSP3RECORD(inSp3NavFile,outPosFile,NavSp3Epo,iEpoSatNum);
printf(" The Total number of SP3/Eph record are %d.\n",iEpoNum);
/*******Sp3 Precise Eph compute the position and clock of the satellite *****/
printf(" Poly the Precise satellite coordinates and clock to each OBS epoch...\n");
SATPOS *satpos=new SATPOS[TOTALOBSNUM*MAXNAVSAT];
long double *dSatClk=new long double[TOTALOBSNUM*MAXNAVSAT];
/********************************************************
***************************/
for(i=0;i<MAXNAVSAT;i++)
{ //UTCCOVGPS(&obsrec[i].epohead.TIME,&gpst[i]);
if(i==8)
int scl=i;
printf(" Process the satellite No %d.\n",i+1);
int PolyOrder;
long double **lCofficient;
for(int j=0;j<iEpoNumofSat[i];j++)
{
int iObsEpoNum=SatObsRec[i].PObsEpoRec[j]->iEpoNum;
if(iObsEpoNum==541)
int csl=j;
satpos[iObsEpoNum*MAXNAVSAT+i].iPRN=i+1;
double dDletaT1;
double dDletaT0;
dDletaT0=0.069;
PolyOrder=11;
// int iOldIni=ININUMFROMTIME(iObsEpoNum-1,PolyOrder);
int iNewIni=ININUMFROMTIME(iObsEpoNum,PolyOrder);
lCofficient=TwoArrayAlloc(PolyOrder,3);
// if((iNewIni>iOldIni)||SatObsRec[i].PObsEpoRec[j]->iCycleSlip>0)
CREATSP3POLYCOEF(iNewIni,PolyOrder,i+1,lCofficient,NavSp3Epo);
do
{ dDletaT1=dDletaT0;
GPSTIME time;
UTCCOVGPS(&obsrec[iObsEpoNum].epohead.TIME,&time);
long double t=time.dWEEKSECONDS-dDletaT1;
GETSATPOSFROMSP3(t,iNewIni,PolyOrder,
NavSp3Epo,lCofficient,satpos[iObsEpoNum*MAXNAVSAT+i]);
double dDis=SSDistance(satpos[iObsEpoNum*MAXNAVSAT+i].xyz.lX,
satpos[iObsEpoNum*MAXNAVSAT+i].xyz.lY,satpos[iObsEpoNum*MAXNAVSAT+i]
.xyz.lZ, head.dapproxmiatexyz.lX,head.dapproxmiatexyz.lY,head.dapproxmiatexyz.lZ);
//dDletaT0=dDis/VC+dSatClkCofVec[0];
dDletaT0=dDis/VC;
} while(fabs(dDletaT0-dDletaT1)>1E-6);
Earthrotation(dDletaT0,&satpos[iObsEpoNum*MAXNAVSAT+i]);
// Theory and experiments proves It need rotation
long double *dSatClkCofVec= new long double[2];
PolyOrder=3;
iNewIni=ININUMFROMTIME(iObsEpoNum,PolyOrder);
lCofficient=TwoArrayAlloc(PolyOrder,1);
CREATCLKPOLYCOEF(iNewIni,PolyOrder,i+1,lCofficient,NavSp3Epo);
GPSTIME time;
UTCCOVGPS(&obsrec[iObsEpoNum].epohead.TIME,&time);
long double t=time.dWEEKSECONDS;
GETSATCLKFROMSP3(t-dDletaT0,iNewIni,PolyOrder,
NavSp3Epo,lCofficient,dSatClkCofVec);
dSatClk[iObsEpoNum*MAXNAVSAT+i]=dSatClkCofVec[0]+dDletaT0
*dSatClkCofVec[1];
delete[]dSatClkCofVec;
// outPosFile<<dSatClk[iObsEpoNum*MAXNAVSAT+i];
}
}
delete[]NavSp3Epo;
/*****************************************************************************/
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);
GPSTIME *gpst=new GPSTIME[TOTALOBSNUM];
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 ITRF frame are : \n "<<head.dapproxmiatexyz.lX
<<" "<<head.dapproxmiatexyz.lY<<" "<<head.dapproxmiatexyz.lZ<<endl;
for(i=0;i<TOTALOBSNUM;i++)
{
int num=obsrec[i].epohead.iSATNUM;
long double **lL0,**lA,**lP,**lDELTA,**lQx;
/*
BOOL *BigResidual = new BOOL[num];
//the residual whether or not go beyond the limit
for(int j=0;j<num;j++)
{ lV[j]=0;
BigResidual[j]=0;
}
*/
if(i>=1010)
int csl=i;
int demention= 1 ;
int row = num*2;
if(!Par.ModelHolding)
demention=4;
int IAmbSmoLen=AmbSmoLen;
int PriJ=0;
//The order with same PRN in the satlist between last epoch and this
for(int j=0;j<num;j++)
{ if(obsrec[i].eporec[j].iCycleFlag)
obsrec[i].eporec[j].iAmbSmoTime=1;
else
{ PriJ=SEAPRN(obsrec[i-1].epohead.SatList,obsrec[i].eporec[j].iPRN);
obsrec[i].eporec[j].iAmbSmoTime=obsrec[i-1].eporec[PriJ].iAmbSmoTime+1;
}
if(obsrec[i].eporec[j].iAmbSmoTime<AmbSmoLen)
demention+=1;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -