// BlackStar.sce
//
// f(x,y,z) = x^2 + y^2 + z^2- 1 = 0 
N = 2000 
theta =  2*%pi*rand(N,1) 
phi   = -%pi/2 + %pi*rand(N,1)
x = cos(theta).*cos(phi) 
y = sin(theta).*cos(phi) 
z =             sin(phi) 
//scatter3d(x, y, z) 
//mtlb_axis("equal")
//
function fout = f1(x,y,z)
    fout = 1 - x.^2 - y.^2 - z.^2;
endfunction

function fout = f2(x,y,z)
    fout = 1 - (x-0.9).^2 - (y-0.9).^2 - (z-0.9).^2;
endfunction

function fout = f(x,y,z)
   fout = min(f1(x,y,z), -f2(x,y,z))
endfunction

// Tirage aleatoire de points (xi0,yi0,zi0)

function Fvec = F(X, xi0,yi0,zi0)
    x = X(1); 
    y = X(2); 
    z = X(3);
    F1 = f(x,y,z);
    F2 = xi0*y - yi0*x;
    F3 = xi0*z - zi0*x;
    Fvec = [F1;F2;F3];
endfunction

Xsamples = [];
Ni = 100
Nj = 50;
for i = 1:Ni
    for j=1:Nj
    //theta =  2*%pi*rand(1,1) 
    //phi   = -%pi/2 + %pi*rand(1,1)
      theta = 2*%pi*(i-1.0)/Ni;
      phi = -%pi/2 + %pi*(j-1.0)/Nj;
      xi0 = 0.5*cos(theta).*cos(phi); 
      yi0 = 0.5*sin(theta).*cos(phi); 
      zi0 = 0.5*            sin(phi);
    //xi0 = rand(1,1,"normal"); 
    //yi0 = rand(1,1,"normal");
   // zi0 = rand(1,1,"normal");
   // norme = norm([xi0,yi0,zi0]);
   // xi0 = xi0 / norme
   // yi0 = yi0 / norme
   // zi0 = zi0 / norme
    X0 = [xi0; yi0; zi0];
    Xout = fsolve(X0, list(F, xi0,yi0,zi0));
    Xsamples = [Xsamples, Xout];
end;end;
clf;
scatter3d(Xsamples(1,:), Xsamples(2,:), Xsamples(3,:))
mtlb_axis("equal")