GPS rinex文件读取的MATLAB代码.rar

function [L1, L2, P1, P2, time, sv, t0] = read_robs_02(obsfil); %[L1, L2, P1, P2, time, sv, t0] = read_robs_02(obsfil) % % debugged on ASE server Matlab v. 6.0. A change was added because % different versions of Matlab process strings differently. % (Thanks to Greg Holt for providing the modification.) % %Read a RINEX file and return the observables %time = column vector of observation times %sv = row vector of satellite PRN numbers %L1(itime,iprn) = L1 observation for time itime, PRN iprn %And similarly for L2, P1, P2. %obsfil is the name of the file enclosed in single quotes, eg. 'data1.o' % %All observables, code and phase, are returned in units of METERS. % % Note -- This function returns NaN values where there is % no observation. This makes it easy to write linear % combinations, but care must be taken with subsequent % averaging of values, taking min or max, etc. % You can use the function nan2zero to convert NaNs to zeros. % Author: Jeff Freymueller, Stanford University % % Based on original code by Chuck Meertens and % Doug Miller, University of Utah % Date: 7 February 1995 disp(' ') % open file obs = fopen(obsfil,'r'); % % tell user if file was successfully opened or not % if obs ~= -1 disp('Observation file successfully opened!') else disp('Unable to open observation file!') return end % % First, read the RINEX file header to get the order of % the observables. % % read the first line of obs file header % frewind(obs) line = fgets(obs); % % determine rinex version % version = sscanf(line,'%d',1) % % For now, we'll hard code the order of observables % % We will order the observables L1, L2, P1, P2, C1 obs_order(3) = 3; obs_order(1) = 1; obs_order(2) = 2; obs_order(4) = 4; % obs_order(5) = 5; % % read next lines of obs file header % if version == 1 % for rinex version 1 exit if no char in line COUNT = 1; while COUNT ~= 0 line = fgets(obs); [s, COUNT, ERRMSG] = sscanf(line,'%s'); end else % for rinex version 2 read until end of header encountered st = 'NOTNOT'; while any(abs(st(1:6)) ~= abs('ENDOFH')), line = fgets(obs); [st, count, err] = sscanf(line,'%s'); end end % % Load our fundamental constants % gpsconst; % % Initialize % n_t = 0; sv_index = zeros(1,37); nsv = 0; % % sv epoch time tag line = fgets(obs); while (line ~= -1) n_t = n_t + 1; % coded for up to a 12 channel GPS receiver ft = sscanf(line,'%f %f %f %f %f %f %f %f G%fG%fG%fG%fG%fG%fG%fG%fG%fG%fG%fG%f'); % convert rinex date/time to int Y int M real D Y = 2000 + ft(1); M = ft(2); D = ft(3) + (ft(4)/24) + (ft(5)/1440) + (ft(6)/86400); time(n_t,1) = date2j(Y,M,D); nsvs_t = ft(8); % read observations for all svs this epoch for n = 1:nsvs_t this_sv = ft(8+n); if (sv_index(this_sv) == 0) nsv = nsv + 1; sv_index(this_sv) = nsv; sv(1,nsv) = this_sv; end % Read the line line = fgets(obs); obsval(1) = read_strval(line( 1:16)); obsval(2) = read_strval(line(17:32)); obsval(3) = read_strval(line(33:48)); if length(line)<63 obsval(4) = 0; else obsval(4) = read_strval(line(49:63)); end % obsval(5) = read_strval(line(65:78)); % NOTE: velocity of light m/sec: c = 299792458.0e+0 % L1 and L2 frequencies : f1 = 1575.42e6, f2 = 1227.60e6 % lam1 = c/f1, lam2 = c/f2 % L1_t = obsval(obs_order(1))*lam1; % L2_t = obsval(obs_order(2))*lam2; L1_t = obsval(obs_order(1)); L2_t = obsval(obs_order(2)); % P1_t = obsval(obs_order(3)); P1_t = 0; P2_t = obsval(obs_order(4)); C1_t = obsval(obs_order(3)); if (P1_t == 0) P1_t = C1_t; end % Plug these into the matrices L1, L2, P1, P2 L1(n_t, sv_index(this_sv)) = L1_t; L2(n_t, sv_index(this_sv)) = L2_t; P1(n_t, sv_index(this_sv)) = P1_t; P2(n_t, sv_index(this_sv)) = P2_t; % line = fgets(obs); end line = fgets(obs); end % % Transform time into seconds past start of file % t0 = time(1); %time = time - ones(n_t,1)*t0; %time = time*86400; % % Quick and dirty -- convert all zeros to NaN % [n, m] = size(L1); for i = 1:n for j = 1:m if L1(i,j) == 0 L1(i,j) = NaN; end if L2(i,j) == 0 L2(i,j) = NaN; end if P1(i,j) == 0 P1(i,j) = NaN; end if P2(i,j) == 0 P2(i,j) = NaN; end end end