part4

0001 clear;
0002 warning('off','all');
0003 warning;
0004 
0005 nsats = 4;
0006 
0007 rho = 26570*1000; % distance from origin in meters
0008 
0009 niter = 10;
0010 
0011 x0 = [ 0 0 6370 0 ]';
0012 
0013 % backward error (i.e. input error)
0014 delta_t = [ 10^-6 10^-7 10^-8 10^-9 10^-10 10^-11 10^-12 ];
0015 delta_t = fliplr(delta_t);
0016 
0017 emf = zeros(niter,1);
0018 err = zeros(niter,1);
0019 maxEMF = zeros( length(delta_t), 1 );
0020 minEMF = zeros( length(delta_t), 1 );
0021 maxERR = zeros( length(delta_t), 1 );
0022 
0023 distEMF = zeros( length(delta_t), niter );
0024 distERR = zeros( length(delta_t), niter );
0025 
0026 for x = 1:niter,
0027     for i = 1:length(delta_t), % for a given satellite clock error
0028         for j = 1:niter,
0029             % generate a set of random clock positions
0030             sat = genRandomSatellites( nsats, rho );
0031             % and solve the sys. of eqns., noting the error mag. factor
0032             % and the relative position error
0033             [emf(j), err(j)] = pertubateClockError( sat, x0, delta_t(i), @solveNavEqns1 );
0034             % RINSE, REPEAT ...
0035         end % END for j
0036         % our condition number of the system is the max EMF
0037         maxEMF(i) = max(emf);
0038         % for posterity
0039         minEMF(i) = min(emf);
0040         % and record the max position error
0041         maxERR(i) = max(err);
0042         % REPEAT for the next clock error
0043     end % END for i
0044     % REPEAT above randomization several times and record the trend
0045     distEMF(1:length(maxEMF),x) = maxEMF;
0046     distERR(1:length(maxERR),x) = maxERR;
0047 end % END for x
0048 
0049 meanEMF = zeros( length(delta_t), 1 );
0050 meanERR = zeros( length(delta_t), 1 );
0051 % take the avg of the randomized trend
0052 for i = 1:length(delta_t),
0053     meanEMF(i) = mean(distEMF(i,1:niter));
0054     meanERR(i) = mean(distERR(i,1:niter));
0055 end
0056 
0057 figure(1);
0058 semilogx(delta_t, meanEMF,'b-',delta_t, meanEMF,'ro');
0059 ylabel('cond(EMF)');
0060 xlabel('clock error');
0061 title('cond vs clock error');
0062 
0063 figure(2);
0064 semilogx(delta_t, meanERR,'b-',delta_t, meanERR,'ro');
0065 ylabel('position error (m)');
0066 xlabel('clock error');
0067 title('position error vs clock error');
0068 
0069 figure(3);
0070 semilogx(delta_t(1:3), meanERR(1:3),'b-', delta_t(1:3), meanERR(1:3),'ro');
0071 set(gca,'XTick', delta_t(1:3));
0072 ylabel('position error (m)');
0073 xlabel('clock error');
0074 title('position error vs clock error');