搜索
GPS卫星位置的计算利用C++编写了一段能计算单一瞬时卫星坐标的程序，在运行程序之前，需做部分准备工作：（1）在F盘下建立一名为“单一卫星广播星历”的txt文件。（2）从“广播星历.txt”文件中拷贝从卫星PRN号开始的8行数据到“单一卫星广播星历.txt”中（3）在编辑选项中，将全部的“D”替换为“E”。下面为我所选取的一个广播星历：1806825600.0-2.472363412380E-04-1.023181539495E-120.000000000000E+001.410000000000E+02-1.721875000000E+014.502687555010E-091.413760604187E+00-7.990747690201E-077.598234573379E-031.118145883083E-055.153709835052E+034.536000000000E+05-1.303851604462E-08-1.095067942661E-011.527369022369E-079.571235745530E-011.640000000000E+02-2.656176299285E+00-8.037477650349E-09-5.193073455211E-101.000000000000E+001.389000000000E+030.000000000000E+002.000000000000E+000.000000000000E+00-1.024454832077E-081.410000000000E+024.464490000000E+054.000000000000E+00程序设计部分：#includestdio.h#includemath.hintmain(){inti=0;doublen[50],n0,nn,t,tk,Mk,Ek,Vk,Yk,Gu,Gr,Gi,uk,rk,ik,xk,yk,zk,X,Y,Z,Lk,UT,yy,mm,JD,gpsz;FILE*fp;fp=fopen(F:\\单一卫星广播星历.txt,r);if(fp==NULL){printf(文件打开失败!\n);return0;}while(!feof(fp)){fscanf(fp,%lf,&n[i]);i++;}n0=(sqrt(3986005E+8))/pow(n[17],3);nn=n0+n[12];/*计算卫星运行的平均角速度*/UT=n[4]+(n[5]/60)+(n[7]/3600);/*民用日的时分秒化为实数时*/if(n[1]=80)/*广播星历中年只有后两位，化为4位，参考1980年1月6日0点*/{if(n[1]==80&&n[2]==1&&n[3]6){n[1]=n[1]+2000;}n[1]=n[1]+1900;}else{n[1]=n[1]+2000;}if(n[2]=2){yy=n[1]-1;mm=n[2]+12;}if(n[2]2){yy=n[1];mm=n[2];}JD=(int)(365.25*yy)+(int)(30.6001*(mm+1))+n[3]+(UT/24)+1720981.5;/*化为儒略日*/gpsz=(int)((JD-2444244.5)/7);/*计算GPS周*/t=(JD-2444244.5-7*gpsz)*24*3600;/*得出GPS秒*/tk=t-n[18];/*tk1为中间值，用以判断tk与正负302400的关系，然后返回到tk上*/while(tk302400||tk-302400){if(tk302400){tk=tk-604800;}else{tk=tk+604800;}}/*计算归化观测时间*/Mk=n[13]+nn*tk;/*观测时刻的卫星平近点角*/Ek=Mk;Ek=Mk+n[15]*sin(Ek);Ek=Mk+n[15]*sin(Ek);/*迭代两次计算观测时刻的偏近点角*/Vk=atan(sqrt(1-n[15]*n[15])*sin(Ek))/(cos(Ek)-n[15]);/*真近点角*/Yk=Vk+n[24];/*升交距角*/Gu=n[14]*cos(2*Yk)+n[16]*sin(2*Yk);Gr=n[23]*cos(2*Yk)+n[11]*sin(2*Yk);Gi=n[19]*cos(2*Yk)+n[21]*sin(2*Yk);/*摄动改正项*/uk=Yk+Gu;rk=n[17]*n[17]*(1-n[15]*cos(Ek))+Gr;ik=n[22]+Gi+n[26]*tk;/*经摄动改正后的升交距角、卫星矢径、轨道倾角*/xk=rk*cos(uk);yk=rk*sin(uk);zk=0;/*卫星在轨道坐标系的坐标*/Lk=n[20]+(n[25]-7.29211515E-5)*tk-7.29211515E-5*n[18];/*观测时刻t的升交点经度*/X=xk*cos(Lk)-yk*cos(ik)*sin(Lk);Y=xk*sin(Lk)+yk*cos(ik)*cos(Lk);Z=yk*sin(ik);/*卫星在WGS-84坐标系的坐标*/printf(该卫星在WGS-84坐标系中的坐标为：\nX=%lfm\nY=%lfm\nZ=%lfm\n,X,Y,Z);fclose(fp);return0;}计算结果：该卫星在WGS-84坐标系中的坐标为：X=9223153.692525mY=24133486.931401mZ=6032585.919385m