📄 preclk.cpp
字号:
#include "stdAfx.h"
#include "PreClk.h"
#include "math.h"
#include "errorcorrection.h"
#include "matrixcsl.h"
#include "iomanip.h"
int READCLKFILE(ifstream inClkFile,char *MarkerName,PRECLK &PreClk)
{ char str[81];
CString s;
char *p;
PreClk.Makername=MarkerName;
int realsatnum;
do
{
inClkFile.getline(str,82,'\n');
if(p=strstr(str,"COMMENT"))
continue;
if(p=strstr(str,"END OF HEADER "))
break;
if(p=strstr(str,MarkerName))
{ s=str;
PreClk.RecPos[0]=atof(s.Mid(24,12))/1000;
PreClk.RecPos[1]=atof(s.Mid(36,12))/1000;
PreClk.RecPos[2]=atof(s.Mid(48,12))/1000;
}
if(p=strstr(str,"SOLN SATS"))
{
s=str;
realsatnum=atoi(s.Mid(0));
}
}while(1);
int i,j;
for (i=0;i<MAXCLKRECORDS;i++)
{ bool EpochSucceed=false;
while(1)
{ inClkFile.getline(str,82,'\n');
if(p=strstr(str,MarkerName))
{ s=str;
PreClk.pRecClk[i].UTCT.iYEAR=atoi(s.Mid(7));
PreClk.pRecClk[i].UTCT.iMONTH=atoi(s.Mid(12));
PreClk.pRecClk[i].UTCT.iDAY=atoi(s.Mid(15));
PreClk.pRecClk[i].UTCT.iHOUR=atoi(s.Mid(18));
PreClk.pRecClk[i].UTCT.iMINUTE=atoi(s.Mid(21));
PreClk.pRecClk[i].UTCT.dSECOND=atof(s.Mid(24));
PreClk.pRecClk[i].lClk[0]=atof(s.Mid(39));
PreClk.pRecClk[i].lClk[1]=atof(s.Mid(59));
}
if(p=strstr(str,"AS G01"))
{ s=str;
int iPrn=atoi(s.Mid(4));
PreClk.pSavClk[i].UTCT[iPrn-1]=PreClk.pRecClk[i].UTCT;
PreClk.pSavClk[i].clk[iPrn-1][0]=atof(s.Mid(39));
PreClk.pSavClk[i].clk[iPrn-1][1]=atof(s.Mid(59));
for(j=1;j<realsatnum;j++)
{ inClkFile.getline(str,82,'\n');
s=str;
int iPrn=atoi(s.Mid(4));
if(!iPrn)
{ printf(" The clk file is damaged for absence of satellite clock record.\n");
break;
}
PreClk.pSavClk[i].UTCT[iPrn-1]=PreClk.pRecClk[i].UTCT;
PreClk.pSavClk[i].clk[iPrn-1][0]=atof(s.Mid(39));
PreClk.pSavClk[i].clk[iPrn-1][1]=atof(s.Mid(59));
}
}
if(p=strstr(str,"AR GPST"))
{
s=str;
PreClk.pDifToReference[i].UTCT=PreClk.pRecClk[i].UTCT;
PreClk.pDifToReference[i].lClk[0]=atof(s.Mid(39));
EpochSucceed=TRUE;
}
if(EpochSucceed)
break;
}
}
return 1;
}
int READCODCLKFILE(ifstream inClkFile,char *MarkerName,PRECLK &PreClk)
{ char str[81];
CString s;
char *p;
PreClk.Makername=NULL;
int realsatnum=0;
do
{
inClkFile.getline(str,82,'\n');
if(p=strstr(str,"COMMENT"))
continue;
if(p=strstr(str,"END OF HEADER "))
break;
if(p=strstr(str,MarkerName))
{ PreClk.Makername= MarkerName;
s=str;
PreClk.RecPos[0]=atof(s.Mid(24,12))/1000;
PreClk.RecPos[1]=atof(s.Mid(36,12))/1000;
PreClk.RecPos[2]=atof(s.Mid(48,12))/1000;
}
if(p=strstr(str,"SOLN SATS"))
{
s=str;
realsatnum=atoi(s.Mid(0));
}
}while(1);
int i,j;
for (i=0;i<MAXCLKRECORDS*10;i++)
{ bool EpochSucceed=false;
while(1)
{ inClkFile.getline(str,82,'\n');
if(p=strstr(str,"AR AMC2"))
{
s=str;
PreClk.pDifToReference[i].UTCT.iYEAR=atoi(s.Mid(7));
PreClk.pDifToReference[i].UTCT.iMONTH=atoi(s.Mid(12));
PreClk.pDifToReference[i].UTCT.iDAY=atoi(s.Mid(15));
PreClk.pDifToReference[i].UTCT.iHOUR=atoi(s.Mid(18));
PreClk.pDifToReference[i].UTCT.iMINUTE=atoi(s.Mid(21));
PreClk.pDifToReference[i].UTCT.dSECOND=atof(s.Mid(24));
PreClk.pDifToReference[i].lClk[0]=atof(s.Mid(39));
continue;
}
if(((i/10)*10==i)&&strstr(str,MarkerName))
{ s=str;
PreClk.pRecClk[i/10].UTCT=PreClk.pDifToReference[i].UTCT;
PreClk.pRecClk[i/10].lClk[0]=atof(s.Mid(39));
PreClk.pRecClk[i/10].lClk[1]=atof(s.Mid(59));
continue;
}
if(p=strstr(str,"AS G01"))
{ s=str;
int iPrn=atoi(s.Mid(4));
PreClk.pSavClk[i].UTCT[iPrn-1]=PreClk.pDifToReference[i].UTCT;
PreClk.pSavClk[i].clk[iPrn-1][0]=atof(s.Mid(39));
PreClk.pSavClk[i].clk[iPrn-1][1]=atof(s.Mid(59));
for(j=1;j<realsatnum;j++)
{ inClkFile.getline(str,82,'\n');
s=str;
int iPrn=atoi(s.Mid(4));
if(!iPrn)
{
printf(" The clk file is damaged for absence of satellite clock record.\n");
}
PreClk.pSavClk[i].UTCT[iPrn-1]=PreClk.pDifToReference[i].UTCT;
PreClk.pSavClk[i].clk[iPrn-1][0]=atof(s.Mid(39));
PreClk.pSavClk[i].clk[iPrn-1][1]=atof(s.Mid(59));
}
EpochSucceed=TRUE;
}
if(EpochSucceed)
break;
}
}
return 1;
}
int GETWCLKRMS(ofstream pFile,long double *lObs,SIGCLK *pRecClk,
SIGCLK *pDifToReference,int n)
{ long double VTV=0;
long double average = 0;
pFile<<" The clock figured out by TRIP-T relative to true clock from the IGS Precise CLK file ："<<endl;
pFile<<" No CLK By TRIP-T CLK From IGS The Difference "<<endl;
for(int i=0;i<n;i++)
{
long double lDif=lObs[i]-(pRecClk[i].lClk[0]-pDifToReference[i].lClk[0]);
if(pRecClk[i].lClk[0]>1)
lDif=0;
// to avoid the case of +indefinitive eg. 4.2e+66
pFile<<setw(4)<<i+1<<" "<<setw(15)<<lObs[i]<<" ";
pFile<<setw(15)<<pRecClk[i].lClk[0]-pDifToReference[i].lClk[0]<<" "<<setw(15)<<lDif<<endl;
VTV+=lDif*lDif;
average+=lDif;
}
VTV/=n;
long double RMS=pow(VTV,0.5);
pFile<<" RMS is："<<RMS<<endl;
average/=n;
for(i=0;i<n;i++)
{ long double lDif=lObs[i]-(pRecClk[i].lClk[0]-pDifToReference[i].lClk[0])-average;
if(pRecClk[i].lClk[0]>1)
lDif=0;
// to avoid the case of +indefinitive eg. 4.2e+66
VTV+=lDif*lDif;
}
VTV/=(n-1);
RMS=pow(VTV,0.5);
pFile<<"偏差均值为："<<average<<" 去均值 RMS is："<<RMS<<endl;
return 1;
}
int GETCLKRMS(ofstream pFile,long double *lObs,int n)
{//Note: lObs is the vector of computed(corrected) observations
int i,j;
long double RMS=0;
long double *lEvalu=new long double[n];
//the average be seemed to real values
// long double lEverge=0;
// for(i=0;i<n;i++)
// lEverge+=lObs[i];
// lEverge/=n;
// for(i=0;i<n;i++)
// lEvalu[i]=lEverge;
pFile.flags(ios::scientific);
int order=6; //the order of polynomial;
long double **plB=TwoArrayAlloc(n,order);
long double **plL=TwoArrayAlloc(n,1);
long double **plX=TwoArrayAlloc(order,1);
for (i=0;i<n;i++)
{ for(j=0;j<order;j++)
plB[i][j]=pow(i*30,j);
plL[i][0]=lObs[i];
}
LEASTSQUAREAJUST(plX,plB,n,order,plL);
pFile<<" The clock fitted by polynomial are："<<endl;
pFile<<" Raw value Fitting value Difference \n";
for (i=0;i<n;i++)
{ lEvalu[i]=0;
for(j=0;j<order;j++)
lEvalu[i]+=plX[j][0]*plB[i][j];
}
for(i=0;i<n;i++)
{
long double ldif=lObs[i]-lEvalu[i];
pFile<<setprecision(12)<<setw(15)<<lObs[i]<<" "<<setprecision(12)<<setw(15)<<lEvalu[i]<<" "<<setprecision(12)<<setw(15)<<ldif<<endl;
RMS+=ldif*ldif;
}
RMS/=n;
RMS=pow(RMS,0.5);
pFile<<" The clock RMS relative to its fitting value\n according to the clock model is ： "<<RMS<<endl;
delete []lEvalu;
return 1;
}
int GETPOSRMS(ofstream pFile,long double *lV,int n)
{ long double RMS,RMSSQUARE=0;
int i;
for(i=0;i<n;i++)
RMSSQUARE+=lV[i]*lV[i];
RMSSQUARE/=n;
RMS=pow(RMSSQUARE,0.5);
pFile<<" RMS value is："<<RMS<<endl;
return 1;
}
int Clk30SecTo5Min(ofstream pFile,long double *lInput,long double *lOutput,
int outputNum,int groups) //least square
{ int i=0,j=0,k=0;
pFile<<" The receiver clock with 5 min interval："<<endl;
for(i=0;i<outputNum;i++)
{ int order=2; //the order of polynomial;
lOutput[i]=0;
long double **plB=TwoArrayAlloc(groups,order);
long double **plL=TwoArrayAlloc(groups,1);
long double **plX=TwoArrayAlloc(order,1);
for(j=0;j<groups;j++)
{
for(k=0;k<order;k++)
plB[j][k]=pow(j*30,k); //units as second
plL[j][0]=lInput[groups*i+j];
}
pFile.flags(ios::fixed);
LEASTSQUAREAJUST(plX,plB,groups,order,plL);
for(k=0;k<order;k++)
lOutput[i]+=plX[k][0]*pow(0*30,k);
pFile<<setw(3)<<i+1<<" "<<setw(15)<<setprecision(12)<<lOutput[i]<<endl;
}
return 1;
}
//int Clk30SecToMin(long double *lInput,long double *lOutput,int outputNum,
// int groups) //average the clk
//{ int i=0,j=0;
// for(i=0;i<outputNum;i++)
// {
// lOutput[i]=0;
// for(j=0;j<groups;j++)
// lOutput[i]+=lInput[groups*i+j];
// lOutput[i]/=groups;
// cout<<setw(12)<<lOutput[i]<<endl;
// }
// return 1;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -