function A2

%clear all variables
clear;

global data;

%initialize 
[data,P,r]=ginit(10);

fprintf('Anfangspunkt (%.4f,%.4f) mit Radius %.4f\n',P(1),P(2),r);

%%% now main loop
%tolerance for Newton iteration
tol=10^(-16);
%set d to enter main loop
d=tol+1;
%iteration counter
iter=1;
itmax=50;
%store successive errors
clear err;
while d>tol & iter<itmax
  x=[P';r];
  xn=x-dfi(P,r)*f(P,r);
  Pn=[xn(1),xn(2)];
  rn=xn(3);
  d=sqrt(sum((Pn-P).^2)+(rn-r)^2);
  err(iter)=d;
  P=Pn; r=rn;
%plot a green circle
  circle(P,r,sprintf('%d',iter),'g');
  drawnow;
  iter=iter+1;
end;
%% end of main loop
%release figure
hold off;
figure;
plot(log10(err),'r-*');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% returns function value
function [q]=f(P,r)

global data;

q=(P(1)-data(:,1)).^2+(P(2)-data(:,2)).^2-r^2;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%returns inverse of Jacobian
function [M]=dfi(P,r)

global data;

M=[2*(P(1)-data(:,1)),2*(P(2)-data(:,2)),-2*r*ones(3,1)]^(-1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%draws a circle
function circle(P,r,s,col)
%usage circle(P,r,s,col) 
%plot a circle around P (row vector) with radius r
%the colorstring has the same syntax as in the plot command 
%the string s is put on the circle

global PHIVAR;

%initialize global variable
if isempty(PHIVAR) PHIVAR=0; end;

%number of points for plotting
N=100;
%angles for plotting
theta=linspace(0,2*pi,N)';
Q=repmat(P,N,1)+r*[cos(theta),sin(theta)];
plot(Q(:,1),Q(:,2),col);
%get random angle
x=P(1)+cos(PHIVAR)*r;
y=P(2)+sin(PHIVAR)*r;
text(x,y,s);
%increase writing angle
PHIVAR=PHIVAR+23/180*pi;
%adjust axes if necessary
Lc=max(max(abs(Q)));
L=max(abs(axis));
if Lc>L
  axis([-Lc,Lc,-Lc,Lc]);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%gets initial values

function [data,P,r]=ginit(L)
% usage: [data,P,r]=ginit(L)
% sets up square [-L,L]x[L,L]
% and returns three point coordinates in the rows of data
% as well as a point in P and a radius r

%clear all variables
clear data;
%set up a new figure and keep all plots
figure;
hold on;
%put axis and make square plot
axis([-L,L,-L,L]);
axis square;
%get three points
title('Geben Sie drei Punkte ein');
for i=1:3 
  P=ginput(1);
  plot(P(1),P(2),'b*');
  data(i,:)=P;
end;
%now coordinates are stored in the rows of variable data
% next: get initial data
title('Markieren Sie den Startmittelpunkt');
P=ginput(1);
plot(P(1),P(2),'r*');
title('Markieren Sie den Startradius'); 
Q=ginput(1);
r=sqrt(sum((P-Q).^2));
%plot as red circle
circle(P,r,'','r');
%and clear title
title('');
