pertubateClockError

 Tests the conditioning of the GPS problem for a given satellite clock
 error
 INPUT: sat - an array of structs representing the satellite positions
 x0 - initial guess of the receiver position
 delta_t - the satellite clock error
 solveNavEqns - function handle to a solver for the sys. of eqns.
 OUTPUT: emf - error magnification factor; fwderr - the forward error of
 the position of the receiver caused by the satellite clock error

0001 % Tests the conditioning of the GPS problem for a given satellite clock
0002 % error
0003 % INPUT: sat - an array of structs representing the satellite positions
0004 % x0 - initial guess of the receiver position
0005 % delta_t - the satellite clock error
0006 % solveNavEqns - function handle to a solver for the sys. of eqns.
0007 % OUTPUT: emf - error magnification factor; fwderr - the forward error of
0008 % the position of the receiver caused by the satellite clock error
0009 function [emf, fwderr] = pertubateClockError( sat, x0, delta_t, solveNavEqns )
0010 c = 299792458; % speed of light in m/s
0011 d = 0.0001; % error of Earth receiver clock
0012 
0013 % Euclidean distances of satellites to X,Y,Z position on Earth
0014 R = zeros( length(sat), 1 );
0015 % R1 = sqrt( sat(1).x^2 + sat(1).y^2 + (sat(1).z - x0(3))^2 );
0016 for i = 1:length(R),
0017     R(i) = sqrt( sat(i).x^2 + sat(i).y^2 + (sat(i).z - x0(3))^2 );
0018 end
0019 
0020 % calculate travel time of signal from satellite to receiver
0021 for i = 1:length(sat),
0022     sat(i).t = d + R(i)/c;
0023 end
0024 
0025 [ x1 y1 z1 d1 ] = solveNavEqns( sat, x0, c );
0026 
0027 % find the ratio of the incorrect answer if we change our clock error by
0028 % a small amount
0029 % update the measured times from each satellite with the clock error
0030 deltaT = zeros( length(sat), 1 );
0031 for i = 1:length(sat),
0032     sgn = 1;
0033     if( mod(i,2) ~= 0 ) sgn = -1;
0034     sat(i).t = sat(i).t + sgn*delta_t;
0035     deltaT(i) = sgn*delta_t;
0036 end
0037 
0038 [ x2 y2 z2 d2 ] = solveNavEqns( sat, x0, c );
0039 
0040 deltaX = (x1 - x2);
0041 deltaY = (y1 - y2);
0042 deltaZ = (z1 - z2);
0043 
0044 infNormXYZ = norm( [ deltaX deltaY deltaZ ], Inf );
0045 infNormT = c * norm( deltaT, Inf );
0046 emf = infNormXYZ / infNormT;
0047 fwderr = infNormXYZ;
0048 end
0049