amflab.m

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

% -----------------------------------------------------------------------------
% GPSLab:    AMFLAB.M
%            Relative positioning s1->s2 by double difference ambiguity mapping function
% -----------------------------------------------------------------------------
%
% Kommentar: Die folgenden Files sind noetig: 
%            1 Eph.file, 1 Codephasenfile (C1) je Station,
%            2 Traegerphasenfiles je Station (L1 und  L2 je Station), 
%            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
%
% AMF="ambiguity mapping function", manchmal auch AFM="ambiguity function method" genannt
%
% Die Variablen "trop_typ" und "el_mask" sollten vor dem Start dieses Programms gesetzt werden;
% Dies erreicht man i.a. durch Aufruf des Skript-Files OPTIONS.M; ansonsten wird als default
% gesetzt: trop_typ='mh' , el_mask=0 , d.h. Troposph鋜enkorrekturmodell Modified Hopfield und 
% Elevationsmaske 0?(also alle Beobachtungen werden verwendet).
%
% Output: Die erzeugten Variablen sind in der Datei var_amf.txt nachzulesen.
%
% Literatur: Anwendung der Ambiguity Function bei GPS:
%               Counselman+Gourevitch(1981), Remondi(1984, S.148 ff.), Mader(1992/1)+(1992/2), 
%               Han+Rizos(1996), Lachapelle(1992), Lu+Lachapelle(1992), W黚bena(1991),
%               Abidin+Wells+Kleusberg(1991), Abidin+Wells+Kleusberg(1992),
%               Leick(1995, S.375 ff.), Hofmann-Wellenhof(1991, S.198 ff.), ...
%            Grundlagen zur Ambiguity Function (aus der Radartechnik):
%               Urkowitz(1966), Wohlleben+Mattes(1973), Ulaby+Moore+Fung(1982)
%
% -----------------------------------------------------------------------------
% (c) iapg 1999 Zebhauser



load s1head;
load s2head;


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

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


% Richtungskosinusse f黵 die Ortsvektoren der Empf鋘gerstandpunkte
s1str=sqrt(xs1^2+ys1^2+zs1^2);
s2str=sqrt(xs2^2+ys2^2+zs2^2);
s1ricox=xs1/s1str;
s1ricoy=ys1/s1str;
s1ricoz=zs1/s1str;
s2ricox=xs2/s2str;
s2ricoy=ys2/s2str;
s2ricoz=zs2/s2str;

% B,L,H auf beiden Punkten
[app1_b,app1_l,app1_h]=xyz2blh(xs1,ys1,zs1);
[app2_b,app2_l,app2_h]=xyz2blh(xs2,ys2,zs2);

% Standard atmosphere bei H鰄e 0 m
ttt=19.85;	% [癈] - z.T. sind Formeln auf 癈, z.T. auf 癒 ausgelegt (T=t+273.15)
ppp=1013;		% [mbar]
hhh=50;		% [Prozent]

% Standardatmosph鋜e f黵 Referenz und Rover
temp1 = ttt-0.0065*app1_h;
druck1= ppp*(1-0.0000226*app1_h)^(5.225);
humid1= hhh*exp(-0.0006396*app1_h);

temp2 = ttt-0.0065*app2_h;
druck2= ppp*(1-0.0000226*app2_h)^(5.225);
humid2= hhh*exp(-0.0006396*app2_h);

% 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] - Achtung: offizieller Wert des WGS84; jede 膎derung dieses c0-Werts bedeutet eine 膎derung des Systemma遱tabs
omega_e=7.2921151467e-5;        % [rad/s] - Achtung: offizieller WGS84-Wert
lambda1=c0/10.23e6/154;         % L1-Wellenl鋘ge
lambda2=c0/10.23e6/120;         % L2-Wellenl鋘ge

if ~exist('el_mask')                  % wenn OPTIONS.M noch nicht aufgerufen wurde
   el_mask=10;
end

if ~exist('raster')                   % wenn OPTIONS.M noch nicht aufgerufen wurde
   raster=.01;
end

if ~exist('comment')                  % wenn OPTIONS.M noch nicht aufgerufen wurde
   comment='Project: , Baseline: , Supervisor: ';
end

if ~exist('trop_typ')                   % wenn OPTIONS.M nicht aufgerufen wurde
   trop_typ='mh';
end

if ~exist('what_trop')                  % wenn OPTIONS.M nicht aufgerufen wurde
   what_trop='Modified Hopfield';
end
% --> Bem.: Die Variable "what_trop" enth鋖t einen String mit der Bezeichnung des Modells

% Nullsetzen mit Initialisierung von Variablen/Matrizen:
ddamfall=zeros(21,21,21); % Nullsetzen der Summationsmatrix 黚er alles                          
used_all=0;               % Beob.z鋒ler f黵 die Anzahl von L1- und L2-Beob. 黚er alles
n_accu=0;                 % Nullsetzen der aufakkumulierten NGlMatrix
obs_is_used=0;            % Nullsetzen des Indexes der tats鋍hlich genutzten Beobachtungen

% Gitter erstellen f黵 dx,dy,dz:
%    --> Multipl. mit lambda kostet zwar doppelten Speicherplatz, 
%        daf黵 weniger Rechenzeit, weil au遝rhalb der Schleife
% #############
%raster=.1;    % in [m] ; in Zukunft au遝rhalb definieren, je nach gew黱schter Rasterweite
% #############
[dxcube1,dycube1,dzcube1]=ndgrid(-10:10,-10:10,-10:10);

% &&&&&&&&&&&&&&&&&&&&&&&&&&&&??????????????####################???????????????
% evtl. erst sp鋞er mit Abfrage, ob mit_l2==1, um Speicherplatz zu sparen !!!
dxcube2=raster./lambda2.*dxcube1;
dycube2=raster./lambda2.*dycube1;
dzcube2=raster./lambda2.*dzcube1;
% &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

dxcube1=raster./lambda1.*dxcube1;
dycube1=raster./lambda1.*dycube1;
dzcube1=raster./lambda1.*dzcube1;

% 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 AMFLAB.M during the reading of the Pseudoranges.',...
             ' For Reference neither C/A- nor P-Codephases ' ...
             ' at L1 exists. For example check the RINEX-Import ...'], ...   
       'GPSLab: Break -  1998 zeb');
       return;
end
load s1l1;                      % L1-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 AMFLAB.M during the reading of the Pseudoranges.',...
             ' For Rover neither C/A- nor P-Codephases ' ...
             ' at L1 exists. For example check the RINEX-Import ...'], ...   
       'GPSLab: Break -  1998 zeb');
       return;
end
load s2l1;                      % L1-Tr鋑erphasen Rov

% L2 vorhanden ? Und falls ja, squaring ? Dann darauf verzichten !

if s1head(14)==2 | s2head(14)==2 | s1head(8)==0 | s2head(8)==0
   % keine L2-Daten ben黷zen, da halbe cycles m鰃lich (squaring)
%   warndlg(['Im diesem Lauf des Programms AMFLAB.M werden keine L2-Beobachtungen verwendet,',...
%             ' da sie entweder nicht vorhanden oder aber quadriert sind, ' ...
%             ' d.h. auch halbe cycles und cycle ambiguities vorkommen k鰊nen. ' ...
%             ' Solche sind jedoch in der Ambiguity Mapping Function nicht verwertbar.'], ...   
%          'GPSLab: Hinweis -  1998 zeb');
       mit_l2=0;
else
       load s1l2;                      % L2-Tr鋑erphasen Referenzpunkt
       load s2l2;                      % L2-Tr鋑erphasen Roverpunkt
       mit_l2=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

% Elimination von nicht plausiblen Beobachtungsgr鲞en (PR) unterbleibt bei AMF ------
% ###### evtl. noch einf黦en (黚ernehmen aus DGPSLab) #####

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

satmask1=sum(s1c1>0,2);  % Anzahl der beob. Satelliten je Epoche auf Referenzpunkt
satmask2=sum(s2c1>0,2);  % Anzahl der beob. Satelliten je Epoche auf Roverpunkt

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),:); % -"-
s1l1=s1l1(logical(ge4mask1),:);   % Verk黵zung der Matrizen auf die Zeilen mit > 3 Sat.s
s2l1=s2l1(logical(ge4mask2),:);   % -"-
if mit_l2==1
   s1l2=s1l2(logical(ge4mask1),:);   % Verk黵zung der Matrizen auf die Zeilen mit > 3 Sat.s
   s2l2=s2l2(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 S1L1, S2L1 as also S1C1, S2C1 contains' ...
             ' no single epoch with at least 4 Satellites. So it is e.g. not possible' ...
             ' to get a CDGPS-solution.'], ...