function [APXYZ DeltaHEN signal_type obsdata] = read_obs(filename)
% %%function [signal_type obsdata] = read_obs(filename)
%程序描述：读取观测数据文件，rinex格式
%---------------------------------------------------------------
%-- Inputs:
%-- filename
%-- Outputs:
%-- obsdata %-- observation message
% ---modified by Zhang 2017-12-1
%---------------------------------------------------------------
%-----Open the file-------%
fid = fopen(filename);
disp(['读取卫星观测数据文件 ...... : ',filename]);
% [filename,pathname]=uigetfile('*.**o','读取卫星观测数据文件');
% fid=fopen(strcat(pathname,filename),'rt');
% if(fid==-1)
% error('error to open the '+ filename);
% end
%% -----读取头部-----%
hang = zeros(1,80);
head_lines = 0;
while 1
head_lines = head_lines + 1;
hang = fgetl(fid);
head(head_lines).str = str2mat(hang);
if(~isempty(findstr(hang,'APPROX POSITION XYZ'))) %obtain approximate marker position (WGS-84)，获得观测点大概位置（WGS-84）
APXYZ=[0,0,0];
APXYZ(1) = str2num(hang(1:14));
APXYZ(2) = str2num(hang(15:28));
APXYZ(3) = str2num(hang(29:42));
end
if(~isempty(findstr(hang,'ANTENNA: DELTA H/E/N'))) %obtain Antenna height、Eccentri
cities of antenna center relative to marker to the east and north (all units in meters)
%天线高度：天线底部相交于标记点的高度，天线中心相对标记点东向和北向距离（单位：米）
DeltaHEN=[0,0,0];
DeltaHEN(1) = str2num(hang(1:14));
DeltaHEN(2) = str2num(hang(15:28));
DeltaHEN(3) = str2num(hang(29:42));
end
if(~isempty(findstr(hang,'WAVELENGTH FACT L1/2'))) %获取L1和L2载波的缺省系数，1:?整周；2:半周；0(L2载波):单频接收机
L1_Wavalength = str2num(hang(1:6));
L2_Wavalength = str2num(hang(7:12));
end
if(~isempty(findstr(hang,'# / TYPES OF OBSERV'))) %文件中观测类型数量
%格式例如：21 L1 P1 C1 L2 P2 D1 D2 S1 S2# / TYPES OF OBSERV
% L5 C5 D5 S5 L7 C7 D7 S7 L8# / TYPES OF OBSERV
% C8 D8 S8 # / TYPES OF OBSERV
type = str2num(hang(1:6));
if(type <= 9)%观测类型数量少于9个
for t = 1:type
signal_type(t).types = hang(6 * t +5 : 6 * t + 6);
end
else%观测类型数量多于9个，分多行显示
for t = 1:9%由于第1行已经读取，先存储第1行
signal_type(t).types = hang(6 * t +5 : 6 * t + 6);
end
for m=1:floor(type/9)%再循环存储其他行
head_lines = head_lines + 1;
hang = fgetl(fid);
head(head_lines).str = str2mat(hang);
for t = 9*m+1:9*m+9
if(t<=type)
signal_type(t).types = hang(6 * (t - 9*m) +5 : 6 * (t - 9*m) + 6);
end
end
end
end
end
%文件中观测类型数量读取完毕，相应说明：
% TYPES OF OBSERV：在本数据文件中所存储不同观测值类型的数量和观测值类型列表，如果超过9种观测值类型, 则使用续行。
% 说明：
% L1,L2：L1和L2上的相位观测值
% C1：采用L1上C/A 码所测定的伪距；
% P1,P2：采用L1 、L2 上的P 码所测定的伪距；
% D1,D2：L1 和L2 上的多普勒频率；
% T1,T2：子午卫星的150(T1)和400 MHz(T2) 信号上的多普勒积分；
% S1,S2：接收机所给出的L1、L2相位观测值的原始信号强度或SNR值。
% 观测值的单位: 载波相位为周,伪距为m,多普勒为Hz,子午卫星为周,SNR等则与接收机有关。
if(~isempty(findstr(hang,'INTERVAL')))%观测时间间隔，单位：秒
interval = str2num(hang(1:6));
end
if(~isempty(findstr(hang,'END OF HEADER')))
break;
end
end
% if head(1).str(21) ~= 'O'
% error('This file is not a navgation message file.Please selete a correct file');
% elseif head(1).str(41) ~= 'G'
% error('This is not a GPS system');
% end
%% --------读取观测数据-------------
data_num = 0;
while 1
%data_num = data_num + 1;
%----- First line of observation data -----
hang = fgetl(fid);
if(hang == -1) % end of file
break;
end
a=blanks(80);
a(1:length(hang)) = hang;
hang = a;
disp(hang);
if(isstrprop(hang(2),'digit')==1)fopen和open区别
data_num = data_num + 1;
%-------- Epoch / GPS time system --------
obsdata(data_num).year = 2000 + str2double(hang(1:3));
obsdata(data_num).month = str2num(hang(4:6));
obsdata(data_num).day = str2num(hang(7:9));
版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。
发表评论