From 4b3bfb99326383be1e482e67b37a23409d7f5bd9 Mon Sep 17 00:00:00 2001
From: Zeno Rogue <zeno@attnam.com>
Date: Fri, 6 Jan 2023 00:09:12 +0100
Subject: [PATCH] euc_in_product

---
 3d-models.cpp    |  7 ++++--
 cell.cpp         |  2 ++
 geometry.cpp     | 37 ++++++++++++++++++++++++++++----
 geometry2.cpp    |  6 +++++-
 graph.cpp        | 10 +++++----
 hyperpoint.cpp   | 55 ++++++++++++++++++++++++++++++++++++++----------
 hypgraph.cpp     | 29 ++++++++++++++++++-------
 nonisotropic.cpp | 14 +++++++++---
 sky.cpp          |  2 ++
 9 files changed, 129 insertions(+), 33 deletions(-)

diff --git a/3d-models.cpp b/3d-models.cpp
index 22eb696f..077cc532 100644
--- a/3d-models.cpp
+++ b/3d-models.cpp
@@ -34,6 +34,7 @@ vector<hyperpoint> geometry_information::get_shape(hpcshape sh) {
 hyperpoint get_center(const vector<hyperpoint>& vh) {
   hyperpoint h = Hypc;
   for(auto h1: vh) h = h + h1;
+  if(geom3::euc_in_product()) return h / isize(vh);
   return normalize_flat(h);
   }
 
@@ -743,8 +744,10 @@ hyperpoint psmin(hyperpoint H) {
 
 void geometry_information::adjust_eye(hpcshape& eye, hpcshape head, ld shift_eye, ld shift_head, int q, ld zoom) {
   hyperpoint center = Hypc;
-  for(int i=eye.s; i<eye.e; i++) if(q == 1 || hpc[i][1] > 0) center += hpc[i];
-  center = normalize_flat(center);
+  int c = 0;
+  for(int i=eye.s; i<eye.e; i++) if(q == 1 || hpc[i][1] > 0) center += hpc[i], c++;
+  if(geom3::euc_in_product()) center /= c;
+  else center = normalize_flat(center);
   // center /= (eye.e - eye.s);
   ld rad = 0;
   for(int i=eye.s; i<eye.e; i++) if(q == 1 || hpc[i][1] > 0) rad += hdist(center, hpc[i]);
diff --git a/cell.cpp b/cell.cpp
index 368a373b..f4bb307e 100644
--- a/cell.cpp
+++ b/cell.cpp
@@ -1269,6 +1269,8 @@ EX int clueless_celldistance(cell *c1, cell *c2) {
 
 EX int celldistance(cell *c1, cell *c2) {
 
+  if(embedded_plane) return IPF(celldistance(c1, c2));
+
   if(fake::in()) return FPIU(celldistance(c1, c2));
 
   if(mhybrid) return hybrid::celldistance(c1, c2);
diff --git a/geometry.cpp b/geometry.cpp
index 8a1a2ae5..23eae267 100644
--- a/geometry.cpp
+++ b/geometry.cpp
@@ -596,7 +596,7 @@ void geometry_information::prepare_lta() {
       lta[1][1] *= geom3::euclid_embed_scale * geom3::euclid_embed_scale_y;
       lta = cspin(0, 1, geom3::euclid_embed_rotate * degree) * lta;
       }
-    if(geom3::euc_in_nil()) lta = cspin90(2, 1) * lta;
+    if(geom3::euc_vertical()) lta = cspin90(2, 1) * lta;
     if(geom3::hyp_in_solnih()) lta = cspin90(0, 1) * cspin90(1, 2) * cspin90(0, 1) * lta;
     }
   actual_to_logical = inverse(lta);
@@ -1108,7 +1108,10 @@ EX namespace geom3 {
     seProduct,
     seNil,
     seSol, seNIH, seSolN,
-    seCliffordTorus
+    seCliffordTorus,
+    seProductH,
+    seProductS,
+    seSL2
     };
   #endif
 
@@ -1123,6 +1126,9 @@ EX namespace geom3 {
     {"stretched hyperbolic",              "Embed into stretched hyperbolic geometry. Works only with Euclidean. You need to set the variation to Pure."},
     {"stretched Sol",              "Embed into stretched Sol geometry. Works only with Euclidean. You need to set the variation to Pure."},
     {"Clifford Torus",              "Embed Euclidean torus into S3. You need to set the variation to Pure."},
+    {"hyperbolic product", "embed Euclidean or hyperbolic plane in the H2xR product space. For E2, set the variation to Pure."},
+    {"spherical product", "embed Euclidean or spherical plane in the H2xR product space. For E2, set the variation to Pure."},
+    {"SL(2,R)",           "Embed Euclidean plane in twisted product geometry. Set the variation to Pure."}
     };
 
   EX eSpatialEmbedding spatial_embedding = seDefault;
@@ -1157,6 +1163,18 @@ EX namespace geom3 {
     return ggclass() == gcNil && mgclass() == gcEuclid;
     }
 
+  EX bool euc_in_product() {
+    return ggclass() == gcProduct && mgclass() == gcEuclid;
+    }
+
+  EX bool euc_in_sl2() {
+    return ggclass() == gcSL2 && mgclass() == gcEuclid;
+    }
+
+  EX bool euc_vertical() {
+    return mgclass() == gcEuclid && among(ggclass(), gcNil, gcProduct, gcSL2);
+    }
+
   EX bool euc_in_solnih() {
     return among(ggclass(), gcSol, gcNIH, gcSolN) && mgclass() == gcEuclid;
     }
@@ -1166,7 +1184,7 @@ EX namespace geom3 {
     }
 
   EX bool euc_in_noniso() {
-    return among(ggclass(), gcNil, gcSol, gcNIH, gcSolN, gcSphere) && mgclass() == gcEuclid;
+    return among(ggclass(), gcNil, gcSol, gcNIH, gcSolN, gcSphere, gcProduct, gcSL2) && mgclass() == gcEuclid;
     }
 
   EX bool sph_in_euc() {
@@ -1234,12 +1252,18 @@ EX namespace geom3 {
           g.sig[3] = g.sig[2];
           g.sig[2] = g.sig[1];
 
-          if(spatial_embedding == seProduct && g.kind != gcEuclid) {
+          if(among(spatial_embedding, seProduct, seProductH, seProductS) && g.kind != gcEuclid) {
             g.kind = gcProduct;
             g.homogeneous_dimension--;
             g.sig[2] = g.sig[3];
             }
 
+          if(among(spatial_embedding, seProductH, seProductS) && g.kind == gcEuclid) {
+            g.kind = gcProduct;
+            g.homogeneous_dimension--;
+            g.sig[2] = spatial_embedding == seProductH ? -1 : 1;
+            }
+
           if(spatial_embedding == seLowerCurvature) {
             if(g.kind == gcEuclid) g = ginf[gSpace534].g;
             if(g.kind == gcSphere) g = ginf[gCubeTiling].g;
@@ -1279,6 +1303,11 @@ EX namespace geom3 {
             g = ginf[gSolN].g;
             g.gameplay_dimension = 2;
             }
+
+          if(spatial_embedding == seSL2 && ieuclid) {
+            g = giSL2;
+            g.gameplay_dimension = 2;
+            }
           }
         }
       }
diff --git a/geometry2.cpp b/geometry2.cpp
index 0e88adac..1d7cbd0c 100644
--- a/geometry2.cpp
+++ b/geometry2.cpp
@@ -239,6 +239,8 @@ void horo_distance::become(hyperpoint h1) {
   #endif
   else if(mhybrid)
     a = 0, b = hdist(h1, C0);
+  else if(geom3::euc_in_product())
+    a = 0, b = hdist(h1, C0);
   else
     a = 0, b = intval(h1, tile_center());
   }
@@ -249,6 +251,8 @@ horo_distance::horo_distance(shiftpoint h1, const shiftmatrix& T) {
   else
 #endif
   if(sn::in() || mhybrid || nil) become(inverse_shift(T, h1));
+  else if(geom3::euc_in_product())
+    a = 0, b = hdist(h1.h, unshift(T * tile_center(), h1.shift));
   else
     a = 0, b = intval(h1.h, unshift(T * tile_center(), h1.shift));
   }
@@ -585,7 +589,7 @@ hyperpoint hrmap_standard::get_corner(cell *c, int cid, ld cf) {
     }
   #endif
   if(PURE) {
-    if(geom3::euc_in_nil()) {
+    if(geom3::euc_in_noniso()) {
       return lspinpush0(spin_angle(c, cid) + M_PI/S7, cgi.hcrossf * 3 / cf);
       }
     return ddspin(c,cid,M_PI/S7) * lxpush0(cgi.hcrossf * 3 / cf);
diff --git a/graph.cpp b/graph.cpp
index 3ca4e5f5..ddd39579 100644
--- a/graph.cpp
+++ b/graph.cpp
@@ -350,7 +350,8 @@ EX transmatrix lpispin() {
   }
 
 EX const transmatrix& lmirror() {
-  if(geom3::euc_in_nil()) return MirrorZ;
+  if(geom3::euc_in_product()) return Id;
+  if(geom3::euc_vertical()) return MirrorZ;
   if(geom3::hyp_in_solnih()) return MirrorZ;
   return Mirror;
   }
@@ -691,7 +692,7 @@ transmatrix otherbodyparts(const shiftmatrix& V, color_t col, eMonster who, doub
 
   shiftmatrix Tright, Tleft;
   
-  if(GDIM == 2 || mhybrid) {
+  if(GDIM == 2 || mhybrid || geom3::euc_in_product()) {
     Tright = VFOOT * xpush(rightfoot);
     Tleft = VFOOT * lmirror() * xpush(-rightfoot);
     }
@@ -3145,7 +3146,8 @@ EX bool drawMonster(const shiftmatrix& Vparam, int ct, cell *c, color_t col, col
     
     if(!nospins) {
       shiftmatrix& where = (c->monst == moMirrorSpirit && inmirrorcount) ? ocwtV : cwtV;
-      if(WDIM == 2 || mproduct) {
+      if(geom3::euc_in_product()) { }
+      else if(WDIM == 2 || mproduct) {
         hyperpoint V0 = inverse_shift(Vs, where * tile_center());
         ld z = 0;
         if(gproduct) {
@@ -3818,7 +3820,7 @@ EX void pushdown(cell *c, int& q, const shiftmatrix &V, double down, bool rezoom
       auto pp = dynamic_cast<dqi_poly*> (&*ptds[q++]);
       if(!pp) continue;
       auto& ptd = *pp;
-      ptd.V = ptd.V * zpush(+down);
+      ptd.V = ptd.V * lzpush(+down);
       }
     return;
     }
diff --git a/hyperpoint.cpp b/hyperpoint.cpp
index 2f024a27..248fceac 100644
--- a/hyperpoint.cpp
+++ b/hyperpoint.cpp
@@ -567,7 +567,15 @@ EX hyperpoint ultra_normalize(hyperpoint H) {
 
 /** normalize, and in product geometry, also flatten */
 EX hyperpoint normalize_flat(hyperpoint h) {
-  if(gproduct) return product_decompose(h).second;  
+  if(gproduct) {
+    if(geom3::euc_in_product()) {
+      ld bz = zlevel(h);
+      auto h1 = h / exp(bz);
+      ld bx = atan_auto(h1[0] / h1[2]);
+      return zpush(bz) * xpush(bx) * C0;
+      }
+    return product_decompose(h).second;
+    }
   if(sl2) h = slr::translate(h) * zpush0(-atan2(h[2], h[3]));
   if(geom3::euc_in_nil()) h[1] = 0;
   if(geom3::euc_in_solnih()) h[2] = 0;
@@ -609,7 +617,7 @@ EX hyperpoint mid(const hyperpoint& H1, const hyperpoint& H2) {
     auto d1 = product_decompose(H1);
     auto d2 = product_decompose(H2);
     hyperpoint res1 = PIU( mid(d1.second, d2.second) );
-    hyperpoint res = orthogonal_move(res1, (d1.first + d2.first) / 2);
+    hyperpoint res = res1 * exp((d1.first + d2.first) / 2);
     return res;
     }
   return normalize(H1 + H2);
@@ -660,7 +668,8 @@ EX transmatrix cspin180(int a, int b) {
 
 /** rotate by alpha degrees in the XY plane */
 EX transmatrix spin(ld alpha) {
-  if(embedded_plane && geom3::euc_in_nil()) return cspin(0, 2, alpha);
+  if(embedded_plane && geom3::euc_in_product()) return Id;
+  if(embedded_plane && geom3::euc_vertical()) return cspin(0, 2, alpha);
   if(embedded_plane && geom3::hyp_in_solnih()) return cspin(1, 2, alpha);
   return cspin(0, 1, alpha);
   }
@@ -671,21 +680,24 @@ EX transmatrix unswap_spin(transmatrix T) {
 
 /** rotate by 90 degrees in the XY plane */
 EX transmatrix spin90() {
-  if(embedded_plane && geom3::euc_in_nil()) return cspin90(0, 2);
+  if(embedded_plane && geom3::euc_in_product()) return Id;
+  if(embedded_plane && geom3::euc_vertical()) return cspin90(0, 2);
   if(embedded_plane && geom3::hyp_in_solnih()) return cspin90(1, 2);
   return cspin90(0, 1);
   }
 
 /** rotate by 180 degrees in the XY plane */
 EX transmatrix spin180() {
-  if(embedded_plane && geom3::euc_in_nil()) return cspin180(0, 2);
+  if(embedded_plane && geom3::euc_in_product()) return Id;
+  if(embedded_plane && geom3::euc_vertical()) return cspin180(0, 2);
   if(embedded_plane && geom3::hyp_in_solnih()) return cspin180(1, 2);
   return cspin180(0, 1);
   }
 
 /** rotate by 270 degrees in the XY plane */
 EX transmatrix spin270() {
-  if(embedded_plane && geom3::euc_in_nil()) return cspin90(2, 0);
+  if(embedded_plane && geom3::euc_in_product()) return Id;
+  if(embedded_plane && geom3::euc_vertical()) return cspin90(2, 0);
   if(embedded_plane && geom3::hyp_in_solnih()) return cspin90(2, 1);
   return cspin90(1, 0);
   }
@@ -775,7 +787,7 @@ EX transmatrix cpush(int cid, ld alpha) {
 
 EX transmatrix lzpush(ld z) {
   if(geom3::hyp_in_solnih()) return cpush(0, z);
-  if(geom3::euc_in_nil()) return cpush(1, z);
+  if(geom3::euc_vertical()) return cpush(1, z);
   return cpush(2, z);
   }
 
@@ -853,6 +865,14 @@ EX hyperpoint orthogonal_move(const hyperpoint& h, ld z) {
     hf[3] = h[3] / ty * sin(z0);
     return hf;
     }
+  if(geom3::euc_in_product()) {
+    ld bz = zlevel(h);
+    auto h1 = h / exp(bz);
+    ld by = asin_auto(h1[1]);
+    ld bx = atan_auto(h1[0] / h1[2]);
+    by += z;
+    return zpush(bz) * xpush(bx) * ypush(by) * C0;
+    }
   if(GDIM == 2) return scale_point(h, geom3::scale_at_lev(z));
   if(gproduct) return scale_point(h, exp(z));
   if(sl2) return slr::translate(h) * cpush0(2, z);
@@ -929,6 +949,11 @@ EX void swapmatrix(transmatrix& T) {
       return;
       }
     }
+  else if(geom3::euc_in_product()) {
+    hyperpoint h1 = cgi.logical_to_actual * get_column(T, 2);
+    T = xpush(h1[0]) * zpush(h1[2]);
+    return;
+    }
   else if(geom3::in_product()) {
     /* just do nothing */
     }
@@ -950,6 +975,11 @@ EX void swapmatrix(transmatrix& T) {
 
 /** Just like swapmatrix but for hyperpoints. */
 EX void swapmatrix(hyperpoint& h) {
+  if(geom3::euc_in_product()) {
+    h = cgi.logical_to_actual * h;
+    h = xpush(h[0]) * zpush(h[2]) * C0;
+    return;
+    }
   if(geom3::in_product()) return;
   if(geom3::sph_in_euc()) { h[3] = 1; return; }
   if(geom3::sph_in_hyp()) { h[0] *= sinh(1); h[1] *= sinh(1); h[2] *= sinh(1); h[3] = cosh(1); return; }
@@ -1098,7 +1128,8 @@ EX transmatrix rspintox(const hyperpoint& H) {
   }
 
 EX transmatrix lspintox(const hyperpoint& H) {
-  if(geom3::euc_in_nil()) return spintoc(H, 0, 2);
+  if(geom3::euc_in_product()) return Id;
+  if(geom3::euc_vertical()) return spintoc(H, 0, 2);
   if(geom3::hyp_in_solnih()) return spintoc(H, 1, 2);
   if(WDIM == 2 || gproduct) return spintoc(H, 0, 1);
   transmatrix T1 = spintoc(H, 0, 1);
@@ -1106,7 +1137,8 @@ EX transmatrix lspintox(const hyperpoint& H) {
   }
 
 EX transmatrix lrspintox(const hyperpoint& H) {
-  if(geom3::euc_in_nil()) return rspintoc(H, 0, 2);
+  if(geom3::euc_in_product()) return Id;
+  if(geom3::euc_vertical()) return rspintoc(H, 0, 2);
   if(geom3::hyp_in_solnih()) return rspintoc(H, 1, 2);
   if(WDIM == 2 || gproduct) return rspintoc(H, 0, 1);
   transmatrix T1 = spintoc(H, 0, 1);
@@ -1562,7 +1594,7 @@ EX hyperpoint tile_center() {
   }
 
 EX transmatrix orthogonal_move(const transmatrix& t, double level) {
-  if(gproduct) return scale_matrix(t, exp(level));
+  if(gproduct && !geom3::euc_in_product()) return scale_matrix(t, exp(level));
   if(GDIM == 3) return t * lzpush(level);
   return scale_matrix(t, geom3::lev_to_factor(level));
   }
@@ -1849,7 +1881,7 @@ EX hyperpoint ztangent(ld z) { return ctangent(2, z); }
 /** tangent vector in logical direction Z */
 EX hyperpoint lztangent(ld z) {
   if(geom3::hyp_in_solnih()) return ctangent(0, z);
-  if(geom3::euc_in_nil()) return ctangent(1, z);
+  if(geom3::euc_vertical()) return ctangent(1, z);
   return ctangent(2, z);
   }
 
@@ -1957,6 +1989,7 @@ EX unsigned bucketer(hyperpoint h) {
     auto d = product_decompose(h);
     h = d.second;
     dx += bucketer(d.first) * 50;
+    if(geom3::euc_in_product() && in_h2xe()) h /= h[2];
     }
   dx += bucketer(h[0]) + 1000 * bucketer(h[1]) + 1000000 * bucketer(h[2]);
   if(MDIM == 4) dx += bucketer(h[3]) * 1000000001;
diff --git a/hypgraph.cpp b/hypgraph.cpp
index c8ce1dfe..b2421276 100644
--- a/hypgraph.cpp
+++ b/hypgraph.cpp
@@ -2207,16 +2207,23 @@ void ballgeometry() {
   queuereset(pmodel, PPR::CIRCLE);
   }
 
+EX transmatrix logical_to_actual_units() {
+  transmatrix T = cgi.logical_to_actual;
+  for(int i=0; i<3; i++) set_column(T, i, get_column(T, i) / hypot_d(3, get_column(T, i)));
+  return T;
+  }
+
 EX void resetview() {
   DEBBI(DF_GRAPH, ("reset view"));
   // EUCLIDEAN
   NLP = Id;
   stretch::mstretch_matrix = Id;
+  auto& vo = get_view_orientation();
   if(cwt.at) {
     centerover = cwt.at;
     View = iddspin(cwt.at, cwt.spin);
-    if(!flipplayer) View = spin180() * View;
-    if(cwt.mirrored) View = lmirror() * View;
+    if(!flipplayer) vo = spin180() * vo;
+    if(cwt.mirrored) vo = lmirror() * vo;
 
     if(centering) {
       hyperpoint vl = View * get_corner_position(cwt.at, cwt.spin);
@@ -2236,12 +2243,11 @@ EX void resetview() {
 
   adjust_eye(View, cwt.at, -1);
 
-  if(WDIM == 2) View = spin(M_PI + vid.fixed_facing_dir * degree) * View;
-  if(WDIM == 3 && !gproduct) View = cspin90(0, 2) * View;
-  if(gproduct) NLP = cspin90(0, 2);
-  View = cgi.actual_to_logical * View;
-  if(embedded_plane) get_view_orientation() = cspin90(1, 2) * get_view_orientation();
-  if(embedded_plane && vid.wall_height < 0) View = cspin180(0, 1) * View;
+  if(WDIM == 2) vo = spin(M_PI + vid.fixed_facing_dir * degree) * vo;
+  if(WDIM == 3) vo = cspin90(0, 2) * vo;
+  vo = inverse(logical_to_actual_units()) * vo;
+  if(embedded_plane) vo = cspin90(1, 2) * vo;
+  if(embedded_plane && vid.wall_height < 0) vo = cspin180(0, 1) * vo;
 
   cwtV = shiftless(View);
   current_display->which_copy = 
@@ -3349,6 +3355,13 @@ void shift_view_by_matrix(const transmatrix T, eShiftMethod sm) {
 /* like rgpushxto0 but keeps the map orientation correct */
 EX transmatrix map_relative_push(hyperpoint h) {
   if(!embedded_plane) return rgpushxto0(h);
+  if(geom3::euc_in_product()) {
+    ld bz = zlevel(h);
+    auto h1 = h / exp(bz);
+    ld by = asin_auto(h1[1]);
+    ld bx = atan_auto(h1[0] / h1[2]);
+    return zpush(bz) * xpush(bx) * ypush(by);
+    }
   if(geom3::same_in_same()) {
     ld z = -asin_auto(h[2]);
     ld u = 1 / cos_auto(z);
diff --git a/nonisotropic.cpp b/nonisotropic.cpp
index 01d3ebee..f18243e4 100644
--- a/nonisotropic.cpp
+++ b/nonisotropic.cpp
@@ -1056,7 +1056,11 @@ EX namespace hybrid {
   EX geometry_information *underlying_cgip;
 
   EX eGeometryClass under_class() {
-    if(embedded_plane) return geom3::ginf_backup[geometry].cclass;
+    if(embedded_plane) {
+      auto c = geom3::ginf_backup[geometry].cclass;
+      if(c == gcEuclid) c = cginf.g.sig[2] > 0 ? gcSphere : gcHyperbolic;
+      return c;
+      }
     return ginf[hybrid::underlying].cclass;
     }
 
@@ -1366,6 +1370,10 @@ EX namespace hybrid {
   template<class T> auto in_underlying_geometry(const T& f) -> decltype(f()) {
     if(!mhybrid && !gproduct) return f();
     if(embedded_plane) {
+      if(geom3::euc_in_product()) {
+        dynamicval<eGeometryClass> dgc(cginf.g.kind, cginf.g.sig[2] < 0 ? gcHyperbolic : gcSphere);
+        return f();
+        }
       geom3::light_flip(true);
       finalizer ff([] { geom3::light_flip(false); });
       return f();
@@ -1645,7 +1653,7 @@ EX namespace product {
   EX hyperpoint inverse_exp(hyperpoint h) {
     hyperpoint res;
     res[2] = zlevel(h);
-    h = orthogonal_move(h, -res[2]);
+    h = h * exp(-res[2]);
     ld r = hypot_d(2, h);
     if(hybrid::under_class() == gcEuclid) {
       res[0] = h[0];
@@ -1671,7 +1679,7 @@ EX namespace product {
     res[0] = h[0] * cd;
     res[1] = h[1] * cd;
     res[2] = cos_auto(d);
-    return orthogonal_move(res, h[2]);
+    return res * exp(h[2]);
     }
 
   EX bool validate_spin() {
diff --git a/sky.cpp b/sky.cpp
index 28b88e89..bd2c34cf 100644
--- a/sky.cpp
+++ b/sky.cpp
@@ -74,6 +74,8 @@ void compute_skyvertices(const vector<sky_item>& sky) {
   if(among(geom3::ggclass(), gcSol, gcSolN)) return;  /* errors */
   if(among(geom3::ggclass(), gcNil)) return;  /* errors sometimes too */
   if(geom3::hyp_in_solnih()) return;
+  if(geom3::euc_in_product()) return;
+  if(geom3::euc_in_sl2()) return;
 
   int sk = get_skybrightness();