# Objective: convex nonlinear
# Constraints: bounds
# Feasible set: convex

param pi := 4*atan(1);

param n0 := 32;
param n  := n0+1;

set S := {0..n};
set T := {0..n};

set S0 := {0..n0};
set T0 := {0..n0};

param ssize := 2;
param tsize := 2;

param hs := ssize/n;
param ht := tsize/n;

var x {s in S, t in T} := (1+s/n)*cos(t*2*pi/n);
var y {s in S, t in T} := (1+s/n)*sin(t*2*pi/n);
var z {s in S, t in T} := log(1+s/n+sqrt((1+s/n)^2-1));

minimize area: 
   sum {s in S0, t in T0}
      (
      sqrt(
           (
	     (y[s+1,t] - y[s,t] + y[s+1,t+1] - y[s,t+1])*
	     (z[s,t+1] - z[s,t] + z[s+1,t+1] - z[s+1,t])
	     -
	     (z[s+1,t] - z[s,t] + z[s+1,t+1] - z[s,t+1])*
	     (y[s,t+1] - y[s,t] + y[s+1,t+1] - y[s+1,t])
	   )^2 + 
           (
	     (z[s+1,t] - z[s,t] + z[s+1,t+1] - z[s,t+1])*
	     (x[s,t+1] - x[s,t] + x[s+1,t+1] - x[s+1,t])
	     -
	     (x[s+1,t] - x[s,t] + x[s+1,t+1] - x[s,t+1])*
	     (z[s,t+1] - z[s,t] + z[s+1,t+1] - z[s+1,t])
	   )^2 + 
	   (
	     (x[s+1,t] - x[s,t] + x[s+1,t+1] - x[s,t+1])*
	     (y[s,t+1] - y[s,t] + y[s+1,t+1] - y[s+1,t])
	     -
	     (y[s+1,t] - y[s,t] + y[s+1,t+1] - y[s,t+1])*
	     (x[s,t+1] - x[s,t] + x[s+1,t+1] - x[s+1,t])
	   )^2
	  )
      );

fix {s in S} x[s,0];
fix {t in T} x[n,t];
fix {s in S} x[s,n];
fix {t in T} x[0,t];

fix {s in S} y[s,0];
fix {t in T} y[n,t];
fix {s in S} y[s,n];
fix {t in T} y[0,t];

fix {s in S} z[s,0];
fix {t in T} z[n,t];
fix {s in S} z[s,n];
fix {t in T} z[0,t];

fix {s in S, t in T} x[s,t];
fix {s in S, t in T} y[s,t];

solve;

printf {s in S, t in T}: "%7.3f %7.3f %7.3f\n", 
	x[s,t], y[s,t], z[s,t] > height3.wrl;

printf {s in S0, t in T0}: 
  "%7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f \n", 
  x[s,t], y[s,t], z[s,t], 
  x[s+1,t], y[s+1,t], z[s+1,t],
  x[s,t+1], y[s,t+1], z[s,t+1] > "faces"; 

printf {s in S0, t in T0}: 
  "%7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f \n", 
  x[s+1,t+1], y[s+1,t+1], z[s+1,t+1], 
  x[s+1,t], y[s+1,t], z[s+1,t],
  x[s,t+1], y[s,t+1], z[s,t+1] > "faces";