From e636d3f40c4e32086a325f77f1e9ac6e6491e8d2 Mon Sep 17 00:00:00 2001
From: Zeno Rogue <zeno@attnam.com>
Date: Sun, 18 Aug 2019 16:58:44 +0200
Subject: [PATCH] preliminary S2xE

---
 control.cpp      |  2 +-
 drawing.cpp      | 31 ++++++++++++++++++++++++++-----
 hyperpoint.cpp   | 13 +++++++++++--
 nonisotropic.cpp | 10 +++++++---
 4 files changed, 45 insertions(+), 11 deletions(-)

diff --git a/control.cpp b/control.cpp
index 3bbdd7a9..49073450 100644
--- a/control.cpp
+++ b/control.cpp
@@ -113,7 +113,7 @@ EX void remission() {
  }
 
 EX hyperpoint move_destination_vec(int d) {
-  hyperpoint Forward = prod ? forward_dir(1) : tC0(pushone());
+  hyperpoint Forward = prod ? forward_dir(.5) : tC0(pushone());
   if(WDIM == 2) return spin(-d * M_PI/4) * tC0(pushone());
   // else if(WDIM == 2 && pmodel == mdPerspective) return cspin(0, 2, d * M_PI/4) * tC0(pushone());
   // else if(WDIM == 2) return spin(-d * M_PI/4) * tC0(pushone());
diff --git a/drawing.cpp b/drawing.cpp
index 144ade54..8770d0c7 100644
--- a/drawing.cpp
+++ b/drawing.cpp
@@ -890,20 +890,41 @@ ld get_width(dqi_poly* p) {
 
 void debug_this() { }
 
+glvertex junk = glhr::makevertex(0,0,1);
+
 void dqi_poly::draw() {
   if(flags & POLY_DEBUG) debug_this();
 
   if(prod && 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(product::inverse_exp(V * glhr::gltopoint( (*tab)[offset+i]))));
-      
     npoly.offset = 0;
     npoly.tab = &glcoords;
     npoly.V = Id;
     set_width(1);
-    npoly.gldraw();
+    if(product::product_sphere()) {
+      for(int gen=-5; gen<=5; gen++) {
+        glcoords.clear();
+        int junks = 0;
+        for(int i=0; i<cnt; i++) {
+          hyperpoint h = V * glhr::gltopoint( (*tab)[offset+i]);
+          if(hypot_d(2, h) < 1e-6 && (gen || h[2] < 0)) {
+            junks = 3 - (i % 3);
+            }
+          if(junks == 3)
+            glcoords.push_back(junk), junks--;
+          else if(junks)
+            glcoords.push_back(glcoords.back()), junks--;
+          else        
+            glcoords.push_back(glhr::pointtogl(product::inverse_exp(h, gen)));
+          }
+        npoly.gldraw();
+        }
+      }
+    else {
+      glcoords.clear();
+      for(int i=0; i<cnt; i++) glcoords.push_back(glhr::pointtogl(product::inverse_exp(V * glhr::gltopoint( (*tab)[offset+i]))));
+      npoly.gldraw();
+      }
     return;
     }
 
diff --git a/hyperpoint.cpp b/hyperpoint.cpp
index b2990eaf..6f42bf18 100644
--- a/hyperpoint.cpp
+++ b/hyperpoint.cpp
@@ -171,7 +171,7 @@ ld inverse_tanh(ld x) { return log((1+x)/(1-x)) / 2; } */
 
 EX ld squar(ld x) { return x*x; }
 
-EX int sig(int z) { return prod ? (z<2?1:-1) : (sphere || sol || z<GDIM)?1:-1; }
+EX int sig(int z) { return prod ? PIU(sig(z)) : (sphere || sol || z<GDIM)?1:-1; }
 
 EX int curvature() {
   switch(cgclass) {
@@ -230,6 +230,15 @@ EX ld asin_auto_clamp(ld x) {
     }
   }
 
+EX ld acos_auto_clamp(ld x) {
+  switch(cgclass) {
+    case gcHyperbolic: return x < 1 ? 0 : acosh(x);
+    case gcSphere: return x > 1 ? 0 : x < -1 ? M_PI : acos(x);
+    case gcProduct: return PIU(acos_auto_clamp(x));
+    default: return x;
+    }
+  }
+
 EX ld cos_auto(ld x) {
   switch(cgclass) {
     case gcEuclid: return 1;
@@ -339,7 +348,7 @@ EX ld hypot_d(int d, const hyperpoint& h) {
 EX ld zlevel(const hyperpoint &h) {
   if(translatable) return h[LDIM];
   else if(sphere) return sqrt(intval(h, Hypc));
-  else if(prod) return log(sqrt(-intval(h, Hypc)));
+  else if(prod) return log(sqrt(abs(intval(h, Hypc)))); /* abs works with both underlying spherical and hyperbolic */
   else return (h[LDIM] < 0 ? -1 : 1) * sqrt(-intval(h, Hypc));
   }
 
diff --git a/nonisotropic.cpp b/nonisotropic.cpp
index c850f3be..2f77d5cd 100644
--- a/nonisotropic.cpp
+++ b/nonisotropic.cpp
@@ -640,18 +640,22 @@ EX namespace product {
     dynamicval<geometry_information*> gc(cgip, underlying_cgip);
     return mscale(get_corner_position(c, i), exp(plevel * z/2));
     }
+
+  EX bool product_sphere() { return ginf[underlying].cclass == gcSphere; }
   
-  EX hyperpoint inverse_exp(hyperpoint h) {
+  EX hyperpoint inverse_exp(hyperpoint h, int gen IS(0)) {
     hyperpoint res;
     res[2] = zlevel(h);
     h = zshift(h, -res[2]);
     ld r = hypot_d(2, h);
-    if(r < 1e-6 || h[2] < 1) {
+    if(r < 1e-6) {
       res[0] = h[0];
       res[1] = h[1];
       }
     else {
-      r = acosh(h[2]) / r;
+      auto c = acos_auto_clamp(h[2]);
+      c += 2 * M_PI * gen;
+      r = c / r;
       res[0] = h[0] * r;
       res[1] = h[1] * r;
       }