c=299792.458;
p=26570;

a=[0,0,6370,.0001]';
achange=[];
s1=[p,1,1.05]';
s2=[p,1.03,1]';
s3=[p,1.02,1.02]';
s4=[p,1.05,1.01];

%print ranges
R1=sqrt((p*cos(s1(2))*cos(s1(3)))^2+(p*cos(s1(2))*sin(s1(3)))^2+(p*sin(s1(2))-6370)^2);
R2=sqrt((p*cos(s2(2))*cos(s2(3)))^2+(p*cos(s2(2))*sin(s2(3)))^2+(p*sin(s2(2))-6370)^2);
R3=sqrt((p*cos(s3(2))*cos(s3(3)))^2+(p*cos(s3(2))*sin(s3(3)))^2+((p*sin(s3(2)))-6370)^2);
R4=sqrt((p*cos(s4(2))*cos(s4(3)))^2+(p*cos(s4(2))*sin(s4(3)))^2+(p*sin(s4(2))-6370)^2);

%original vectors in rectangular coordinates
sr1=[p*cos(s1(2))*cos(s1(3)), p*cos(s1(2))*sin(s1(3)),p*sin(s1(2)), a(4)+R1/c]';
sr2=[p*cos(s2(2))*cos(s2(3)), p*cos(s2(2))*sin(s2(3)),p*sin(s2(2)), a(4)+R2/c]';
sr3=[p*cos(s3(2))*cos(s3(3)), p*cos(s3(2))*sin(s3(3)),p*sin(s3(2)), a(4)+R3/c]';
sr4=[p*cos(s4(2))*cos(s4(3)), p*cos(s4(2))*sin(s4(3)),p*sin(s4(2)), a(4)+R4/c]';

%variables we are going to change
src1=sr1;
src2=sr2;
src3=sr3;
src4=sr4;



%t and the t we want to change
t=[sr1(4),sr2(4),sr3(4),sr4(4)]';
changet=t;

%emf vector
emfmat=[];
%this loop makes all errors -1e-8 and makes one other 1e-8
for i=1:4
    for j=1:4
        if j==i
            changet(j)=t(j)+1e-8;
        else
            changet(j)=t(j)-1e-8;
        end
    end
    
    %new t values
    src1(4)=changet(1);
    src2(4)=changet(2);
    src3(4)=changet(3);
    src4(4)=changet(4);
    
    %new result and achange is difference in position
    %last result of a change is the timing delay, we made it 0 for our
    %purposes
    achange=a-newtonmgps1(a,src1,src2,src3,src4);
    achange(1)=achange(1);
    achange(2)=achange(2);
    achange(3)=achange(3);
    achange(4)=0;
    %report max error and emf
    maxposition(i)=norm(achange, inf)*1000; %report it in meters
    emf(i)=(norm(achange, inf))/(c*(1e-8));
end
%does same thing as lat loop but changes error in t differently
for i=1:4
    for j=1:4
        if j==i
            changet(j)=t(j)-1e-8;
        else
            changet(j)=t(j)+1e-8;
        end
    end
    
    %new t values
    src1(4)=changet(1);
    src2(4)=changet(2);
    src3(4)=changet(3);
    src4(4)=changet(4);
    
    %new result and achange is difference in position
    %last result of a change is the timing delay, we made it 0 for our
    %purposes
    achange=a-newtonmgps1(a,src1,src2,src3,src4);
    achange(1)=achange(1);
    achange(2)=achange(2);
    achange(3)=achange(3);
    achange(4)=0;
    maxposition(i+4)=norm(achange, inf)*1000; %report it in meters
    emf(i+4)=(norm(achange, inf))/(c*(1e-8));
end
max_position_error=max(maxposition)

condition=max(emf)