function ex29();
% find the circle through 3 given points
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
clear all;
close all;

global x y

% 3 given points;
x = [0 2 2]';
y = [0 0 1]';


TOL = 0.001;
differ = 10*TOL;

xc0 = 0;
yc0 = 0;
r0 = 1;

while differ > TOL;
  UPD = [xc0;yc0;r0] - ...
      inv(DF(xc0,yc0,r0))*F(xc0,yc0,r0);
  
  xc1 = UPD(1);
  yc1 = UPD(2);
  r1  = UPD(3);
  
  differ = max(abs([xc0-xc1; yc0-yc1; r0-r1]));
  
  xc0 = xc1;
  yc0 = yc1;
  r0 = r1;
end

r = r0
xc = xc0
yc = yc0

% plot the result:
t = linspace(0,2*pi,100);

X = r*cos(t)+xc;
Y = r*sin(t)+yc;

figure('name','Circle through 3 given points')
hold on

plot(x,y,'ro')
plot(X,Y,'b')
xlabel('x')
ylabel('y')
grid on
axis equal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Y = F(xc,yc,r);

global x y

Y = (x-xc).^2 + (y-yc).^2 - r^2;


function Y = DF(xc,yc,r);
% Jacobian
  
global x y

Y = [-2*(x-xc),-2*(y-yc),-2*r*ones(3,1)];