% function SATsimulation

clear

clc

% GNSS星座参数设计

% 北斗(BDS)，其中IGSO自定义

We = 7.292115e-5;                                                        % 地球自转速度(rad)

Geo = 3.986004418e14;                                                    % 地心引力常数

A = 1e3*[42164*ones(1,5),42164*ones(1,3),27878*ones(1,27)];              %轨道长半轴

e = [0.002*ones(1,5),0.002*ones(1,3),0.002*ones(1,27)];                  %轨道偏心率

I = pi/180*[0*ones(1,5),55*ones(1,3),55*ones(1,27)];                     %轨道倾角

W = pi/180*[0*ones(1,5),0*ones(1,3),0*ones(1,27)];                       %近地点幅角

M = pi/180*[0*ones(1,5),[0 120 240],[[0:40:320],40/3+[0:40:320],80/3+[0:40:320]]]; %平近点角

Omega = pi/180*[[58.75 80 110.5 140 160],[118 335 261],[0*ones(1,9) 120*ones(1,9) 240*ones(1,9)]]; %升交点赤经

N = length(A);

TL = 24*3600; %second

for i=1:N

   aa=A(i);

   ee=e(i);

   ii=I(i);

   mm=M(i);

   ww=W(i);

   om=Omega(i);

   % 轨道真近心角f(rad)计算公式

   nn = sqrt(Geo/aa^3);   % 卫星平运动速度(rad/s)

   T = 2*pi / nn;         % 卫星运行周期

   Ts = 5;                % 采样时间间隔（s）

   t = [0:Ts:TL]';        % 采样数

   Mt = mm + nn * t;      % 卫星平近点角（rad）

   Et = Mt + (ee-ee^3/8+ee^5/192-ee^7/9216)*sin(Mt)+ (ee^2/2-ee^4/6+ee^6/98)*sin(2*Mt)...

        + (ee^3*3/8-ee^5*27/128+ee^7*243/5120)*sin(3*Mt)+ (ee^4/3-ee^6*4/15)*sin(4*Mt)...

        + (ee^5*125/384-ee^7*3125/921)*sin(5*Mt)+ ee^6*27/80*sin(6*Mt) + ee^7*16807/46080*sin(7*Mt);  % 偏近心角(rad)

   f = 2*atan( sqrt((1+ee)/(1-ee))*tan(Et/2) );  % 轨道真近心角(rad)

   R = aa*(1-ee^2)./(1+ee*cos(f));  % 卫星距地心的距离

   u=f+ww;

   x=R.*cos(u);

   y=R.*sin(u);

   omNEW = om - We*t;

   Xs(:,i)=x.*cos(omNEW)-y.*(cos(ii)*sin(omNEW));

   Ys(:,i)=x.*sin(omNEW)+y.*(cos(ii)*cos(omNEW));

   Zs(:,i)=y*sin(ii);

end

L=atan2(Ys,Xs)*180/pi;  

L= (sign(L)<0)*360 + L;

S=sqrt(Xs.^2+Ys.^2);

B=atan(Zs./S)*180/pi; %近似纬度

figure,plot(L(:,1:5),B(:,1:5),'.')     %GEO卫星

figure,plot(L(:,6:8),B(:,6:8),'.')     %IGSO卫星（自定义）

figure,plot(L(:,9:17),B(:,9:17),'.')   %第一轨道面MEO卫星

figure,plot(L(:,18:26),B(:,18:26),'.') %第二轨道面MEO卫星

figure,plot(L(:,27:35),B(:,27:35),'.') %第三轨道面MEO卫星

(北斗GEO星下点轨迹)

(北斗IGSO星下点轨迹)

(北斗MEO星下点轨迹-第一轨道面)