function [R_vector,GDOP,PDOP,TDOP] = gpspos(SV_matrix,P_vector,R_est_vector)

n = size(SV_matrix,1);
D=zeros(4,n);
R_calc_vector = zeros(n,1);
R_vector = zeros(n,1);
L_vector = zeros(n,1);
A = zeros(4,n);


for ctr1 = 1:n
    disp(['Calculating range from receiver position estimate to SV number ',num2str(ctr1)]);
    R_calc_vector(ctr1) = ((SV_matrix(ctr1,1)-R_est_vector(1))^2 + (SV_matrix(ctr1,2)-R_est_vector(2))^2 + (SV_matrix(ctr1,3)-R_est_vector(3))^2)^.5;
	
	disp(['Calculating Observed - Predicted range for SV number ',num2str(ctr1)]);
	L_vector(ctr1) = rem(R_calc_vector(ctr1),299792.458) - P_vector(ctr1);

	disp(['Computing directional derivatives (XYZT) for SV number ',num2str(ctr1)]);
	D(ctr1,1) = (SV_matrix(ctr1,1)-R_est_vector(1))/R_calc_vector(ctr1);
	D(ctr1,2) = (SV_matrix(ctr1,2)-R_est_vector(2))/R_calc_vector(ctr1);
	D(ctr1,3) = (SV_matrix(ctr1,3)-R_est_vector(3))/R_calc_vector(ctr1);
	D(ctr1,4) = -1;

	for ctr2 = 1:4
	    A(ctr1,ctr2) =  D(ctr1,ctr2);
	end
end

disp('Finding DOPs...');
P = inv(A'*A);
GDOP = sqrt(trace(P))
PDOP = sqrt(trace(P(1:3,1:3)))
TDOP = sqrt(P(4,4))

disp('Solving for correction to receiver position estimate...');
dR = P*A'*L_vector;

disp('Applying corrections to receiver XYZ, and calculating time bias...');
for ctr1 = 1:3
    R_vector(ctr1) = R_est_vector(ctr1) + dR(ctr1);
end

R_vector(4) = dR(4);

disp('Done.');
return;