spplab.m

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

% -----------------------------------------------------------------------------
% SPP Lab:   Absolute Positioning epoch by epoch performing Least Squares Estimation 
%            (SPP, Pseudoranging)
%            regarding aberration (earth rotation correction)
% -----------------------------------------------------------------------------
%
% comment:   The following files are needed: 
%            1 Ephemeris file, 1 Code phase file (C1) for each Station,
%            1 Epoch file (EPO) for each Station, 
%            1 Header file for each Station
%
% Output: a description of the generated variables can be found in var_spp.txt 
%
% -----------------------------------------------------------------------------
% (c) iapg 1998 Zebhauser

% 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]

% 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 SPPLAB.M while reading the Pseudoranges.',...
             ' For the Reference neither C/A- nor P-Code pseudoranges ' ...
             ' on L1 are available. Check e.g. the RINEX import.'], ...   
       'GPSLab: break -  1998 zeb');
       return;
end

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

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 SPPLAB.M while reading the Pseudoranges.',...
             ' For the Rover neither C/A- nor P-Code pseudoranges ' ...
             ' on L1 are available. Check e.g. the RINEX import.'], ...   
       'GPSLab: break -  1998 zeb');
end

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

% 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 SPPLAB.M while reading the ephemeris.',...
             ' The length of one of the files S1EPH and S2EPH is no multiple of 24' ...
             ' and so seems corrupted. Check e.g. the RINEX import.'], ...   
       'GPSLab: break -  1998 zeb');
       return;
end

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

if nsats1<4 | nsats2<4
       errordlg(['error in SPPLAB.M while reading the ephemeris.',...
             ' At least one of the files S1EPH,S2EPH contains only ephemeris for' ...
             ' less than 4 satellites. Due to this fact a position solution is unpossible.'], ...   
       'GPSLab: break -  1998 zeb');
    return;
end


% 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);


% 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
s2c1=s2c1(logical(ge4mask2),:);   % -"-
s1epo=s1epo(logical(ge4mask1),:); % -"-
s2epo=s2epo(logical(ge4mask2),:); % -"-

[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(['Break in SPPLAB.M while checking the pseudoranges.',...
             ' At least one of the files S1C1,S2C1 contains no epoche with' ...
             ' at least 4 Satelliten. So on this station no' ...
             ' position solution and so also no DGPS is possible.'], ...   
       'GPSLab: Break -  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); % -"-



% ACHTUNG f黵 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 Empfaengerpositionen fuer 2 Stationen ueber alle Epochen  ----------------
% ---> wie gross ist die Variation  gg.ueber den Naeherungskoord. (RINEX Header)
% Berechnung der Koord.unterschiede ueber alle Epochen  -----------------------------------
% ---> wie gross ist die Variation gg.ueber Naeherungskoord.differenz (RINEX-Header)

h=waitbar(0,'calculating the positions of the reference station');

for t=1:nobs1          % Schleife 黚er alle Beobachtungen
   
clear as1 ls1z ls1 vvs1;
   
versatz=0;   
for isat=1:satmask1(t) % Schleife 黚er alle beobachteten Satelliten
while s1c1(isat+versatz)==0,
   versatz=versatz+1;         % fehlende Beobachtungen innerhalb einer Zeile 黚erspringen
end

ephsatz=1;
while s1eph(ephsatz)~=s1epo(t,isat+1+versatz)
   ephsatz=ephsatz+24;        % PRN raussuchen (einfach die erste passende)
end

[s1x(t,isat),s1y(t,isat),s1z(t,isat),s1dts(t,isat),s1rel(t,isat)]=eph2xyzn(s1epo(t,1),s1c1(t,isat+versatz),s1eph(ephsatz:(ephsatz+23)));
[s1b(t,isat),s1l(t,isat),s1h(t,isat)]=xyz2blh(s1x(t,isat),s1y(t,isat),s1z(t,isat));
% Erdrot.korr. Empf.koord: Ss20,dOmega_e,Xe0,Ye0,Ze0 je (t,isat)(nicht indizieren, aber je s2 und s1)