loadrinexobs.m

利用gps多天线载波相位技术

%%========================================
%%     Toolbox for attitude determination
%%     Zhen Dai
%%     [e-mail]
%%     ZESS, University of Siegen, Germany
%%     Last Modified  : 1.Sep.2008
%%========================================
%% Functions:
%%      Read  the RINEX observation file and save the measurement
%% Input parameters:
%%      filename --> RINEX observation file name
%%      id -> This indicates the file from which antenna. 1- master antena,
%%              2,3... slave antennas
%%      totalepoch -> total epochs to be processed. -999 implies processing
%%               all the data in this file 
%% Output:
%%      mSatID -> Satellite visibility of each epoch
%%      vTime-> Time of each epoch scaled to "second of week"
%%      vTimeStr -> Time expressed in STRING
%%      interval -> Sampling interval of the data
%%      mCode -> C/A code on  L1 signal
%%      mPhase -> Carrier phase data on L1 signal
%%  Remarks:
%%      1. Only read the data from the GPS satellites
%%      2. Only consider the L1 C/A and L1 phase data
%%      3. Outputs are saved in a data file
%%      4. It may provide some unexpected result when processing some
%%      specially formatted RINEX observation file.
 
function datafilename=LoadRinexOBS(filename,id,totalepoch) 

%% GUI Initialization
f1=figure;
close(f1);
figure;
title=sprintf('Processing RINEX observation file #%d',id);
set(gcf,'MenuBar','none','NumberTitle','off','Name',title,...
    'Position',[400 500 450 150],'Visible','off');
 movegui(gcf,'center');
hlistbox = uicontrol(gcf,'Style', 'listbox', 'String', 'Clear',...
    'Position', [0 0 450 150],'String','','Foregroundcolor','b');
set(gcf,'Visible','on');

%% Start reading the RINEX file
totalGPSsatellite=32;
vSatID=zeros(totalGPSsatellite,1);
interval=0;
%% Initialize the data matrix.
if totalepoch==-999, %% When the total number of epochs are unknown
    vTime=[];  mCode=[]; mPhase=[]; vTimeStr=[];mDoppler=[];
elseif   totalepoch>0,%% With the known number of epochs  
    vTime=zeros(totalepoch,1);
    vTimeStr=cell(totalepoch,1);
    mCode=zeros(totalGPSsatellite,totalepoch);
    mPhase=zeros(totalGPSsatellite,totalepoch);
    mDoppler=zeros(totalGPSsatellite,totalepoch);
else
    errordlg('Invalid total number of epochs')
end

%% Open the file
fid = fopen(filename);
if fid==-1,
    errmsg=sprintf('Can not open this file%s',filename);
    errordlg(errmsg);
    return;
end
 
line80num=0;
while (1)
    line80num = line80num + 1;
    %% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    %%              Read the header section
    %% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    line80 = GetLine80(fid);
    if findstr(line80,'END OF HEADER'),
        break;
    end
    %% Antenna Height
    if findstr(line80,'ANTENNA: DELTA H/E/N'),
        vAntennaDelta(1) = str2num(line80(1:14));
        vAntennaDelta(2) = str2num(line80(15:28));
        vAntennaDelta(3) = str2num(line80(29:42));
    end
    %% Sampling rate
    if findstr(line80,'INTERVAL'),
        interval = str2num(line80(1:10));
    end
    %% Check the order of the observation
    if findstr(line80,'TYPES OF OBSERV'),
        obsC1pos=0; %% Indicates which observation is C/A code
        obsL1pos=0; %% Indicates which observation is phase on L1
        pos=0;
        totalobs=str2num(line80(6:8));
        for j=11:6:11+totalobs*6,
            pos=pos+1;
            %% For C/A code  
            if (line80(j:j+1)=='C1')  ,
                obsC1pos=pos;
            end
            %% For phase data L1
            if (line80(j:j+1)=='L1')
                obsL1pos=pos;
            end
            if (line80(j:j+1)=='D1')
                obsD1pos=pos;
            end
            
        end 
    end  
end %% End of reading header

%  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%            Read the data epoch by epoch
%  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
epoch=0;
line80 = GetLine80(fid);
while(1)
    vSatID=zeros(totalGPSsatellite,1);
    %% check whether the pointer of the file reaches the end
    if line80==-1, 
        break;
    end
    %% Exit of the file reading
    epoch=epoch+1;
    if (totalepoch~=-999), %%  -999= no limit
        if epoch>totalepoch, break; end
    end
    %% Get the time
    year = str2double(line80(1:3));
    month = str2double(line80(4:6));
    day = str2double(line80(7:9));
    hour = str2double(line80(10:12));
    minute = str2double(line80(13:15));
    second = str2double(line80(16:26));
    %% Get GPS week and seconds of week
    [gps_week,week_seconds]=GetWeekSeconds(year, month ,day, hour, minute, second);
    %% Now format the receipt epoch expressed in STRING
    vTime(epoch)=week_seconds;
    if year<80,
        vTimeStr{epoch}=sprintf('%4d.%02d.%02d  %02d:%02d:%05.2f',2000+year,month,day,hour,minute,second);
    else
        vTimeStr{epoch}=sprintf('%4d.%02d %02d  %02d:%02d:%05.2f',1900+year,month,day,hour,minute,second); 
    end
    %% Read the satellite id
    sat_num = str2double(line80(30:32));
    %% Start checking the satellite types and PRN
    %% Only record the GPS satellites, ignore the GLONASS,  GALILEO. etc
    pos=33;
    j=1;
    sat_gps=0;
    while (j<=sat_num),
        switch line80(pos)
            case 'G' %% GPS , PRN 1~50
                satid=str2num(line80(pos+1:pos+2));
                sat_gps=sat_gps+1;
                vSatID(sat_gps)=satid;
            case 'R' %% GLONASS,not to be considered
            case 'S' %% SBAS,not to be considered
            case 'E' %% GALILEO, not to be considered
                ;
            otherwise %% satellite information might be in the next line 
                pos=30; %% for next line, make a offset because it will be added with 3
                line80 = GetLine80(fid);
                j=j-1;
        end
        %% for next satellite
        pos=pos+3;
        j=j+1;
    end
    %% Now read the observation data
    %% Note that the loop is implemented for GPS satellites
    %% so we use sat_gps as the final value of the loop
    %% Cannot use sat_num which represents the total number of all GNSS satellites
    %% But sat_num is still useful to filter out the observation data of
    %% the other GNSS satellite.
    line80 = GetLine80(fid);
    for j=1:1:sat_gps,
        pos=1;
        for k=1:1:totalobs,
            %% Get the numerical value
            obs=str2num(line80(pos:pos+16-1));
            if length(obs)>1, obs=obs(1); end
            %% Record C/A data
            if k==obsC1pos,
                %% When no code data avaliable, we then use 0 instead
                if isempty(obs),
                    mCode(vSatID(j),epoch)=0;
                    str=sprintf('At epoch %d %d %d %d %d %2.2f, sat= %d  No corresponding  code data avaliable, use 0 instead', year,month,...
                        day,hour,minute,second,vSatID(j));
                     WriteListbox(str,hlistbox);
                else
                    mCode(vSatID(j),epoch)=obs;
                end
            end
            %% Record L1 phase observation
            if k==obsL1pos,
                 % When no phase data avaliable, we then use 0 instead
                if isempty(obs),
                    mPhase(vSatID(j),epoch)=0;
                    str=sprintf('At epoch %2d %2d %2d %2d %2d %2.2f, sat= %d  No corresponding  phase data avaliable, use 0 instead', year,month,...
                        day,hour,minute,second,vSatID(j)); 
                    WriteListbox(str,hlistbox);
                else
                    mPhase(vSatID(j),epoch)=obs;
                end
            end
            %% Record D1 phase observation
            if k==obsD1pos,
                 % When no phase data avaliable, we then use 0 instead
                if isempty(obs),
                    mDoppler(vSatID(j),epoch)=0;
                    str=sprintf('At epoch %2d %2d %2d %2d %2d %2.2f, sat= %d  No corresponding  phase data avaliable, use 0 instead', year,month,...
                        day,hour,minute,second,vSatID(j)); 
                    WriteListbox(str,hlistbox);
                else
                    mDoppler(vSatID(j),epoch)=obs;
                end
            end
            
            %% Next observation
            pos=pos+16;
            %% Next line80 may includ also the observation of this satellite
            if (pos>80) && (k<totalobs),
                line80 = GetLine80(fid);
                pos=1;
            end
        end %% finished reading all the observation of this GPS satellite
        line80 = GetLine80(fid);
    end %% all GPS satellites of this epoch
    %% Now there might be two cases. In one case, the RINEX file
    %% contains only GPS data, then the next paragraph indicates the
    %% observation data of the next epoch. In the other case, some non-GPS
    %% satellite data are still presented, we then filter them out by 
    %% keep searching for 'G' in each line. Finding it implies that
    %% a paragraph for the next satellite starts.
    lp=1;
    while ( lp==1)
        if findstr(line80,'G'),
            lp=0;
        else
            line80 = GetLine80(fid);
            if line80==-1, 
                break;
            end
        end
    end
    %% Show progress
    if rem(epoch,100)==0,
        str=sprintf('%d epochs processed. \n', epoch);
        WriteListbox(str,hlistbox);
    end
end
%% If the actual total number of epochs is less than the expected number
if epoch<totalepoch,
    vTime(epoch+1:totalepoch)=[];
end

%% For some receivers, RINEX observation file does not contain 
%% the term "INTERVAL". In this case, we have to calculate it.
%% By checking the time span of the first 10 epochs, we could specify the
%% interval. 
totalepoch=epoch-1;
if interval==0, 
    test_epoch=10;
    if epoch<10, test_epoch=totalepoch; end    
    interval1=vTime(2)-vTime(1);
    for i=2:1:test_epoch,
        interval_temp=vTime(i+1)-vTime(i);
        if interval1~=interval_temp,
            errordlg('Can not specify the sampling interval of the observation data!!')
            return
        end        
    end
    interval=interval1;
end

%% Check the data interruption
for i=1:1:totalepoch-1,
    if abs(vTime(i+1)-vTime(i)-interval)>interval,
        str=sprintf('Data interruption at %s',vTimeStr{i});
        errordlg(str);
        return
    end
end

%% Save it into a data file. The data file starts with a common name 
%% "DataRecord" and specified by its ID. 
datafilename=sprintf('DataRecord%1d',id);
save(datafilename,'vTime','interval','mCode','mPhase','mDoppler','vTimeStr');
WriteListbox('Processing finished.',hlistbox)
pause(.5)
 


%% Function for writting information on the ListBox
function WriteListbox(str,hlistbox)
 
mStrs=get(hlistbox,'String');
len=length(mStrs);
mStrs{len+1}=str;
set(hlistbox,'String',mStrs);
set(hlistbox,'Value',len+1);
pause(0.01)