MATLAB绘制GPS星下点轨迹图

轨迹计算部分参考链接：

https://wenku.baidu.com/view/45bd098d4a7302768e9939cf.html

本文对上述matlab代码进行了整理与修改：

增加了底图并进行了相关的图形美化。

将轨道六参数设置为GPS相关参数。

从原有的1颗卫星轨迹绘制增加至24颗。

可进一步修改的部分：

通过修改底图绘制部分相关参数可以进行图片的基础修改。

通过修改轨道六根数的大小以及循环次数可以实现其他卫星导航系统的星下点轨迹绘制。

代码如下：

clc

close all

clear all

%% % % % % % % % % % % % % % % % % 底图绘制 % % % % % % % % % % % % % % % % % % % % % %

h = geoshow('landareas.shp', 'FaceColor', 'c');

grid on

hold on

xlabel('Longitude');

ylabel('Latitude');%坐标轴标题

set(gca,'Ylim',[-90,90],'ytick',[-90:30:90]);

set(gca,'yticklabel',{'90°S','60°','30°','0°','30°','60°','90°N'});

set(gca,'Xlim',[-180,180],'xtick',[-180:30:180]);

set(gca,'xticklabel',{'180°W','150°','120°','90°','60°','30°','0°','30°','60°','90°','120°','150°','180°E'});

%坐标轴范围及刻度分划，坐标轴文字替代

set(gca,'Box','on');%坐标轴是否为四面

set(gca,'FontSize',10,'Fontname', 'Times New Roman','Fontweight', 'bold');%字号、字体、是否加粗

set(gca,'GridAlpha',1,'GridLineStyle','--');%格网透明度(0-1)及线型

title('Track of GPS Satellite Point','FontSize',14,'Fontweight', 'bold');%图标题

%% % % % % % % % % % % % % % % % % 变量定义 % % % % % % % % % % % % % % % % % % % % % %

PI = 3.1415926;

W_EARTH = 7.29e-5;% 地球自转速度(rad)

GGC = 3.986e5;% 地心引力常数

N_T =2;% 实验观测周期数，可延长观察时间

for ii = 1 : 24

% 轨道六根数

a = 26578;% 轨道长半轴(Km)

e = 0.009047195664607;% 轨道偏心率

i = 55 * PI/180;% 轨道倾角(rad)

w = 60 * PI/180;% 轨道近地点幅角(rad)

RAAN = 15 * ii * PI/180;% 升交点赤经

T = 2 * PI * sqrt((a^3) / GGC);

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

Ts = 30;% 采样时间间隔(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)-180;

%plot(alpha1,delta,'.','Color','m');%将轨迹设置为同一颜色

plot(alpha1,delta,'.');

end

print -f1 -r1200 -dpng Track_of_GPS_Satellite_Point;%图片保存

标签：cos,180,so,MATLAB,星下点,90,PI,sin,GPS

来源： https://blog.csdn.net/sinat_39238867/article/details/97892803