卫星坐标计算

#includestdafx.h#includeiostream#includestdio.h#includemath.h#includeiomanip#includefstream#includestringusingnamespacestd;doublepower(doublex,doublen){doubleval=1.0;while(n--)val*=x;returnval;}constdoublePI=3.1415926;constdoublewe=7.29211567/power(10,5);constdoubleGM=3.986005*power(10,14);classGPS{public://Mo为0时刻的平近点角，e为GPS卫星轨道偏心率，an为导航电文给出的摄动改正数dn，no为平均角速度，toe参考时刻，Ek偏近点角，Mk平近点角，//Vk真近点角，Qk升交距角，w卫星电文给出的近地点角距doubledn,Mo,e,w,Cuc,Cus,Crc,Crs,Cis,Cic,Wo,W,n0;doublen,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;intw_n,Y,M,D,H,min;doublesec,a1,a2,a0,aode,L2,GPSn_week,L2P,wxjd,MSB,TGD,IODC,sendtime,wm,ek1,ek2,cosVk,sinVk;doubleTIME(intnYear,intnMonth,intnDay,intnHour,intnMinute,doubledSecond);voidmult(GPS&p);};doubleGPS::TIME(intnYear,intnMonth,intnDay,intnHour,intnMinute,doubledSecond){//计算时间，有GPS接收时间转化成儒略日doublem_dWeekSecond;intDayofMonth=0,DayofYear=0,weekno=0,dayofWeek,m;if(nYear1980||nMonth1||nMonth12||nDay31)return-1;for(m=1980;mnYear;m++){if((m%4==0&&m%100!=0)||m%400==0)DayofYear+=366;elseDayofYear+=365;}//计算从1980年到当前的前一年的天数//for(m=0;mnMonth;m++){if(m==1||m==3||m==5||m==7||m==8||m==10||m==12)DayofMonth+=31;elseif(m==4||m==6||m==9||m==11)DayofMonth+=30;elseif(m==2){if((nYear%4==0&&nYear%100!=0)||(nYear%400==0))DayofMonth+=29;elseDayofMonth+=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;returnm_dWeekSecond;}voidGPS::mult(GPS&p){GPShhd[175];fstreamin;in.open(b3062240.11n);ofstreamfout;fout.open(结果文件.txt);stringstr;while(getline(in,str)){string::size_typepos=str.find(ENEOFHEAEER);if(pos!=string::npos)break;}intj=0;for(j=0;j175;j++){//do//{inhhd[j].w_nhhd[j].Yhhd[j].Mhhd[j].Dhhd[j].Hhhd[j].minhhd[j].sechhd[j].a0hhd[j].a1hhd[j].a2hhd[j].aodehhd[j].Crshhd[j].dnhhd[j].Mohhd[j].Cuchhd[j].ehhd[j].Cushhd[j].ahhd[j].toehhd[j].Cichhd[j].Wohhd[j].Cishhd[j].Ihhd[j].Crchhd[j].whhd[j].Whhd[j].Idothhd[j].L2hhd[j].GPSn_weekhhd[j].L2Phhd[j].wxjdhhd[j].MSBhhd[j].TGDhhd[j].IODChhd[j].sendtime;if(hhd[j].Y=80)hhd[j].Y=1900+hhd[j].Y;elseif(Y80&&Y=0)hhd[j].Y=2000+hhd[j].Y;hhd[j].t=hhd[j].TIME(hhd[j].Y,hhd[j].M,hhd[j].D,hhd[j].H,hhd[j].min,hhd[j].sec);hhd[j].tk=hhd[j].t-hhd[j].toe;//计算归化时间if(hhd[j].tk302400)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);//计算观测时刻的偏近点角Ekhhd[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].cosVk0&&hhd[j].sinVk0)hhd[j].Vk=atan(hhd[j].sinVk/hhd[j].cosVk);elseif(hhd[j].sinVk0&&hhd[j].cosVk0)hhd[j].Vk=acos(hhd[j].cosVk);elseif(hhd[j].sinVk0&&hhd[j].cosVk0)hhd[j].Vk=PI-asin(hhd[j].sinVk);elsehhd[j].Vk=2*PI-acos(hhd[j].cosVk);/*Vk=atan((sqrt(1-e*e)*sin(Ek))/(cos(Ek)-e));*///计算真近点角Vkhhd[j].Qk=hhd[j].Vk+hhd[j].w;//计算升交距角Qkhhd[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,rihhd[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和轨道倾角ikhhd[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时刻升交点经度Lkhhd[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卫星hhd[j].w_n在hhd[j].Y年hhd[j].M月hhd[j].D日hhd[j].H+8时hhd[j].min分hhd[j].sec秒的相关参数：endl;fout卫星平均角速度n0=hhd[j].n0fixedendl;fout卫星改化平均角速度n=hhd[j].nendl;fout归化时间tk=hhd[j].tkendl;fout摄动改正项：endl;foutdu=hhd[j].du;dr=hhd[j].dr;di=hhd[j].diendl;fout卫星在轨道平面坐标系的坐标：endl;foutxk=hhd[j].xk;yk=hhd[j].ykendl;fout观测时刻升交点经度Lk：hhd[j].Lkendl;fout观测时刻卫星平近点角Mk：hhd[j].Mkendl;fout观测时刻卫星偏近点角Ek：hhd[j].Ekendl;fout观测时刻卫星真近点角Vk：hhd[j].Vkendl;fout观测时刻卫星升交距角Qk：hhd[j].Qkendl;fout摄动改正后升交距角uk，卫星矢径rk，轨道倾角ik：endl;foutuk=hhd[j].uk;rk=hhd[j].rk;ik=hhd[j].ikendl;fout卫星坐标为：endl;foutXk=hhd[j].Xk;Yk=hhd[j].Yk;Zk=hhd[j].Zkendl;fout.unsetf(ios::fixed);fout==============================================================================================endl;foutendl;foutendl;fout==========================