fixed graphics

hyperpoint.cpp

@@ -169,38 +169,33 @@ ld hypot_auto(ld x, ld y) {
     }
   }
 
-// get the center of the line segment from H1 to H2
-hyperpoint mid(const hyperpoint& H1, const hyperpoint& H2) {
-
-  hyperpoint H3;
-  H3[0] = H1[0] + H2[0];
-  H3[1] = H1[1] + H2[1];
-  H3[2] = H1[2] + H2[2];
-  
-  ld Z = 2;
-  
-  if(sphere) Z = sqrt(intval(H3, Hypc));
+// move H back to the sphere/hyperboloid/plane
+hyperpoint normalize(hyperpoint H) {
+  ld Z;
+  if(sphere) Z = sqrt(intval(H, Hypc));
   else if(!euclid) { 
-    Z = intval(H3, Hypc); 
+    Z = intval(H, Hypc); 
     Z = sqrt(-Z);
     }
+  else Z = H[2];
+  for(int c=0; c<3; c++) H[c] /= Z;
+  return H;
+  }
 
-  for(int c=0; c<3; c++) H3[c] /= Z;
-  
-  return H3;
+// get the center of the line segment from H1 to H2
+hyperpoint mid(const hyperpoint& H1, const hyperpoint& H2) {
+  using namespace hyperpoint_vec;
+  return normalize(H1 + H2);
   }
 
 // like mid, but take 3D into account
 hyperpoint midz(const hyperpoint& H1, const hyperpoint& H2) {
-
-  hyperpoint H3;
-  H3[0] = H1[0] + H2[0];
-  H3[1] = H1[1] + H2[1];
-  H3[2] = H1[2] + H2[2];
+  using namespace hyperpoint_vec;
+  hyperpoint H3 = H1 + H2;
   
   ld Z = 2;
   
-  if(sphere || !euclid) Z = zlevel(H3) * 2 / (zlevel(H1) + zlevel(H2));
+  if(!euclid) Z = zlevel(H3) * 2 / (zlevel(H1) + zlevel(H2));
   for(int c=0; c<3; c++) H3[c] /= Z;
   
   return H3;
@@ -327,6 +322,15 @@ void set_column(transmatrix& T, int i, const hyperpoint& H) {
     T[j][i] = H[j];
   }
 
+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;
+  }
+
 // reverse of spintox(H)
 transmatrix rspintox(const hyperpoint& H) {
   transmatrix T = Id;

hypgraph.cpp

@@ -472,31 +472,16 @@ void drawrec(cell *c, const transmatrix& V) {
     }
   } */
 
-  hyperpoint atz(const transmatrix& T, loc at) {
-    int sp = 0;
-    while(at.first < 0 || at.second < 0)
-      at = at * eudir(1), sp++;
-    if(sp>3) sp -= 6;
-
-    hyperpoint h = spin(2*M_PI*sp/S7) * T * hpxyz(at.first, at.second, 1);
-    h = mid(h,h);
-    return h;
-    }
-  
-  void drawrec(cell *c, const transmatrix& V, const transmatrix& T, whirl::loc at, int dir) {
+  void drawrec(cell *c, const transmatrix& V, whirl::loc at, int dir, int maindir) {
     if(dodrawcell(c)) {
-      hyperpoint h = atz(T, at);
-      hyperpoint hl = atz(T, at + eudir(dir));
-
-      transmatrix T1 = V * rgpushxto0(h) * rspintox(gpushxto0(h) * hl) * spin(M_PI);
-      drawcell(c, T1, 0, false);
+      drawcell(c, V * Tf[maindir][at.first&31][at.second&31][fix6(dir)], 0, false);
       }
     for(int i=0; i<c->type; i++) {
       cell *c2 = c->mov[i];
       if(!c2) continue;
       if(c2->mov[0] != c) continue;
       if(c2 == c2->master->c7) continue;
-      drawrec(c2, V, T, at + eudir(dir+i), dir + i + 3);
+      drawrec(c2, V, at + eudir(dir+i), dir + i + 3, maindir);
       }
     }
   
@@ -508,19 +493,7 @@ void drawrec(cell *c, const transmatrix& V) {
       if(!c2) continue;
       if(c2->mov[0] != c) continue;
       if(c2 == c2->master->c7) continue;
-      transmatrix T;
-      set_column(T, 0, C0);
-      set_column(T, 1, ddspin(c, i) * xpush(tessf) * C0);
-      set_column(T, 2, ddspin(c, i+1) * xpush(tessf) * C0);
-      transmatrix corners;
-      set_column(corners, 0, hpxyz(0, 0, 1));
-      set_column(corners, 1, hpxyz(whirl::param.first, whirl::param.second, 1));
-      loc nx = param * loc(0,1);
-      set_column(corners, 2, hpxyz(nx.first, nx.second, 1));
-      // corners * e[i] = corner[i]
-      T = T * inverse(corners);
-      
-      drawrec(c2, V, T, whirl::loc(1,0), 3);
+      drawrec(c2, V, whirl::loc(1,0), 3, i);
       }
     }
   }

surface.cpp

@@ -40,15 +40,6 @@ ld integral(ld 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) };
   

whirl.cpp

@@ -253,6 +253,53 @@ namespace whirl {
     WHD( printf("DONE\n\n"); )
     }
   
+  hyperpoint loctoh_ort(loc at) {
+    return hpxyz(at.first, at.second, 1);
+    }
+  
+  hyperpoint atz(const transmatrix& T, const transmatrix& corners, loc at) {
+    int sp = 0;
+    again:
+    auto corner = corners * loctoh_ort(at);
+    if(corner[1] < -1e-6 || corner[2] < -1e-6) {
+      at = at * eudir(1);
+      sp++;
+      goto again;
+      }
+    if(sp>3) sp -= 6;
+
+    return normalize(spin(2*M_PI*sp/S7) * T * corner);
+    }
+  
+  transmatrix Tf[8][32][32][6];
+  
+  void prepare_matrices() {
+    transmatrix corners = inverse(build_matrix(
+      loctoh_ort(loc(0,0)),
+      loctoh_ort(param),
+      loctoh_ort(param * loc(0,1))
+      ));
+    printf("corners = \n");
+    display(corners);
+    for(int i=0; i<S7; i++) {
+      cell cc; cc.type = S7;
+      transmatrix T = spin(-alpha) * build_matrix(
+        C0, 
+        ddspin(&cc, i) * xpush(tessf) * C0,
+        ddspin(&cc, i+1) * xpush(tessf) * C0
+        );
+      for(int x=-10; x<10; x++)
+      for(int y=-10; y<10; y++)
+      for(int d=0; d<6; d++) {
+        loc at = loc(x, y);
+        
+        hyperpoint h = atz(T, corners, at);
+        hyperpoint hl = atz(T, corners, at + eudir(d));
+        Tf[i][x&31][y&31][d] = rgpushxto0(h) * rspintox(gpushxto0(h) * hl) * spin(M_PI);
+        }
+      }       
+    }
+    
   void compute_geometry() {
     if(whirl) {
       int x = param.first;
@@ -269,6 +316,7 @@ namespace whirl {
       base_distlimit = (base_distlimit + log(scale) / log(2.618)) / scale;
       if(base_distlimit > 30)
         base_distlimit = 30;
+      prepare_matrices();
       }
     else {
       scale = 1;