showlvg_etrack.m

建立GEO和LEO的轨道

function var=ShowLVG_ETrack(GEO_Radius,LEO_LongAxes,LEO_e,LEO_Angle)
       a=LEO_LongAxes;
       e=LEO_e;
       b=(a^2-a*e)^(1/2);
       u=3.986013*10^5;

LEO_T=2*pi*(a)^(3/2)/u^(1/2)/3600;
GEO_T=2*pi*GEO_Radius^(3/2)/u^(1/2)/3600;
%--------------------------------------------------------------
%计算LEO和GEO的坐标，以GEO轨道平面为参照（为了画图）
        t=linspace(0,GEO_T,160);
        for(i=1:length(t))
            LEO_x(i)=(a*cos(2*pi/LEO_T*t(i))-a*e)*cos(LEO_Angle/180*pi);
            LEO_y(i)=a*(1-e^2)^(1/2)*sin(2*pi/LEO_T*t(i));
            LEO_z(i)=(a*cos(2*pi/LEO_T*t(i))-a*e)*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;
       end
%-------------------------------------------------------------
%显示地球
[Earth_x,Earth_y,Earth_z]=sphere(30);
figure;
mesh(6378*Earth_x,6378*Earth_y,6378*Earth_z);hold on 
%--------------------------------------------------------------------------
%显示LEO和GEO轨迹
for(i=1:4: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
       %轨迹显示完毕 
%--------------------------------------------------------------------------
%求方位和俯仰角速度，链路距离和盲区
m=1;
t=linspace(0,GEO_T,5000);
for(i=1:length(t))
    %----------------------------------------------------------------------
    %以LEO轨道平面为参考，重新计算坐标（为了求角速度和角加速度）
    LEO_x(i)=(a*cos(2*pi/LEO_T*t(i))-a*e);
    LEO_y(i)=a*(1-e^2)^(1/2)*sin(2*pi/LEO_T*t(i));
    LEO_z(i)=0;
    GEO_x(i)=GEO_Radius*cos(2*pi/GEO_T*t(i))*cos(LEO_Angle/180*pi);
    GEO_y(i)=GEO_Radius*sin(2*pi/GEO_T*t(i));
    GEO_z(i)=GEO_Radius*cos(2*pi/GEO_T*t(i))*cos(pi/2-LEO_Angle/180*pi);
    %---------------------------------------------------------------------
    %求距离
    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(GEO_z(i)/D(i));%求俯仰角
    %----------------------------------------------------------------------
    %求方位角，GEO与LEO连线与X轴的夹角
    if((GEO_y(i)-LEO_y(i))>0&&(GEO_x(i)-LEO_x(i))>=0)
        cita(i)=atan((GEO_y(i)-LEO_y(i))/(GEO_x(i)-LEO_x(i)))/pi*180;
    else if((GEO_y(i)-LEO_y(i))>0&&(GEO_x(i)-LEO_x(i))<0)
            cita(i)=(pi+atan((GEO_y(i)-LEO_y(i))/(GEO_x(i)-LEO_x(i))))/pi*180;
        else if((GEO_y(i)-LEO_y(i))<=0&&(GEO_x(i)-LEO_x(i))<=0)
                cita(i)=(pi+atan((GEO_y(i)-LEO_y(i))/(GEO_x(i)-LEO_x(i))))/pi*180;
            else if((GEO_y(i)-LEO_y(i))<=0&&(GEO_x(i)-LEO_x(i))>0)
                    cita(i)=(2*pi+atan((GEO_y(i)-LEO_y(i))/(GEO_x(i)-LEO_x(i))))/pi*180;
                end
            end
        end
    end
    %----------------------------------------------------------------------
    %求方位角，GEO与LEO连线和GEO与地心连线的夹角
    if(LEO_y(i)<=0&&LEO_x(i)<=0)
        cita1(i)=(atan(LEO_y(i)/LEO_x(i)))/pi*180;
    else if(LEO_y(i)<=0&&LEO_x(i)>0)
            cita1(i)=(pi+atan(LEO_y(i)/LEO_x(i)))/pi*180;
        else if(LEO_y(i)>0&&LEO_x(i)>=0)
                cita1(i)=(pi+atan(LEO_y(i)/LEO_x(i)))/pi*180;
            else if(LEO_y(i)>0&&LEO_x(i)<0)
                    cita1(i)=(2*pi+atan(LEO_y(i)/LEO_x(i)))/pi*180;
                end
            end
        end
    end 
    if(cita(i)-cita1(i)>180)
        cita_xy(i)=cita(i)-cita1(i)-180;
    else if(cita(i)-cita1(i)<-180)
            cita_xy(i)=cita(i)-cita1(i)+180;
        else
            cita_xy(i)=cita(i)-cita1(i);
        end
    end
    %===================================================================
    %求多普勒频移
    if(LEO_y(i)<=0&&LEO_x(i)<=0)
        cita2(i)=(2*pi-(pi/2-atan(LEO_y(i)/LEO_x(i))))/pi*180;
    else if(LEO_y(i)<=0&&LEO_x(i)>0)
            cita2(i)=(pi/2+atan(LEO_y(i)/LEO_x(i)))/pi*180;
        else if(LEO_y(i)>0&&LEO_x(i)>=0)
                cita2(i)=(pi/2+atan(LEO_y(i)/LEO_x(i)))/pi*180;
            else if(LEO_y(i)>0&&LEO_x(i)<0)
                    cita2(i)=(pi+pi/2+atan(LEO_y(i)/LEO_x(i)))/pi*180;
                end
            end
        end
    end 
    
        cita_duopule(i)=cita(i)-cita2(i);
        r(i)=(LEO_x(i)^2+LEO_y(i)^2)^(1/2);
        LEO_xSpeed(i)=(u*(2/r(i)-1/a))^(1/2);
        
        f(i)=LEO_xSpeed(i)*1000*1/0.8*10^6*cos(cita_duopule(i)/180*pi);
%======================================================================
    %求盲区
    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
%=========================================================================
%显示方位和俯仰角速度，死区

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;
    f(k(i))=0;
end
title(['方位角','[',num2str(min(cita_xy(2:end))),',',num2str(max(cita_xy(2:end))),']',...
        '俯仰角','[',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  on;
%==========================================================================
%显示方位和俯仰角加速度和死区
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('卫星距离','死区');
%=================================================================
figure;
plot(t,f);
xlabel('归一化周期 单位：小时');
ylabel('多普勒频移');
title('LEO光端机的多普勒频移');