这是做卫星通信的课程设计时涉及的卫星星下点的作图程序,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('卫星纬度');