rinexnepoch.cpp
来自「计算摄动力模型中关于海潮固体潮极潮以及利用行星星历表计算太阳月亮位置的程序」· C++ 代码 · 共 516 行
CPP
516 行
// RinexNEpoch.cpp: implementation of the CRinexNEpoch class.
//
//////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include "compute.h"
#include "RinexNEpoch.h"
#include "math.h"
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
//--------读文件行--------------//
void fReadLn(FILE *fpa,char *s)
{ UINT i=0;
char ch;
while(!feof(fpa))
{ ch=fgetc(fpa);
if(ch==0xa||ch==EOF)
break;
s[i++]=ch;
}
s[i]=0;
}
//------判断行是否为空--------//
BOOL IsBlank(CString s)
{ BOOL b=TRUE;
int len;
int i;
len=s.GetLength();
for(i=0;i<len;i++)
{ if(s[i]!=' '&&s[i]!='\t'&&s[i]!='\n')
return FALSE;
}
return TRUE;
}
//-------------------定义卫星导航数据结构-------------------------//
CRinexNEpoch::CRinexNEpoch()
{
int i;
yr=0;
mon=0;
dy=0;
hr=0;
mi=0;
prn=0;
sec=0.0;
for (i=0;i<28;i++)
{
eph[i]=0;
}
}
//-------------------读卫星导航数据文件头-------------------------//
int fReadHd(FILE *fhd)
{ char s[81],ss[13];
int i;
int line=1;
fReadLn(fhd,s);
for(i=60;i<73;i++)
{
ss[i-60]=s[i];
}
while(strstr(ss,"END OF HEADER")==NULL)
{
line++;
fReadLn(fhd,s);
for(i=60;i<73;i++)
{
ss[i-60]=s[i];
}
}
return (line);
}
//-------------------转换时间为儒略日-------------------------//
void TimeJd(int year,int month,int day,int hour,
int minute,double second,double& jd,double& mjd)
{
int i,j,k,m;
float const c1=24.0,c2=60.0;
float t;
t=float(day);
j=2000+year;
m=month;
if(m<=2)
{
j-=1;
m+=12;
}
i=j/100;
k=2-i+i/4;
mjd=(365.25*j-fmod(365.25*j,1.0))-679006.0;
mjd+=int(30.6001*(m+1))+t+k/1.0;
mjd+=hour/c1+minute/c1/c2+second/c2/c2/c1;
jd=mjd+2400000.5;
}
//-------------------读导航数据-------------------------//
void CRinexNEpoch::ReadEph(FILE *fp1)
{
int j,n;
char buff[81];
fscanf(fp1,"%2d %2d %2d %2d %2d %2d %f",&prn,&yr,&mon,
&dy,&hr,&mi,&sec);
for (j=0;j<28;j++)
{
fscanf(fp1,"%lf%*c%d",&eph[j],&n);
eph[j]*=pow(10,n);
}
// fscanf(fp1,"%*lf%*c%*d%*lf%*c%*d%*lf%*c%*d");
fReadLn(fp1,buff);
}
int Com_Eph(CString* fName,int nFile,CString ComFile,CRinexNEpoch* epoch,
int& n_prn,CString& File_status)
{
char Filebuff[81],buff[81];
CRinexNEpoch epo,epo_sel;
int m,i,j,k,p,nlin,hlin;
int year,month,day,hour,minute;
double second,jd;//float second,jd;
FILE *fp2,*fp1;
double mjd,tim[neph],tim0;
bool FOUND=false;
p=0;m=0;
if((fp2=fopen(ComFile,"w"))==NULL)
{
sprintf(Filebuff,"打开 %s 文件错！",ComFile);
File_status=Filebuff;
return 1;
}
for(i=0;i<nFile;i++)
{
if((fp1=fopen(fName[i],"r"))==NULL)
{
sprintf(Filebuff,"打开 %s 文件错！",fName[i]);
File_status=Filebuff;
return 1;
}
hlin=fReadHd(fp1);
// printf("\nhead line=%d\n",hlin);
nlin=0;
while (!feof(fp1))
{
fReadLn(fp1,buff);
if (!IsBlank(buff)) nlin++;
}
nlin/=8;
// printf("%d\n",nlin);
rewind(fp1);
for (j=0;j<hlin;j++)
{
fReadLn(fp1,buff);
if (i==0) fprintf(fp2,"%s\n",buff);
}
if(p==0)
{
epoch[p].ReadEph(fp1);
epoch[p].GetTime(year,month,day,hour,minute,second);
TimeJd(year,month,day,0,0,0.0,jd,mjd);
tim[p]=mjd-(int(mjd/1000)*1000);
}
for (j=1;j<nlin;j++)
{
epo.ReadEph(fp1);
epo.GetTime(year,month,day,hour,minute,second);
TimeJd(year,month,day,0,0,0.0,jd,mjd);
tim0=mjd-(int(mjd/1000)*1000);
for(k=0;k<j;k++)
{
FOUND=EphComp(epo,tim0,epoch[k],tim[k]);
if(FOUND)
{
m++;
break;
}
}
if(!FOUND)
{
p++;
epoch[p]=epo;
epoch[p].GetTime(year,month,day,hour,minute,second);
TimeJd(year,month,day,0,0,0.0,jd,mjd);
tim[p]=mjd-(int(mjd/1000)*1000);
}
}
}
n_prn=p+1;
fclose(fp1);
for (i=0;i<n_prn;i++)
{
epoch[i].Output(fp2);
fprintf(fp2,"\n");
}
fclose(fp2);
if((fp2=fopen(ComFile,"r"))==NULL)
{
sprintf(Filebuff," %s 文件合并失败",ComFile);
File_status=Filebuff;
return 1;
}
return 0;
}
void Getstr(char* str)
{
int i,n;
char* str2;
char str1[]={"E"};
for (i=0;i<2;i++)
{
str[17+i]=str[18+i];
}
str[19]='\0';
str2=strstr(str,str1);
n=str2-str;
str[n]='D';
}
//--------输出文件----------------//
void CRinexNEpoch::Output(FILE *fp2)
{
int i,j;
float a=0.0;
char str[20];
fprintf(fp2,"%2d %2d %2d %2d %2d %2d %4.1f",prn,yr,mon,
dy,hr,mi,sec);
for (i=0;i<3;i++)
{
sprintf(str,"%15.12E",eph[i]);
Getstr(str);
fprintf(fp2,"%19s",str);
}
for(i=0;i<6;i++)
{
fprintf(fp2," \n ");
for(j=0;j<4;j++)
{
sprintf(str,"%15.12E",eph[i*4+j+3]);
Getstr(str);
fprintf(fp2,"%19s",str);
}
}
fprintf(fp2," \n ");
sprintf(str,"%15.12E",eph[27]);
Getstr(str);
fprintf(fp2,"%19s",str);
for (i=0;i<3;i++)
{
sprintf(str,"%15.12E",a);
Getstr(str);
fprintf(fp2,"%19s",str);
}
}
//-------获取时间------//
void CRinexNEpoch::GetTime(int& year,int& month,int& day,
int& hour,int& minute,double& second)
{
year=yr;
month=mon;
day=dy;
hour=hr;
minute=mi;
second=sec;
}
//-------获取toe时间-------//
double CRinexNEpoch::GetToe()
{
return (eph[11]);
}
bool EphSel(int prn_sel,int n_eph,CRinexNEpoch* epoch,double tmjd,
int& n1,int thr,int tmi,double tsec)
{
int i,n0,year,month,day,hour,minute;
double jd,second;//float jd,second;
double dt0,dt1,mjd,tmjd_i;
bool sat_match=false;
// if(fmod((thr+(tmi+tsec/60.0)/60.0),2)>1.0)
// thr=thr+1;
// if(fmod((thr+(tmi+tsec/60.0)/60.0),2)<1.0 && fmod((thr+(tmi+tsec/60.0)/60.0),2)!=0)
// thr=thr-1;
for (i=0;i<n_eph;i++)
{
if(prn_sel==epoch[i].prn)
{
epoch[i].GetTime(year,month,day,hour,minute,second);
TimeJd(year,month,day,0,0,0.0,jd,mjd);
tmjd_i=mjd-(int(mjd/1000)*1000);
dt0=fabs(tmjd*86400-(tmjd_i*86400+epoch[i].eph[11]));
n0=i;
n1=i;
sat_match=true;
//break;
//modify by liuhx
//if(thr==hour && tmi==minute && tsec==second) break;
//if(thr==hour) break;
if(fabs((thr*3600+tmi*60+tsec)-(hour*3600+minute*60+second))<=3600) break;
}
}
for (i=n0+1;i<n_eph;i++)
{
if(prn_sel==epoch[i].prn)
{
epoch[i].GetTime(year,month,day,hour,minute,second);
TimeJd(year,month,day,0,0,0.0,jd,mjd);
tmjd_i=mjd-(int(mjd/1000)*1000);
dt1=fabs(tmjd*86400-(tmjd_i*86400+epoch[i].eph[11]));
if (dt1<dt0)
{
dt0=dt1;
n1=i;
}
}
}
return (sat_match);
}
bool EphComp(CRinexNEpoch a,double tim1, CRinexNEpoch b,double tim2)
{
bool FOUND=false;
double fifmin;
tim1=tim1*86400+a.eph[11];
tim2=tim2*86400+b.eph[11];
fifmin=15.0*60.0+fZero;
if((a.prn==b.prn)&&((fabs(tim1-tim2))<fifmin))
FOUND=true;
return (FOUND);
}
//-------------------求反正切-------------------------//
void juqu(double &a,double ca,double sa)
{
double const pi=3.1415926535897932;
a=atan2(sa,ca);
}
//-------------------获取GPS周-------------------------//
void mjdgps(double mjd,int& week,double& trgps)
{
double dt,xx;
dt=mjd-44244.0;
xx=int(dt/7.0);
trgps=(dt-xx*7.0)*86400.0;
week=int(xx);
}
//-------------------计算卫星时间-------------------------//
void GetSat_Time(double& trgps,double& ts,double& fRange,CRinexNEpoch a)
{
double dt,delta_t;
dt=trgps-a.eph[11];
delta_t=a.eph[0]+a.eph[1]*dt+a.eph[2]*dt*dt;
fRange+=c*delta_t;
ts=trgps-delta_t-fRange/c;
}
CRinexNEpoch GetEph(double toe,int prn_sel,int n_prn,CRinexNEpoch* epoch)
{
int i,nsel;
for (i=0;i<n_prn;i++)
{
if (prn_sel==epoch[i].prn&&toe==epoch[i].eph[11])
{
nsel=i;
break;
}
}
return (epoch[nsel]);
}
//-------------------用卫星导航数据计算卫星的地心坐标-------------------------//
void svp(double trgps,CRinexNEpoch E,double* xs)
{
int k;
double i0,m0,ik,mk,a11;
double Crs,dn,Cuc,ec,Cus,sa,toe,Cic,Crc,w,omgd,di;
double a,n0,n,ts,tk,ek,q,v1,v2,v,omg0;
double fa,a1,a2,dta_uk,dta_rk,dta_ik,rk,uk,omgk;
double b1,b2,b3,b4,b5,Cis;
double const we=7.2921151467e-5,u0=3.986005e+14;
//将导航数据赋值进来
a11=E.eph[0];//lhx12.8
Crs=E.eph[4];
dn=E.eph[5];
m0=E.eph[6];
Cuc=E.eph[7];
ec=E.eph[8];
Cus=E.eph[9];
sa=E.eph[10];
toe=E.eph[11];
Cic=E.eph[12];
omg0=E.eph[13];
Cis=E.eph[14];
i0=E.eph[15];
Crc=E.eph[16];
w=E.eph[17];
omgd=E.eph[18];
di=E.eph[19];
//计算平均角速度n
a=sa*sa;
n0=sqrt(u0)/a/sa;
n=n0+dn;
//计算t时刻的平近点角mk和偏近点角ek
ts=trgps;
tk=ts-toe;//-0.075 why not ?
if (tk>302400.0) tk=tk-604800.0;
if(tk<-302400.0) tk=tk+604800.0;
mk=m0+n*tk;
ek=mk;
for (k=0;k<9;k++)
{
ek=mk+ec*sin(ek);
}
q=1.0-ec*cos(ek);
//观测时刻t的真近点角v和升交距角fa
v1=(cos(ek)-ec)/q;
v2=sin(ek)*sqrt(1.0-ec*ec)/q;
v=fabs(asin(v2));
juqu(v,v1,v2);
fa=w+v;
//摄动改正项
a1=cos(2.0*fa);
a2=sin(2.0*fa);
dta_uk=Cuc*a1+Cus*a2;
dta_rk=Crc*a1+Crs*a2;
dta_ik=Cic*a1+Cis*a2;
//观测时刻的升交距角uk、卫星地心距rk、轨道倾角ik、
rk=a*q+dta_rk;
uk=fa+dta_uk;
ik=i0+di*tk+dta_ik;
//升交点经度omgk
omgk=omg0+(omgd-we)*tk-we*toe;
//卫星在协议地球坐标系中的空间直角坐标
b1=cos(omgk);
b2=sin(omgk);
b3=cos(uk);
b4=sin(uk);
b5=cos(ik);
xs[0]=(b1*b3-b2*b4*b5)*rk;
xs[1]=(b2*b3+b1*b4*b5)*rk;
xs[2]=b4*sin(ik)*rk;
}
//------------主程序调用计算程序-----------------//
void svp_Calc(CString* fName,int nFile,CString ComFile,
int tyr,int tmon,int tdy,int thr,int tmi,
double tsec,int sv_sel,double fRange,
double* xs,CString& File_status,double& toe,
bool base)
{
CRinexNEpoch epoch[neph],epo,epo_sel;
int eph_i,n_prn,week;
char buff[81];
double jd;//float jd;
double mjd,mjd0,trgps,ts;
bool FOUND=false;
bool sat_match=false;
Com_Eph(fName,nFile,ComFile,epoch,n_prn,File_status);
TimeJd(tyr,tmon,tdy,thr,tmi,tsec,jd,mjd);
mjd0=mjd-(int(mjd/1000)*1000);
sat_match=EphSel(sv_sel,n_prn,epoch,mjd0,eph_i,thr,tmi,tsec);
if (sat_match)
{
if (base)
epo_sel=epoch[eph_i];
else
epo_sel=GetEph(toe,sv_sel,n_prn,epoch);
mjdgps(mjd,week,trgps);
GetSat_Time(trgps,ts,fRange,epo_sel);
svp(ts,epo_sel,xs);
toe=epo_sel.GetToe();
}
else
{
sprintf(buff,"没找到匹配的卫星！");
File_status=buff;
}
}
//void main()
//{
// int nFile;
// int tyr,tmon,tdy,sv_sel,thr,tmi;
// float tsec;
// double xs[3],fRange,toe;
// bool base=true;
// CString File_status;
//// CString fName[]={"e:\\cheng\\gpscalc\\shao3260.98n","e:\\cheng\\gpscalc\\taej3260.98n"};
// CString fName[]={"c:\\00091081.00N"};
// CString ComFile="c:\\NAV1";
// tyr=99,tmon=8,tdy=24,sv_sel=13,thr=9,tmi=0,tsec=0.0;
// nFile=1;
// svp_Calc(fName,nFile,ComFile,tyr,tmon,tdy,thr,tmi,
// tsec,sv_sel,fRange,xs,File_status,toe,base);
//}
CRinexNEpoch::~CRinexNEpoch()
{
}
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?