amflab.m

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

       '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
s1l1=s1l1(logical(epomask1),:);   % ...durch Herausstreichen der anderen Beobachtungen
s1epo=s1epo(logical(epomask1),:); % -"-
s2c1=s2c1(logical(epomask2),:);   % -"-
s2l1=s2l1(logical(epomask2),:);   % -"-
s2epo=s2epo(logical(epomask2),:); % -"-

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

if mit_l2==1
   s1l2=s1l2(logical(epomask1),:);   % ...durch Herausstreichen der anderen Beobachtungen
   s2l2=s2l2(logical(epomask2),:);   % ...durch Herausstreichen der anderen Beobachtungen
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,S1L1,S2L1  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); % -"-

% =========================================================================================
% Berechnung f黵 S2 (=ROVER) bzgl. S1 (=REF.)
% Berechnung der Orbits aller Sat. fuer 2 Stat. und fuer alle Epochen  --------------------
% Korrektion der Empfaengerkoordinaten wegen Erdrotation  ---------------------------------
% Berechnung der double diff. zwischen den 2 Stationen und allen Satelliten ueber alle Epochen  ----------------
% etc. ...
% ---> Abweichung gg.ueber Naeherungskoord.differenz (RINEX-Haeder) ???

h=waitbar(0,['L1L2-AMF, 21x21x21 (',num2str(raster),' m-) Grid, all epochs']);

for t=1:nobs1          % Schleife 黚er alle Beobachtungen #################################
   
clear as2 ;   % jede Epoche nicht nur neu, sd. bei wechselnder Sat.anz. auch verschieden gro? Daher vor jeder Epoche neu l鰏chen.
   
versatz=0;   

isat=0;             % number of double diff.: simultan auf beiden Stationen beob. Satelliten


ddamfobs=zeros(21,21,21); % Nullsetzen der Summationsmatrizen 黚er die Satelliten, die erst
                          %    zur jeweiligen Gesamtmatrix addiert werden, wenn alle Sats
                          %    der Epoche abgearbeitet und es 黚er 4 Sats sind. 
                          % Bei jeder Epoche wird wieder bei Null begonnen. [74088 KB]
                          
used_obs1=0;               % Beob.z鋒ler f黵 die Anzahl von L1-Beob. w鋒rend 1 Epoche
used_obs2=0;               % Beob.z鋒ler f黵 die Anzahl von L2-Beob. w鋒rend 1 Epoche
                          

for ksat=1:satmask1(t) % Schleife 黚er alle beobachteten Satelliten #####################
   
while s1c1(ksat+versatz)==0,
   versatz=versatz+1;         % fehlende Beobachtungen innerhalb einer Zeile 黚erspringen
end

% 蹷ERPR蹻EN, OB AUF S1 EMPF. SAT AUCH AUF S2 DA:

% Zu Beob. in S1L1 passende Beob. in S2L1:
rovmask=(s2epo(t,:)==s1epo(t,ksat+1+versatz));  % Maske auf Zeile in Rover (dieselbe epoche): Suche nach demselben Sat. dort

if sum(rovmask)>0   % ->>> nur wenn Sat auch auf Rover vorhanden, weitermachen ...

   
rovindx=find(rovmask)-1; % an welcher Stelle in der Zeile bei S2L1 sich dieselbe Sat.nr. befindet

isat=isat+1;       % number of double diff.: simultan auf beiden Stationen beob. Satelliten

% richtige PRN aus Eph raussuchen (einfach die erste passende aus S1EPH)
ephsatz=1;           
while s1eph(ephsatz)~=s1epo(t,ksat+1+versatz)
   ephsatz=ephsatz+24;        
end

% Koordinaten auch f黵 S2, d.h. mit denselben eph, aber mit s2c1 !!!
[s1x(t,isat),s1y(t,isat),s1z(t,isat),s1dts(t,isat),s1rel(t,isat)]=eph2xyzn(s1epoms(t,1),s1c1(t,ksat+versatz),s1eph(ephsatz:(ephsatz+23)));
[s2x(t,isat),s2y(t,isat),s2z(t,isat),s2dts(t,isat),s2rel(t,isat)]=eph2xyzn(s2epoms(t,1),s2c1(t,rovindx),s1eph(ephsatz:(ephsatz+23)));

[s1b(t,isat),s1l(t,isat),s1h(t,isat)]=xyz2blh(s1x(t,isat),s1y(t,isat),s1z(t,isat));
[s2b(t,isat),s2l(t,isat),s2h(t,isat)]=xyz2blh(s2x(t,isat),s2y(t,isat),s2z(t,isat));

% Erdrot.korr. Empf.koord: sS1_0,dOmega_e,Xe0,Ye0,Ze0 je (t,isat)(nicht indizieren, aber je s2 und s1)
%
% ber. range je Satellit fuer beide Punkte
rs1=sqrt((s1x(t,isat)-xs1)^2+(s1y(t,isat)-ys1)^2+(s1z(t,isat)-zs1)^2);
% Winkel aus Erdrotation
dos1=rs1*omega_e/c0;
% um Erdrot. und um Antennenh鰄e korr. Empf.koord
xs10=xs1*cos(dos1)-ys1*sin(dos1)+ants1dx;
ys10=ys1*cos(dos1)+xs1*sin(dos1)+ants1dy;
zs10=zs1+ants1dz;

% Erdrot.korr. Empf.koord: sS2_0,dOmega_e,Xe0,Ye0,Ze0 je (t,isat)(nicht indizieren, aber je s2 und s1)
%
% ber. range je Satellit fuer beide Punkte
rs2=sqrt((s2x(t,isat)-xs2)^2+(s2y(t,isat)-ys2)^2+(s2z(t,isat)-zs2)^2); % Sat.koord-Rec.koord
% Winkel aus Erdrotation
dos2=rs2*omega_e/c0;                                                   % [m*rad/s/(m/s)=rad]
% um Erdrot. und um Antennenh鰄e korr. Empf.koord
xs20=xs2*cos(dos2)-ys2*sin(dos2)+ants2dx;
ys20=ys2*cos(dos2)+xs2*sin(dos2)+ants2dy;
zs20=zs2+ants2dz;

% #############################################################################################
% AMF-Ansatz:
% Funktion berechnen:  (t,isat)
% (nur A bzgl. Sat indizieren, Rest wird nicht wieder gebraucht, allerdings je s2 und s1)
% #############################################################################################

% Strecken, um Erdrotation und Antennenh鰄e korrigiert
ss1=sqrt((s1x(t,isat)-xs10)^2+(s1y(t,isat)-ys10)^2+(s1z(t,isat)-zs10)^2);
ss2=sqrt((s2x(t,isat)-xs20)^2+(s2y(t,isat)-ys20)^2+(s2z(t,isat)-zs20)^2); % Sat.koord-Rec.koord

% Richtungskosinusse f黵 die Vektoren Empf鋘ger --> Satellit (f黵 die Epoche, je Station, f黵 den Satellit)
rcs1(isat,1)=(s1x(t,isat)-xs10)/ss1;            
rcs1(isat,2)=(s1y(t,isat)-ys10)/ss1;
rcs1(isat,3)=(s1z(t,isat)-zs10)/ss1;
rcs2(isat,1)=(s2x(t,isat)-xs20)/ss2;
rcs2(isat,2)=(s2y(t,isat)-ys20)/ss2;
rcs2(isat,3)=(s2z(t,isat)-zs20)/ss2;

% Berechnung der Zenitdistanz und der troposph鋜ischen Korrektur (f黵 die Epoche, je Station, f黵 den Satellit)
s1ko_z=s1ricox*rcs1(isat,1)+s1ricoy*rcs1(isat,2)+s1ricoz*rcs1(isat,3);
s2ko_z=s2ricox*rcs2(isat,1)+s2ricoy*rcs2(isat,2)+s2ricoz*rcs2(isat,3);


% Abfrage wg. Elevationsmaske:
% d.h. Sat. in dieser Epoche 黚erspringen, wenn Beob. unterhalb der Elevationsmaske

if abs(s1ko_z)>cos((90-el_mask)/180*pi)&abs(s2ko_z)>cos((90-el_mask)/180*pi)    
   % je gr鲞er z, desto kleiner E und desto kleiner cos(z)
   % Ende dieses if-Blocks bei Ende der for-Schleife 黚er ksat (alle Sat.)
   
   
% Auswahl und Anbringen eines Troposph鋜enmodells: Berechnung der Korrektur wegen tropsph鋜ischer Refraktion

if trop_typ=='mh'
s1tropo=tropo_mh(s1ko_z,temp1,druck1,humid1);
s2tropo=tropo_mh(s2ko_z,temp2,druck2,humid2);
elseif trop_typ=='sa'
s1tropo=tropo_sa(s1ko_z,temp1,druck1,humid1);
s2tropo=tropo_sa(s2ko_z,temp2,druck2,humid2);
else
    warndlg(['Error in program AMFLAB.M when applying the troposphere model.',...
             ' No valid model was selected.' ...
             ' The processing will be stopped.'], ...   
             'GPSLab: Break -  1999 zeb');
    return;
end

lss1(isat)=ss1+s1tropo; % Gerechnete Strecke Sat.-->S1, korr. wg. Aberration und wg. troposph. Refraktion
lss2(isat)=ss2+s2tropo; % Gerechnete Strecke Sat.-->S2, korr. wg. Aberration und wg. troposph. Refraktion


% --------------------- phi(obs)-phi(calc) f黵 double difference afm -----------------------------------------------

% Referenzsatellit: Hier m黶ste theoretisch der Auswahlalgorithmus stehen
%                   Hier wird aber einfach der je Epoche erstbeste Satellit als Ref.sat. verwendet

% wenn isat=1, Bildung von REFSAT-Variablen, 
% wenn isat>1, Bildung von Richtungscos. bzgl. X,Y,Z (A-Matrix-Zeile),
%              Multiplikation dieser mit der Versatzmatrix (4-dim.) dividiert durch lambda1/2
%              ergibt richtungstransformierte Versatzmatrix (3-dim.) in cycles,
%              diese wird vom Widerspruchsvektor (in cycles) abgezogen und mit pi multipliziert,
%              --> cos davon ergibt den afm-Wert f黵 aktuelle Epoche und Satelliten
%              Nach allen Satelliten die Summe berechnen (--> wirklich nur bei >3, d.h. >2 dd ???), 
%              Anz. der Beob. (bei >3, d.h. >2 dd ???) zu Gesamtsumme addieren (f黵 Normierung & Statistik)
%              dasselbe parallel f黵 L1 und L2 ausf黨ren, 
%              in der afm-Summationsmatrix aufaddieren ... --> n鋍hste Epoche bearbeiten
%              Zum Schlu?Matrix normieren, auswerten und Statistik erstellen.

% Checks: Auch wenn C1 existiert, hei遲 das noch lange nicht, da?auch L1 oder gar L2 da sind !!!

if isat==1
refs=lss1(1)-lss2(1);   % Differenz der 2 Strecken von den Bodenpunkten zum Ref.sat.
refx=rcs2(1,1);         % Richtungscosinus bzgl. X-Komponente des Roverpunkts und Ref.sat.
refy=rcs2(1,2);         % Richtungscosinus bzgl. Y-Komponente des Roverpunkts und Ref.sat.
refz=rcs2(1,3);         % Richtungscosinus bzgl. Z-Komponente des Roverpunkts und Ref.sat.
refs1indx=ksat+versatz;        % Wo Refsat in s1c1 steht
refs2indx=rovindx;             % Wo Refsat in s2c1 steht
end     % Ende des IF-Blocks [if isat==1]


% Die 黚rigen Satelliten ...

if isat>1
as2(isat-1,1)=-rcs2(isat,1)+refx;     % Double Diff. Richtungscosinus bzgl. X f黵 Rover-Koord.鋘derung
as2(isat-1,2)=-rcs2(isat,2)+refy;     % Double Diff. Richtungscosinus bzgl. Y f黵 Rover-Koord.鋘derung
as2(isat-1,3)=-rcs2(isat,3)+refz;     % Double Diff. Richtungscosinus bzgl. Z f黵 Rover-Koord.鋘derung

% F黵 L1
if s2l1(t,rovindx)~=0 & s2l1(t,refs2indx)~=0 & s1l1(t,ksat+versatz)~=0 & s1l1(t,refs1indx)~=0