clc; clear;
xc=input('Xc= ');yc=input('yc= ');
r=input('r= ');
hold on;
i=0; 
trace=zeros(1,8);
axis([-100 100 -100 100])
x=0; y=r;
while x<=y
 plot(xc+x,yc+y,'g.');
 plot(xc+x,yc-y,'b.');
 plot(xc-x,yc+y,'k.');
 plot(xc-x,yc-y,'r.');
 plot(xc+y,yc+x,'b.');
 plot(xc+y,yc-x,'k.');
 plot(xc-y,yc+x,'g.');
 plot(xc-y,yc-x,'b.');
  da=(x+1)^2+y^2-r^2;
 db=(x+1)^2+(y-1)^2-r^2;
  i=i+1;
  disp('     x    y    da    x+1     y    db x+1    y-1'    )
  trace(i,:)=[x y da x+1 y db x+1 y-1]
 if abs(db)<abs(da)
    y=y-1;
 end   
 x=x+1;
end