function [xs,ys,zs]=gps_sat(t);

%*********************************
% ** gps 星座(地心地固坐標) 
% **********************************

%得出時刻t的衛星坐標
a = 26559800*ones(6,4); % a 軌道面半長軸
% ne=0.01*ones(6,4); % e 偏心率
ne=0.0*ones(6,4); % e 偏心率
i=55*ones(6,4); % i 軌道面與赤道面傾角
bomigar=[296,296,296,296;  % omiga 升交點赤經
    356,356,356,356;
    56,56,56,56;
    116,116,116,116;
    176,176,176,176;
    236,236,236,236];
w=0*ones(6,4); % w 近地點角
v0=[225,315,45,135;  % v0 真近地點角初始值
    255,345,75,165;
    195,285,15,105;
    225,315,45,135;
    255,345,75,165;
    195,285,15,105];
% 衛星運行周期 11小時58分43080。近似12小時43200
v = 360/43080*t*ones(6,4)+v0;
we = 0;%7.2921151467e-5*360/(2*pi);   %0;  % we 地球自轉角速度7.2921151467e-5弧度
bomigarG0=0;
bomigarG=bomigarG0+we*t; % omiga_G 地心地固坐標系旋轉
trans=2*pi/360;   
r=a.*(ones(6,4)-ne.^2)./(ones(6,4)+ne.*cos(trans*v));
xs=r.*cos(trans*(v+w)).*cos(trans*(bomigar-bomigarG))-r.*sin(trans*(v+w)).*sin(trans*(bomigar-bomigarG)).*cos(trans*i);
ys=r.*cos(trans*(v+w)).*sin(trans*(bomigar-bomigarG))+r.*sin(trans*(v+w)).*cos(trans*(bomigar-bomigarG)).*cos(trans*i);
zs=r.*sin(trans*(v+w)).*sin(trans*i);

% a=26559800*ones(1,24);%a是軌道半徑
% i=55.0*ones(6,4); % i 軌道面與赤道面傾角
% ne=0.0*ones(1,24); % ne是軌道離心率，用圓軌道，所以為0
% bomigar=[272.85,272.85,272.85,272.85;  % omiga 升交點赤經
%     332.85,332.85,332.85,332.85;
%     32.85,32.85,32.85,32.85;
%     92.85,92.85,92.85,92.85;
%     152.85,152.85,152.85,152.85;
%     212.85,212.85,212.85,212.85];
% w=0*ones(6,4); % w 近地點角
% v0=[11.68,41.81,161.79,268.13;  % v0 真近地點角初始值
%     80.96,173.34,204.38,309.98;
%     111.88,241.57,339.67,11.80;
%     135.27,167.36,265.45,35.16;
%     197.05,302.60,333.69,66.07;
%     238.89,345.23,105.21,135.35];
% v=360/43080*t*ones(6,4)+v0;
% we=(7.2921151467e-5)*360/(2*pi);                      %地球自轉角速度    
% bomigarG0=0;
% bomigarG=bomigarG0+we*t; % omiga_G 地心地固坐標系旋轉
% trans=2*pi/360;   
% r=a.*(ones(6,4)-ne.^2)./(ones(6,4)+ne.*cos(trans*v));
% xs=r.*cos(trans*(v+w)).*cos(trans*(bomigar-bomigarG))-r.*sin(trans*(v+w)).*sin(trans*(bomigar-bomigarG)).*cos(trans*i);
% ys=r.*cos(trans*(v+w)).*sin(trans*(bomigar-bomigarG))+r.*sin(trans*(v+w)).*cos(trans*(bomigar-bomigarG)).*cos(trans*i);
% zs=r.*sin(trans*(v+w)).*sin(trans*i);