dgpslab.m

这是国外关于卫星导航方面一书的源代码

% -----------------------------------------------------------------------------
% GPSLab:    DGPSLAB.M
%            Relative positioning s1->s2 by pseudorange double differencing
% -----------------------------------------------------------------------------
%            
%            *** erweitert um eine Zweifrequenz- (=ionosph鋜enfreie) L鰏ung ***
%            Vor dem Aufruf dieses Programms sollte der Parameter "what_modl" definiert werden:
%            M鰃liche Werte sind:   what_modl = 'ionosph鋜enfreie L鰏ung'
%                             oder  what_modl = 'C/A-Code-L鰏ung'
%
% Kommentar: Die folgenden Files sind noetig: 
%            1 Eph.file, mind. 1 C/A-Codephasenfile (C1) je Station,
%            1 P2-Codephasenfile je Station bei Zweifrequenzl鰏ung
%            1 Epochenfile (EPO) je Station, 
%              welcher in der 1. Spalte die Epochen und in den darauffolgenden 
%              die Satellitennummern, wie die Beobachtungen je Epoche enth鋖t.
%            1 Headerfiles je Station
%
%            Achtung: Zweifrequenzl鰏ung nur bei P/W-Tracking-Receivern 
%                     besser als reine C/A-Codel鰏ung
%
% Output: Die erzeugten Variablen sind in der Datei var_dgps.txt nachzulesen.
%
% -----------------------------------------------------------------------------
% (c) iapg 1998, 1999 Zebhauser


if ~exist('what_modl')                   % wenn der Parameter what_modl noch nicht definiert wurde
what_modl=questdlg('Which processing model should be used for Pseudorange-DGPS       (if P2 exists)?', ...
                         'GPS Lab: processing model for DGPS -  1999 zeb', ...
                         'ionosphere free solution','C/A-Code solution','ionosphere free solution (default)');
end
   
     switch what_modl,
       case 'ionosphere free solution', 
          ionofree=1;
       case 'C/A-Code solution',
          ionofree=0;
       case 'ionosphere free solution (default)',
          errordlg(['Break at model selection (either ionosphere free solution or not). ' ...
                    'default) is ionosphere free (if P2 exists).'], ...   
                    'GPSLab: Break -  1999 zeb');
          ionofree=1;
      return;
     end % switch
     
% "what_modl" enth鋖t String f黵 das Auswertemodell


% Naeherungskoordinaten -------------------------

xs1=4231107.740; % K黨alm
ys1= 839574.578;
zs1=4684782.848;

xs2=4233906.848; % Kr黱
ys2= 843981.554;
zs2=4680440.960;

load s1head;
load s2head;

xs1=s1head(1);
ys1=s1head(2);
zs1=s1head(3);

xs2=s2head(1);
ys2=s2head(2);
zs2=s2head(3);

[app1_b,app1_l,app1_h]=xyz2blh(xs1,ys1,zs1);
[app2_b,app2_l,app2_h]=xyz2blh(xs2,ys2,zs2);

% Antennenh鰄e auf den beiden Punkten
ants1=s1head(4); ants2=s2head(4);   
deg2rad=pi/180;

% Verbesserungen wg. Antennenh鰄en (*** ACHTUNG: XYZ2BLH liefert DEG, nicht RAD ! ***)
ants2dx=ants2*cos(app2_b*deg2rad)*cos(app2_l*deg2rad);
ants2dy=ants2*cos(app2_b*deg2rad)*sin(app2_l*deg2rad);
ants2dz=ants2*sin(app2_b*deg2rad);

ants1dx=ants1*cos(app1_b*deg2rad)*cos(app1_l*deg2rad);
ants1dy=ants1*cos(app1_b*deg2rad)*sin(app1_l*deg2rad);
ants1dz=ants1*sin(app1_b*deg2rad);


% WGS84 Konstanten  ------------------------------------------------------------------
c0=299792458;                   % [m/s]
omega_e=7.2921151467e-5;        % [rad/s]
f1q=(10.23e6*154)^2;         % quadrierte L1-Frequenz
f2q=(10.23e6*120)^2;         % quadrierte L2-Frequenz
f2pr1=f1q/(f1q-f2q);         % Term aus quadr. Frequ. zu Pseudorange auf L1 f黵 iono.free
f2pr2=f2q/(f1q-f2q);         % Term aus quadr. Frequ. zu Pseudorange auf L1 f黵 iono.free

% Data loading  ----------------------------------------------------------------------
load s1eph;                     % Ephemeriden Ref
load s2eph;                     % Ephemeriden Rov

load s1epo;                     % Epochen/Satelliten Ref

if s1head(9)==1
   load s1c1;                      % C/A-Codephasen (PR's) Ref
elseif s1head(9)==0 & s1head(10)==1
   load s1p1;                      % P-Codephasen (PR's) Ref, falls C/A nicht vorhanden
   s1c1=s1p1;
   clear s1p1;
else
       errordlg(['Error in program DGPSLAB.M during the reading of the Pseudoranges.',...
             ' For Reference neither C/A- nor P-Codephase measurements on L1 are available.' ...
             ' Check the RINEX-Import.'], ...   
       'GPSLab: Break -  1998 zeb');
       return;
end

% load s1l1;                      % L1-Tr鋑erphasen Ref
% load s1l2;                      % L2-Tr鋑erphasen Ref

load s2epo;                     % Epochen/Satelliten Rov

if s2head(9)==1
   load s2c1;                      % C/A-Codephasen (PR's) Rover
elseif s2head(9)==0 & s2head(10)==1
   load s2p1;                      % P-Codephasen (PR's) Rover, falls C/A nicht vorhanden
   s2c1=s2p1;
   clear s2p1;
else
       errordlg(['Error in program DGPSLAB.M during the reading of the Pseudoranges.',...
             ' For Rover neither C/A- nor P-Codephase measurements on L1 are available.' ...
             ' Check the RINEX-Import.'], ...   
       'GPSLab: Break -  1998 zeb');
end

% load s2l1;                      % L1-Tr鋑erphasen Rov
% load s2l2;                      % L2-Tr鋑erphasen Rov

% !!!!!!!  P2 vorhanden ? Wenn ja, verwenden f黵 ionosph鋜enfreie L鰏ung !!!!!!!
%          -------------------------------------------------------------

if s1head(11)==0 | s2head(11)==0 | ionofree==0
       mit_p2=0;
else
       load s1p2;                      % Pseudoranges auf L2 am Referenzpunkt
       load s2p2;                      % Pseudoranges auf L2 am Roverpunkt
       mit_p2=1;
end

% Dimensionen (Anz. Satelliten, Anz. Epochen)  ---------------------------------------
[lephs1,dummy]=size(s1eph);
[lephs2,dummy]=size(s2eph);
[nobs1,dummy]=size(s1epo);  % muss identisch mit s1c1 sein
[nobs2,dummy]=size(s2epo);  % -"-

% Checks  ----------------------------------------------------------------------------
if rem(lephs1,24)~=0 | rem(lephs2,24)~=0 
       errordlg(['Error in program AMFLAB.M during the reading of the Ephemeris.',...
             ' The length iof one of the files S1EPH,S2EPH is not equal to the multiple of 24' ...
             ' and so seems corrupted. For example check the RINEX-Import ...'], ...   
       'GPSLab: Break -  1998 zeb');
       return;
end

nsats1=lephs1/24;
nsats2=lephs2/24;

if nsats1<4 | nsats2<4
       errordlg(['Error in program AMFLAB.M during the reading of the Ephemeris.',...
             ' At least one of the files S1EPH,S2EPH contains Ephemeris for' ...
             ' less than 4 Satellites. So no position calculation is possible.'], ...   
       'GPSLab: Break -  1998 zeb');
    return;
end

obs_is_used=0;            % Nullsetzen des Indexes der tats鋍hlich genutzten Beobachtungen


% Elimination von nicht plausiblen Beobachtungsgr鲞en --------------------------
% Nur Pseudoranges >10000 km und <30000 km zugelassen

s1c1=s1c1.*(s1c1>10000000 & s1c1<30000000); % Null gesetzt
s2c1=s2c1.*(s2c1>10000000 & s2c1<30000000);

if mit_p2==1
s1p2=s1p2.*(s1p2>10000000 & s1p2<30000000); % Null gesetzt
s2p2=s2p2.*(s2p2>10000000 & s2p2<30000000);
end

% Masken -----------------------------------------------------------------------

satmask1=sum(s1c1>0,2);
satmask2=sum(s2c1>0,2);

ge4mask1=satmask1>3;   % Bei welcher abgespeicherten Epoche mehr als 3 Sat. ? (S1)
ge4mask2=satmask2>3;   % Bei welcher abgespeicherten Epoche mehr als 3 Sat. ? (S2)

s1c1=s1c1(logical(ge4mask1),:);   % Verk黵zung der Matrizen auf die Zeilen mit > 3 Sat.s auf C1
s2c1=s2c1(logical(ge4mask2),:);   % -"-
s1epo=s1epo(logical(ge4mask1),:); % -"-
s2epo=s2epo(logical(ge4mask2),:); % -"-

if mit_p2==1
s1p2=s1p2(logical(ge4mask1),:);   % -"-
s2p2=s2p2(logical(ge4mask2),:);   % -"-
end

[nobs1,dummy]=size(s1epo);  % Neuberechnung der L鋘ge des Beobachtungsvektors (enth鋖t nur noch Epochen mit mind. 4 Sat.s)
[nobs2,dummy]=size(s2epo);  % -"-

if nobs1==0 | nobs2==0
       warndlg(['Error in program AMFLAB.M during the check of the observation files.',...
             ' At least one of the files S1C1, S2C1 contains' ...
             ' no single epoch with at least 4 Satellites. So it is e.g. not possible' ...
             ' to get a DGPS-solution.'], ...   
       'GPSLab: Abbruch -  1998 zeb');
    return;
end

satmask1=sum(s1c1>0,2); % Neuberechnung der Maske mit der Anz. der beob. Sat.s, passend zu nobs1 bzw. nobs2
satmask2=sum(s2c1>0,2); % -"-

% --- Runden der Epochen auf n鋍hste Integer, Warnung, wenn diff>10 msec !!
%     Verdacht auf unganzzahlige Epochen wie halbe Sekunden o.? ...
if max(s1epo(:,1)-round(s1epo(:,1)))>.01 | max(s2epo(:,1)-round(s2epo(:,1)))>.01
   warndlg(['Error in program AMFLAB.M during the check of the epochs (time tags).',...
             ' At least one of the files S1EPO,S2EPO contains epochs with' ...
             ' offsets from integer seconds > 10 msec or a registration rate' ...
             ' < 1 sec. Such data were not computed by this module.'], ...   
       'GPSLab: Break -  1999 zeb');
   return;
end

% Vorher: Sichern der Original-Epochen (nur diese, keine PRN-Nummern) mit etwaigen 
% Millisekunden-Spr黱gen f黵 die richtige Berechnung der Satellitenkoordinaten; 
% die gerundeten Epochen jedoch vereinfachen die Abfrage nach "simultanen" Epochen,
% um dort double differences bilden zu k鰊nen:

s1epoms=s1epo(:,1);
s2epoms=s2epo(:,1);

s1epo(:,1)=round(s1epo(:,1));
s2epo(:,1)=round(s2epo(:,1));

% --- Suchen der simultanen Epochen auf S1 und S2 -----------------------------------
epomask1=zeros(nobs1,1);
epomask2=zeros(nobs2,1);
   
such1=1;such2=1;

while such1<=nobs1 & such2<=nobs2
   if s1epo(such1,1)<s2epo(such2,1)
      such1=such1+1;
   elseif s1epo(such1,1)>s2epo(such2,1)
      such2=such2+1;
   else
      epomask1(such1)=1;
      epomask2(such2)=1;
      such1=such1+1;
      such2=such2+1;
   end
end

% --- Reduzierung der Felder auf die Elemente mit simultanen Epochen -----------------
s1c1=s1c1(logical(epomask1),:);   % ...durch Herausstreichen der anderen Beobachtungen
s1epo=s1epo(logical(epomask1),:); % -"-
s2c1=s2c1(logical(epomask2),:);   % -"-
s2epo=s2epo(logical(epomask2),:); % -"-

s1epoms=s1epoms(logical(epomask1),:);  % -"-
s2epoms=s2epoms(logical(epomask2),:);  % -"-

if mit_p2==1
s1p2=s1p2(logical(epomask1),:);   % -"-
s2p2=s2p2(logical(epomask2),:);   % -"-
end


[nobs1,dummy]=size(s1epo);  % Neuberechnung der L鋘ge des Beobachtungsvektors (enth鋖t nur noch Epochen mit simultanen Epochen auf s1 und s2)
[nobs2,dummy]=size(s2epo);  % -"-

if nobs1==0 | nobs2==0
       warndlg(['Error in program AMFLAB.M during the check of the observation files.',...
             ' The files S1C1,S2C1 contain no simultaneous observations' ...
             ' with at least 4 satellites.' ...
             ' So no DGPS is possible.'], ...   
       'GPSLab: Break -  1999 zeb');
    return;
end

satmask1=sum(s1c1>0,2); % Neuberechnung der Maske mit der Anz. der beob. Sat.s, passend zu nobs1 bzw. nobs2
satmask2=sum(s2c1>0,2); % -"-


% bisher nicht, bei Gl鋞tung n鰐ig: s1l1,s2l1

% ACHTUNG f黵 S2 (=ROVER) bzgl. S1 (=REF.)
% Berechnung der Position aller Sat. fuer 2 Stationen und alle Epochen  ------------------&
% Berechnung der Orbithoehen ueber Grund aller Sat. fuer 2 Stat. fuer alle Epochen  -------
% Korrektion der Empfaengerkoordinaten wegen Erdrotation  ---------------------------------
% Berechnung der Vektoren zwischen den 2 Stationen ueber alle Epochen  ----------------
% ---> wie gross ist die Variation gg.ueber Naeherungskoord.differenz (RINEX-Haeder)