这是做卫星通信的课程设计时涉及的卫星星下点的作图程序,matlab编写。

% 轨道真近心角f(rad)，由以下公式计算

Ts = 5;% 采样时间间隔（s）

t = [0:Ts:fix(N_T*T)];% 采样时间点

tp = 5400;

n = sqrt(GGC/a^3);% 卫星平运动速度(rad)

M = n * (t - tp);% 卫星平近点角（rad）

f = M + (2*e - e^3/4)*sin(M) + 1.25*e^2 * sin(2*M) +13/12 * e^3 * sin(3*M);% 轨道真近心角(rad)

%% % % % % % % % % % % % % % % % % 计算星下点轨迹 % % % % % % % % % % % % % % Rz_RAAN = [cos(RAAN) -sin(RAAN) 0 ; sin(RAAN) cos(RAAN) 0 ; 0 0 1];

Rz_w = [cos(w) -sin(w) 0 ; sin(w) cos(w) 0 ; 0 0 1];

Rx_i = [1 0 0 ; 0 cos(i) -sin(i) ; 0 sin(i) cos(i)];

R = a*(1-e^2)./(1+e*cos(f));% 卫星距地心的距离，考虑f离散值

r_sv = [R;R;R] .* [cos(f);sin(f);zeros(1,size(f,2))];% 卫星在轨道坐标系中的坐标

r_so = Rz_RAAN * Rx_i * Rz_w * r_sv;% 卫星在地心惯性坐标系中的坐标

x_so = r_so(1,:) ; y_so = r_so(2,:) ;z_so = r_so(3,:) ;% 卫星在地心惯性坐标系的分量

% 卫星纬度delta，大于零为北纬

delta = atan( z_so ./sqrt((x_so .^2+y_so .^2)) ) * 180/PI;

% 卫星经度alpha，大于180为西经

for m = 1:1:size(r_so,2)

if (r_so(1,m) < 0)

alpha(m) = 180 + atan(r_so(2,m)/r_so(1,m)) * 180/PI;

else

if (r_so(2,m)>0)

alpha(m) = atan(r_so(2,m)/r_so(1,m)) * 180/PI;

else

alpha(m) = 360 + atan(r_so(2,m)/r_so(1,m)) * 180/PI;

end

end

end

% 计算地心经度alpha1

alpha1 = rem((alpha - W_EARTH*t*180/PI + N_T*360),360);

plot(alpha1,delta,'.');grid; xlabel('卫星经度');ylabel('卫星纬度');