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];

%calculate 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 rectangular vectors
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;

emfmat=[];
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)=norm(achange, inf)*1000; %report it in meters
    emf(i)=(norm(achange, inf))/(c*(1e-8));
end
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)