clear;clf;
pas=0.01;

// Cercle de base
// sur un quart du plan (x>=0, y>=0)

Y=[0:pas:1]';
NY=size(Y); ny=NY(1,1);

for i=1:ny
   X(i)=sqrt(1-Y(i)^2);
end

// Cercles

for j=1:8
   r1=cos(j*%pi/18); // décalage de 10°
   r2=sin(j*%pi/18);
   for i=1:ny
      X1(i,j)=X(i)*r1; // cooredonnées des poitns du cercle incliné
      Z(i,j)=X(i)*r2;
   end
end

// Projection stéréographique

// Pôle Sud S(0,0,-1), pôle P(X1,Y,Z)
// Droite : S + m*SP = (m*X1)i + (m*Y)j + (-1+m*(Z+1))k

for j=1:8
   for i=1:ny
      m = 1/(1+Z(i,j)); // z=0
      X2(i,j)=m*X1(i,j);
      Y2(i,j)=m*Y(i);
   end
end

// Tracé

angle=['10','20','30','40','50','60','70','80'];

for j=1:8
   plot2d(X2(:,j),Y2(:,j),style=2)
   plot2d(-X2(:,j),Y2(:,j),style=2)
   plot2d(X2(:,j),-Y2(:,j),style=2)
   plot2d(-X2(:,j),-Y2(:,j),style=2)
   xstring(X2(1,j)+0.01,0.01,angle(j))
end

plot2d(X,Y)
plot2d(-X,Y)
plot2d(X,-Y)
plot2d(-X,-Y)

xpoly([0,0],[-1,1])
xpoly([-1,1],[0,0])

axe=get('current_axes');
axe.axes_visible=['off','off'];
axe.isoview='on';
axe.box='off';