卫星坐标计算

#include "stdafx.h"
#include
#include
#include
#include
#include
#include
using namespace std;
double power(double x,double n)
{
double val=1.0;
while(n--)
val*=x;
return val;
}
const double PI=3.1415926;
const double we=7.29211567/power(10,5);
const double GM=3.986005*power(10,14);

class GPS
{
public:
//Mo为0时刻的平近点角，e为GPS卫星轨道偏心率，an为导航电文给出的摄动改正数dn，no为平均角速度，toe参考时刻，Ek偏近点角，Mk平近点角，
//Vk真近点角，Qk升交距角，w卫星电文给出的近地点角距
double dn,Mo,e,w,Cuc,Cus,Crc,Crs,Cis,Cic,Wo,W,n0;
double n,Ek,Ek0,Mk,Vk,Ok,du,di,dr,uk,rk,ik,xk,yk,zk,Lk,Xk,Yk,Zk,Qk,a,UT,io,wo,I,JD,tk,toe,Idot,t;
int w_n ,Y,M,D,H,min;
double sec,a1,a2,a0,aode,L2,GPSn_week,L2P,wxjd,MSB,TGD,IODC,sendtime,wm,ek1,ek2,cosVk,sinVk;
double TIME(int nYear,int nMonth,int nDay,int nHour,int nMinute,double dSecond);
void mult(GPS &p );
};
double GPS::TIME(int nYear,int nMonth,int nDay,int nHour,int nMinute,double dSecond)
{
//计算时间，有GPS接收时间转化成儒略日
double m_dWeekSecond;
int DayofMonth=0,DayofYear=0,weekno=0,dayofWeek,m;
if(nYear<1980||nMonth<1||nMonth>12||nDay>31)return -1;
for(m=1980;m{
if((m%4==0 && m%100!=0) || m%400==0)
DayofYear+=366;
else
DayofYear+=365;
}//计算从1980年到当前的前一年的天数
//
for(m=0;m{
if(m==1||m==3||m==5||m==7||m==8||m==10||m==12)
DayofMonth+=31;
else if(m==4||m==6||m==9||m==11)
DayofMonth+=30;
else if(m==2)
{
if((nYear%4==0&&nYear%100!=0)||(nYear%400==0))
DayofMonth+=29;
else
DayofMonth+=28;
}
}//计算当前一年元月到当前前一月的天数
DayofMonth=DayofMonth+nDay-6;//加上当月天数减去1980年元月6日
weekno=(DayofMonth+DayofYear)/7;//计算GPS周
dayofWeek=(DayofMonth+DayofYear)%7;//计算GPS周秒时间
m_dWeekSecond=dayofWeek*86400+nHour*3600+nMinute*60+dSecond;
return m_dWeekSecond;
}
void GPS::mult(GPS &p)
{
GPS hhd[175];
fstream in;
in.open("b3062240.11n");
ofstream fout;
fout.open("结果文件.txt");
string str;
while(getline(in,str))
{
string::size_type pos=str.find("ENE OF HEAEER");
if(pos!=string::npos)break;
}
int j=0;
for(j=0;j<175;j++)
{
//do
//{
in>>hhd[j].w_n>>hhd[j].Y>>hhd[j].M>>hhd[j].D>>hhd[j].H>>hhd[j].min>>hhd[j].sec
>>hhd[j].a0>>hhd[j].a1>>hhd[j].a2>>hhd[j].aode>>hhd[j].Crs>>hhd[j].dn>>hhd[j].Mo>>hhd[j].Cuc>>hhd[j].e>>hhd[j].Cus>>hhd[j].a>>
hhd[j].toe>>hhd[j].Cic>>hhd[j].Wo>>hhd[j].Cis>>hhd[j].I>>hhd[j].Crc>>hhd[j].w>>hhd[j].W>>hhd[j].Idot>>hhd[j].L2>>hhd[j].GPSn_week>>
hhd[j].L2P>>hhd[j].wxjd>>hhd[j].MSB>>hhd[j].TGD>>hhd[j].IODC>>hhd[j].sendtime;
if(hhd[j].Y>=80)
hhd[j].Y=1900+hhd[j].Y;
else if(Y<80&&Y>=0)
hhd[j].Y=2000+hhd[j].Y;
hhd[j].t=hhd[j].TIME(hhd[j].Y,hhd[j].M,hhd[j].D,h

hd[j].H,hhd[j].min,hhd[j].sec);

hhd[j].tk=hhd[j].t-hhd[j].toe;//计算归化时间
if(hhd[j].tk>302400)
hhd[j].tk-=604800;
if(hhd[j].tk<-302400)
hhd[j].tk+=604800;
hhd[j].n0=sqrt(GM)/pow(hhd[j].a,3);
hhd[j].n=hhd[j].n0+hhd[j].dn;
hhd[j].Mk=hhd[j].Mo+hhd[j].n*hhd[j].tk;
hhd[j].Ek0=hhd[j].Mk;
hhd[j].ek1,hhd[j].ek2;
hhd[j].ek1=hhd[j].Mk+hhd[j].e*sin(hhd[j].Ek0);
hhd[j].ek2=hhd[j].Mk+hhd[j].e*sin(hhd[j].ek1);
hhd[j].Ek=hhd[j].Mk+hhd[j].e*sin(hhd[j].ek2);

//计算观测时刻的偏近点角Ek
hhd[j].cosVk,hhd[j].sinVk;
hhd[j].cosVk=(cos(hhd[j].Ek)-hhd[j].e)/(1-hhd[j].e*cos(hhd[j].Ek));
hhd[j].sinVk=sqrt(1-hhd[j].e*hhd[j].e)*sin(hhd[j].Ek)/(1-hhd[j].e*cos(hhd[j].Ek));
if(hhd[j].cosVk>0 && hhd[j].sinVk>0 )
hhd[j].Vk=atan(hhd[j].sinVk/hhd[j].cosVk);
else if(hhd[j].sinVk>0 && hhd[j].cosVk<0)
hhd[j].Vk=acos(hhd[j].cosVk);
else if(hhd[j].sinVk<0 && hhd[j].cosVk<0)
hhd[j].Vk=PI-asin(hhd[j].sinVk);
else
hhd[j].Vk=2*PI-acos(hhd[j].cosVk);
/*Vk=atan((sqrt(1-e*e)*sin(Ek))/(cos(Ek)-e));*///计算真近点角Vk
hhd[j].Qk=hhd[j].Vk+hhd[j].w; //计算升交距角Qk
hhd[j].du=hhd[j].Cuc*cos(2*hhd[j].Qk)+hhd[j].Cus*sin(2*hhd[j].Qk);
hhd[j].dr=hhd[j].Crc*cos(2*hhd[j].Qk)+hhd[j].Crs*sin(2*hhd[j].Qk);
hhd[j].di=hhd[j].Cic*cos(2*hhd[j].Qk)+hhd[j].Cis*sin(2*hhd[j].Qk);//计算摄动改正项ru,rr,ri
hhd[j].uk=hhd[j].Qk+hhd[j].du;
hhd[j].rk=(1-hhd[j].e*cos(hhd[j].Ek))*hhd[j].a*hhd[j].a+hhd[j].dr;
hhd[j].ik=hhd[j].I+hhd[j].di+hhd[j].Idot*hhd[j].tk; //计算摄动改正的升交距角uk，卫星向径rk和轨道倾角ik
hhd[j].xk=hhd[j].rk*cos(hhd[j].uk);
hhd[j].yk=hhd[j].rk*sin(hhd[j].uk); //计算卫星的轨道坐标系的位置
hhd[j].Lk=hhd[j].Wo+(hhd[j].W-we)*hhd[j].tk-we*hhd[j].toe; //计算t时刻升交点经度Lk

hhd[j].Xk=hhd[j].xk*cos(hhd[j].Lk)-hhd[j].yk*cos(hhd[j].ik)*sin(hhd[j].Lk);
hhd[j].Yk=hhd[j].xk*sin(hhd[j].Lk)+(hhd[j].Lk)*cos(hhd[j].ik)*cos(hhd[j].Lk);
hhd[j].Zk=hhd[j].yk*sin(hhd[j].ik); //计算卫星在WGS-84坐标系的位置

fout<<"卫星"<fout<<"卫星平均角速度n0="<fout<<"卫星改化平均角速度n="<fout<<"归化时间tk="<fout<<"摄动改正项："<fout<<"du="<fout<<"xk="<fout<<"观测时刻卫星平近点角Mk："<fout<<"观测时刻卫星偏近点角Ek："<fout<<"观测时刻卫星

真近点角Vk："<fout<<"观测时刻卫星升交距角Qk："<fout<<"摄动改正后升交距角uk，卫星矢径rk，轨道倾角ik："<fout<<"uk="<fout<<"Xk="<fout<fout<fout<<"=============================================================================================="<}
}
/*}while(j<175);*/
void main()
{
fstream fin("b3062240.11n", ios::in|ios::out);
string str;
string strKey = "D";
string strReplace = "E";
string text;
string line;
int p;
while(getline(fin, line))
{
while((p=line.find(strKey.c_str()))!=-1)
{
line.replace(p, strKey.length(), strReplace.c_str());
}
text+=line+"\n";
}
fin.clear();
fin.sync();
fin.seekp(0);
fin<
GPS G;
GPS hhd[175];
G.mult(*hhd);

}