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);

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');

toc   % over
return