$title maximum area of an N-gon 

option limrow=0, limcol=0;

set vertices/1*20/;
alias (vertices, i, j);

scalar N "number of vertices (max 20)" /8/;

variable rho(vertices) "polar radius"
	 theta(vertices) "polar angle"
	 area;

scalar pi;
pi = 4.0 * arctan(1.0);

rho.lo(vertices) = 0.0;
rho.up(vertices) = 1.0;
theta.lo(vertices) = 0.0;
theta.up(vertices) = pi;

equations
	inscribe(i,j)       "ensure that each vertex pair no more than 1 apart"
	monotonic(vertices) "enforce monotonic increase in polar angles"
	objective           "area of polygon";

inscribe(i,j)$(ord(i) < ord(j) and ord(j) <= N)..
        rho(i)**2 + rho(j)**2 - 2*rho(i)*rho(j)*cos(theta(i) - theta(j)) =l= 1;
monotonic(vertices)$(ord(vertices) < N)..
        theta(vertices) =l= theta(vertices+1);
objective..
	area =e= 
	  .5*sum(i$(ord(i) <N), rho(i)*rho(i+1)*sin(theta(i+1)-theta(i)));


* set initial values
loop(i$(ord(i)<N),
   rho.l(i)   = 4*ord(i)*(N+1-ord(i))/(N+1)**2;
   theta.l(i) = pi*ord(i)/N;
);

* fix the final vertex position at the origin
rho.fx(i)$(ord(i) = N)   = 0.0;
theta.fx(i)$(ord(i) = N) = pi;

model ngon/inscribe,monotonic,objective/;

option nlp=minos;
solve ngon using nlp maximizing area;

option N:0; option area:8; 
display N, area.l;
display rho.l, theta.l;

* convert to Cartesian coordinates
set components/x,y/;
parameter cartesian(vertices,components);
loop(i$(ord(i) <= N),
  cartesian(i,'x') = rho.l(i) * cos(theta.l(i));
  cartesian(i,'y') = rho.l(i) * sin(theta.l(i));
);
option cartesian:8; display cartesian;

* write to an output file for plotting
file plotit /ngon_coords.dat/;
put plotit

loop(i$(ord(i) <= N),
  put cartesian(i,'x'):14:8, cartesian(i,'y'):14:8;
  put /;
);
putclose plotit;