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

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

%得出時刻t的衛星坐標
a=29993707*ones(3,9); % a 軌道面半長軸
% a=23616000*ones(3,9); % a 軌道面半長軸
ne=zeros(3,9); % e 偏心率
i=56*ones(3,9); % i 軌道面與赤道面傾角
bomigar=[0,0,0,0,0,0,0,0,0;  % omiga 升交點赤經
    120,120,120,120,120,120,120,120,120;
    240,240,240,240,240,240,240,240,240];
w=0*ones(3,9); % w 近地點角
v0=[0,40,80,120,160,200,240,280,320;  % v0 真近地點角初始值
    25,65,105,145,185,225,265,305,345;
    50,90,130,170,210,250,290,330,10];
% v=360/51720*t*ones(3,9)+v0; % 運行周期 14 h 21min//14.4h 51840
v=360/51840*t*ones(3,9)+v0;
we=0;    %7.2921151467e-5*360/(2*pi); % we 地球自轉角速度7.2921151467e-5弧度
bomigarG0=0;
bomigarG=bomigarG0+we*t; % omiga_G 地心地固坐標系旋轉
trans=2*pi/360;   
r=a.*(ones(3,9)-ne.^2)./(ones(3,9)+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);