calpos.asv

包含三个M文件

function[x,y,z]= CalPos(PRN,time)
global EphDat
t1=TimetoJD(time(1),time(2),time(3),time(4),time(5),time(6));
t2=TimetoJD(time(1),time(2),time(3),0,0,0);
tic     %start
EphDat=ReadGpsData;
sz=size(EphDat,2);
gg=0;
for i=1:sz
    if(EphDat(i).PRN==PRN & abs(EphDat(i).toc-t1)<(1/24))
        G=EphDat(i);
        gg=1;
        break;
    end
end
if (gg==0)
    msgbox('未读到符合要求的星历文件','warning','warn');
    return;
end
t3=t2-G.weekno*7;
% 减整周 数
t=t3*24*3600+time(4)*3600+time(5)*60+time(6);
%化为GPS秒
u=3.986004418E+14;   %地球引力常数
we=7.292115E-5;      %地球自转速度
a=G.sqa^2;           %地球长半轴
n0=sqrt(u/a^3);      %计算的平运动  
tk=t-G.toe;          %从参考历元开始的时间 
if tk>=302400
    tk=tk-604800;
elseif tk<-302400
    tk=tk + 604800;
end
n=n0+G.dn;           %改正后的运动 
Mk=G.M0+n*tk;
Ek=Mk;                                         %偏近点角  弧度
for i=1:2
    Ek=Mk+G.e*sin(Ek);
end
fk=2*atan((sqrt((1+G.e)/(1-G.e))*tan(Ek/2)));   %真近点角
if(fk<0)
    fk=fk+2*pi;
end
cosf=(cos(Ek)-G.e)/(1-G.e*cos(Ek));
sinf=(sqrt(1-G.e^2)*sin(Ek))/(1-G.e*cos(Ek));
uk=fk+G.w;                                       %升交角矩及轨道摄动改正项
iuk=G.Cuc*cos(2*uk)+G.Cus*sin(2*uk);
irk=G.Crc*cos(2*uk)+G.Crs*sin(2*uk);
iik=G.Cic*cos(2*uk)+G.Cis*sin(2*uk);
uk=uk+iuk ;           %改正后的升角距
r=a*(1-G.e*cos(Ek))+irk;            %改正后的地心相径
i=G.i0+iik+G.iDot*tk;               %改正后的倾角 
xx=r*cos(uk);                        %在轨道面中的位置
yy=r*sin(uk);
omg=G.omg0+(G.omg-we)*tk-we*G.toe;   %改正后的升交点经度
x=xx*cos(omg)-yy*cos(i)*sin(omg);
y=xx*sin(omg)+yy*cos(i)*cos(omg);
z=yy*sin(i);
CalPos=[x,y,z];
x=num2str(x,'%12.6f');
y=num2str(y,'%12.6f');
z=num2str(z,'%12.6f');
fprintf('卫星号PRN:%2.0f',PRN); fprintf('\n');   
fprintf('时间time:%4.0f%4.0f%4.0f%4.0f%4.0f%4.1f',time); fprintf('\n');
fprintf('所求卫星在地心坐标系中的空间坐标为：\n');
fprintf('X = ');fprintf(x,'\n');fprintf('m\n');
fprintf('Y = ');fprintf(y,'\n');fprintf('m\n');
fprintf('Z = ');fprintf(z,'\n');fprintf('m\n');
return
toc   % over