gui.m

本人编写的天际反卫星场景仿真代码。希望大家喜欢

clear
clc
rou=[400000 405000];
A=[pi*0.31 pi*0.32];
E=[pi*0.21 pi*0.22];
t=[20 24];
s=[pi/3 pi/3+pi/10800];                                %时角序列%
zemuda=pi/3;
fai=pi/3;
n=length(A);
tt=t-t(1)*ones(1,n);
%计算两点数据%
Re=6378000;
u=3.986e14;
c1=[-sin(zemuda) cos(zemuda) 0;-cos(zemuda)*sin(fai) -sin(fai)*sin(zemuda) cos(zemuda);cos(zemuda)*cos(fai) sin(zemuda)*cos(fai) sin(fai)];    %测站坐标到地球坐标系转换矩阵%
R=[Re*cos(fai)*cos(zemuda);Re*cos(fai)*sin(zemuda);Re*sin(fai)];                        %测量站到地心矢量在地球坐标系中分量%
rc1=[rou(1)*cos(E(1))*sin(A(1));rou(1)*cos(E(1))*cos(A(1));rou(1)*sin(E(1))];           %卫星在测站坐标系中的矢量%
rc2=[rou(2)*cos(E(2))*sin(A(2));rou(2)*cos(E(2))*cos(A(2));rou(2)*sin(E(2))];    
c21=[1 0 0;0 cos(s(1)) sin(s(1));0 -sin(s(1)) cos(s(1))];                               %地球坐标系与地心惯性坐标系转换矩阵%
c22=[1 0 0;0 cos(s(2)) sin(s(2));0 -sin(s(2)) cos(s(2))];   
r1=(c21)'*(R+(c1)'*(rc1));                                                              %卫星在地心惯性坐标系中表达%
r2=(c22)'*(R+(c1)'*(rc2)); 
%用两点数据计算出r0 v0%
tao=t(2)-t(1);
r00=r1;                                                       %初始位置%
F0=1-tao^2/(2*(norm(r00))^3);
G0=tao-tao^3/(6*(norm(r00))^3);
v00=(r2-F0*r1)/G0;                                          %初始速度%
%循环初始化%
r0=r00;
v0=v00;
r1=[0;0;0];
v1=[0;0;0];
a1=r0;
b1=v0;
e11=0.01;
e22=0.01;
while(norm(r1-r0)>e11|norm(v1-v0)>e22)
    r0=a1;v0=b1;
    %由r0 v0计算轨道六根数%
    a=u*norm(r0)/(2*u-norm(r0)*(norm(v0))^2);                       %半长轴%
    h=cross(r0,v0);
    p=(norm(h))^2/u;
    e=(1-p/a)^0.5;                                                %偏心率%
    z=[0;0;1];
    N=cross(z,h)/norm(cross(z,h));
    if (N(1)>=0)&(N(2)>=0) W=asin(N(2));
    elseif (N(1)>=0)&(N(2)<0) W=asin(N(2))+2*pi;
    elseif (N(1)<0)&(N(2)>=0) W=acos(N(1));
    else W=acos(-N(1))+pi; 
    end                                                              %升交点赤经%
    i=acos(dot(z,h)/norm(h));                                      %轨道倾角%
    if dot(r0,v0)>=0 theta=acos((1-p/norm(r0))/e);
    else theta=acos((1-p/norm(r0))/e)+pi;
    end                                                              %真近角%
    if (r0(3)>0) U=acos((norm(dot(r0,N)))/(norm(r0)));
    else U=acos((norm(dot(r0,N)))/(norm(r0)))+pi;
    end                                                           %纬度幅角%
    w=U-theta;                                                    %近心点角距%
    psy=2*atan(((1-e)/(1+e))^0.5*(tan(theta/2)));                  %偏近角%
    m=psy-e*sin(psy);                                              %  m  %
    psy
    %初始值%
    A1=0;B=0;C=0;D=0;L=[0;0;0];M=[0;0;0];
    %内循环%
    for k=1:n
        rc=[rou(k)*cos(E(k))*sin(A(k));rou(k)*cos(E(k))*cos(A(k));rou(k)*sin(E(k))];       %卫星在测站坐标系中的矢量%
        c1=[-sin(zemuda) cos(zemuda) 0;-cos(zemuda)*sin(fai) -sin(fai)*sin(zemuda) cos(zemuda);cos(zemuda)*cos(fai) sin(zemuda)*cos(fai) sin(fai)];    %测站坐标到地球坐标系转换矩阵%
        R=[Re*cos(fai)*cos(zemuda);Re*cos(fai)*sin(zemuda);Re*sin(fai)];                       %测量站到地心矢量在地球坐标系中分量%
        c2=[1 0 0;0 cos(s(k)) sin(s(k));0 -sin(s(k)) cos(s(k))];                                 %地球坐标系与地心惯性坐标系转换矩阵%
        r=(c2)'*(R+(c1)'*rc);  
        %用牛顿迭代计算t(k)时刻的delapsy(k)%
        a
        tao1=tt(k);
        syms dpsy f;
        f=a*tao1+(1-norm(r0)/a)*sin(dpsy)-norm(r0)*norm(v0)*(1-cos(dpsy))/(a)^0.5-dpsy;
        df=diff(f);
        df
        dpsy0=0.1; 
        e1=0.01;
        e2=0.01;
        dpsy=dpsy0;
        abs(eval(f))
        while(abs(eval(f))>e1)
            dpsy1=dpsy0-eval(f)/eval(df);
            if abs(dpsy1-dpsy0)<e2
                   break
            else
               dpsy0=dpsy1;
               dpsy=dpsy1;
               dpsy
            end
        end
        %计算t(k)时刻的F(k),G(k)%
        F=1-a/norm(r0)*(1-cos(dpsy));
        G=tao1-(dpsy-sin(dpsy))/(a^(-1.5));
        A1=A1+F^2;
        B=F*G;
        C=C+G^2;
        D=D+A1*C-B^2;
        L=L+F*r;
        M=M+G*r;
    end
r1=(C*L-B*M)/D;
v1=(A1*M-B*L)/D;
a1=r1;
b1=v1;
end
out=[a e i W w m];
text = sprintf('轨道六根数为：%f %f %f %f %f %f\n',[a e i W w m]);
disp(text);