function E = energy(R,B,phi,theta,alpha)
N=length(phi);
E1 = 0; E2 = 0;
for i = 1:N-1
   for j = i+1:N
      rr = R(:,i)-R(:,j);
      d  = dot(rr,rr);
      E1 = E1 - 1/d;
   end
end
for i = 2:3
   E2 = E2 + cos(3*phi(i));
end
   E  = E1 + alpha*E2;
