function [Q, iter] = newton_raphson(x, y, z, Q)

    iter = 0;
    while(1)
        a = Q(1);
        b = Q(2);
        c = Q(3);
        r = Q(4);
    
        % Jacobian
        J = [ -2*(x(1)-a)   -2*(y(1)-b)   -2*(z(1)-c)   -2*r;
              -2*(x(2)-a)   -2*(y(2)-b)   -2*(z(2)-c)   -2*r;
              -2*(x(3)-a)   -2*(y(3)-b)   -2*(z(3)-c)   -2*r;
              -2*(x(4)-a)   -2*(y(4)-b)   -2*(z(4)-c)   -2*r ];

        % function value
        f = [ (x(1)-a)^2 + (y(1)-b)^2 + (z(1)-c)^2 - r^2;
              (x(2)-a)^2 + (y(2)-b)^2 + (z(2)-c)^2 - r^2;
              (x(3)-a)^2 + (y(3)-b)^2 + (z(3)-c)^2 - r^2;
              (x(4)-a)^2 + (y(4)-b)^2 + (z(4)-c)^2 - r^2];

        % calculate {dq}
        dq = J \ (-f);
        
        %update
        Q = Q + dq';
        iter = iter + 1;

        % check convergence
        if(max(abs(dq)) < 1.e-6), break, end 
        if(iter > 20), break, end
    end
end