Added surfaces of constant negative curvature (surface.cpp)

compileunits.h

@@ -84,6 +84,7 @@ bool inv::activating;
 #if CAP_TOUR
 #include "tour.cpp"
 #endif
+#include "surface.cpp"
 #if ISMOBILE==0
 #include <unistd.h>
 #endif

hyper.h

@@ -2870,5 +2870,16 @@ ld hypot2(const hyperpoint& h);
 transmatrix pushxto0(const hyperpoint& H);
 transmatrix rpushxto0(const hyperpoint& H);
 transmatrix spintox(const hyperpoint& H);
+transmatrix ypush(ld alpha);
 
 extern int ewhichscreen;
+
+#if CAP_SURFACE
+namespace surface {
+
+  enum eShape { dsNone, dsTractricoid, dsDini, dsKuen, dsHyperlike, dsHyperboloid, dsHemisphere };
+  extern eShape sh;
+  void show_surfaces();
+  
+  }
+#endif

init.cpp

@@ -125,6 +125,10 @@
 #define CAP_RUG (!ISMINI && CAP_GL)
 #endif
 
+#ifndef CAP_SURFACE
+#define CAP_SURFACE CAP_RUG
+#endif
+
 #ifndef CAP_EDIT
 #define CAP_EDIT (CAP_FILES && !ISWEB && !ISMINI)
 #endif

rug.cpp

@@ -1702,6 +1702,14 @@ void show() {
   dialog::addSelItem(XLAT("automatic move speed"), fts(ruggo), 'G');
   dialog::addSelItem(XLAT("anti-crossing"), fts(anticusp_factor), 'A');
 
+#if CAP_SURFACE  
+  if(hyperbolic) {
+    if(gwhere == gEuclid)
+      dialog::addItem(XLAT("smooth surfaces"), 'c');
+    else dialog::addBreak(100);
+    }
+#endif
+
   dialog::display();
   keyhandler = [] (int sym, int uni) {
     dialog::handleNavigation(sym, uni);
@@ -1709,7 +1717,12 @@ void show() {
     if(uni == 'h') gotoHelp(makehelp());
     else if(uni == 'u') {
       if(rug::rugged) rug::close();
-      else rug::init();
+      else { 
+        #if CAP_SURFACE
+        surface::sh = surface::dsNone; 
+        #endif
+        rug::init(); 
+        }
       }
     else if(uni == 'R')
       dialog::editNumber(finger_range, 0, 1, .01, .1, XLAT("finger range"),
@@ -1719,7 +1732,7 @@ void show() {
       dialog::editNumber(finger_force, 0, 1, .01, .1, XLAT("finger force"),
         XLAT("Press 1 to enable the finger force.")
         );
-    else if(uni == 'o' && !rug::rugged)
+    else if(uni == 'o')
       renderonce = !renderonce;
     else if(uni == 'G') {
       dialog::editNumber(ruggo, -1, 1, .1, 0, XLAT("automatic move speed"),
@@ -1793,6 +1806,10 @@ void show() {
       texturesize *= 2;
       if(texturesize == 8192) texturesize = 64;
       }
+#if CAP_SURFACE
+    else if(uni == 'c') 
+      pushScreen(surface::show_surfaces);
+#endif
     else if(handlekeys(sym, uni)) ;
     else if(doexiton(sym, uni)) popScreen();
     };

surface.cpp

@@ -0,0 +1,826 @@
+// Hyperbolic Rogue
+// Copyright (C) 2011-2018 Zeno Rogue, see 'hyper.cpp' for details
+// This file implements the surfaces of constant negative curvature
+// See http://webmath2.unito.it/paginepersonali/sergio.console/CurveSuperfici/AG15.pdf for a nice reference
+
+#if CAP_SURFACE
+namespace surface {
+
+using namespace hyperpoint_vec;
+
+ld sech(ld d) { return 1 / cosh(d); }
+
+string shape_name[] = { "hypersian rug", "tractricoid", "Dini's surface", "Kuen surface", "concave barrel",
+  "hyperboloid", "hemisphere" };
+  
+eShape sh;
+
+hyperpoint unit_vector[3] = {hpxyz(1,0,0), hpxyz(0,1,0), hpxyz(0,0,1)};
+
+ld last_int_of = 0, last_int = 0;
+
+ld dini_b = .15;
+
+ld hyper_b = 1;
+
+ld f(ld x) {
+  return sqrt(1 - pow(hyper_b * sinh(x), 2));
+  }
+
+int kuen_branch(ld v, ld u);
+
+ld integral(ld x) {
+  if(x == 0) {
+    last_int = last_int_of = 0;
+    }
+  else {
+    last_int += (x - last_int) * f((x + last_int)/2);
+    last_int_of = x;
+    }
+  return last_int;
+  }
+
+transmatrix build_matrix(hyperpoint h1, hyperpoint h2, hyperpoint h3) {
+  transmatrix T;
+  for(int i=0; i<3; i++)
+    T[i][0] = h1[i],
+    T[i][1] = h2[i],
+    T[i][2] = h3[i];
+  return T;
+  }
+
+hyperpoint coord(hyperpoint h) {
+  // return { cos(u)*sin(v), cos(u)*cos(v), sin(u) };
+  
+  /*
+  ld t = h[0];
+  ld v = h[1];
+  ld r = 1 / cosh(t);
+  ld x = t - tanh(t);
+  return { x + v * .1, r * sin(v), r * cos(v) };
+  */
+
+  switch(sh) {
+    case dsTractricoid: {
+      ld t = h[0];
+      ld v = h[1];
+      ld r = 1 / cosh(t);
+      ld x = t - tanh(t);
+
+      return { r * sin(v), r * cos(v), x };
+      break;
+      }
+    
+    case dsDini: {
+      ld t = h[0]; // atan(h[0])/2 + M_PI * 3/ 4;
+      ld v = h[1];
+      
+      ld a = sqrt(1-dini_b*dini_b);
+      
+      return { a * sin(v) * sin(t), a * cos(v) * sin(t),  a * (cos(t) + log(tan(t/2))) + dini_b * v };
+      break;
+      }
+    
+    case dsKuen: {
+      ld v = h[0];
+      ld u = h[1];
+      
+      ld deno = 1 / (1 + u * u * sin(v) * sin(v));
+      
+      return {
+        2 * (cos(u) + u * sin(u)) * sin(v) * deno,
+        2 * (sin(u) - u * cos(u)) * sin(v) * deno,
+        log(tan(v/2)) + 2 * cos(v) * deno
+        };
+      }
+    
+    case dsHyperlike: {
+      ld u = h[0];
+      ld v = h[1];
+
+      ld phi = hyper_b * cosh(v);
+      ld psi = integral(v);
+      
+      return { phi * cos(u), phi * sin(u), psi};
+      }
+    
+    default:
+      return h;
+    }
+  }
+
+ld det(hyperpoint h1, hyperpoint h2, hyperpoint h3) {
+  return det(build_matrix(h1, h2, h3));
+  }
+
+ld epsd = 1e-5;
+
+hyperpoint coord_derivative(hyperpoint h, int cc) {
+  switch(sh) {
+    case dsHyperlike: {
+      ld u = h[0];
+      ld v = h[1];
+      if(cc == 0) {
+        ld phi = hyper_b * cosh(v);
+        return { phi * -sin(u), phi * cos(u), 0 };
+        }
+      else {
+        return { hyper_b * sinh(v) * cos(u), hyper_b * sinh(v) * sin(u), f(v) };
+        }
+      }
+    case dsKuen: {
+      ld v = h[0];
+      ld u = h[1];      
+      ld denom = pow(sin(v),2)*(u*u)+1;
+      ld denom2 = denom * denom;
+      if(cc == 1) 
+        return {
+          2*sin(v)/denom*u*cos(u)+-4*(sin(u)*u+cos(u))*pow(sin(v),3)/denom2*u,
+          -4*pow(sin(v),3)*(sin(u)-u*cos(u))/denom2*u+2*sin(u)*sin(v)/denom*u,
+          -4*pow(sin(v),2)/denom2*u*cos(v)
+          };
+      else return {
+          2*(sin(u)*u+cos(u))/denom*cos(v)+-4*(sin(u)*u+cos(u))*pow(sin(v),2)/denom2*(u*u)*cos(v),
+          2*(sin(u)-u*cos(u))/denom*cos(v)+-4*pow(sin(v),2)*(sin(u)-u*cos(u))/denom2*(u*u)*cos(v),
+         -4*sin(v)/denom2*(u*u)*pow(cos(v),2)+1/tan(v/2)*(pow(tan(v/2),2)+1)/2+-2*sin(v)/denom
+         };
+      break;
+      }
+    default:
+      return (coord(h + unit_vector[cc] * epsd) - coord(h)) / epsd;
+    }
+  }
+
+ld compute_curvature(hyperpoint at) {
+  hyperpoint xu = coord_derivative(at, 0);
+  hyperpoint xv = coord_derivative(at, 1);
+  hyperpoint xuu = (coord_derivative(at + unit_vector[0] * epsd, 0) - xu) / epsd;
+  hyperpoint xuv = (coord_derivative(at + unit_vector[1] * epsd, 0) - xu) / epsd;
+  hyperpoint xvv = (coord_derivative(at + unit_vector[1] * epsd, 1) - xv) / epsd;
+  return 
+    (det(xuu, xu, xv) * det(xvv, xu, xv) - pow(det(xuv, xu, xv), 2)) /
+    pow((xu|xu) * (xv|xv) - pow((xu|xv), 2), 2);
+  }
+
+hyperpoint shape_origin() {
+  switch(sh) {
+    case dsDini:
+      return hpxyz(M_PI * .82, 0, 0);
+    case dsTractricoid:
+      return hpxyz(1, 0, 0);
+    case dsKuen:
+      return hpxyz(M_PI * .500001, M_PI * 1, 0);
+    case dsHyperlike:
+      return hpxyz(0,0,0);
+    default:
+      return Hypc;
+    }
+  }
+
+ld hyperlike_bound() { return asinh(1 / hyper_b); }
+
+bool flag_clamp_min(ld& coord, ld minv) {
+  if(coord < minv) { coord = minv; return false; }
+  return true;
+  }
+
+bool flag_clamp_max(ld& coord, ld maxv) {
+  if(coord > maxv) { coord = maxv; return false; }
+  return true;
+  }
+
+bool flag_clamp(ld& coord, ld minv, ld maxv) {
+  return flag_clamp_min(coord, minv) & flag_clamp_max(coord, maxv);
+  };
+
+bool flag_clamp_sym(ld& coord, ld v) {
+  return flag_clamp(coord, -v, v);
+  }
+
+int surface_branch(hyperpoint p) {
+  if(sh == dsKuen) return kuen_branch(p[0], p[1]);
+  return 0;
+  }
+
+bool inbound(ld& x, ld& y) {
+  switch(sh) {
+    case dsDini:
+      return flag_clamp(x, M_PI/2, M_PI);
+    
+    case dsTractricoid:
+      return flag_clamp_min(x, 0) & flag_clamp_sym(y, M_PI);
+    
+    case dsKuen: 
+      return flag_clamp(x, 0, M_PI) & flag_clamp(y, 0, 2*M_PI);
+    
+    case dsHyperlike:
+      return flag_clamp_sym(x, M_PI) & flag_clamp_sym(y, hyperlike_bound());
+    
+    default: 
+      return true;
+    }
+  }
+
+bool is_inbound(hyperpoint h) {
+  return inbound(h[0], h[1]);
+  }
+
+int precision = 100;
+
+using rug::dexp_data;
+
+struct dexp_origin {
+  transmatrix H;     // isometry of H2 moving zero to C0
+  transmatrix M;     // local coordinates on H2 to local coordinates on surface
+  hyperpoint zero;   // parameters of the zero point
+  };
+
+dexp_data dexp(hyperpoint p, hyperpoint t) {
+  ld eps = 1. / precision;
+  int b = surface_branch(p);
+  
+  for(ld u=0; u<1; u += eps) {
+    
+    transmatrix T = build_matrix(coord_derivative(p, 0), coord_derivative(p, 1), Hypc);
+    
+    // printf("Tt = %lf\n", hypot3(T * t));
+
+    p += t * eps;
+
+    if(!is_inbound(p) || surface_branch(p) != b)
+      return { p - t * eps, t, hypot3(t) * (1-u) / precision };
+      
+    auto v0 = coord_derivative(p, 0);
+    auto v1 = coord_derivative(p, 1);
+
+    transmatrix T2 = build_matrix(v0, v1, v0 ^ v1);
+      
+    t = inverse(T2) * T * t;
+    t[2] = 0;
+    }
+  return { p, t, 0 };
+  }
+
+dexp_data map_to_surface(hyperpoint p, const dexp_origin& dor) {
+  hyperpoint h = dor.H * p;
+  ld rad = hypot2(h);
+  if(rad == 0) rad = 1;
+  ld d = hdist0(h);
+
+  hyperpoint direction;
+  direction[0] = d * h[0] / rad;
+  direction[1] = d * h[1] / rad;
+  direction[2] = 0;
+
+  return dexp(dor.zero, dor.M * direction);
+  }
+
+transmatrix create_M_matrix(hyperpoint zero, hyperpoint v1) {
+  hyperpoint Te0 = coord_derivative(zero, 0);
+  hyperpoint Te1 = coord_derivative(zero, 1);
+
+  transmatrix T = build_matrix(Te0, Te1, Hypc);
+
+  v1 = v1 / hypot3(T*v1);
+  hyperpoint v2 = hpxyz(1e-3, 1e-4, 0);
+  v2 = v2 - v1 * ((T*v1) | (T*v2)) / hypot3(T*v1);
+  v2 = v2 / hypot3(T*v2);
+  
+  if((((T*v1) ^ (T*v2)) | ((T*unit_vector[0]) ^ (T*unit_vector[1]))) < 0)
+    v2 = v2 * -1;
+  
+  transmatrix M = build_matrix(v1, v2, Hypc);
+  
+  display(M);
+
+  printf("M matrix test: %lf %lf %lf %lf\n", hypot3(T*M*unit_vector[0]), hypot3(T*M*unit_vector[1]), hypot3(T*M*(unit_vector[0]+unit_vector[1])),
+    ((T*M*unit_vector[0]) | (T*M*unit_vector[1])));
+    
+  return M;
+  }
+
+dexp_origin at_zero(hyperpoint zero, transmatrix start) {
+
+  printf("zero = %s\n", display(zero));
+
+  printf("curvature at zero = %lf\n", compute_curvature(zero));
+  printf("curvature at X1 = %lf\n", compute_curvature(zero + hpxyz(.3, 0, 0)));
+  printf("curvature at X2 = %lf\n", compute_curvature(zero + hpxyz(0, .3, 0)));
+  printf("curvature at X3 = %lf\n", compute_curvature(zero + hpxyz(.4, .3, 0)));
+  
+  return {start, create_M_matrix(zero, unit_vector[0]), zero};
+  }
+
+dexp_origin at_other(dexp_origin& o1, hyperpoint h) {
+  printf("\n\nmapping %s...\n", display(h));
+  display(o1.H); display(o1.M);
+  auto dd = map_to_surface(h, o1);
+  
+  hyperpoint newzero = dd.params;
+  printf("error = %lf\n", dd.remaining_distance);
+
+  transmatrix Spin = spintox(o1.H * h);
+  transmatrix T = pushxto0(Spin * o1.H * h) * Spin;
+  
+  printf("h is = %s\n", display(h));
+  printf("T*c0 is = %s\n", display(T * C0));
+  printf("T*h is = %s\n", display(T * o1.H * h));
+  
+  return {T * o1.H, create_M_matrix(newzero, dd.cont), newzero};
+  }
+
+void addTriangleV(rug::rugpoint *t1, rug::rugpoint *t2, rug::rugpoint *t3, ld len = 1) {
+  if(t1 && t2 && t3)
+    rug::addTriangle(t1, t2, t3, len);
+  }
+
+hyperpoint kuen_cross(ld v, ld u) {
+  auto du = coord_derivative(hpxyz(v,u,0), 0);
+  auto dv = coord_derivative(hpxyz(v,u,0), 1);
+  return du^dv;
+  }
+
+ld kuen_hypot(ld v, ld u) {
+  auto du = coord_derivative(hpxyz(v,u,0), 0);
+  auto dv = coord_derivative(hpxyz(v,u,0), 1);
+  auto n = hypot3(du^dv);
+  return n;
+  }
+
+int kuen_branch(ld v, ld u) {
+  if(v > M_PI/2)
+  return kuen_cross(v, u)[2] > 0 ? 1 : 2;
+  else
+  return kuen_cross(v, u)[2] < 0 ? 1 : 2;
+  }
+
+int dexp_colors[16] = { 
+  0xFF0000, 0x00FF00, 0x0000FF, 0xFFFF00 };
+
+int dexp_comb_colors[16] = { 
+  0x000000, 0x0000FF, 0x00FF00, 0x00FFFF,
+  0xFF0000, 0xFF00FF, 0xFFFF00, 0xFFFFFF,
+  0xFFD500, 
+  0x123456, 0x123456, 0x123456, 0x123456, 0x123456, 0x123456, 0x123456
+   };
+
+int coverage_style;
+vector<pair<hyperpoint, int> > coverage;
+
+void draw_kuen_map() {
+  SDL_Surface *kuen_map = SDL_CreateRGBSurface(SDL_SWSURFACE,512,512,32,0,0,0,0);
+  
+  ld nmax = 0;
+
+  for(int i=0; i<2; i++) {
+    for(int r=0; r<512; r++)
+    for(int h=0; h<512; h++) {
+      ld v = M_PI * (r+.5) / 512;
+      ld u = 2 * M_PI * (h+.5) / 512;
+      auto du = coord_derivative(hpxyz(v,u,0), 0);
+      auto dv = coord_derivative(hpxyz(v,u,0), 1);
+      auto n = hypot3(du^dv);
+
+      if(n > nmax) nmax = n;
+      if(i == 1) {
+        auto vv = kuen_cross(v, u);
+        auto& px = qpixel(kuen_map, r, h);
+        px |= 0xFF000000;
+        for(int k=0; k<3; k++)
+          part(px, k) = (vv[k] > 0 ? 0xFF : 0);
+        px = 0xFF000000 + (((int)(n*255/nmax)) * (kuen_branch(v,u) == 1 ? 0x10101 : 0x10001));
+        }
+      }
+    printf("nmax = %lf\n", nmax);
+    }
+
+  for(auto p: rug::points) {
+    auto hp = p->surface_point.params;
+    int x = int(512 * hp[0] / M_PI);
+    int y = int(512 * hp[1] / 2 / M_PI);
+    qpixel(kuen_map, x, y) = 0xFF000000 | dexp_colors[p->dexp_id];
+    }
+
+  IMAGESAVE(kuen_map, "kuen.png");    
+  }
+
+void full_mesh() {
+  rug::clear_model();
+  rug::buildRug();
+  rug::qvalid = 0;
+  for(auto p: rug::points) p->valid = true, rug::qvalid++;
+  while(rug::subdivide_further()) rug::subdivide();
+  rug::sort_rug_points();
+  for(auto p: rug::points) p->valid = false;
+  rug::good_shape = true;
+  }
+
+char rchar(int id) {
+  return 33 + id % 94;
+  }
+  
+void run_hyperlike() {
+  if(!rug::rugged) rug::reopen();
+  rug::clear_model();
+  
+  int lim = (int) sqrt(rug::vertex_limit);
+  for(int r=0; r<lim; r++)
+  for(int h=0; h<lim; h++) 
+    rug::addRugpoint(xpush(2 * M_PI * hyper_b * (2*r-lim) / lim) * ypush(hyperlike_bound() * (2*h-lim) / lim) * C0, -1);
+  for(int r=0; r<lim-1; r++)
+  for(int h=0; h<lim-1; h++) {
+    addTriangle(rug::points[lim*r+h], rug::points[lim*r+h+1], rug::points[lim*r+h+lim]);
+    addTriangle(rug::points[lim*r+h+lim+1], rug::points[lim*r+h+lim], rug::points[lim*r+h+1]);
+    }
+
+  vector<ld> integral_table;
+  for(int i=0; i<=precision; i++)
+    integral_table.push_back(integral(hyperlike_bound() * i / precision));
+
+  int id = 0;
+  for(auto p: rug::points) {
+    auto h = p->h;
+    coverage.emplace_back(h, rchar(id++) + 7 * 256);
+    
+    ld y = asinh(h[1]);
+    ld x = asinh(h[0] / cosh(y)) / hyper_b;
+    p->surface_point.remaining_distance = !inbound(x, y);
+    p->surface_point.params = hpxyz(x, y, 0);
+
+    int sgn = y > 0 ? 1 : -1;
+    ld phi = hyper_b * cosh(y);
+    int pt = y * precision * sgn / hyperlike_bound();
+    p->flat = hpxyz(phi * cos(x), phi * sin(x), sgn * integral_table[pt]);
+    p->valid = true;
+    }
+  }
+
+void run_kuen() {
+  full_mesh();
+
+  auto H = Id; // spin(-M_PI / 4) * xpush(2);
+  auto Hi = inverse(H);
+  
+  auto frontal_map = at_zero(hpxyz(M_PI * .500001, M_PI * 1, 0), Id);
+  auto back0 = at_zero(hpxyz(M_PI * .500001, .67, 0), H);
+  auto back1 = at_other(back0, Hi * spin(-M_PI/2) * hpxy(0.511, -0.5323));
+  auto back2 = at_other(back0, Hi * spin(-M_PI/2) * hpxy(0.511, 0.5323));
+  
+  frontal_map.H = frontal_map.H * ypush(2.6);
+  back0.H = back0.H * ypush(.4);
+  back1.H = back1.H * ypush(.4);
+  back2.H = back2.H * ypush(.4);
+  
+  int it = 0;
+  for(auto p: rug::points) p->dexp_id = it++;
+  
+  vector<rug::rugpoint*> mesh = move(rug::points);
+  vector<rug::triangle> old_triangles = move(rug::triangles);
+  
+  rug::clear_model();
+  
+  int part = 0;
+  
+  vector<int> coverages(size(mesh), 0);
+  
+  for(auto& m: { frontal_map, back0, back1, back2 } ) {
+    part++;
+    int pid[5] = {0, 8, 1, 2, 4};
+    string captions[5] = {"", "the upper component", "the lower center", "the lower left", "the lower right"};
+    
+    vector<rug::rugpoint*> newmesh(size(mesh), NULL);
+    for(auto p: mesh) {
+      auto px = map_to_surface(p->h, m);
+      p->surface_point = px;
+      conformal::progress(XLAT("solving the geodesics on: %1, %2/%3", XLAT(captions[part]), its(p->dexp_id), its(size(mesh))));
+      }
+    for(auto p: mesh) {
+      // make it a bit nicer by including the edges where only one endpoint is valid
+
+      auto& px = p->surface_point;
+      p->valid = px.remaining_distance == 0;
+      for(auto e: p->edges) if(e.target->surface_point.remaining_distance == 0)
+        p->valid = true;
+      if(p->valid) {
+        rug::rugpoint *np = new rug::rugpoint;
+        newmesh[p->dexp_id] = np;
+        rug::points.push_back(np);
+        np->x1 = p->x1;
+        np->y1 = p->y1;
+        np->valid = true;
+        np->inqueue = false;
+        np->dist = 0;
+        np->h = p->h;
+        np->flat = coord(px.params);
+        np->surface_point = px;
+        np->dexp_id = p->dexp_id;
+        coverages[p->dexp_id] |= pid[part];
+        }
+      }
+    for(auto& t: old_triangles) {
+      rug::rugpoint* r[3];
+      for(int i=0; i<3; i++)
+        r[i] = newmesh[t.m[i]->dexp_id];
+      bool looks_good = true;
+      for(int i=0; i<3; i++)
+        if(!r[i]) looks_good = false;
+      if(!looks_good) continue;
+      for(int i=0; i<3; i++)
+        if(hypot3(r[i]->flat - r[(i+1)%3]->flat) > .2)
+          looks_good = false;
+      if(looks_good)
+        addTriangleV(r[0], r[1], r[2]);
+      }
+    }
+  
+  for(auto t: mesh) {
+    int c = coverages[t->dexp_id];
+    coverage.emplace_back(t->h, rchar(t->dexp_id) + 256 * c);
+    }
+
+  // delete the old mesh
+  for(auto t: mesh) delete t;
+
+  draw_kuen_map();
+  }
+
+template<class T> void run_function(T f) {
+  full_mesh();
+  for(auto p: rug::points)
+    p->flat = f(p->h),
+    p->valid = true;
+  }
+
+void run_other() {
+  full_mesh();
+  auto dp = at_zero(shape_origin(),  spin(M_PI/2));
+
+  int it = 0;
+  for(auto p: rug::points) {
+    it++;
+    auto h = p->h;
+
+    p->surface_point = map_to_surface(h, dp);
+    p->flat = coord(p->surface_point.params);
+    conformal::progress(XLAT("solving the geodesics on: %1, %2/%3", shape_name[sh], its(it), its(size(rug::points))));
+    if(p->surface_point.remaining_distance == 0)
+      coverage.emplace_back(h, rchar(it) + 256 * 7);
+    }
+  
+  clearMessages();
+
+  for(auto p: rug::points) {
+    // make it a bit nicer by including the edges where only one endpoint is valid
+    p->valid = p->surface_point.remaining_distance == 0;
+    if(sh != dsKuen) {
+      for(auto e: p->edges) if(e.target->surface_point.remaining_distance == 0)
+        p->valid = true;
+      }
+    }
+  }
+
+void run_shape(eShape s) {
+  coverage.clear();
+  need_mouseh = true;
+  sh = s;
+  transmatrix crot = rug::rugged ? rug::currentrot : Id;
+  rug::apply_rotation(inverse(crot));
+  
+  if(rug::rugged) rug::close();
+
+  rug::init();
+  // if(!rug::rugged) rug::reopen();
+  
+  pushScreen(conformal::progress_screen);
+  if(sh != dsNone) rug::good_shape = true;
+
+  switch(sh) {
+    case dsNone:
+      break;
+      
+    case dsKuen:
+      run_kuen();
+      break;
+    
+    case dsHyperlike:
+      run_hyperlike();
+      break;
+    
+    default:
+      run_other();
+      break;
+
+    case dsHyperboloid:
+      run_function([] (hyperpoint h) { return h; });
+      break;
+
+    case dsHemisphere:
+      run_function([] (hyperpoint h) { h = h / h[2]; h[2] = sqrt(1 - sqhypot2(h)); return h; });
+      break;      
+    }
+
+  if(sh != dsNone) rug::good_shape = true;
+  
+  rug::qvalid = 0;
+  
+  popScreen();
+  
+  if(sh != dsHyperboloid && sh != dsHemisphere && sh != dsNone) {
+
+    ld minx = 1e9, maxx = -1e9;
+  
+    for(auto p: rug::points) if(p->valid) {
+      minx = min(p->flat[2], minx);
+      maxx = max(p->flat[2], maxx);
+      rug::qvalid++;
+      }
+    
+    printf("minx = %lf maxx = %lf\n", minx, maxx);
+  
+    ld shift = -(minx + maxx) / 2;
+    for(auto p: rug::points) if(p->valid)
+      p->flat[2] += shift;
+    }
+
+  rug::apply_rotation(crot);
+  if(rug::rug_perspective)
+    rug::push_all_points(2, -rug::model_distance);
+  }
+
+void cancel_shape() {
+  if(sh) {
+    sh = dsNone;
+    rug::good_shape = false;
+    rug::qvalid = 0;
+    for(auto p: rug::points) 
+      p->valid = p->surface_point.remaining_distance == 0;
+    for(auto p: rug::points) 
+      if(p->valid) 
+        rug::qvalid++, rug::enqueue(p);
+    }
+  }
+
+cell *coverage_center;
+transmatrix coverage_matrix;
+
+void show_surfaces() {
+  cmode = sm::SIDE;
+  gamescreen(0);
+  dialog::init(XLAT("constant curvature surfaces"), iinf[itPalace].color, 150, 0);
+  
+  bool b = rug::rugged || coverage_style;
+
+  dialog::addBoolItem(XLAT("tractricoid"), b && sh == dsTractricoid, '1');
+  dialog::addBoolItem(XLAT("concave barrel"), b && sh == dsHyperlike, '2');
+  dialog::addSelItem(XLAT("     parameter"), fts(hyper_b), '@');
+  dialog::addBoolItem(XLAT("Dini's surface"), b && sh == dsDini, '3');
+  dialog::addSelItem(XLAT("     parameter") + " ", fts(dini_b), '#');
+  dialog::addBoolItem(XLAT("Kuen surface"), b && sh == dsKuen, '4');
+  
+  dialog::addBreak(50);
+  
+  dialog::addTitle(XLAT("other 3D models"), iinf[itPalace].color, 150);
+
+  dialog::addBoolItem(XLAT("Hypersian Rug"), b && sh == dsNone, '5');
+  dialog::addBoolItem(XLAT("Minkowski hyperboloid"), b && sh == dsHyperboloid, '6');
+  dialog::addBoolItem(XLAT("hemisphere model"), b && sh == dsHemisphere, '7');
+
+  dialog::addBreak(100);
+
+  dialog::addItem(XLAT("what's this?"), 'h');
+  dialog::addSelItem(XLAT("precision"), its(precision), 'p');
+  string cstyles[4] = { "OFF", "on surface", "on H² (static)", "on H² (dynamic)" };
+  if((rug::rugged && sh && sh != dsHyperboloid && sh != dsHemisphere) || coverage_style)
+    dialog::addSelItem(XLAT("display coverage"), cstyles[coverage_style], 'c');
+  else dialog::addBreak(100);
+  dialog::addItem(XLAT("take me back"), 'q');
+  
+  dialog::display();
+  keyhandler = [] (int sym, int uni) {
+    dialog::handleNavigation(sym, uni);
+
+    if(uni == 'h') 
+      gotoHelp(
+        "In this menu you can choose from several known smooth surfaces of constant negative curvature. "
+        "Since the curvature of hyperbolic plane is also constant negative, these surfaces "
+        "are great to draw hyperbolic tesselations on. While they look great, only a small part "
+        "of the hyperbolic plane can be visibly represented in this way, so they are not "
+        "good for playing HyperRogue; however, the coverage extends far away in specific directions, "
+        "so first increasing sight range in graphicss config and texture size in hypersian rug settings may improve the effect. "
+        "For convenience, you can also choose other 3D models from this menu."
+        );
+
+    else if(uni == '1')
+      run_shape(dsTractricoid);
+    else if(uni == '2')
+      run_shape(dsHyperlike);
+    else if(uni == '3')
+      run_shape(dsDini);
+    else if(uni == '4')
+      run_shape(dsKuen);
+    else if(uni == '5')
+      run_shape(dsNone);
+    else if(uni == '6')
+      run_shape(dsHyperboloid);
+    else if(uni == '7')
+      run_shape(dsHemisphere);
+
+    else if(uni == '@') {
+      dialog::editNumber(hyper_b, -1, 1, .05, 1, XLAT("parameter"),
+        XLAT("Controls the inner radius.")
+        );
+      dialog::reaction = [] () {
+        if(sh == dsHyperlike) run_shape(sh);
+        };
+      }
+    else if(uni == 'x')
+      for(auto p: rug::points)
+        p->flat = p->surface_point.params;
+    else if(uni == '#')
+      dialog::editNumber(dini_b, -1, 1, .05, .15, XLAT("parameter"),
+        XLAT("The larger the number, the less twisted it is.")
+        );
+    else if(uni == 'p') {
+      dialog::editNumber(precision, 1, 10000, .1, 100, XLAT("precision"),
+        XLAT("Computing these models involves integrals and differential equations, which are currently solved numerically. This controls the precision.")
+        );
+      dialog::scaleLog();
+      }
+    else if(uni == 'c') {
+      coverage_style = (1 + coverage_style) % 4;
+      if(coverage_style == 0) {
+        rug::reopen();
+        }
+      if(coverage_style == 2) {
+        if(rug::rugged) rug::close();
+        }
+      coverage_matrix = inverse(shmup::ggmatrix(coverage_center = cwt.c));
+      }
+    else if(rug::handlekeys(sym, uni)) ;
+    else if(doexiton(sym, uni)) popScreen();
+    };
+  }
+
+#if CAP_COMMANDLINE
+int surface_args() {
+  using namespace arg;
+           
+  if(0) ;
+  else if(argis("-kuen")) {
+    calcparam();
+    run_shape(dsKuen);
+    }
+
+  else if(argis("-dini")) {
+    calcparam();
+    shift();
+    dini_b = argf();
+    run_shape(dsDini);
+    }
+
+  else if(argis("-barrel")) {
+    calcparam();
+    shift();
+    hyper_b = argf();
+    run_shape(dsHyperlike);
+    }
+
+  else if(argis("-tractricoid")) {
+    calcparam();
+    run_shape(dsTractricoid);
+    }
+
+  else if(argis("-hemi")) {
+    calcparam();
+    run_shape(dsHemisphere);
+    }
+
+  else if(argis("-hyperb")) {
+    calcparam();
+    run_shape(dsHyperboloid);
+    }
+
+  else return 1;
+  return 0;
+  }
+
+auto surface_hook = addHook(hooks_args, 100, surface_args);
+#endif
+
+void display_coverage() {
+
+  transmatrix M = 
+    coverage_style == 3 ? shmup::ggmatrix(coverage_center) * coverage_matrix
+    : Id;
+
+  if(coverage_style)
+    for(auto p : coverage)
+      queuechr(M * p.first, 10, char(p.second), dexp_comb_colors[(p.second >> 8) & 15]);
+/*      if(p->valid && p->surface_point.remaining_distance == 0)
+        queuechr(p->h, 10, 'x', dexp_colors[p->dexp_id]); */
+  }
+
+auto surface_hook2 = addHook(hooks_frame, 0, display_coverage);
+
+}
+#endif