From 4f27b12ca2f6c2fb5b36cfb7f8c55faffbc86cdd Mon Sep 17 00:00:00 2001
From: Zeno Rogue <zeno@attnam.com>
Date: Tue, 6 Aug 2019 12:00:46 +0200
Subject: [PATCH] nilv:: preliminary implementation

---
 basegraph.cpp    |   8 +-
 bigstuff.cpp     |   2 +-
 cell.cpp         |   3 +-
 classes.cpp      |   1 +
 classes.h        |   4 +-
 config.cpp       |   2 +-
 conformal.cpp    |   4 +-
 control.cpp      |   2 +-
 drawing.cpp      |  14 ++
 geom-exp.cpp     |   8 +-
 geometry.cpp     |   2 +-
 graph.cpp        |  25 +--
 hyper.h          |  31 ++-
 hyperpoint.cpp   |  44 ++--
 hypgraph.cpp     |  21 +-
 nonisotropic.cpp | 521 +++++++++++++++++++++++++++++++++--------------
 pattern2.cpp     |   2 +-
 polygons.cpp     |   6 +-
 rug.cpp          |   2 +-
 19 files changed, 476 insertions(+), 226 deletions(-)

diff --git a/basegraph.cpp b/basegraph.cpp
index f9512844..ff617b0a 100644
--- a/basegraph.cpp
+++ b/basegraph.cpp
@@ -233,7 +233,7 @@ void display_data::set_projection(int ed) {
       shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::halfplane3;
     if(DIM == 3 && hyperbolic && apply_models && pmodel == mdPerspective)
       shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardH3, pers3 = true;
-    if(DIM == 3 && (euclid || sol) && apply_models && pmodel == mdPerspective)
+    if(DIM == 3 && (euclid || sol || nil) && apply_models && pmodel == mdPerspective)
       shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardR3, pers3 = true;
     if(DIM == 3 && apply_models && pmodel == mdSolPerspective)
       shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardSolv, pers3 = true;
@@ -334,8 +334,8 @@ void display_data::set_projection(int ed) {
     if(pers3) {
       glhr::projection_multiply(glhr::frustum(current_display->tanfov, current_display->tanfov * cd->ysize / cd->xsize));
       glhr::projection_multiply(glhr::scale(1, -1, -1));
-      if(solv::local_perspective_used())
-        glhr::projection_multiply(glhr::tmtogl_transpose(solv::local_perspective));
+      if(nisot::local_perspective_used())
+        glhr::projection_multiply(glhr::tmtogl_transpose(nisot::local_perspective));
       }
     else if(DIM == 3) {
       glhr::glmatrix M = glhr::ortho(cd->xsize/current_display->radius/2, -cd->ysize/current_display->radius/2, 1);
@@ -344,6 +344,8 @@ void display_data::set_projection(int ed) {
       M[2][2] = 2 / (clip_max - clip_min);
       M[3][2] = (clip_min + clip_max) / (clip_max - clip_min);
       glhr::projection_multiply(M);
+      if(nisot::local_perspective_used())
+        glhr::projection_multiply(glhr::tmtogl_transpose(nisot::local_perspective));
       }
     else {
       glhr::projection_multiply(glhr::frustum(cd->xsize / cd->ysize, 1));
diff --git a/bigstuff.cpp b/bigstuff.cpp
index 4baaef7b..1a09c7c2 100644
--- a/bigstuff.cpp
+++ b/bigstuff.cpp
@@ -1240,7 +1240,7 @@ bool good_for_wall(cell *c) {
   }
   
 void buildBigStuff(cell *c, cell *from) {
-  if(sphere || quotient || sol || (penrose && !binarytiling)) return;
+  if(sphere || quotient || nonisotropic || (penrose && !binarytiling)) return;
   if(chaosmode > 1) return;
   bool deepOcean = deep_ocean_at(c, from);
   
diff --git a/cell.cpp b/cell.cpp
index 89fc9f7d..797fd12b 100644
--- a/cell.cpp
+++ b/cell.cpp
@@ -220,6 +220,7 @@ void initcells() {
   
   hrmap* res = callhandlers((hrmap*)nullptr, hooks_newmap);
   if(res) currentmap = res;  
+  else if(nonisotropic) currentmap = nisot::new_map();
   #if CAP_CRYSTAL
   else if(geometry == gCrystal) currentmap = crystal::new_map();
   #endif
@@ -241,7 +242,6 @@ void initcells() {
   else if(sphere) currentmap = new hrmap_spherical;
   else if(quotient) currentmap = new quotientspace::hrmap_quotient;
   #if CAP_BT
-  else if(sol) currentmap = solv::new_map();
   else if(binarytiling) currentmap = binary::new_map();
   #endif
   else currentmap = new hrmap_hyperbolic;
@@ -407,6 +407,7 @@ int compdist(int dx[]) {
 int celldist(cell *c) {
   if(fulltorus && WDIM == 2) 
     return torusmap()->dists[decodeId(c->master)];
+  if(nil) return DISTANCE_UNKNOWN;
   if(euwrap)
     return torusconfig::cyldist(decodeId(c->master), 0);
   if(masterless)
diff --git a/classes.cpp b/classes.cpp
index 9b664f20..4a9c510a 100644
--- a/classes.cpp
+++ b/classes.cpp
@@ -563,6 +563,7 @@ vector<geometryinfo> ginf = {
   {"sol",    "none",    "Solv geometry",                              "sol",      8, 3, qBINARY,   gcSol,        0x41600, {{7, 5}}, eVariation::pure},
   {"kd2",    "none",    "kite-and-dart",                              "kd2",      4, 3, qPENROSE,  gcEuclid,     0x48000, {{7, 7}}, eVariation::pure},
   {"kd3",    "none",    "kite-and-dart on horospheres",               "kd3",     12, 3, qsBP,      gcHyperbolic, 0x48200, {{7, 3}}, eVariation::pure},
+  {"nil",    "none",    "Nil geometry",                               "nil",     22, 3, 0,         gcNil,        0x48600, {{7, 5}}, eVariation::pure},
 //{"product","none",    "product space",                              "product",  7, 3, 0,         gcProduct,    0x48400, {{7, 3}}, eVariation::pure},
   };
 
diff --git a/classes.h b/classes.h
index 17ae50eb..5b11eba1 100644
--- a/classes.h
+++ b/classes.h
@@ -201,10 +201,10 @@ enum eGeometry {
   gHoroTris, gHoroRec, gHoroHex,
   gField435, gField534,
   gBinary4, gSol,
-  gKiteDart2, gKiteDart3, PROD2(gProduct,)
+  gKiteDart2, gKiteDart3, gNil, PROD2(gProduct,)
   gGUARD};
 
-enum eGeometryClass { gcHyperbolic, gcEuclid, gcSphere, gcSol, PROD(gcProduct) };
+enum eGeometryClass { gcHyperbolic, gcEuclid, gcSphere, gcSol, gcNil, PROD(gcProduct) };
 
 enum class eVariation { bitruncated, pure, goldberg, irregular, dual };  
 
diff --git a/config.cpp b/config.cpp
index 2b65eb05..759a2282 100644
--- a/config.cpp
+++ b/config.cpp
@@ -465,7 +465,7 @@ void initConfig() {
 
   addsaver(bounded_mine_percentage, "bounded_mine_percentage");
 
-  addsaver(solv::geodesic_movement, "solv_geodesic_movement", true);
+  addsaver(nisot::geodesic_movement, "solv_geodesic_movement", true);
 
 #if CAP_SHMUP  
   shmup::initConfig();
diff --git a/conformal.cpp b/conformal.cpp
index 5a948ba0..f6638273 100644
--- a/conformal.cpp
+++ b/conformal.cpp
@@ -780,8 +780,8 @@ namespace conformal {
     
     dialog::addBreak(100);
     
-    if(sol)
-      dialog::addBoolItem_action(XLAT("geodesic movement in Sol"), solv::geodesic_movement, 'G');
+    if(nonisotropic)
+      dialog::addBoolItem_action(XLAT("geodesic movement in Sol/Nil"), nisot::geodesic_movement, 'G');
 
     dialog::addBoolItem(XLAT("rotation"), do_rotate == 2, 'r');
     if(do_rotate == 0) dialog::lastItem().value = XLAT("NEVER");
diff --git a/control.cpp b/control.cpp
index f2d58a8e..aae6d3f9 100644
--- a/control.cpp
+++ b/control.cpp
@@ -64,7 +64,7 @@ movedir vectodir(const hyperpoint& P) {
 
   transmatrix U = ggmatrix(cwt.at);
   if(GDIM == 3 && WDIM == 2)  U = radar_transform * U;
-  if(solv::local_perspective_used()) U = solv::local_perspective * U;
+  if(nisot::local_perspective_used()) U = nisot::local_perspective * U;
 
   hyperpoint H = sphereflip * tC0(U);
   transmatrix Centered = sphereflip * rgpushxto0(H);
diff --git a/drawing.cpp b/drawing.cpp
index 1dea8813..4771762c 100644
--- a/drawing.cpp
+++ b/drawing.cpp
@@ -793,6 +793,20 @@ void debug_this() { }
 void dqi_poly::draw() {
   if(flags & POLY_DEBUG) debug_this();
 
+  if(nil && vid.usingGL && pmodel == mdPerspective && (current_display->set_all(global_projection), shaderside_projection)) {
+    auto npoly = *this;
+    glcoords.clear();
+    for(int i=0; i<cnt; i++)
+      glcoords.push_back(glhr::pointtogl(nisot::inverse_exp(V * glhr::gltopoint( (*tab)[offset+i]), nisot::iTable)));
+      
+    npoly.offset = 0;
+    npoly.tab = &glcoords;
+    npoly.V = Id;
+    set_width(1);
+    npoly.gldraw();
+    return;
+    }
+
   dynamicval<ld> bs(hr::band_shift, band_shift);
   if(!hyperbolic && among(pmodel, mdPolygonal, mdPolynomial)) {
     bool any = false;
diff --git a/geom-exp.cpp b/geom-exp.cpp
index f0305d04..3876a22c 100644
--- a/geom-exp.cpp
+++ b/geom-exp.cpp
@@ -725,8 +725,8 @@ void showEuclideanMenu() {
     }
   dialog::addSelItem(XLAT("projection"), current_proj_name(), '1');
   dialog::add_action_push(conformal::model_menu);
-  if(sol)
-    dialog::addBoolItem_action(XLAT("geodesic movement in Sol"), solv::geodesic_movement, 'G');
+  if(nonisotropic)
+    dialog::addBoolItem_action(XLAT("geodesic movement in Sol/Nil"), nisot::geodesic_movement, 'G');
   #if CAP_CRYSTAL && MAXMDIM >= 4
   crystal::add_crystal_transform('x');  
   #endif
@@ -803,6 +803,10 @@ void showEuclideanMenu() {
       dialog::addSelItem(XLAT("Curvature"), XLAT("Sol"), 0);
       break;
 
+    case gcNil:
+      dialog::addSelItem(XLAT("Curvature"), XLAT("Nil"), 0);
+      break;
+
     PROD( case gcProduct:
       dialog::addSelItem(XLAT("Curvature"), XLAT("Product"), 0);
       break; )
diff --git a/geometry.cpp b/geometry.cpp
index 05be75a1..63c9a01d 100644
--- a/geometry.cpp
+++ b/geometry.cpp
@@ -135,7 +135,7 @@ void geometry_information::prepare_basics() {
   #endif
   #if CAP_BT
   if(binarytiling) hexvdist = rhexf = 1, tessf = 1, scalefactor = 1, crossf = hcrossf7;
-  if(geometry == gHoroRec || penrose || sol) hexvdist = rhexf = .5, tessf = .5, scalefactor = .5, crossf = hcrossf7/2;
+  if(geometry == gHoroRec || penrose || sol || nil) hexvdist = rhexf = .5, tessf = .5, scalefactor = .5, crossf = hcrossf7/2;
   #endif
   #if CAP_BT && MAXMDIM >= 4
   if(binarytiling) binary::build_tmatrix();
diff --git a/graph.cpp b/graph.cpp
index f1d5a5b9..8b686002 100644
--- a/graph.cpp
+++ b/graph.cpp
@@ -769,13 +769,13 @@ pair<bool, hyperpoint> makeradar(hyperpoint h) {
   using namespace hyperpoint_vec;
   ld d = hdist0(h);
 
-  if(sol && solv::geodesic_movement) {
-    h = solv::inverse_exp(h, true);
+  if(sol && nisot::geodesic_movement) {
+    h = nisot::inverse_exp(h, nisot::iLazy);
     ld r = hypot_d(3, h);
     if(r < 1) h = h * (atanh(r) / r);
     else return {false, h};
     }
-  if(solv::local_perspective_used()) h = solv::local_perspective * h;
+  if(nisot::local_perspective_used()) h = nisot::local_perspective * h;
   
   if(WDIM == 3) {
     if(d >= vid.radarrange) return {false, h};
@@ -820,7 +820,7 @@ color_t kind_outline(eItem it) {
 
 transmatrix face_the_player(const transmatrix V) {
   if(DIM == 2) return V;
-  if(sol) return V * cspin(0, 2, ptick(618, 0));
+  if(nonisotropic) return V * cspin(0, 2, ptick(618, 0));
   return rgpushxto0(tC0(V));
   }
 
@@ -4429,6 +4429,7 @@ int get_darkval(int d) {
   if(binarytiling) return darkval_hbt[d];
   if(hyperbolic && S7 == 6) return darkval_e6[d];
   if(hyperbolic && S7 == 12) return darkval_s12[d];
+  if(nil) return ((d % 11) * 3) % 7;
   return 0;
   }
 
@@ -4504,7 +4505,7 @@ void make_clipping_planes() {
     hyperpoint sx = point3(y1 * z2 - y2 * z1, z1 * x2 - z2 * x1, x1 * y2 - x2 * y1);
     sx /= hypot_d(3, sx);
     sx[3] = 0;
-    if(solv::local_perspective_used()) sx = inverse(solv::local_perspective) * sx;
+    if(nisot::local_perspective_used()) sx = inverse(nisot::local_perspective) * sx;
     clipping_planes.push_back(sx);
     };
   ld tx = current_display->tanfov;
@@ -5047,7 +5048,7 @@ void drawcell(cell *c, transmatrix V, int spinv, bool mirrored) {
     }
   if(just_gmatrix) return;
 #if MAXMDIM >= 4
-  if(WDIM == 3 && pmodel == mdPerspective) {
+  if(WDIM == 3 && pmodel == mdPerspective && !nil) {
     using namespace hyperpoint_vec;
     hyperpoint H = tC0(V);
     for(hyperpoint& cpoint: clipping_planes) if((H|cpoint) < -sin_auto(cgi.corner_bonus)) {
@@ -5061,7 +5062,7 @@ void drawcell(cell *c, transmatrix V, int spinv, bool mirrored) {
     hyperpoint H = tC0(V);
     if(abs(H[0]) <= 2 && abs(H[1]) <= 2 && abs(H[2]) <= 2) ;
     else {
-      hyperpoint H2 = solv::inverse_exp(H, true);
+      hyperpoint H2 = nisot::inverse_exp(H, nisot::iLazy);
       for(hyperpoint& cpoint: clipping_planes) if((H2|cpoint) < -.2) return;
       }
     noclipped++;
@@ -5986,7 +5987,7 @@ void drawcell(cell *c, transmatrix V, int spinv, bool mirrored) {
           
           for(int a=0; a<c->type; a++)
             if(c->move(a) && !isWall3(c->move(a), dummy)) {
-              if(pmodel == mdPerspective && !sphere && !quotient && !penrose && !sol) {
+              if(pmodel == mdPerspective && !sphere && !quotient && !penrose && !nonisotropic) {
                 if(a < 4 && among(geometry, gHoroTris, gBinary3) && celldistAlt(c) >= celldistAlt(viewctr.at->c7)) continue;
                 else if(a < 2 && among(geometry, gHoroRec) && celldistAlt(c) >= celldistAlt(viewctr.at->c7)) continue;
                 else if(c->move(a)->master->distance > c->master->distance && c->master->distance > viewctr.at->distance && !quotient) continue;
@@ -7144,11 +7145,11 @@ void make_actual_view() {
       actual_view_transform = solmul(zpush(wall_radar((masterless ? centerover.at : viewctr.at->c7), inverse(View), max)), actual_view_transform * View) * inverse(View); 
     camera_level = asin_auto(tC0(inverse(actual_view_transform * View))[2]);
     }
-  if(sol) {
+  if(nonisotropic) {
     transmatrix T = actual_view_transform * View;
     transmatrix T2 = eupush( tC0(inverse(T)) );
-    solv::local_perspective = T * T2;
-    actual_view_transform = inverse(solv::local_perspective) * actual_view_transform;
+    nisot::local_perspective = T * T2;
+    actual_view_transform = inverse(nisot::local_perspective) * actual_view_transform;
     }
   #endif
   }
@@ -7938,7 +7939,7 @@ cell *viewcenter() {
 bool inscreenrange(cell *c) {
   if(sphere) return true;
   if(euclid) return celldistance(viewcenter(), c) <= get_sightrange_ambush();
-  if(sol) return gmatrix.count(c);
+  if(nonisotropic) return gmatrix.count(c);
   return heptdistance(viewcenter(), c) <= 8;
   }
 
diff --git a/hyper.h b/hyper.h
index 52894ea5..16db5ba2 100644
--- a/hyper.h
+++ b/hyper.h
@@ -112,8 +112,10 @@ void addMessage(string s, char spamtype = 0);
 #define euclid (cgclass == gcEuclid)
 #define sphere (cgclass == gcSphere)
 #define sol (cgclass == gcSol)
+#define nil (cgclass == gcNil)
 #define hyperbolic (cgclass == gcHyperbolic)
-#define eusol (euclid || sol)
+#define nonisotropic (sol || nil)
+#define translatable (euclid || nonisotropic)
 #define nonorientable (ginf[geometry].flags & qNONORIENTABLE)
 #define elliptic (ginf[geometry].flags & qELLIPTIC)
 #define quotient (ginf[geometry].flags & qANYQ)
@@ -552,7 +554,7 @@ struct gcell {
 #define NOBARRIERS 127
 #define MODFIXER (2*10090080*17)
 
-#define MAX_EDGE 18
+#define MAX_EDGE 22
 
 template<class T> struct walker;
 
@@ -5639,19 +5641,32 @@ namespace kite {
   }
 #endif
 
-namespace solv {
+/* nonisotropic */
+namespace nisot {
   extern transmatrix local_perspective;
-  inline bool local_perspective_used() { return sol; }
-  
+  inline bool local_perspective_used() { return nonisotropic; }
   hrmap *new_map();
-  hyperpoint inverse_exp(hyperpoint h, bool lazy);
+  transmatrix translate(const hyperpoint h);
+  bool in_table_range(hyperpoint h);
+  
+  enum iePrecision { iLazy, iTable };
 
-  transmatrix get_solmul(const transmatrix T, const transmatrix V);
-  transmatrix get_solmul_pt(const transmatrix Position, const transmatrix T);
+  transmatrix parallel_transport(const transmatrix Position, const transmatrix T);
+  transmatrix transport_view(const transmatrix T, const transmatrix V);
   transmatrix spin_towards(const transmatrix Position, const hyperpoint goal);
+  hyperpoint inverse_exp(const hyperpoint h, iePrecision p);
+  }
+
+namespace solv {  
   extern string solshader;
   }
 
+
+namespace nilv {
+  extern array<vector<hyperpoint>, 22> facevertices;
+  void software_renderer(dqi_poly *p);
+  }
+
 bool in_perspective();
 
 extern int noclipped;
diff --git a/hyperpoint.cpp b/hyperpoint.cpp
index 0bc9dec1..45879638 100644
--- a/hyperpoint.cpp
+++ b/hyperpoint.cpp
@@ -142,11 +142,11 @@ ld edge_of_triangle_with_angles(ld alpha, ld beta, ld gamma) {
 // (this is analogous to representing a sphere with points such that x^2+y^2+z^2 == 1)
 
 hyperpoint hpxy(ld x, ld y) { 
-  return hpxyz(x,y, eusol ? 1 : sphere ? sqrt(1-x*x-y*y) : sqrt(1+x*x+y*y));
+  return hpxyz(x,y, translatable ? 1 : sphere ? sqrt(1-x*x-y*y) : sqrt(1+x*x+y*y));
   }
 
 hyperpoint hpxy3(ld x, ld y, ld z) { 
-  return hpxyz3(x,y,z, eusol ? 1 : sphere ? sqrt(1-x*x-y*y-z*z) : sqrt(1+x*x+y*y+z*z));
+  return hpxyz3(x,y,z, translatable ? 1 : sphere ? sqrt(1-x*x-y*y-z*z) : sqrt(1+x*x+y*y+z*z));
   }
 
 // origin of the hyperbolic plane
@@ -201,7 +201,7 @@ ld dhypot_d(int d, const hyperpoint& a, const hyperpoint& b) {
   }
   
 ld zlevel(const hyperpoint &h) {
-  if(eusol) return h[GDIM];
+  if(translatable) return h[GDIM];
   else if(sphere) return sqrt(intval(h, Hypc));
   else return (h[GDIM] < 0 ? -1 : 1) * sqrt(-intval(h, Hypc));
   }
@@ -278,20 +278,15 @@ transmatrix eupush(ld x, ld y) {
   return T;
   }
 
-transmatrix eupush3(ld x, ld y, ld z) {
+transmatrix eupush(hyperpoint h) {
+  if(nonisotropic) return nisot::translate(h);
   transmatrix T = Id;
-  T[0][GDIM] = x;
-  T[1][GDIM] = y;
-  if(GDIM == 3) T[2][GDIM] = z;
-  if(sol) T[0][0] = exp(-z), T[1][1] = exp(+z);
+  for(int i=0; i<GDIM; i++) T[i][GDIM] = h[i];
   return T;
   }
 
-transmatrix eupush(hyperpoint h) {
-  transmatrix T = Id;
-  for(int i=0; i<GDIM; i++) T[i][GDIM] = h[i];
-  if(sol) T[0][0] = exp(-h[2]), T[1][1] = exp(+h[2]);
-  return T;
+transmatrix eupush3(ld x, ld y, ld z) {
+  return eupush(point3(x, y, z));
   }
 
 transmatrix euscalezoom(hyperpoint h) {
@@ -312,12 +307,8 @@ transmatrix euaffine(hyperpoint h) {
 
 transmatrix cpush(int cid, ld alpha) {
   transmatrix T = Id;
-  if(sol) {
-    T[cid][GDIM] = alpha;
-    if(cid == 2)
-      T[0][0] = exp(-alpha), T[1][1] = exp(+alpha);
-    return T;
-    }
+  if(sol || nil) 
+    return eupush3(cid == 0 ? alpha : 0, cid == 1 ? alpha : 0, cid == 2 ? alpha : 0);
   T[GDIM][GDIM] = T[cid][cid] = cos_auto(alpha);
   T[cid][GDIM] = sin_auto(alpha);
   T[GDIM][cid] = -curvature() * sin_auto(alpha);
@@ -339,7 +330,8 @@ bool eqmatrix(transmatrix A, transmatrix B, ld eps) {
 // in the 3D space, move the point h orthogonally to the (x,y) plane by z units
 hyperpoint orthogonal_move(const hyperpoint& h, ld z) {
   if(!hyperbolic) return rgpushxto0(h) * cpush(2, z) * C0;
-  if(eusol) return hpxy3(h[0], h[1], h[2] + z);
+  if(nil) return nisot::translate(h) * cpush0(2, z);
+  if(translatable) return hpxy3(h[0], h[1], h[2] + z);
   ld u = 1;
   if(h[2]) z += asin_auto(h[2]), u /= acos_auto(z);
   u *= cos_auto(z);
@@ -473,7 +465,7 @@ transmatrix rpushxto0(const hyperpoint& H) {
   }
 
 transmatrix ggpushxto0(const hyperpoint& H, ld co) {
-  if(eusol) {
+  if(translatable) {
     using namespace hyperpoint_vec;
     return eupush(co * H);
     }
@@ -505,7 +497,7 @@ transmatrix rgpushxto0(const hyperpoint& H) {
 // (without using this, imprecision could accumulate)
 
 void fixmatrix(transmatrix& T) {
-  if(sol) ; // fixed inside solmul
+  if(nonisotropic) ; // T may be inverse... do not do that
   else if(euclid) {
     for(int x=0; x<GDIM; x++) for(int y=0; y<=x; y++) {
       ld dp = 0;
@@ -784,17 +776,19 @@ bool asign(ld y1, ld y2) { return signum(y1) != signum(y2); }
 ld xcross(ld x1, ld y1, ld x2, ld y2) { return x1 + (x2 - x1) * y1 / (y1 - y2); }
 
 transmatrix solmul(const transmatrix T, const transmatrix V) {
-  if(sol) return solv::get_solmul(T, V);
+  if(nonisotropic) return nisot::transport_view(T, V);
   else return T * V;
   }
 
 transmatrix solmul_pt(const transmatrix Position, const transmatrix T) {
-  if(sol) return solv::get_solmul_pt(Position, T);
+  if(nonisotropic) return nisot::parallel_transport(Position, T);
   else return Position * T;
   }
 
 transmatrix spin_towards(const transmatrix Position, const hyperpoint goal, int dir, int back) {
-  transmatrix T = sol ? solv::spin_towards(Position, goal) : rspintox(inverse(Position) * goal);
+  transmatrix T = 
+    nonisotropic ? nisot::spin_towards(Position, goal) : 
+    rspintox(inverse(Position) * goal);
   if(back < 0) T = T * spin(M_PI);
   if(dir) T = T * cspin(dir, 0, -M_PI/2);
   T = Position * T;
diff --git a/hypgraph.cpp b/hypgraph.cpp
index 13232e7d..4777cc30 100644
--- a/hypgraph.cpp
+++ b/hypgraph.cpp
@@ -313,7 +313,7 @@ void applymodel(hyperpoint H, hyperpoint& ret) {
   switch(pmodel) {
     case mdPerspective: {
       ld ratio = vid.xres / current_display->tanfov / current_display->radius / 2;
-      if(solv::local_perspective_used()) H = solv::local_perspective * H;
+      if(nisot::local_perspective_used()) H = nisot::local_perspective * H;
       ret[0] = H[0]/H[2] * ratio;
       ret[1] = H[1]/H[2] * ratio;
       ret[2] = 1;
@@ -321,8 +321,8 @@ void applymodel(hyperpoint H, hyperpoint& ret) {
       }
 
     case mdSolPerspective: {
-      auto S = solv::inverse_exp(H, false);
-      if(solv::local_perspective_used()) S = solv::local_perspective * S;
+      auto S = nisot::inverse_exp(H, nisot::iTable);
+      if(nisot::local_perspective_used()) S = nisot::local_perspective * S;
       ld ratio = vid.xres / current_display->tanfov / current_display->radius / 2;
       ret[0] = S[0]/S[2] * ratio;
       ret[1] = S[1]/S[2] * ratio;
@@ -345,7 +345,7 @@ void applymodel(hyperpoint H, hyperpoint& ret) {
       }
     
     case mdDisk: {
-      if(solv::local_perspective_used()) H = solv::local_perspective * H;
+      if(nisot::local_perspective_used()) H = nisot::local_perspective * H;
       ld tz = get_tz(H);
       if(!vid.camera_angle) {
         ret[0] = H[0] / tz;
@@ -886,8 +886,8 @@ transmatrix actualV(const heptspin& hs, const transmatrix& V) {
 bool point_behind(hyperpoint h) {
   if(sphere) return false;
   if(!in_perspective()) return false;
-  if(pmodel == mdSolPerspective) h = solv::inverse_exp(h, true);
-  if(solv::local_perspective_used()) h = solv::local_perspective * h;
+  if(pmodel == mdSolPerspective) h = nisot::inverse_exp(h, nisot::iLazy);
+  if(nisot::local_perspective_used()) h = nisot::local_perspective * h;
   return h[2] < 0;
   }
 
@@ -910,6 +910,7 @@ bool invalid_point(const hyperpoint h) {
 bool in_smart_range(const transmatrix& T) {
   hyperpoint h = tC0(T);
   if(invalid_point(h)) return false;
+  if(nil) return cells_drawn < 2000;
   if(pmodel == mdSolPerspective) return solv::in_table_range(h);
   hyperpoint h1;
   applymodel(h, h1);
@@ -1610,7 +1611,7 @@ void queuestraight(hyperpoint X, int style, color_t lc, color_t fc, PPR p) {
 void draw_boundary(int w) {
 
   if(w == 1) return;
-  if(sol) return;
+  if(nonisotropic || euclid) return;
 
   dynamicval<ld> lw(vid.linewidth, vid.linewidth * vid.multiplier_ring);
 
@@ -1880,7 +1881,11 @@ bool do_draw(cell *c, const transmatrix& T) {
   if(WDIM == 3) {
     if(cells_drawn > vid.cells_drawn_limit) return false;
     if(pmodel == mdSolPerspective) {
-      if(!solv::in_table_range(tC0(T))) return false;
+      if(!nisot::in_table_range(tC0(T))) return false;
+      if(!limited_generation(c)) return false;
+      }
+    else if(nil) {
+      if(!nisot::in_table_range(tC0(T))) return false;
       if(!limited_generation(c)) return false;
       }
     else if(vid.use_smart_range) {
diff --git a/nonisotropic.cpp b/nonisotropic.cpp
index c0839f84..7a65e650 100644
--- a/nonisotropic.cpp
+++ b/nonisotropic.cpp
@@ -5,20 +5,37 @@
 
 namespace hr {
 
+namespace nisot {
+  typedef array<float, 3> ptlow;
+
+  transmatrix local_perspective;
+  
+  bool geodesic_movement = true;
+
+  transmatrix translate(hyperpoint h) {
+    transmatrix T = Id;
+    for(int i=0; i<GDIM; i++) T[i][GDIM] = h[i];
+    if(sol) {
+      T[0][0] = exp(-h[2]);
+      T[1][1] = exp(+h[2]);
+      }
+    if(nil) 
+      T[2][1] = h[0];
+    return T;
+    }
+
+  }
+
 namespace solv {
   
   int PRECX, PRECY, PRECZ;
   
-  typedef array<float, 3> ptlow;
-
-  vector<ptlow> inverse_exp_table;
+  vector<nisot::ptlow> inverse_exp_table;
   
   bool table_loaded;
   
   string solfname = "solv-geodesics.dat";
   
-  hyperpoint inverse_exp(hyperpoint h);
-
   void load_table() {
     if(table_loaded) return;
     FILE *f = fopen(solfname.c_str(), "rb");
@@ -28,75 +45,20 @@ namespace solv {
     fread(&PRECY, 4, 1, f);
     fread(&PRECZ, 4, 1, f);
     inverse_exp_table.resize(PRECX * PRECY * PRECZ);
-    fread(&inverse_exp_table[0], sizeof(ptlow) * PRECX * PRECY * PRECZ, 1, f);
+    fread(&inverse_exp_table[0], sizeof(nisot::ptlow) * PRECX * PRECY * PRECZ, 1, f);
     fclose(f);
     table_loaded = true;    
     }
   
-  void step2(ld& x, ld &y, ld &z, ld &vx, ld &vy, ld &vz) {
-  
-    ld ax = -2 * vx * vz;
-    ld ay = 2 * vy * vz;
-    ld az = exp(2*z) * vx * vx - exp(-2*z) * vy * vy;
-  
-    // ld x2 = x + vx/2;
-    // ld y2 = y + vy/2;
-    ld z2 = z + vz/2;
-    
-    ld vx2 = vx + ax/2;
-    ld vy2 = vy + ay/2;
-    ld vz2 = vz + az/2;
-    
-    ld ax2 = -2 * vx2 * vz2;
-    ld ay2 = 2 * vy2 * vz2;
-    ld az2 = exp(2*z2) * vx2 * vx2 - exp(-2*z2) * vy2 * vy2;
-    
-    x += vx + ax/2;
-    y += vy + ay/2;
-    z += vz + az/2;
-    
-    vx += ax2;
-    vy += ay2;
-    vz += az2;
+  hyperpoint christoffel(const hyperpoint at, const hyperpoint velocity, const hyperpoint transported) {
+    return hpxyz3(
+      -velocity[2] * transported[0] - velocity[0] * transported[2],
+       velocity[2] * transported[1] + velocity[1] * transported[2],
+       velocity[0] * transported[0] * exp(2*at[2]) - velocity[1] * transported[1] * exp(-2*at[2]),
+       0
+       );
     }
   
-  hyperpoint direct_exp(hyperpoint v, int steps) {
-    hyperpoint at;
-    at[0] = 0;
-    at[1] = 0;
-    at[2] = 0;
-    at[3] = 1;
-    v[0] /= steps;
-    v[1] /= steps;
-    v[2] /= steps;
-    for(int i=0; i<steps; i++) step2(at[0], at[1], at[2], v[0],v[1],v[2]);
-    return at;
-    }
-
-  void parallel(ld x, ld y, ld z, ld vx, ld vy, ld vz, ld& wx, ld& wy, ld& wz, ld t) {
-    // ld dax = -dz * gxyz * ay;
-    
-    ld dwx = -vz * wx - vx * wz;
-    ld dwy =  vz * wy + vy * wz;
-    ld dwz = vx * wx * exp(2*z) - vy * wy * exp(-2*z);
-      
-    wx += dwx * t;
-    wy += dwy * t;
-    wz += dwz * t;
-    }
-
-  hyperpoint direct_exp(hyperpoint v, int steps, vector<hyperpoint> transported) {
-    ld x = 0, y = 0, z = 0;
-    v[0] /= steps;
-    v[1] /= steps;
-    v[2] /= steps;
-    for(int i=0; i<steps; i++) {
-      for(auto t: transported) parallel(x,y,z, v[0],v[1],v[2], t[0], t[1], t[2], 1);
-      step2(x,y,z, v[0],v[1],v[2]);
-      }
-    return hpxy3(x, y, z);
-    }
-
   ld x_to_ix(ld u) {
     if(u == 0.) return 0.;
     ld diag = u*u/2.;
@@ -111,7 +73,7 @@ namespace solv {
     return 0.5 - atan((0.5-x) / y) / M_PI;
     }
 
-  hyperpoint inverse_exp(hyperpoint h, bool lazy) {
+  hyperpoint get_inverse_exp(hyperpoint h, bool lazy) {
     load_table();
     
     ld ix = h[0] >= 0. ? x_to_ix(h[0]) : x_to_ix(-h[0]);
@@ -160,10 +122,6 @@ namespace solv {
     return res * atanh(r) / r; */
     }
 
-  transmatrix local_perspective;
-  
-  bool geodesic_movement = true;
-
   struct hrmap_sol : hrmap {
     hrmap *binary_map;
     unordered_map<pair<heptagon*, heptagon*>, heptagon*> at;
@@ -299,7 +257,15 @@ namespace solv {
 
     };
 
-  hrmap *new_map() { return new hrmap_sol; }
+  pair<heptagon*,heptagon*> getcoord(heptagon *h) {
+    return ((hrmap_sol*)currentmap)->coords[h];
+    }
+
+  heptagon *get_at(heptagon *h1, heptagon *h2, bool gen) {
+    auto m = ((hrmap_sol*)currentmap);
+    if(!gen && !m->at.count(make_pair(h1, h2))) return nullptr;
+    return m->get_at(h1, h2);
+    }
   
   ld solrange_xy = 15;
   ld solrange_z = 4;
@@ -311,62 +277,6 @@ namespace solv {
     return abs(h[0]) < solrange_xy && abs(h[1]) < solrange_xy && abs(h[2]) < solrange_z;
     }
 
-  void fixmatrix(transmatrix& T) {
-    transmatrix push = eupush( tC0(T) );
-    transmatrix push_back = inverse(push);
-    transmatrix gtl = push_back * T;
-    { dynamicval<eGeometry> g(geometry, gSphere); hr::fixmatrix(gtl); }
-    T = push * gtl;
-    }
-  
-  transmatrix spt(transmatrix Pos, transmatrix T) {
-  
-    solv::fixmatrix(Pos);
-  
-    hyperpoint h = tC0(T);
-    h[3] = 0;
-  
-    h = Pos * h;
-    
-    int steps = 100;
-    using namespace hyperpoint_vec;
-    h /= steps;
-  
-    for(int i=0; i<steps; i++) {
-      for(int j=0; j<3; j++) 
-        parallel(Pos[0][3], Pos[1][3], Pos[2][3], 
-          h[0], h[1], h[2],
-          Pos[0][j], Pos[1][j], Pos[2][j],
-          +1);
-      step2(Pos[0][3], Pos[1][3], Pos[2][3], h[0], h[1], h[2]);
-      }
-                                                                                  
-    return Pos;
-    }
-
-  transmatrix get_solmul(const transmatrix T, const transmatrix V) {
-    if(!geodesic_movement) return V * eupush(inverse(V) * T * V * C0);
-    
-    return inverse(spt(inverse(V), inverse(T)));
-    }
-
-  transmatrix get_solmul_pt(const transmatrix Position, const transmatrix T) {
-    if(!geodesic_movement) return Position * T;
-    return spt(Position, T);
-    }
-  
-  transmatrix spin_towards(const transmatrix Position, const hyperpoint goal) {
-
-    hyperpoint at = tC0(Position);
-    transmatrix push_back = inverse(eupush(at));
-    hyperpoint back_goal = push_back * goal;
-    back_goal = inverse_exp(back_goal, false);
-    
-    transmatrix back_Position = push_back * Position;
-
-    return rspintox(inverse(back_Position) * back_goal);
-    }
-  
   int approx_distance(heptagon *h1, heptagon *h2) {
     auto m = (hrmap_sol*) currentmap;
     dynamicval<eGeometry> g(geometry, gBinary4); 
@@ -412,36 +322,335 @@ namespace solv {
     
     "return res;"
     "}";
+  }
 
-auto sol_config = addHook(hooks_args, 0, [] () {
-  using namespace arg;
-  if(argis("-solrange")) {
-    shift_arg_formula(solrange_xy);
-    shift_arg_formula(solrange_z);
-    return 0;
+namespace nilv {
+
+  hyperpoint christoffel(const hyperpoint Position, const hyperpoint Velocity, const hyperpoint Transported) {
+    ld x = Position[0];
+    return point3(
+      x * Velocity[1] * Transported[1] - 0.5 * (Velocity[1] * Transported[2] + Velocity[2] * Transported[1]),
+      -.5 * x * (Velocity[1] * Transported[0] + Velocity[0] * Transported[1]) + .5 * (Velocity[2] * Transported[0] + Velocity[0] * Transported[2]),
+      -.5 * (x*x-1) * (Velocity[1] * Transported[0] + Velocity[0] * Transported[1]) + .5 * x * (Velocity[2] * Transported[0] + Velocity[0] * Transported[2])
+      );
     }
-  else if(argis("-fsol")) {
-    shift(); solfname = args();
-    return 0;
+
+  hyperpoint formula_exp(hyperpoint v) {
+    // copying Modelling Nil-geometry in Euclidean Space with Software Presentation
+    // v[0] = c cos alpha
+    // v[1] = c sin alpha
+    // v[2] = w
+    
+    if(v[0] == 0 && v[1] == 0) return v;
+    
+    if(v[2] == 0) return point31(v[0], v[1], v[0] * v[1] / 2);
+    
+    ld alpha = atan2(v[1], v[0]);
+    ld w = v[2];
+    ld c = hypot(v[0], v[1]) / v[2];
+    
+    return point31(
+      2 * c * sin(w/2) * cos(w/2 + alpha), 
+      2 * c * sin(w/2) * sin(w/2 + alpha), 
+      w * (1 + (c*c/2) * ((1 - sin(w)/w) + (1-cos(w))/w * sin(w + 2 * alpha)))
+      );
     }
-  else if(argis("-solglitch")) {
-    shift_arg_formula(glitch_xy);
-    shift_arg_formula(glitch_z);
-    return 0;
-    }
-  else if(argis("-solgeo")) {
-    geodesic_movement = true;
-    pmodel = mdSolPerspective;
-    return 0;
-    }
-  else if(argis("-solnogeo")) {
-    geodesic_movement = false;
-    pmodel = mdPerspective;
-    return 0;
-    }
-  return 1;
-  });
   
+  hyperpoint get_inverse_exp(hyperpoint h, int iterations) {
+    ld wmin, wmax;
+    
+    ld side = h[2] - h[0] * h[1] / 2;
+    
+    if(hypot_d(2, h) < 1e-6) return h;
+    else if(side > 1e-6) {
+      wmin = 0, wmax = 2 * M_PI;
+      }
+    else if(side < -1e-6) {
+      wmin = - 2 * M_PI, wmax = 0;
+      }
+    else return point3(h[0], h[1], 0);
+    
+    ld alpha_total = atan2(h[1], h[0]);
+  
+    ld cmul;
+    if(abs(h[0]) > abs(h[1]))
+      cmul = h[0] / 2 / cos(alpha_total);
+    else
+      cmul = h[1] / 2 / sin(alpha_total);  
+    
+    ld sa = sin(2 * alpha_total);
+    
+    for(int it=0;; it++) {
+      ld w = (wmin + wmax) / 2;    
+      ld c = cmul / sin(w/2);
+      ld z = w * (1 + (c*c/2) * ((1 - sin(w)/w) + (1-cos(w))/w * sa));
+      
+      if(it == iterations) {
+        ld alpha = alpha_total - w/2;
+        return point3(c * w * cos(alpha), c * w * sin(alpha), w);
+        }
+      if(h[2] > z) wmin = w;
+      else wmax = w;    
+      }
+    }
+  
+  struct mvec : array<int, 3> {
+    
+    mvec() { }
+  
+    mvec(int x, int y, int z) { 
+      auto& a = *this;
+      a[0] = x; a[1] = y; a[2] = z;
+      }
+    mvec inverse() {  
+      auto& a = *this;
+      return mvec(-a[0], -a[1], -a[2]+a[1] * a[0]); 
+      }
+    mvec operator * (const mvec b) {
+      auto& a = *this;
+      return mvec(a[0] + b[0], a[1] + b[1], a[2] + b[2] + a[0] * b[1]);
+      }
+    };
+
+  static const mvec mvec_zero = mvec(0, 0, 0);
+
+  hyperpoint mvec_to_point(mvec m) { return hpxy3(m[0], m[1], m[2]); }
+
+  array<mvec, 22> movevectors = {mvec(0,-1,-1), mvec(1,-1,-1), mvec(-1,0,-1), mvec(0,0,-1), mvec(1,0,-1), mvec(-1,1,-1), mvec(0,1,-1), mvec(-1,-1,0), mvec(0,-1,0), mvec(-1,0,0), mvec(1,1,0), mvec(0,1,1), mvec(-1,1,0), mvec(1,0,1), mvec(0,0,1), mvec(-1,0,1), mvec(1,-1,0), mvec(0,-1,1), mvec(1,1,1), mvec(0,1,0), mvec(1,0,0), mvec(-1,-1,1), };
+
+  array<vector<hyperpoint>, 22> facevertices = {{
+    {point31(0.434375,-0.5,-0.489063), point31(0.414062,-0.417969,-0.578125), point31(0.003125,-0.49375,-0.49375), point31(0.414062,-0.574219,-0.414062), },
+    {point31(0.434375,-0.5,-0.489063), point31(0.490451,-0.560764,-0.292535), point31(0.414062,-0.574219,-0.414062), },
+    {point31(-0.49375,0,-0.49375), point31(-0.575,-0.415625,-0.170313), point31(-0.500977,-0.433594,-0.270508), point31(-0.414062,-0.414062,-0.414062), },
+    {point31(0.414062,-0.417969,-0.578125), point31(0.498047,0.00585938,-0.492188), point31(0.414062,0.414062,-0.414062), point31(-0.0178571,0.5,-0.497768), point31(-0.414062,0.417969,-0.578125), point31(-0.49375,0,-0.49375), point31(-0.414062,-0.414062,-0.414062), point31(0.003125,-0.49375,-0.49375), },
+    {point31(0.498047,0.00585938,-0.492188), point31(0.575,0.415625,-0.170313), point31(0.503906,0.434896,-0.265625), point31(0.414062,0.414062,-0.414062), },
+    {point31(-0.414062,0.574219,-0.414062), point31(-0.489258,0.560547,-0.293945), point31(-0.435268,0.504464,-0.485491), },
+    {point31(-0.0178571,0.5,-0.497768), point31(-0.414062,0.574219,-0.414062), point31(-0.435268,0.504464,-0.485491), point31(-0.414062,0.417969,-0.578125), },
+    {point31(-0.575,-0.415625,-0.170313), point31(-0.558594,-0.491211,-0.0136719), point31(-0.500977,-0.433594,-0.270508), },
+    {point31(0.003125,-0.49375,-0.49375), point31(-0.414062,-0.414062,-0.414062), point31(-0.500977,-0.433594,-0.270508), point31(-0.558594,-0.491211,-0.0136719), point31(-0.489258,-0.560547,0.293945), point31(-0.414062,-0.574219,0.414062), point31(-0.00260417,-0.492188,0.497396), point31(0.414062,-0.414062,0.414062), point31(0.503906,-0.434896,0.265625), point31(0.558594,-0.491211,0.0136719), point31(0.490451,-0.560764,-0.292535), point31(0.414062,-0.574219,-0.414062), },
+    {point31(-0.575,-0.415625,-0.170313), point31(-0.558594,-0.491211,-0.0136719), point31(-0.489258,-0.560547,0.293945), point31(-0.434375,-0.5,0.489063), point31(-0.414062,-0.417969,0.578125), point31(-0.49375,0.0015625,0.495312), point31(-0.575,0.415625,0.170313), point31(-0.558594,0.492969,0.0109375), point31(-0.489258,0.560547,-0.293945), point31(-0.435268,0.504464,-0.485491), point31(-0.414062,0.417969,-0.578125), point31(-0.49375,0,-0.49375), },
+    {point31(0.575,0.415625,-0.170313), point31(0.558594,0.492969,-0.0109375), point31(0.503906,0.434896,-0.265625), },
+    {point31(0.435268,0.504464,0.485491), point31(0.414062,0.574219,0.414062), point31(0.0146484,0.499023,0.499023), point31(0.414062,0.417969,0.578125), },
+    {point31(-0.558594,0.492969,0.0109375), point31(-0.500977,0.433594,0.270508), point31(-0.575,0.415625,0.170313), },
+    {point31(0.503906,-0.434896,0.265625), point31(0.575,-0.415625,0.170313), point31(0.497396,-0.00390625,0.494792), point31(0.414062,-0.414062,0.414062), },
+    {point31(-0.00260417,-0.492188,0.497396), point31(-0.414062,-0.417969,0.578125), point31(-0.49375,0.0015625,0.495312), point31(-0.414062,0.414062,0.414062), point31(0.0146484,0.499023,0.499023), point31(0.414062,0.417969,0.578125), point31(0.497396,-0.00390625,0.494792), point31(0.414062,-0.414062,0.414062), },
+    {point31(-0.49375,0.0015625,0.495312), point31(-0.414062,0.414062,0.414062), point31(-0.500977,0.433594,0.270508), point31(-0.575,0.415625,0.170313), },
+    {point31(0.558594,-0.491211,0.0136719), point31(0.575,-0.415625,0.170313), point31(0.503906,-0.434896,0.265625), },
+    {point31(-0.414062,-0.574219,0.414062), point31(-0.434375,-0.5,0.489063), point31(-0.414062,-0.417969,0.578125), point31(-0.00260417,-0.492188,0.497396), },
+    {point31(0.490451,0.560764,0.292535), point31(0.414062,0.574219,0.414062), point31(0.435268,0.504464,0.485491), },
+    {point31(-0.0178571,0.5,-0.497768), point31(-0.414062,0.574219,-0.414062), point31(-0.489258,0.560547,-0.293945), point31(-0.558594,0.492969,0.0109375), point31(-0.500977,0.433594,0.270508), point31(-0.414062,0.414062,0.414062), point31(0.0146484,0.499023,0.499023), point31(0.414062,0.574219,0.414062), point31(0.490451,0.560764,0.292535), point31(0.558594,0.492969,-0.0109375), point31(0.503906,0.434896,-0.265625), point31(0.414062,0.414062,-0.414062), },
+    {point31(0.414062,-0.417969,-0.578125), point31(0.498047,0.00585938,-0.492188), point31(0.575,0.415625,-0.170313), point31(0.558594,0.492969,-0.0109375), point31(0.490451,0.560764,0.292535), point31(0.435268,0.504464,0.485491), point31(0.414062,0.417969,0.578125), point31(0.497396,-0.00390625,0.494792), point31(0.575,-0.415625,0.170313), point31(0.558594,-0.491211,0.0136719), point31(0.490451,-0.560764,-0.292535), point31(0.434375,-0.5,-0.489063), },
+    {point31(-0.414062,-0.574219,0.414062), point31(-0.434375,-0.5,0.489063), point31(-0.489258,-0.560547,0.293945), },
+    }};
+  
+  struct hrmap_nil : hrmap {
+    unordered_map<mvec, heptagon*> at;
+    unordered_map<heptagon*, mvec> coords;
+    
+    heptagon *getOrigin() override { return get_at(mvec_zero); }
+    
+    heptagon *get_at(mvec c) {
+      auto& h = at[c];
+      if(h) return h;
+      h = tailored_alloc<heptagon> (S7);
+      h->c7 = newCell(S7, h);
+      coords[h] = c;
+      h->dm4 = 0;
+      h->zebraval = c[0];
+      h->emeraldval = c[1];
+      h->fieldval = c[2];
+      h->cdata = NULL;
+      h->alt = NULL;
+      return h;      
+      }
+
+    heptagon *create_step(heptagon *parent, int d) override {
+      auto p = coords[parent];
+      auto q = p * movevectors[d];
+      auto child = get_at(q);
+      parent->c.connect(d, child, (d + 11) % 22, false);
+      return child;
+      }
+
+    transmatrix adjmatrix(int i) {
+      return nisot::translate(mvec_to_point(movevectors[i]));
+      }
+    
+    virtual transmatrix relative_matrix(heptagon *h2, heptagon *h1) override { 
+      return nisot::translate(mvec_to_point(coords[h1].inverse() * coords[h2]));
+      }
+
+    void draw() override {
+      dq::visited.clear();
+
+      dq::enqueue(viewctr.at, cview());
+      
+      while(!dq::drawqueue.empty()) {
+        auto& p = dq::drawqueue.front();
+        heptagon *h = get<0>(p);
+        transmatrix V = get<1>(p);
+        dq::drawqueue.pop();
+        
+        cell *c = h->c7;
+        if(!do_draw(c, V)) continue;
+        drawcell(c, V, 0, false);
+
+        if(0) for(int t=0; t<c->type; t++) {
+          if(!c->move(t)) continue;
+          dynamicval<color_t> g(poly_outline, darkena((0x142968*t) & 0xFFFFFF, 0, 255) );
+          queuepoly(V, cgi.shWireframe3D[t], 0);
+          }
+  
+        for(int i=0; i<S7; i++) {
+          // note: need do cmove before c.spin
+          heptagon *h1 = h->cmove(i);
+          dq::enqueue(h1, V * adjmatrix(i));
+          }
+        }
+      }
+
+    };
+  }
+
+namespace nisot {
+
+  hyperpoint christoffel(const hyperpoint at, const hyperpoint velocity, const hyperpoint transported) {
+    if(sol) return solv::christoffel(at, velocity, transported);
+    else if(nil) return nilv::christoffel(at, velocity, transported);
+    else return point3(0, 0, 0);
+    }
+
+  bool in_table_range(hyperpoint h) {
+    if(sol) return solv::in_table_range(h);
+    return true;
+    }
+  
+  hyperpoint inverse_exp(const hyperpoint h, iePrecision p) {
+    if(sol) return solv::get_inverse_exp(h, p == iLazy);
+    if(nil) return nilv::get_inverse_exp(h, 10);
+    return point3(h[0], h[1], h[2]);
+    }
+
+  void geodesic_step(hyperpoint& at, hyperpoint& velocity) {
+    using namespace hyperpoint_vec;
+    auto acc = christoffel(at, velocity, velocity);
+    
+    auto at2 = at + velocity / 2;
+    auto velocity2 = velocity + acc / 2;
+    
+    auto acc2 = christoffel(at2, velocity2, velocity2);
+    
+    at = at + velocity + acc2 / 2;
+    
+    velocity = velocity + acc;
+    }
+  
+  hyperpoint direct_exp(hyperpoint v, int steps) {
+    using namespace hyperpoint_vec;
+    hyperpoint at = point31(0, 0, 0);
+    v /= steps;
+    v[3] = 0;
+    for(int i=0; i<steps; i++) geodesic_step(at, v);
+    return at;
+    }
+
+  transmatrix transpose(transmatrix T) {
+    transmatrix result;
+    for(int i=0; i<MDIM; i++)
+      for(int j=0; j<MDIM; j++)
+        result[j][i] = T[i][j];
+    return result;
+    }
+  
+  transmatrix parallel_transport_bare(transmatrix Pos, transmatrix T) {
+  
+    hyperpoint h = tC0(T);
+    h[3] = 0;
+  
+    h = Pos * h;
+    
+    int steps = 100;
+    using namespace hyperpoint_vec;
+    h /= steps;
+  
+    auto tPos = transpose(Pos);
+  
+    for(int i=0; i<steps; i++) {
+      for(int j=0; j<3; j++)
+        tPos[j] += christoffel(tPos[3], h, tPos[j]);
+      geodesic_step(tPos[3], h);
+      }
+                                                                                  
+    return transpose(tPos);
+    }
+
+  void fixmatrix(transmatrix& T) {
+    transmatrix push = eupush( tC0(T) );
+    transmatrix push_back = inverse(push);
+    transmatrix gtl = push_back * T;
+    { dynamicval<eGeometry> g(geometry, gSphere); hr::fixmatrix(gtl); }
+    T = push * gtl;
+    }
+  
+  transmatrix parallel_transport(const transmatrix Position, const transmatrix T) {
+    if(!geodesic_movement) return Position * T;
+    auto P = Position;
+    nisot::fixmatrix(P);  
+    return parallel_transport_bare(P, T);
+    }
+  
+  transmatrix transport_view(const transmatrix T, const transmatrix V) {
+    if(!geodesic_movement) return V * eupush(inverse(V) * T * V * C0);
+    return inverse(parallel_transport(inverse(V), inverse(T)));
+    }
+
+  transmatrix spin_towards(const transmatrix Position, const hyperpoint goal) {
+
+    hyperpoint at = tC0(Position);
+    transmatrix push_back = inverse(translate(at));
+    hyperpoint back_goal = push_back * goal;
+    back_goal = inverse_exp(back_goal, iTable);
+    
+    transmatrix back_Position = push_back * Position;
+
+    return rspintox(inverse(back_Position) * back_goal);
+    }
+  
+  hrmap *new_map() { 
+    if(sol) return new solv::hrmap_sol; 
+    if(nil) return new nilv::hrmap_nil;
+    return NULL;
+    }
+  
+  auto config = addHook(hooks_args, 0, [] () {
+    using namespace arg;
+    if(argis("-solrange")) {
+      shift_arg_formula(solv::solrange_xy);
+      shift_arg_formula(solv::solrange_z);
+      return 0;
+      }
+    else if(argis("-fsol")) {
+      shift(); solv::solfname = args();
+      return 0;
+      }
+    else if(argis("-solglitch")) {
+      shift_arg_formula(solv::glitch_xy);
+      shift_arg_formula(solv::glitch_z);
+      return 0;
+      }
+    else if(argis("-solgeo")) {
+      geodesic_movement = true;
+      pmodel = mdSolPerspective;
+      return 0;
+      }
+    else if(argis("-solnogeo")) {
+      geodesic_movement = false;
+      pmodel = mdPerspective;
+      return 0;
+      }
+    return 1;
+    });
+    
   }
 
 }
diff --git a/pattern2.cpp b/pattern2.cpp
index 5f7bd93a..526d47ee 100644
--- a/pattern2.cpp
+++ b/pattern2.cpp
@@ -367,7 +367,7 @@ int fieldval_uniq(cell *c) {
     auto p = cell_to_pair(c);
     return gmod(p.first * torusconfig::dx + p.second * torusconfig::dy, torusconfig::qty);
     }
-  else if(binarytiling || archimedean) return 0;
+  else if(binarytiling || archimedean || nil) return 0;
   else if(&currfp == &fp_invalid) return 0;
   else if(WDIM == 3) return c->master->fieldval;
   else if(ctof(c) || NONSTDVAR) return c->master->fieldval/S7;
diff --git a/polygons.cpp b/polygons.cpp
index 3003e59e..cde85e41 100644
--- a/polygons.cpp
+++ b/polygons.cpp
@@ -909,7 +909,7 @@ void geometry_information::create_wall3d() {
       }
     }
 
-  if(DIM == 3 && !euclid && !binarytiling) {
+  if(DIM == 3 && !euclid && !binarytiling && !nil) {
     reg3::generate();
     int facesize = isize(reg3::cellshape) / S7;
     for(int w=0; w<S7; w++) {
@@ -935,6 +935,10 @@ void geometry_information::create_wall3d() {
     make_wall(7, {pt(-1,00,+1), pt(+1,00,+1), pt(+1,-1,+1), pt(-1,-1,+1)});
     }
   
+  if(geometry == gNil) {
+    for(int i=0; i<S7; i++) make_wall(i, nilv::facevertices[i]);
+    }
+  
   if(penrose) {
     auto kv = kite::make_walls();
     for(auto& v: kv.first) for(auto& h: v) {
diff --git a/rug.cpp b/rug.cpp
index 6771c696..1599d6f8 100644
--- a/rug.cpp
+++ b/rug.cpp
@@ -948,7 +948,7 @@ bincode acd_bin(ld x) {
 
 bincode get_bincode(hyperpoint h) {
   switch(ginf[gwhere].cclass) {
-    case gcEuclid: case gcSol: PROD( case gcProduct: )
+    case gcEuclid: case gcSol: case gcNil: PROD( case gcProduct: )
       return acd_bin(h[0]) + acd_bin(h[1]) * sY + acd_bin(h[2]) * sZ;
     case gcHyperbolic:
       return acd_bin(hypot_d(3, h));