let NB=300;;
let billes=make_matrix NB 3 0.;;
let dPI=6.2831853071795864770;;

let fmod_dPI x=
	let i=int_of_float(x/.dPI) in
	x -. (dPI*.float_of_int(i));;

let VecteurCarthesien i phi theta =
	  billes.(i).(0)<-cos(phi)*.sin(theta);
	  billes.(i).(1)<-sin(theta)*.sin(phi);
	  billes.(i).(2)<-cos(theta);;

let distance i j =
	let t = ref 0. in begin
		for k = 0 to 2 do
			let h = billes.(i).(k) -. billes.(j).(k) in
			t:= !t +. (h *. h)
		done;
		sqrt(!t)
	end;;

let EP () =
	let ep=ref 0. in begin
		for i = 0 to NB-2 do
			for j = i+1 to NB-1 do
				ep:= !ep +. (1. /. (distance i j))
			done
		done;
		!ep
	end;;

let PlacerSpirale hw =
	let K = hw /. sqrt(float_of_int(NB)) and phi = ref 0. in
	for i = 0 to NB-1 do
		let h = float_of_int(2*i) /. float_of_int(NB-1) -. 1. in
		let theta=acos h in begin
			if i=0 or i=NB-1 then
				phi:= 0.
			else
				phi:= fmod_dPI (!phi +. K/.sqrt(1.-.h*.h));
			VecteurCarthesien i !phi theta
		end
	done;
	EP ();;

PlacerSpirale 3.6;;