 clear all, close all
 
 p=0;
 q=0.85;
 r=0.5; 

x(1)=rand() 
y(1)=rand() 
z(1)=rand() 

for j=1:600000 
    rr=rand();
    if rr>0 & rr<=0.25
    x(j+1) = +x(j)/(x(j)^2+y(j)^2+z(j)^2) + p;
    y(j+1) = -y(j)/(x(j)^2+y(j)^2+z(j)^2) + q;
    z(j+1) = +z(j)/(x(j)^2+y(j)^2+z(j)^2) + r;
    end
    if rr>0.25 & rr<=0.5
    w=((x(j)-p)^2+(y(j)-q)^2+(z(j)-r)^2);
    x(j+1) = (x(j)-p)/w; 
    y(j+1) = (-y(j)+q)/w; 
    z(j+1) = (z(j)-r)/w; 
    end 
    if rr>0.5 & rr<=0.75
    x(j+1)=x(j)+2;
    y(j+1)=y(j);
    z(j+1)=z(j);
    end 
    if rr>0.75 & rr<=1
    x(j+1)=x(j)-2;
    y(j+1)=y(j);
    z(j+1)=z(j);
    end 
end
    
 plot3(x,y,z,'.','MarkerSize',4)
 axis([-2 2 -2 2 -2 2])

 

 rr=sqrt(q^2+r^2);
 x=x;
 y=y+2*r/rr; 
 z=z-2*q/rr; 
 
Xp=4*rr.*x./(x.^2+y.^2+z.^2);
Yp=4*rr.*(y./(x.^2+y.^2+z.^2)-r/(4*rr));
Zp=4*rr.*(z./(x.^2+y.^2+z.^2)+q/(4*rr));

figure
 plot3(Xp,Yp,Zp,'.','MarkerSize',4) , hold on
 axis equal
