📄 ppt.cpp
字号:
}
lA=TwoArrayAlloc(row,demention);
lP=TwoArrayAlloc(row,row);
lL0=TwoArrayAlloc(row,1);
lDELTA=TwoArrayAlloc(demention,1);
lQx =TwoArrayAlloc(demention,demention);
RECPOS *recpos=new RECPOS[2];
long double *ldRecClk=new long double[2];
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;
int p=1; // input the value in Delta to IAmbiguity in each epoch
for(j=0;j<num;j++)
{
if(j==4)
int csl=j;
int k=obsrec[i].eporec[j].iPRN-1;
//compute the weight of the different satellite
lP[2*j][2*j]=1.0;
lP[2*j+1][2*j+1]=Par.PseudoSTD/Par.PhaseSTD;
SATXYZ_AE(satpos[i*MAXNAVSAT+k],recpos[0]);
if(satpos[i*MAXNAVSAT+k].dE<30)
{ lP[2*j][2*j] *=sin(satpos[i*MAXNAVSAT+k].dE/180*PI);
lP[2*j+1][2*j+1]*=sin(satpos[i*MAXNAVSAT+k].dE/180*PI);
}
if(!obsrec[i].eporec[j].lOBS[obstype[P1]]||!obsrec[i]
.eporec[j].lOBS[obstype[P2]])
{ lP[2*j][2*j]=0;
//lP[2*j+1][2*j+1]=0;
}
if(!obsrec[i].eporec[j].lOBS[obstype[L1]]||!obsrec[i]
.eporec[j].lOBS[obstype[L2]])
{ //lP[2*j][2*j]=0;
lP[2*j+1][2*j+1]=0;
}
// if(BigResidual[j])
// lP[j][j]=0;
// if(obsrec[i].eporec[j].iCycleSlip==-1)
// lP[i][j]=0;
//just for phase smooth psudorange
//compute the various error correction
if(!obsrec[i].eporec[j].lOBS[obstype[P1]]&&obsrec[i].eporec[j].lOBS[obstype[C1]])
obsrec[i].eporec[j].lOBS[obstype[P1]]=obsrec[i].eporec[j].lOBS[obstype[C1]];
long double ldiono;
if(obsrec[i].eporec[j].lOBS[obstype[P1]]&&obsrec[i].eporec[j].lOBS[obstype[P2]])
ldiono=-1.54573*(obsrec[i].eporec[j]
.lOBS[obstype[P1]]-obsrec[i].eporec[j].lOBS[obstype[P2]]);// just for pseuorange
else ldiono=0;
long double ldtrop=SIMPLEHOPFIELD(Par.TropModel,satpos[i*MAXNAVSAT+k].dE,
PS,TS,RHS,recpos->blh);
long double ldRelCor=0;
ldRelCor=RELATIVECORRECTION(satpos[i*MAXNAVSAT+k]);
//relative correction
long double ldSatCor=dSatClk[i*MAXNAVSAT+k]*VC;
if((i>0)&&(obsrec[i].eporec[j].iAmbSmoTime>=IAmbSmoLen))
{ PriJ=SEAPRN(obsrec[i-1].epohead.SatList,obsrec[i].eporec[j].iPRN);
obsrec[i].eporec[j].lAmbiguity=obsrec[i-1].eporec[PriJ].lAmbiguity;
}
// if(!obsrec[i].eporec[j].lOBS[obstype[P1]])
// obsrec[i].eporec[j].lOBS[obstype[P1]]=obsrec[i].eporec[j].lOBS[obstype[P2]];
//
long double ldPseduRow=obsrec[i].eporec[j].lOBS[P1]
+ldSatCor-ldRecClk[0]*VC-ldiono-ldtrop+ldRelCor;
long double ldPhaseRow=obsrec[i].eporec[j].lOBS[L1]
*VC/FREQ_L1+ldSatCor-ldRecClk[0]*VC+ldiono-ldtrop-ldRelCor;
ldPhaseRow-=obsrec[i].eporec[j].lAmbiguity*VC/FREQ_L1;
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];
if(Par.ModelHolding)
{ lA[2*j][0]=1;
lA[2*j+1][0]=1;}
if(!Par.ModelHolding)
{
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 ;}
for(int t=1;t<demention;t++)
{ lA[2*j][t]=0;
lA[2*j+1][t]=0;
}
//
if(obsrec[i].eporec[j].iAmbSmoTime<IAmbSmoLen)
{ lA[2*j+1][p]=VC/FREQ_L1;p++;}
lL0[2*j][0] =ldPseduRow-ldDistance;
lL0[2*j+1][0]=ldPhaseRow-ldDistance;
obsrec[i].eporec[j].lPseudoRes=0;
obsrec[i].eporec[j].lPhaseRes=0;
}
LEASTSQUAREAJUST(lDELTA,lA,row,demention,1,lL0,lP,lQx);
p=1;
for(j=0;j<num;j++)
{
for(int k=0;k<demention;k++)
{ obsrec[i].eporec[j].lPseudoRes += lA[2*j][k]*lDELTA[k][0];
obsrec[i].eporec[j].lPhaseRes += lA[2*j+1][k]*lDELTA[k][0];
}
if(obsrec[i].eporec[j].iAmbSmoTime<IAmbSmoLen)
{ obsrec[i].eporec[j].lAmbiguity+=lDELTA[p][0];
p++;}
obsrec[i].eporec[j].lPseudoRes -= lL0[2*j][0];
obsrec[i].eporec[j].lPhaseRes -= lL0[2*j+1][0];
// if(fabs(lV[j])>=3*1)
// BigResidual[j]=TRUE;
}
// int recompute=0;
// for(j=0;j<num;j++)
// recompute+=BigResidual[j];
//
// if(recompute)
// { for(j=0;j<num;j++)
// lV[j]=0;
// continue;
// }
if(!Par.ModelHolding)
{
recpos[0].xyz.lX+=lDELTA[0][0];
recpos[0].xyz.lY+=lDELTA[1][0];
recpos[0].xyz.lZ+=lDELTA[2][0];
}
ldRecClk[0]+=lDELTA[0][0]/VC;
if((fabs(lDELTA[0][0])<0.1))
break;
}
dRawRecClk[i] =ldRecClk[0];
for(j=0;j<num;j++)
{
if(!obsrec[i].eporec[j].iCycleFlag&&(obsrec[i].eporec[j].iAmbSmoTime<IAmbSmoLen))
{
int time=obsrec[i].eporec[j].iAmbSmoTime;
PriJ=SEAPRN(obsrec[i-1].epohead.SatList,obsrec[i].eporec[j].iPRN);
obsrec[i].eporec[j].lAmbiguity=(1.0-1.0/time)*obsrec[i-1].eporec[PriJ].lAmbiguity+obsrec[i].eporec[j].lAmbiguity/time;
}
}
//store the epoch receiver clock error to the dRawRecClk
if(!Par.ModelHolding)
{
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;
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;
}
outPosFile<<" "<<setprecision(12)<<setw(16)<<ldRecClk[0]<<endl;
//for(j=0;j<num;j++)
// outPosFile<<" "<<obsrec[i].eporec[j].iPRN<<" "<<obsrec[i].eporec[j].lPseudoRes<<" "<<obsrec[i].eporec[j].lPhaseRes<<" "<<obsrec[i].eporec[j].lAmbiguity<<endl;
TwoArrayFree(lDELTA);
TwoArrayFree(lA) ;
TwoArrayFree(lL0);
TwoArrayFree(lQx);
TwoArrayFree(lP) ;
delete[]recpos ;
delete[]ldRecClk ;
}
if(!Par.ModelHolding)
{ 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(PreClk.Makername)
{ int EpoPerGroup=10; //clock devided to groups
int ClkNum=TOTALOBSNUM/EpoPerGroup;
long double * lGroupClk=new long double[ClkNum];
Clk30SecTo5Min(outClkFile,dRawRecClk,lGroupClk,ClkNum,EpoPerGroup);
GETWCLKRMS(outClkFile,dRawRecClk,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;
}
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);
/************************************************************************/
delete[]iEpoNumofSat;
SatObsRecDesTroy(SatObsRec);
delete[]dRawRecClk;
delete[]dCorrecedRecclk;
delete[]satpos;
delete[]dSatClk;
delete[]gpst;
// delete[]SatObsRec; //have deleted before in SatObsRecDesTroy
delete[]obsrec;
delete[]PreClk.pDifToReference;
delete[]PreClk.pRecClk;
delete[]PreClk.pSavClk;
TwoArrayFree(pPosXYZDif);
inObsFile.close();
outPosFile.close();
outPosFile.close();
outClkFile.close();
return 1;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -