ALOS卫星轨道插值并绘制轨道

已知ALOS卫星为用户提供的轨道数据，总共有28个轨道位置点，包括轨道点的时间、空间位置(x,y,z)及其在该点的速度(dx,dy,dz)，任意两个轨道点的获取时间相隔为60秒。
利用各种数据插值方法计算数据获取时间内任意时刻的卫星所在的位置。
并用三维绘图命令画出地球及相应轨道。

%% 用字符格式方式读取数据
FID=fopen(‘I:\专业课程\MATLAB\实验二-三-卫星参数文件读取及基线计算\LED-ALPSRP063392880-P1.1__D.orbit.txt’);
FID_num=textscan(FID,’%s%d’,1);
FID_head=textscan(FID,’%s%d%s’,2);
FID_data=textscan(FID,’%s%f%f%f%s%s%s’,56);
temp_1=FID_data{1,2};
temp_2=FID_data{1,3};
temp_3=FID_data{1,4};
orbit_data=[temp_1 temp_2 temp_3];
% 初始时间和时间间隔
Time_first=double(FID_head{1, 2}(1,1));
Time_intervel=double(FID_head{1, 2}(2,1));

%% 轨道位置数据插值,以时间为基准
myfrontsize=14;
xyz_data=orbit_data(1:2:56,:);
plot3(xyz_data(:,1),xyz_data(:,2),xyz_data(:,3),‘r’,‘linewidth’,2);
title(‘卫星轨道’);
xlabel(‘X( m)’,‘Fontsize’,18);
ylabel(‘Y( m)’,‘Fontsize’,18);
zlabel(‘Z( m)’,‘Fontsize’,18);
%% 插值
orignal_time=Time_first:Time_intervel:57660;
interp_time=Time_first:1:57660;
xyz_data_inter=zeros(size(interp_time,2),3);
% spline 球面线性插值 pchip三次多项式插值 nearest最近邻插值 lagrange插值
xyz_data_inter(:,1)=interp1(orignal_time,xyz_data(:,1),interp_time,‘nearest’);
xyz_data_inter(:,2)=interp1(orignal_time,xyz_data(:,2),interp_time,‘nearest’);
xyz_data_inter(:,3)=interp1(orignal_time,xyz_data(:,3),interp_time,‘nearest’);
xyz_data_inter(:,1)=lagrange(orignal_time,xyz_data(:,1),interp_time);
xyz_data_inter(:,2)=lagrange(orignal_time,xyz_data(:,2),interp_time);
xyz_data_inter(:,3)=lagrange(orignal_time,xyz_data(:,3),interp_time);
%% 绘制卫星轨道图
plot3(xyz_data_inter(:,1),xyz_data_inter(:,2),xyz_data_inter(:,3),‘r’,‘linewidth’,2);hold on;
plot3(xyz_data(:,1),xyz_data(:,2),xyz_data(:,3),‘g’);
title(‘拉格朗日邻插值卫星轨道’,‘Fontsize’,20);
xlabel(‘X( m)’,‘Fontsize’,18);
ylabel(‘Y( m)’,‘Fontsize’,18);
zlabel(‘Z( m)’,‘Fontsize’,18);
% Plotting the Earth 画地球
equat_rad=6378137.00;
polar_rad=6356752.3142;
[xx yy zz]=ellipsoid (0,0,0,equat_rad, equat_rad, polar_rad);
load(‘topo.mat’,‘topo’,‘topomap1’);
topo2 = [topo(:,181:360) topo(:,1:180)];
pro.FaceColor= ‘texture’;
pro.EdgeColor = ‘none’;
pro.FaceLighting = ‘phong’;
pro.Cdata = topo2;
colormap(topomap1)
omega_earth = 7.292115855377074e-005; % (rad/sec)
J2000_days=103752/24;
Go = 1.727564365843028; % (rad) http://www.amsat.org/amsat/articles/g3ruh/106.html
GMST = Go + omega_earth86400*(J2000_days + 0.5);
GMST = GMST - 2pifloor(GMST/(2pi));
GMST_deg=GMST(180/pi);
earth= surface(xx,yy,zz,pro);
%% 写文本文件
Pos_data=[interp_time;xyz_data_inter(:,1)’;xyz_data_inter(:,2)’;xyz_data_inter(:,3)’];
fid=fopen(‘I:\专业课程\MATLAB\实验四\插值轨道数据.txt’,‘w’);
fprintf(fid,’%s\t\t%s\t\t%s\t\t%s\r\n’,‘t’,‘x’,‘y’,‘z’);
fprintf(fid,’%f\t%6.2f\t%6.2f\t%6.2f\r\n’,Pos_data);
fclose(fid);
%% 利用重复轨道计算轨道之间的基线
FID_2=fopen(‘I:\专业课程\MATLAB\实验二-三-卫星参数文件读取及基线计算\LED-ALPSRP070102880-P1.1__D.orbit.txt’);
FID_num_2=textscan(FID_2,’%s%d’,1);
FID_head_2=textscan(FID_2,’%s%d%s’,2);
FID_data_2=textscan(FID_2,’%s%f%f%f%s%s%s’,56);
temp_1=FID_data_2{1,2};
temp_2=FID_data_2{1,3};
temp_3=FID_data_2{1,4};
orbit_data_2=[temp_1 temp_2 temp_3];
% 初始时间和时间间隔
Time_first_2=double(FID_head_2{1, 2}(1,1));
Time_intervel_2=double(FID_head_2{1, 2}(2,1));
% 样条插值
xyz_data_2=orbit_data_2(1:2:56,:);
orignal_time_2=Time_first_2:Time_intervel_2:57660;
interp_time_2=Time_first_2:1:57660;
xyz_data_inter_2=zeros(size(interp_time_2,2),3);
xyz_data_inter_2(:,1)=interp1(orignal_time_2,xyz_data_2(:,1),interp_time_2,‘spline’);
xyz_data_inter_2(:,2)=interp1(orignal_time_2,xyz_data_2(:,2),interp_time_2,‘spline’);
xyz_data_inter_2(:,3)=interp1(orignal_time_2,xyz_data_2(:,3),interp_time_2,‘spline’);
% 原始轨道数据
plot3(xyz_data(:,1),xyz_data(:,2),xyz_data(:,3),’*r’,‘linewidth’,2);hold on;
plot3(xyz_data_2(:,1),xyz_data_2(:,2),xyz_data_2(:,3),‘db’,‘linewidth’,2);hold on;
title(‘拉格朗日邻插值卫星轨道’,‘Fontsize’,20);
xlabel(‘X( m)’,‘Fontsize’,18);
ylabel(‘Y( m)’,‘Fontsize’,18);
zlabel(‘Z( m)’,‘Fontsize’,18);
% 插值后轨道数据
plot3(xyz_data_inter(:,1),xyz_data_inter(:,2),xyz_data_inter(:,3),’*r’,‘linewidth’,2);hold on;
plot3(xyz_data_inter_2(:,1),xyz_data_inter_2(:,2),xyz_data_inter_2(:,3),‘ob’,‘linewidth’,2);hold on;
title(‘拉格朗日邻插值卫星轨道’,‘Fontsize’,20);
xlabel(‘X( m)’,‘Fontsize’,18);
ylabel(‘Y( m)’,‘Fontsize’,18);
zlabel(‘Z( m)’,‘Fontsize’,18);