切比雪夫轨道拟合MATLAB代码

function [Head_data,SP3_data]=Read_SP3 %读取精密星历文件中的数据 % global SP3_data % global Head_data [filename,filepath]=uigetfile('*.sp3','读取GPS精密星历文件'); fid=fopen(strcat(filepath,filename),'rt'); Head_data=struct; SP3_data =struct; if(fid==-1) msgbox('输入的文件或者路径不正确,无法正确打开sp3文件','警告信息'); return; end % 读取文件头，并将其存入Head_data结构数组中 Head=[]; for i=1:22 %制度文件头的前16行 line=fgetl(fid); Head=[Head;line];%将sp3文件的文件头以字符的形式显示出来，切记此处要保证sp3文件中的头文件的每行字符数是相同的，否则读取将失败 end %提取文件头第一行的数据 Head_data.Version=Head(1,1:2); Head_data.P_or_V=Head(1,3); Head_data.start_year=str2num(Head(1,4:7)); Head_data.start_month=str2num(Head(1,9:10)); Head_data.start_day=str2num(Head(1,12:13)); Head_data.start_hour=str2num(Head(1,15:16)); Head_data.start_minute=str2num(Head(1,18:19)); Head_data.start_second=str2num(Head(1,21:31)); Head_data.EpochNum=str2num(Head(1,33:39)); Head_data.data_type=Head(1,41:45); Head_data.Coordinate_system=Head(1,47:51); Head_data.Orbit_type=Head(1,53:55); Head_data.Institution=Head(1,57:60); %第二行 Head_data.GPS_week=str2num(Head(2,4:7)); Head_data.FirstEpoch_weeksecond=str2num(Head(2,9:23)); Head_data.Epoch_interval=str2num(Head(2,25:38)); Head_data.MJD_zhengshu=str2num(Head(2,40:44)); Head_data.MJD_xiaoshu=str2num(Head(2,46:60)); %第三、四行-------提取观测卫星的PRN编号 Head_data.SatNum=str2num(Head(3,5:6)); PRN=strcat(Head(3,10:60),Head(4,10:60)); L=length(PRN)/3; SP3_PRN=[]; for i=1:L Sp3_PRN=str2num(PRN((3*i-1):3*i)); if Sp3_PRN==0 continue; end SP3_PRN=[SP3_PRN Sp3_PRN]; end Head_data.SP3_PRN=SP3_PRN; %第八、九行-----提取相应卫星的精度 SatAccu=strcat(Head(8,10:60),Head(9,10:60)); L=length(SatAccu)/3; SP3_SatAccu=[]; for i=1:L Sp3_SatAccu=str2num(SatAccu((3*i-1):3*i)); if Sp3_SatAccu==0 continue; end SP3_SatAccu=[SP3_SatAccu Sp3_SatAccu]; end Head_data.SP3_SatAccu=SP3_SatAccu; %第十五行-----读取标准差的浮点基数 Head_data.pos_jishu=str2num(Head(15,4:13)); Head_data.Satclock_jishu=str2num(Head(15,15:26)); %--------------------End of the reading of the Fileheader------------------ %--------------------读取卫星坐标及卫星钟钟差-------------------------------- %读取的卫星的三维坐标的单位为km，卫星钟差的单位为us t=1; while 1 line=fgetl(fid); result=findstr(line,'EOF'); if(isempty(result)) SP3_data(t).year=str2num(line(4:7)); SP3_data(t).month=str2num(line(9:10)); SP3_data(t).day=str2num(line(12:13)); SP3_data(t).hour=str2num(line(15:16)); SP3_data(t).minute=str2num(line(18:19)); SP3_data(t).second=str2num(line(21:22)); SP3_data(t).mm=SP3_data(t).hour*60+SP3_data(t).minute+ SP3_data(t).second/60;%将该历元的时分秒化为分钟，方便后面和观测历元比较 SP3_data(t).TimeOEp=TimetoJD(SP3_data(t).year,SP3_data(t).month,SP3_data(t).day,SP3_data(t).hour,SP3_data(t).minute,SP3_data(t).second);%转换为儒略日 for j=1:Head_data.SatNum line=fgetl(fid); SP3_data(t).PRN(j)=str2num(line(3:4));%将每个历元的卫星编号保存在一个数组中 SP3_data(t).x(j)=str2num(line(5:18)); %提取卫星的三维坐标，单位为km SP3_data(t).y(j)=str2num(line(19:32)); SP3_data(t).z(j)=str2num(line(33:46)); SP3_data(t).Clock_corr(j)=str2num(line(47:60));%提取卫星钟差，单位为us %判断是否有坏的或空缺的位置值 % if(SP3_data(t).x(j)==0||SP3_data(t).y(j)==0||SP3_data(t).z(j)==0) % fprintf('注意：星历文件中的第%g个历元的%g号卫星的位置值是坏的或空缺的！\n',t,SP3_data(t).PRN(j)); % end % %判断是否为坏的或空缺的钟差值 % if(SP3_data(t).Clock_corr(j)==999999.999999) % fprintf('注意：星历文件中的第%g个历元的%g号卫星的钟差值是坏的或空缺的！\n',t,SP3_data(t).PRN(j)); % end %判断X的标准差是否为未知，若是将其置为0 if(length(line)>=73 && isempty(str2num(line(62:63)))) SP3_data(t).X_mi(j)=0; elseif(length(line)>=73 && ~isempty(str2num(line(62:63)))) SP3_data(t).X_mi(j)=str2num(line(62:63)); end %判断Y的标准差是否为未知，若是将其置为0 if(length(line)>=73 && isempty(str2num(line(65:66)))) SP3_data(t).Y_mi(j)=0; elseif(length(line)>=73 && ~isempty(str2num(line(65:66)))) SP3_data(t).Y_mi(j)=str2num(line(65:66)); end %判断Z的标准差是否为未知，若是将其置为0 if(length(line)>=73 && isempty(str2num(line(68:69)))); SP3_data(t).Z_mi(j)=0; elseif(length(line)>=73 && ~isempty(str2num(line(68:69)))) SP3_data(t).Z_mi(j)=str2num(line(68:69)); end %判断卫星钟差的标准差是否为未知，若是将其置为0 if(length(line)>=73 && isempty(str2num(line(71:73)))) SP3_data(t).Clock_mi(j)=0; elseif(length(line)>=73 && ~isempty(str2num(line(71:73)))) SP3_data(t).Clock_mi(j)=str2num(line(71:73)); end end elseif(~isempty(result)) break;%判断读到文件的末尾，到达文件尾则退出 end t=t+1;%记录文件体的行数 end fclose(fid);