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);