showgvl_rtrack.m

建立GEO和LEO的轨道

function var=ShowGVL_RTrack(GEO_Radius,LEO_Radius,LEO_Angle)

u=3.986013*10^5;

LEO_T=2*pi*(LEO_Radius)^(3/2)/u^(1/2)/3600;
GEO_T=2*pi*GEO_Radius^(3/2)/u^(1/2)/3600;
m=1;%求死区时用的
t=linspace(0,GEO_T,5000);
for(i=1:length(t))
    %---------------------------------------------------------------
    %计算LEO和GEO的坐标
    LEO_x(i)=LEO_Radius*cos(2*pi/LEO_T*t(i))*cos(LEO_Angle/180*pi);
    LEO_y(i)=LEO_Radius*sin(2*pi/LEO_T*t(i));
    LEO_z(i)=LEO_Radius*cos(2*pi/LEO_T*t(i))*cos(pi/2-LEO_Angle/180*pi);
    GEO_x(i)=GEO_Radius*cos(2*pi/GEO_T*t(i));
    GEO_y(i)=GEO_Radius*sin(2*pi/GEO_T*t(i));
    GEO_z(i)=0;
    %----------------------------------------------------------------------
    %求距离
    D(i)=((GEO_z(i)-LEO_z(i))^2+(GEO_y(i)-LEO_y(i))^2+(GEO_x(i)-LEO_x(i))^2)^(1/2);
    %----------------------------------------------------------------------
    %求俯仰角
    omiga(i)=asin(LEO_z(i)/D(i));
    %-------------------------------------------------------------------
    %求方位角，GEO与LEO连线与X轴的夹角
    if((LEO_x(i)-GEO_x(i))>0&&(LEO_y(i)-GEO_y(i))>=0)
        cita(i)=atan((LEO_y(i)-GEO_y(i))/(LEO_x(i)-GEO_x(i)))/pi*180;    
    else if((LEO_x(i)-GEO_x(i))<=0&&(LEO_y(i)-GEO_y(i))>0)
            cita(i)=(pi+atan((LEO_y(i)-GEO_y(i))/(LEO_x(i)-GEO_x(i))))/pi*180;
        else if((LEO_x(i)-GEO_x(i))<=0&&(LEO_y(i)-GEO_y(i))<=0)
                cita(i)=(pi+atan((LEO_y(i)-GEO_y(i))/(LEO_x(i)-GEO_x(i))))/pi*180;
            else if((LEO_x(i)-GEO_x(i))>0&&(LEO_y(i)-GEO_y(i))<0)
                    cita(i)=(2*pi+atan((LEO_y(i)-GEO_y(i))/(LEO_x(i)-GEO_x(i))))/pi*180;
                end
            end
        end
    end
    %----------------------------------------------------------------------
    %求方位角，GEO与LEO连线和GEO与地心连线的夹角
    if(GEO_y(i)<=0&&GEO_x(i)<=0)
        cita1(i)=(atan(GEO_y(i)/GEO_x(i)))/pi*180;
    else if(GEO_y(i)<=0&&GEO_x(i)>0)
            cita1(i)=(pi+atan(GEO_y(i)/GEO_x(i)))/pi*180;
        else if(GEO_y(i)>0&&GEO_x(i)>=0)
                cita1(i)=(pi+atan(GEO_y(i)/GEO_x(i)))/pi*180;
            else if(GEO_y(i)>0&&GEO_x(i)<=0)
                    cita1(i)=(2*pi+atan(GEO_y(i)/GEO_x(i)))/pi*180;
                end
            end
        end
    end
    cita_xy(i)=cita(i)-cita1(i);
%======================================================================
    %求盲区
    x=linspace(min(GEO_x(i),LEO_x(i)),max(GEO_x(i),LEO_x(i)),30);
   
    for(j=1:length(x))
        y(j)=(GEO_y(i)-LEO_y(i))/(GEO_x(i)-LEO_x(i))*(x(j)-GEO_x(i))+GEO_y(i);
        z(j)=(GEO_z(i)-LEO_z(i))/(GEO_x(i)-LEO_x(i))*(x(j)-GEO_x(i))+GEO_z(i);
        if((x(j)^2+y(j)^2+z(j)^2)<6378^2)
            if (i<length(t)-1)
                k(m)=i;
                m=m+1;
                break;
            end
        end
    end
    %======================================================================
end

%-----------------------------------------------------------------------
%求角速度
for(i=1:length(t)-1)
    omiga_speed(i)=(omiga(i+1)-omiga(i))/((t(i+1)-t(i))*3600);%俯仰速度
    %-----------------------------------------------------------------
    %方位角速度
    if(cita(i)<=360&&cita(i)>=270&&cita(i+1)<=90&&cita(i+1)>=0);
        cita_speed(i)=((cita(i+1)-cita(i)+360))/((t(i+1)-t(i))*3600);
    else if(cita(i)<=90&&cita(i)>=0&&cita(i+1)<=360&&cita(i+1)>=270);
            cita_speed(i)=((cita(i+1)-cita(i)-360))/((t(i+1)-t(i))*3600);
        else
            cita_speed(i)=((cita(i+1)-cita(i)))/((t(i+1)-t(i))*3600);
        end
    end
end

%--------------------------------------------------------------------------
%求加速度
for(i=1:length(t)-2)
    omiga_ac(i)=(omiga_speed(i+1)-omiga_speed(i))/((t(i+1)-t(i))*3600);
    cita_ac(i)=(cita_speed(i+1)-cita_speed(i))/((t(i+1)-t(i))*3600);
end
%--------------------------------------------------------------
%显示地球
[Earth_x,Earth_y,Earth_z]=sphere(30);
figure;
mesh(6378*Earth_x,6378*Earth_y,6378*Earth_z);hold on 
%--------------------------------------------------------------------------
%显示卫星轨迹
%--------------------------------------------------------------------------

       for(i=1:100:length(t)-1)
           plot3([GEO_x(i),GEO_x(i+1)],[GEO_y(i),GEO_y(i+1)],[GEO_z(i),GEO_z(i+1)],'.');
           axis equal;hold on;
           plot3([0,max(max(GEO_x))],[0,0],[0,0]);
           text(max(max(GEO_x)),0,0,'X轴');
           plot3([0,0],[0,max(max(GEO_y))],[0,0]);
           text(0,max(max(GEO_y)),0,'Y轴');
           plot3([0,0],[0,0],[0,5*max(max(max(LEO_z)),max(max(Earth_z)))]);
           text(0,0,5*max(max(max(LEO_z)),max(max(Earth_z))),'Z轴');
           pause(0);
           for(j=0:round(GEO_T/LEO_T))
               plot3([LEO_x(i+j-fix((i+j)/length(t))*(length(t)-1)),LEO_x(i+j+1-fix((i+j)/length(t))*(length(t)-1))],...
                       [LEO_y(i+j-fix((i+j)/length(t))*(length(t)-1)),LEO_y(i+j+1-fix((i+j)/length(t))*(length(t)-1))],...
                       [LEO_z(i+j-fix((i+j)/length(t))*(length(t)-1)),LEO_z(i+j+1-fix((i+j)/length(t))*(length(t)-1))]);
           end
       end
 %=========================================================================
 %显示方位和俯仰角速度，死区
figure;
plot(t(2:end-1),cita_speed(2:end));hold on;
xlabel('归一化周期 单位：小时');
ylabel('方位和俯仰角速度 单位：度/秒');
%------------------------------------------------------------------------
for(i=1:length(k))
    cita_xy(k(i))=0;
    omiga(k(i))=0;
end
title(['方位角','[',num2str(min(cita_xy)),',',num2str(max(cita_xy)),']',...
        '俯仰角','[',num2str(min(omiga)),',',num2str(max(omiga)),']']);

%---------------------------------------------

plot(t(2:end-1),omiga_speed(2:end),'r');
plot(t(k),cita_speed(k),'r.');
title('方位、俯仰角速度和死区');
legend('方位角速度','俯仰角速度','死区');
hold off;
%=========================================================================
%显示方位和俯仰角加速度和死区
figure;
plot(t(2:end-1),cita_ac);hold on;
xlabel('归一化周期 单位：小时');
ylabel('方位和俯仰角加速度 单位：度/平方秒')
plot(t(2:end-1),omiga_ac,'r');
plot(t(k),cita_ac(k),'r.');
title('方位、俯仰角加速度和死区');
legend('方位角加速度','俯仰角加速度','死区');
%=========================================================================
%显示距离和死区
figure;
plot(t,D);hold on;
xlabel('归一化周期 单位：小时');
ylabel('卫星距离 单位：千米')
plot(t(k),D(k),'r.');
title('链路距离和死区');
legend('卫星距离','死区');