param N := 5;  # number of masses
param n := 32;  # number of terms in Fourier series representation
param m := 60;  # number of terms in numerical approx to integral
		# must be a multiple of 3?

param pi := 4*atan(1);
param theta {i in 0..N-1, j in 0..m-1} := j*2*pi/(m*N) + i*2*pi/N;

param tau := 2*pi;  # period of the orbit
param omega := 2*pi/tau;

var as {ko in 1..2*n-1 by 2} := if (ko = 3) then 1 else 0;
var ac {ke in 2..2*n by 2} := if (ke = 2) then 0 else 0;

var bs {ke in 2..2*n by 2} := if (ke = 2) then 1 else 0;
var bc {ko in 1..2*n-1 by 2} := if (ko = 1) then 0 else 0;

var x {i in 0..N-1, j in 0..m-1} 
	= sum {k in 1..n} (
		as[2*k-1]*sin((2*k-1)*theta[i,j]) +
		ac[2*k]*cos((2*k)*theta[i,j]) 
		);

var y {i in 0..N-1, j in 0..m-1} 
	= sum {k in 1..n} (
		bs[2*k]*sin((2*k)*theta[i,j]) +
		bc[2*k-1]*cos((2*k-1)*theta[i,j]) 
		);

var xdot {i in 0..N-1, j in 0..m-1} 
	= sum {k in 1..n} (
		as[2*k-1]*(2*k-1)*cos((2*k-1)*theta[i,j]) -
		ac[2*k]*(2*k)*sin((2*k)*theta[i,j]) 
		);

var ydot {i in 0..N-1, j in 0..m-1} 
	= sum {k in 1..n} (
		bs[2*k]*(2*k)*cos((2*k)*theta[i,j]) -
		bc[2*k-1]*(2*k-1)*sin((2*k-1)*theta[i,j]) 
		);

var K {j in 0..m-1} 
	= (omega^2 /2)*sum {i in 0..N-1} (xdot[i,j]^2 + ydot[i,j]^2);

var P {j in 0..m-1}
	= - sum {i in 0..N-1, ii in 0..N-1: ii>i} 
	    1/sqrt((x[i,j]-x[ii,j])^2 + (y[i,j]-y[ii,j])^2);
	
#check {j in 0..m-1, i in 0..N-1, ii in 0..N-1: ii>i}: 
#    sqrt((x[i,j]-x[ii,j])^2 + (y[i,j]-y[ii,j])^2) > 0;

minimize A: (2*pi/m)*sum {j in 0..m-1} (K[j] - P[j]);

printf {i in 0..N-1, j in 0..m-1}: 
	"%10.5f %10.5f \n", x[i,j], y[i,j] > "before.out";
printf  "%10.5f %10.5f \n", x[0,0], y[0,0] > "before.out";

solve;

printf {i in 0..N-1, j in 0..m-1}: 
	"%10.5f %10.5f \n", x[i,j], y[i,j] > "after.out";
printf  "%10.5f %10.5f \n", x[0,0], y[0,0] > "after.out";

printf {i in 0..N-1, j in 0..m-1}: "pathx[%d]=%10.5f; pathy[%d]=%10.5f; \n", 
    j+i*m, x[i,j], j+i*m, y[i,j] > "pathxy.out";
printf  "pathx[%d]=%10.5f; pathy[%d]=%10.5f; \n", 
    N*m, x[0,0], N*m, y[0,0] > "pathxy.out";

display as, ac, bs, bc;

printf {i in 0..N-1}: "x[%d]=%10f; y[%d]=%10f; \n",  
	i, x[i,0], i, y[i,0];

printf {i in 0..N-1}: "vx[%d]=%10f; vy[%d]=%10f; \n",  
	i, xdot[i,0], i, ydot[i,0];

printf "atan2(ydot/xdot), sqrt(xdot^2+ydot^2) \n";
printf {i in 0..N-1}: "%10f %10f \n",  
	(180/pi)*atan2(xdot[i,0],ydot[i,0]), sqrt(xdot[i,0]^2+ydot[i,0]^2);