mnwxyd.m

标准GPS星座的仿真程序。注：采用24颗卫星

%function 模拟卫星运动
%function mnwxyd
clear;clc;close all;
DtoR=2*pi/360;  %转化为弧度单位
JiaoStep=360/24*DtoR;
for time=0:11
    figure(time+1);
    a=26560;    %卫星椭圆轨道长半轴
    e=0.02;
    E=[0:0.1:2*pi];
    x=a*(cos(E)-e);
    y=a*sqrt((1-e^2))*sin(E);
    z=0*E;
 
    DrawEarth(time);
    hold on;

    A1=[32.8 92.8 152.8 212.6 272.8 332.8];
    for k=1:6;
        A=A1(k)*DtoR;   %升交点赤经
        B=55*DtoR; %轨道倾角
        C=pi/100;
        %六个轨道平面
        R3=[cos(A) -sin(A) 0;
            sin(A) cos(A) 0;
            0      0      1];
        R1=[1   0   0;
            0   cos(B)  -sin(B);
            0   sin(B)  cos(B)];
        R2=[cos(C)  -sin(C) 0;
            sin(C)  cos(C)  0;
            0   0   1];
        L1=length(E);
        R312=R3*R1*R2;
        Ans=R312*[x;y;z];  %将卫星轨道坐标数据转换到地球赤道坐标
        x1=Ans(1,:);
        y1=Ans(2,:);
        z1=Ans(3,:);
        plot3c(x1,y1,z1,k);
        hold on;
        axis equal;
        axis off;
        grid on;
    end
    
    Ctable=[10  50  160 260;
        80  180 220 320;
        10  130 250 340;
        50  150 170 300;
        100 210 310 340;
        140 150 240 350];       %平均近地点角
    Wx=ones(1,1);
    Wy=ones(1,1);
    Wz=ones(1,1);
    for k=1:6;
        A=A1(k)*DtoR;
        B=55*DtoR;
        for m=1:4
            C=Ctable(k,m)*DtoR+JiaoStep*time;
            
            x=a*(cos(C)-e);
            y=a*sqrt((1-e^2))*sin(C);
            z=0*C;
            %六个轨道平面
            R3=[cos(A) -sin(A) 0;
                sin(A) cos(A) 0;
                0      0      1];
            R1=[1   0   0;
                0   cos(B)  -sin(B);
                0   sin(B)  cos(B)];
            R2=[cos(C)  -sin(C) 0;
                sin(C)  cos(C)  0;
                0   0   1];
            L1=length(E);
            R312=R3*R1*R2;
            Ans=R312*[x;y;z];  %将卫星轨道坐标数据转换到地球赤道坐标
            Wx=[Wx Ans(1,:)];
            Wy=[Wy Ans(2,:)];
            Wz=[Wz Ans(3,:)];
            x1=Ans(1,:);
            y1=Ans(2,:);
            z1=Ans(3,:);
            DrawSatellite(x1,y1,z1,k);
            hold on;
        end
    end
    drawnow;        %刷新屏幕
    M(time+1)=getframe; %获得影像动画的帧并将结果送给M
end
figure(111);
axis equal;
axis off;
movie(M,12,2);