xdists = [];
emf = [];

for iter=1:500
    coef = (0.99)^iter;
    phis = rand(1,4) * pi * coef;
    thetas = rand(1,4) * 2 * pi * coef;
    
    rho = 26570;        % distance (in km) of satellites from center of earth
    sat = zeros(4, 3);  % matrix of satellite positions (known)
    for n=1:4
        sat(n,1) = rho * cos(phis(n)) * cos(thetas(n)); % sat #n’s x coord
        sat(n,2) = rho * cos(phis(n)) * sin(thetas(n)); % sat #n’s y coord
        sat(n,3) = rho * sin(phis(n));                  % sat #n’s z coord
    end
    
    %calculate max dist between satellites
    maxdist = 0;
    for a=1:3
        for b=(a+1):4
            %maxdist = max(maxdist, max(abs(sat(a) - sat(b))));
            maxdist = max(maxdist, sqrt((phis(a) - phis(b))^2 + (thetas(a) - thetas(b))^2));
        end
    end
    
    xdists(iter) = maxdist;
    emf(iter) = part4(10-8, phis, thetas); 
end

close all force
figure;
scatter(xdists, emf);
set(gca, 'xscale', 'log');
set(gca, 'yscale', 'log');
xlabel('Max angular between satellites (radians)');
ylabel('Error Magnification Factor');