initial commit

Author: scienceopen

anheader.m

@@ -0,0 +1,43 @@
+function [Obs_types, ant_delta,ifound_types,eof] = anheader(file)
+%ANHEADER Analyzes the header of a RINEX file and outputs
+%	       the list of observation types and antenna offset.
+%	       End of file is flagged 1, else 0. Likewise for the types.
+%	       Typical call: anheader('pta.96o')
+
+%Kai Borre 09-12-96
+%Copyright (c) by Kai Borre
+%$Revision: 1.0 $   $Date: 1997/09/23  $
+
+fid = fopen(file,'rt');
+eof = 0;
+ifound_types = 0;
+Obs_types = [];
+ant_delta = [];
+
+while 1			   % Gobbling the header
+   line = fgetl(fid);
+   answer = findstr(line,'END OF HEADER');
+   if  ~isempty(answer), break; end;
+   if (line == -1), eof = 1; break; end;
+   answer = findstr(line,'ANTENNA: DELTA H/E/N');
+   if ~isempty(answer)
+      for k = 1:3
+         [delta, line] = strtok(line);
+         del = str2num(delta);
+         ant_delta = [ant_delta del];
+      end;
+   end
+   answer = findstr(line,'# / TYPES OF OBSERV');
+   if ~isempty(answer)
+      [NObs, line] = strtok(line);
+      NoObs = str2num(NObs);
+      for k = 1:NoObs
+         [ot, line] = strtok(line);
+         Obs_types = [Obs_types ot];
+      end;
+      ifound_types = 1;
+   end;
+end;
+
+%fclose(fid);
+%%%%%%%%% end anheader.m %%%%%%%%%

b_point.m

@@ -0,0 +1,59 @@
+function pos = b_point(Eph,pr,sv,time) %%%%%
+%B_POINT Prepares input to the Bancroft algorithm for finding
+%	      a preliminary position of a receiver. The input is
+%	      name of ephemeris file, four or more pseudoranges
+%        and the coordinates of the satellites.
+
+%Kai Borre and C.C. Goad 11-24-96
+%Copyright (c) by Kai Borre
+%$Revision: 1.0 $ $Date: 1997/09/26 $
+
+v_light = 299792458;
+dtr = pi/180;
+
+m = size(pr,1);
+%%%%Eph = get_eph(efile); 
+for t = 1:m
+   col_Eph(t) = find_eph(Eph,sv(t),time);
+end
+
+pos = zeros(4,1);	      % First guess Earth center
+for l = 1:2
+   B = [];
+   for jsat = 1:m
+      k = col_Eph(jsat);
+      tx_RAW = time - pr(jsat)/v_light;
+      toc = Eph(21,k);
+      dt = check_t(tx_RAW-toc);
+      tcorr = (Eph(2,k)*dt + Eph(20,k))*dt + Eph(19,k);
+      tx_GPS = tx_RAW-tcorr;
+      dt = check_t(tx_GPS-toc);
+      tcorr = (Eph(2,k)*dt + Eph(20,k))*dt + Eph(19,k);
+      tx_GPS = tx_RAW-tcorr;
+      X = satpos(tx_GPS, Eph(:,k));
+      if l ~= 1
+         rhosq = (X(1)-pos(1))^2+(X(2)-pos(2))^2+(X(3)-pos(3))^2;
+         traveltime = sqrt(rhosq)/v_light;
+         Rot_X = e_r_corr(traveltime,X);
+         rhosq = (Rot_X(1)-pos(1))^2 + ...
+            (Rot_X(2)-pos(2))^2 + ...
+            (Rot_X(3)-pos(3))^2;
+         traveltime = sqrt(rhosq)/v_light;
+         [az,el,dist] = topocent(Rot_X, Rot_X-pos(1:3,:));
+         trop = tropo(sin(el*dtr),h/1000,1013.0,293.0,50.0,...
+            0.0,0.0,0.0);
+      else trop = 0;
+      end
+      corrected_pseudorange = pr(jsat)+v_light*tcorr-trop;
+      B(jsat,1) = X(1);
+      B(jsat,2) = X(2);
+      B(jsat,3) = X(3);
+      B(jsat,4) = corrected_pseudorange;
+   end; % jsat-loop
+   pos = bancroft(B);
+   [phi,lambda,h] = togeod(6378137, 298.257223563, ...
+      pos(1), pos(2), pos(3));
+end; % l-loop
+clock_offset = pos(4)/v_light;
+pos(4) = clock_offset*1.e+9;	  % ns
+%%%%%%%%%%% b_point.m  %%%%%%%%%%%%%%%%%%%%%%%%%

baseline.m

@@ -0,0 +1,39 @@
+function [omc,bas] = baseline(X_i,obs1,obs2,sats,time,Eph)
+% BASELINE Computation of baseline between master and rover 
+%          from pseudoranges alone
+
+%Kai Borre 31-10-2001
+%Copyright (c) by Kai Borre
+%$Revision: 1.1 $  $Date: 2002/11/24  $
+
+m = size(obs1,1);  % number of svs
+% identify ephemerides columns in Eph
+for t = 1:m
+    col_Eph(t) = find_eph(Eph,sats(t),time);
+end
+% preliminary guess for receiver position 
+X_j = X_i(1:3,1);
+% Computation of weight matrix
+D = [ones(m-1,1) -eye(m-1) -ones(m-1,1) eye(m-1)];
+C = inv(D*D');
+
+for iter = 1:8
+    % k is the reference satellite. We select the first one
+    [tcorr,rhok_j,Xk_ECF] = get_rho(time, obs2(1), Eph(:,col_Eph(1)),X_j);
+    [tcorr,rhok_i,Xk_ECF] = get_rho(time, obs1(1), Eph(:,col_Eph(1)),X_i);
+    for t = 2:m % t runs over PRNs given in sats; ref.sat. is number 1
+        [tcorr,rhol_j,Xl_ECF] = get_rho(time, obs2(t), Eph(:,col_Eph(t)), X_j);
+        [tcorr,rhol_i,Xl_ECF] = get_rho(time, obs1(t), Eph(:,col_Eph(t)), X_i); 
+        A(t-1,:) = [(Xk_ECF(1)-X_j(1))/rhok_j - (Xl_ECF(1)-X_j(1))/rhol_j,  ...
+                    (Xk_ECF(2)-X_j(2))/rhok_j - (Xl_ECF(2)-X_j(2))/rhol_j,  ...
+                    (Xk_ECF(3)-X_j(3))/rhok_j - (Xl_ECF(3)-X_j(3))/rhol_j];
+        observed = obs1(1)-obs2(1)-obs1(t)+obs2(t);
+        calculated = rhok_i-rhok_j-rhol_i+rhol_j;
+        omc(t-1,1) = observed - calculated;
+    end; % t 
+    x = inv(A'*C*A)*A'*C*omc; 
+    X_j = X_j+x; 
+end % iter    
+bas = X_i(1:3,1)-X_j;
+%%%%%%%%%%%%%%%%%%%%%  baseline.m  %%%%%%%%%%%%%%%%%%%%%
+

bdata.dat

@@ -0,0 +1,13 @@
+function tt = check_t(t);
+% CHECK_T repairs over- and underflow of GPS time
+
+%  Written by Kai Borre
+%  April 1, 1996
+
+ half_week = 302400;
+ tt = t;
+
+ if t >  half_week, tt = t-2*half_week; end
+ if t < -half_week, tt = t+2*half_week; end
+
+%%%%%%% end check_t.m  %%%%%%%%%%%%%%%%%

check_t.m

@@ -0,0 +1,117 @@
+function Chi2 = chistart (D,L,a,ncands,factor)
+%CHISTART: Computes the initial size of the search ellipsoid
+%
+% This routine computes or approximates the initial size of the
+% search ellipsoid. If the requested number of candidates is not
+% more than the dimension + 1, this is done by computing the squared
+% distances of partially rounded float vectors to the float vector in
+% the metric of the covariance matrix. Otherwise an approximation is used.
+%
+% Input arguments
+%    L,D   : LtDL-decomposition of the variance-covariance matrix of
+%            the float ambiguities (preferably decorrelated)
+%    a     : float ambiguites (preferably decorrelated)
+%    ncands: Requested number of candidates (default = 2)
+%    factor: Multiplication factor for the volume of the resulting
+%            search ellipsoid (default = 1.5)
+%
+% Output arguments:
+%    Chi2  : Size of the search ellipsoid
+
+% ----------------------------------------------------------------------
+% File.....: chistart.m
+% Date.....: 19-MAY-1999
+% Author...: Peter Joosten
+%            Mathematical Geodesy and Positioning
+%            Delft University of Technology
+% ----------------------------------------------------------------------
+
+% ------------------
+% --- Initialize ---
+% ------------------
+
+if nargin < 4; ncands = 2  ; end;
+if nargin < 5; factor = 1.5; end;
+
+n = max(size(a));
+
+% ----------------------------------------------------------------------
+% --- Computation depends on the number of candidates to be computed ---
+% ----------------------------------------------------------------------
+
+if ncands == 1;
+  
+  % ---------------------------------------------------
+  % --- The squared norm, based on the bootstrapped ---
+  % --- solution will be computed                   ---
+  % ---------------------------------------------------
+
+  afloat = a;
+  afixed = a;
+
+  for i = 1:n;
+
+    dw = 0;
+    for j = 1:n-1;
+      dw = dw + L(i,j) * (a(j) - afixed(j));
+    end;
+
+    a(i) = a(i) - dw;
+    afixed(i) = round (a(i));
+   
+  end;
+
+  Chi2 = (afloat-afixed)' * inv(L'*diag(D)*L) * (afloat-afixed) + 1d-6;
+  
+elseif ncands <= n+1;
+
+  % ----------------------------------------------
+  % --- The right squared norm can be computed ---
+  % ----------------------------------------------
+   
+   Linv = inv(L);
+   Dinv = 1./D;
+   
+   dist = round(a) - a;
+   e    = Linv'*dist;
+   Chi  = [zeros(1,n) sum(Dinv' .* e .* e)];
+
+   % ------------------------------------------
+   % --- Compute the real squared distances ---
+   % ------------------------------------------
+
+   for i = 1:n;
+
+      Chi(i) = 0;
+      for j = 1:i; 
+         Chi(i) = Chi(i) + Dinv(j) * Linv(i,j) * (2*e(j)+Linv(i,j));
+      end;
+      Chi(i) = Chi(n+1) + abs(Chi(i));
+   
+   end;
+
+   % ---------------------------------------------------------------
+   % --- Sort the results, and return the appropriate number     ---
+   % --- Add an "eps", to make sure there is no boundary problem ---
+   % ---------------------------------------------------------------
+
+   Chi  = sort(Chi);
+   Chi2 = Chi(ncands) + 1d-6;
+   
+else
+
+   % -----------------------------------------------------
+   % An approximation for the squared norm is computed ---
+   % -----------------------------------------------------
+
+   Linv = inv(L);
+   Dinv = 1./D;
+   
+   Vn   = (2/n) * (pi ^ (n/2) / gamma(n/2));
+   Chi2 = factor * (ncands / sqrt((prod(1 ./ Dinv)) * Vn)) ^ (2/n);
+
+end;
+
+% ----------------------------------------------------------------------
+% End of routine: chistart
+% ----------------------------------------------------------------------

chistart.m

@@ -0,0 +1,117 @@
+function [Q,Z,L,D,z] = decorrel (Q,a);
+%DECORREL: Decorrelate a (co)variance matrix of ambiguities
+%
+% This routine creates a decorrelated Q-matrix, by finding the
+% Z-matrix and performing the corresponding transformation.
+%
+% The method is described in:
+% The routine is based on Fortran routines written by Paul de Jonge (TUD)
+% and on Matlab-routines written by Kai Borre.
+% The resulting Z-matrix can be used as follows:
+% z = Zt * a; \hat(z) = Zt * \hat(a);
+% Q_\hat(z) = Zt * Q_\hat(a) * Z
+%
+% Input arguments:
+%   a: Original ambiguities
+%   Q: Variance/covariance matrux of ambiguities (original)
+%
+% Output arguments:
+%   Q: Decorrelated variance/covariance matrix
+%   Z: Z-transformation matrix
+%   L: L matrix (from LtDL-decomposition of Q, optional)
+%   D: D matrix (from LtDL-decomposition of Q, optional)
+%   z: Transformed ambiguities (optional)
+
+% ----------------------------------------------------------------------
+% Function.: decorrel
+% Date.....: 19-MAY-1999
+% Author...: Peter Joosten
+%            Mathematical Geodesy and Positioning
+%            Delft University of Technology
+% ----------------------------------------------------------------------
+
+% -----------------------
+% --- Initialisations ---
+% -----------------------
+
+n    = size(Q,1);
+Zti  = eye(n);
+i1   = n - 1;
+sw   = 1;
+
+% --------------------------
+% --- LtDL decomposition ---
+% --------------------------
+
+[L,D] = ldldecom (Q);
+
+% ------------------------------------------
+% --- The actual decorrelation procedure ---
+% ------------------------------------------
+
+
+while sw;
+
+   i  = n;
+   sw = 0;
+
+   while ( ~sw ) & (i > 1)
+
+      i = i - 1;
+      if (i <= i1); 
+      
+         for j = i+1:n
+            mu = round(L(j,i));
+            if mu ~= 0
+               L(j:n,i) = L(j:n,i) - mu * L(j:n,j);
+               Zti(1:n,j) = Zti(1:n,j) + mu * Zti(1:n,i);
+            end
+         end
+
+      end;
+
+      delta = D(i) + L(i+1,i)^2 * D(i+1);
+      if (delta < D(i+1))
+
+         lambda(3)    = D(i+1) * L(i+1,i) / delta;
+         eta          = D(i) / delta;
+         D(i)         = eta * D(i+1);
+         D(i+1)       = delta;
+
+         Help         = L(i+1,1:i-1) - L(i+1,i) .* L(i,1:i-1);
+         L(i+1,1:i-1) = lambda(3) * L(i+1,1:i-1) + eta * L(i,1:i-1);
+         L(i,1:i-1)   = Help;
+         L(i+1,i)     = lambda(3);
+
+         Help         = L(i+2:n,i);
+         L(i+2:n,i)   = L(i+2:n,i+1);
+         L(i+2:n,i+1) = Help;
+         
+         Help         = Zti(1:n,i);
+         Zti(1:n,i)   = Zti(1:n,i+1);
+         Zti(1:n,i+1) = Help;
+
+         i1           = i;
+         sw           = 1;
+
+      end;
+
+   end;
+
+end;
+
+% ---------------------------------------------------------------------
+% --- Return the transformed Q-matrix and the transformation-matrix ---
+% --- Return the decorrelated ambiguities, if they were supplied    ---
+% ---------------------------------------------------------------------
+
+Z = inv(Zti');
+Q = Z' * Q * Z;
+
+if nargin == 2 & nargout >= 5;
+  z = Z' * a; %
+end;
+
+% ----------------------------------------------------------------------
+% End of routine: decorrel
+% ----------------------------------------------------------------------