c=299792.458;
p=26570;

a=[0,0,6370,.0001]';
achange=[];
s1=[p,0,pi]';
s2=[p,pi/4,pi/2]';
s3=[p,pi/6,5*pi/3]';
s4=[p,pi/3,2*pi];
s5=[p,pi/6,pi]';
s6=[p,pi/5,20*pi/12]';
s7=[p,pi/10,0]';
s8=[p,pi/2,pi/3];

%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);
R5=sqrt((p*cos(s5(2))*cos(s5(3)))^2+(p*cos(s5(2))*sin(s5(3)))^2+(p*sin(s5(2))-6370)^2);
R6=sqrt((p*cos(s6(2))*cos(s6(3)))^2+(p*cos(s6(2))*sin(s6(3)))^2+(p*sin(s6(2))-6370)^2);
R7=sqrt((p*cos(s7(2))*cos(s7(3)))^2+(p*cos(s7(2))*sin(s7(3)))^2+(p*sin(s7(2))-6370)^2);
R8=sqrt((p*cos(s8(2))*cos(s8(3)))^2+(p*cos(s8(2))*sin(s8(3)))^2+(p*sin(s8(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]';
sr5=[p*cos(s5(2))*cos(s5(3)), p*cos(s5(2))*sin(s5(3)),p*sin(s5(2)), a(4)+R5/c]';
sr6=[p*cos(s6(2))*cos(s6(3)), p*cos(s6(2))*sin(s6(3)),p*sin(s6(2)), a(4)+R6/c]';
sr7=[p*cos(s7(2))*cos(s7(3)), p*cos(s7(2))*sin(s7(3)),p*sin(s7(2)), a(4)+R7/c]';
sr8=[p*cos(s8(2))*cos(s8(3)), p*cos(s8(2))*sin(s8(3)),p*sin(s8(2)), a(4)+R8/c]';

%variables we are going to change
src1=sr1;
src2=sr2;
src3=sr3;
src4=sr4;
src5=sr5;
src6=sr6;
src7=sr7;
src8=sr8;



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

emfmat=[];
for i=1:8
    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);
    src5(4)=changet(5);
    src6(4)=changet(6);
    src7(4)=changet(7);
    src8(4)=changet(8);
    
    %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-newtonmgps6(a,src1,src2,src3,src4,src5,src6,src7,src8);
    achange(1)=achange(1);
    achange(2)=achange(2);
    achange(3)=achange(3);
    achange(4)=0;
    maxposition(i)=norm(achange, inf)*1000; %report it in meters
    emf(i)=(norm(achange, inf))/(c*(1e-8));
end

for i=1:8
    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);
    src5(4)=changet(5);
    src6(4)=changet(6);
    src7(4)=changet(7);
    src8(4)=changet(8);
    
    %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-newtonmgps6(a,src1,src2,src3,src4,src5,src6,src7,src8);
    achange(1)=achange(1);
    achange(2)=achange(2);
    achange(3)=achange(3);
    achange(4)=0;
    maxposition(i+8)=norm(achange, inf)*1000; %report it in meters
    emf(i+8)=(norm(achange, inf))/(c*(1e-8));
end
maxposition'
max_position_error=max(maxposition)
emf'
condition=max(emf)