function SatPos = SatellitePos(Satellite,GASTw,t)
%  計算t時刻衛星Satellite的空間位置,不考慮攝動項
%  Satellite = [ i0 Omiga0 a e w M0 ]; 　　　　其中6個分量為初始時刻t0時的初始值
%　　6個分量分別是軌道面傾角，升交點赤經，長半軸，軌道偏心率，近地點角距，平近地角
%  GASTw : 星期六午夜至星期天子夜的交換時刻格林尼治視恒星時
%  t = tk-t0;
%  i_gps = 55 * 2*pi/360; i_Galileo = 56 * 2*pi/360;
%  e = 0.02; % assumed
%  a_gps = 26562km; a_Galileo = 29584km; 
%  Ts_gps = 11h58m; Ts_Galileo =14h04m;
DtoR = 2*pi/360;
i0 = Satellite(1) * DtoR;
Omiga0 = Satellite(2) * DtoR;
a = Satellite(3);
e = Satellite(4);
w = Satellite(5) * DtoR;
M0 = Satellite(6) * DtoR;

% 計算平均角速度n
Miu = 3.986005e014;
n0 = sqrt( Miu / a^3 );
Delta_n = 0;
n = n0 + Delta_n;
% 計算歸化時間tk
t0 = 0;    % 初始化時刻
tk = t - t0;
% 計算觀測時刻衛星平近點角Mk
Mk = M0 + n * tk;
% 計算偏近點角Ek
% Ek=Mk+esin(Ek),故用迭代算法求出，e通常較小，所以2次迭代即可求出
Ek = Mk;
for j=1:5
    Ek = Mk + e * sin(Ek);
end
% 計算衛星矢徑r
r = a * (1- e * cos(Ek));
% 計算真近點角fk
%%%%%% temp = ( sqrt(1-e^2) * sin(Ek) ) / ( cos(Ek) - e ) ;
%%%%%% fk = atan(temp);   % 該方法求得的fk只有-90度-90度
temp = sqrt((1+e)/1-e);
fk = 2 * atan( temp * tan(Ek/2) );
% 計算升交點角距Thita_k
Thita_k = fk + w;
% 計算攝動改正項Delta_u,Delta_r,Delta_i
Delta_u=0;
Delta_r=0;
Delta_i=0;
% 計算經過攝動改正的升交點角距uk，衛星矢徑rk，軌道傾角ik
uk = Thita_k + Delta_u;
rk = r + Delta_r;
    Dot_i = 0;
ik = i0 + Delta_i + Dot_i*tk;
% 計算衛星在軌道平面直角坐標系中的位置
% 其中軌道平面直角坐標系X軸指向升交點
xk = rk * cos(uk);
yk = rk * sin(uk);
zk = 0;
SatPos = [xk yk zk]';  % 衛星在軌道平面直角坐標系中的位置
% 計算觀測時刻升交點經度Omiga_k
Omiga = Omiga0;
we = 7.29211567e-005;        % 地球自轉角速率
GAST = GASTw + we * tk;
Omiga_k = Omiga - GASTw;
% 計算衛星在地心固定坐標系統中的位置
R3 = [cos(Omiga) -sin(Omiga) 0;
      sin(Omiga) cos(Omiga) 0;
      0          0          1;];
R1 = [1          0          0;
      0          cos(ik)    -sin(ik);
      0          sin(ik)    cos(ik);];
R2 = [cos(w) -sin(w) 0;
      sin(w) cos(w)  0;
      0      0       1;];
SatPos = R3 * R1 * [xk yk 0]';  %  SatPos =[Xk Yk Zk]'