param n := 10;
param pi := 4*atan(1);

option randseed '861276191';

param theta {i in 1..n} := 2*pi*Uniform01();
param phi   {i in 1..n} :=   pi*Uniform01();

var x {1..n};
var y {1..n};
var z {1..n};

minimize potential:
    sum {i in 1..n} sum {j in i+1..n} 
	1/sqrt( (x[i] - x[j])^2 + (y[i] - y[j])^2 + (z[i] - z[j])^2 );

subject to on_the_ball {i in 1..n}:
    x[i]^2 + y[i]^2 + z[i]^2 <= 1;

let {i in 1..n} x[i] := cos(theta[i])*sin(phi[i]);
let {i in 1..n} y[i] := sin(theta[i])*sin(phi[i]);
let {i in 1..n} z[i] :=               cos(phi[i]);

option solver loqo;
option loqo_options "pred_corr=0 mufactor=0.0 steplen=0.5 \
	iterlim=400 verbose=2 sigfig=12 inftol=1.0e-6 bndpush=1000";

#option solver minos;   ## fails!!!

#fix x[1];
#fix y[1];
#fix z[1];

display potential;

solve;

display x, y, z;
#display {i in 1..n, j in i+1..n} x[i]*x[j] + y[i]*y[j] + z[i]*z[j];
display max {i in 1..n} (x[i]^2 + y[i]^2 + z[i]^2);
display min {i in 1..n} (x[i]^2 + y[i]^2 + z[i]^2);

printf: '#VRML V2.0 utf8 \
NavigationInfo { type "EXAMINE" } \n' > fekete.wrl;

printf {i in 1..n}:
"Transform { \
     translation %7.4f %7.4f %7.4f \
     children [ \
         Shape { \
             appearance Appearance { \
                 material Material { \
		     diffuseColor  1 0.8 0 \
		     specularColor 1 1 1 \
		     shininess     0.7 \
		 } \
             } \
             geometry Sphere { radius 0.04 } \
         } \
     ] \
 } \n", 
x[i], y[i], z[i] > fekete.wrl;

printf:
"Shape { \
     appearance Appearance { \
	 material Material { \
	     diffuseColor  0.5 0 0.5 \
#	     specularColor 1 1 1 \
#	     shininess     0.7 \
	 } \
     } \
     geometry Sphere { radius 1 } \
 } \n"
 > fekete.wrl;