analog_move_satellite.m

GPS 演示程序,是教科书上的,个人认为比较有参考价值

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);    %调用程序5
    hold on;
    A1=[32.8    92.8    152.8   212.6   272.8   332.8];     %依据图B-2
    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];    %平均近地角，参见卫星星座（图B-2）
        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;   %近地点的幅角，总共有6个卫星轨道平面，24卫星轨道椭圆，由于e很小，无法区分，因此画出6个椭圆表示
                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);
    %播放12次影像动画
    %每秒显示两祯