model "Polygon"
 uses "mmxslp"

parameters
 N=5
end-parameters
 
declarations
 area: gexp
 rho : array(1..N) of mpvar
 theta : array(1..N) of mpvar
 objdef: mpvar
 D: array(1..N,1..N) of genctr
 rTheta: array (1..50) of xvitem
 BIGARRAY: array(1..100) of real
end-declarations

! Objective - sum of areas
SLPDATA("UF","MoselArea","","DOUBLE,INTEGER","MOSEL","Polygon","BIGARRAY")
SLPDATA("XV",rTheta)

forall (i in 1..N-1) do
	rTheta(i*2-1) := XVitem(rho(i))
	rTheta(i*2) := XVitem(theta(i))
end-do

area := Func("MoselArea",rTheta)
objdef = area
objdef is_free

! Initialisation of SLP variables, plus bounds
forall (i in 1..N-1) do
 	rho(i) >= 0.1
 	rho(i) <= 1
 	SLPDATA("IV",rho(i),4*i*(N + 1 - i)/((N+1)^2))
 	SLPDATA("IV",theta(i),M_PI*i/N)
end-do
 
! Third side of all triangles <= 1
forall (i in 1..N-2) do
 	forall (j in i+1..N-1) do
 		D(i,j) := rho(i)^2 + rho(j)^2 - rho(i)*rho(j)*2*cos(theta(j)-theta(i)) <= 1
 	end-do
end-do
 
! Vertices in increasing order 
forall (i in 2..N-1) do
	theta(i) >= theta(i-1) +.01
end-do
 
! Boundary conditions
theta(N-1) <= M_PI

SLPloadprob(objdef)
SLPmaximise

writeln("Area = ",getobjval)
forall (i in 1..N-1) do
 	writeln("V",i,": r=",getsol(rho(i))," theta=",getsol(theta(i)))
end-do

function MoselArea (A:array(RA:range) of real, N:integer) : real
 declarations
 	n: integer
 	r1,r2: real		! distances
 	t1,t2: real		! angles
 	area: real
 end-declarations
 
 n := 4
 area := 0
 while (n <= N) do
 	r1 := A(n-3)
 	r2 := A(n-1)
 	t1 := A(n-2)
 	t2 := A(n)
 	n := n+2
 	area := area + 0.5*r1*r2*sin(t2-t1)
 end-do
 returned := area
end-function

end-model