fixed raycaster in twisted geometries; weirdly twisted spherical geometry

hyperpoint.cpp

@@ -1300,7 +1300,7 @@ EX hyperpoint tangent_length(hyperpoint dir, ld length) {
 EX hyperpoint direct_exp(hyperpoint v) {
   if(sn::in()) return nisot::numerical_exp(v);
   if(nil) return nilv::formula_exp(v);
-  if(sl2 || stretch::in()) return rots::formula_exp(v);
+  if(sl2 || stretch::in()) return stretch::mstretch ? nisot::numerical_exp(v) : rots::formula_exp(v);
   if(prod) return product::direct_exp(v);
   ld d = hypot_d(GDIM, v);
   if(d > 0) for(int i=0; i<GDIM; i++) v[i] = v[i] * sin_auto(d) / d;

nonisotropic.cpp

@@ -2282,12 +2282,99 @@ EX namespace stretch {
 
   EX ld factor;
   
+  EX bool mstretch;
+  
+  EX transmatrix m_itoa, m_atoi, m_pd;
+  EX ld ms_christoffel[3][3][3];
+  
+  EX void enable_mstretch() {
+    mstretch = true;
+
+    for(int a=0; a<4; a++)
+    for(int b=0; b<4; b++)
+      if(a==3 || b==3) m_atoi[a][b] = (a==b);
+
+    m_itoa = inverse3(m_atoi);
+    
+    for(int a=0; a<4; a++)
+    for(int b=0; b<4; b++)
+      if(a==3 || b==3)
+        m_itoa[a][b] = m_atoi[a][b] = 0;
+
+    for(int j=0; j<3; j++)
+    for(int k=0; k<3; k++) {
+      m_pd[j][k] = 0;
+      for(int i=0; i<3; i++)
+        m_pd[j][k] += m_atoi[i][j] * m_atoi[i][k];
+      }
+
+    auto& c = ms_christoffel;
+
+    ld A00 = m_pd[0][0];
+    ld A11 = m_pd[1][1];
+    ld A22 = m_pd[2][2];
+    ld A01 = m_pd[0][1] + m_pd[1][0];
+    ld A02 = m_pd[0][2] + m_pd[2][0];
+    ld A12 = m_pd[2][1] + m_pd[1][2];
+    ld B01 = A01 * A01;
+    ld B02 = A02 * A02;
+    ld B12 = A12 * A12;
+    ld B00 = A00 * A00;
+    ld B11 = A11 * A11;
+    ld B22 = A22 * A22;
+    
+    ld den = (-4*A00*A11*A22 + A00*B12 + B01*A22 - A01*A02*A12 + B02*A11);
+    
+    /* to do: these formulas are not correct for SL; also the engine currently depends on rotational symmetry, which may not hold */
+
+    c[ 0 ][ 0 ][ 0 ] =  (A01*(A01*A12 - 2*A02*A11) + A02*(2*A01*A22 - A02*A12))/den ;
+    c[ 0 ][ 0 ][ 1 ] =  (A02*(4*A11*A22 - B12)/2 + A12*(2*A01*A22 - A02*A12)/2 - (A00 - A11)*(A01*A12 - 2*A02*A11))/den ;
+    c[ 0 ][ 0 ][ 2 ] =  (-A01*(4*A11*A22 - B12)/2 + A12*(A01*A12 - 2*A02*A11)/2 - (A00 - A22)*(2*A01*A22 - A02*A12))/den ;
+    c[ 0 ][ 1 ][ 0 ] =  (A02*(4*A11*A22 - B12)/2 + A12*(2*A01*A22 - A02*A12)/2 - (A00 - A11)*(A01*A12 - 2*A02*A11))/den ;
+    c[ 0 ][ 1 ][ 1 ] =  (-A01*(A01*A12 - 2*A02*A11) + A12*(4*A11*A22 - B12))/den ;
+    c[ 0 ][ 1 ][ 2 ] =  (B01*A22 - B02*A11 + 4*B11*A22 - A11*B12 - 4*A11*B22 + B12*A22)/(4*A00*A11*A22 - A00*B12 - B01*A22 + A01*A02*A12 - B02*A11) ;
+    c[ 0 ][ 2 ][ 0 ] =  (-A01*(4*A11*A22 - B12)/2 + A12*(A01*A12 - 2*A02*A11)/2 - (A00 - A22)*(2*A01*A22 - A02*A12))/den ;
+    c[ 0 ][ 2 ][ 1 ] =  (B01*A22 - B02*A11 + 4*B11*A22 - A11*B12 - 4*A11*B22 + B12*A22)/(4*A00*A11*A22 - A00*B12 - B01*A22 + A01*A02*A12 - B02*A11) ;
+    c[ 0 ][ 2 ][ 2 ] =  -(A02*(2*A01*A22 - A02*A12) + A12*(4*A11*A22 - B12))/den ;
+    c[ 1 ][ 0 ][ 0 ] =  -(A01*(2*A00*A12 - A01*A02) + A02*(4*A00*A22 - B02))/den ;
+    c[ 1 ][ 0 ][ 1 ] =  (-A02*(2*A01*A22 - A02*A12)/2 - A12*(4*A00*A22 - B02)/2 + (A00 - A11)*(2*A00*A12 - A01*A02))/den ;
+    c[ 1 ][ 0 ][ 2 ] =  (-4*B00*A22 + A00*B02 + A00*B12 + 4*A00*B22 - B01*A22 - B02*A22)/(4*A00*A11*A22 - A00*B12 - B01*A22 + A01*A02*A12 - B02*A11) ;
+    c[ 1 ][ 1 ][ 0 ] =  (-A02*(2*A01*A22 - A02*A12)/2 - A12*(4*A00*A22 - B02)/2 + (A00 - A11)*(2*A00*A12 - A01*A02))/den ;
+    c[ 1 ][ 1 ][ 1 ] =  (A01*(2*A00*A12 - A01*A02) - A12*(2*A01*A22 - A02*A12))/den ;
+    c[ 1 ][ 1 ][ 2 ] =  (A01*(4*A00*A22 - B02)/2 + A02*(2*A00*A12 - A01*A02)/2 + (A11 - A22)*(2*A01*A22 - A02*A12))/den ;
+    c[ 1 ][ 2 ][ 0 ] =  (-4*B00*A22 + A00*B02 + A00*B12 + 4*A00*B22 - B01*A22 - B02*A22)/(4*A00*A11*A22 - A00*B12 - B01*A22 + A01*A02*A12 - B02*A11) ;
+    c[ 1 ][ 2 ][ 1 ] =  (A01*(4*A00*A22 - B02)/2 + A02*(2*A00*A12 - A01*A02)/2 + (A11 - A22)*(2*A01*A22 - A02*A12))/den ;
+    c[ 1 ][ 2 ][ 2 ] =  (A02*(4*A00*A22 - B02) + A12*(2*A01*A22 - A02*A12))/den ;
+    c[ 2 ][ 0 ][ 0 ] =  (A01*(4*A00*A11 - B01) + A02*(2*A00*A12 - A01*A02))/den ;
+    c[ 2 ][ 0 ][ 1 ] =  (4*B00*A11 - A00*B01 - 4*A00*B11 - A00*B12 + B01*A11 + B02*A11)/(4*A00*A11*A22 - A00*B12 - B01*A22 + A01*A02*A12 - B02*A11) ;
+    c[ 2 ][ 0 ][ 2 ] =  (-A01*(A01*A12 - 2*A02*A11)/2 + A12*(4*A00*A11 - B01)/2 - (A00 - A22)*(2*A00*A12 - A01*A02))/den ;
+    c[ 2 ][ 1 ][ 0 ] =  (4*B00*A11 - A00*B01 - 4*A00*B11 - A00*B12 + B01*A11 + B02*A11)/(4*A00*A11*A22 - A00*B12 - B01*A22 + A01*A02*A12 - B02*A11) ;
+    c[ 2 ][ 1 ][ 1 ] =  (-A01*(4*A00*A11 - B01) + A12*(A01*A12 - 2*A02*A11))/den ;
+    c[ 2 ][ 1 ][ 2 ] =  (A00*A01*A12 + 2*A00*A02*A11 - B01*A02 + A01*A11*A12 - A01*A12*A22 - 2*A02*B11 + 2*A02*A11*A22)/(4*A00*A11*A22 - A00*B12 - B01*A22 + A01*A02*A12 - B02*A11) ;
+    c[ 2 ][ 2 ][ 0 ] =  (-A01*(A01*A12 - 2*A02*A11)/2 + A12*(4*A00*A11 - B01)/2 - (A00 - A22)*(2*A00*A12 - A01*A02))/den ;
+    c[ 2 ][ 2 ][ 1 ] =  (A00*A01*A12 + 2*A00*A02*A11 - B01*A02 + A01*A11*A12 - A01*A12*A22 - 2*A02*B11 + 2*A02*A11*A22)/(4*A00*A11*A22 - A00*B12 - B01*A22 + A01*A02*A12 - B02*A11) ;
+    c[ 2 ][ 2 ][ 2 ] =  -(A02*(2*A00*A12 - A01*A02) + A12*(A01*A12 - 2*A02*A11))/den ;
+    
+    for(int i=0; i<3; i++)
+    for(int j=0; j<3; j++)
+    for(int k=0; k<3; k++)
+      if(c[i][j][k])
+        println(hlog, tie(i,j,k), " : ", c[i][j][k]);
+      
+    
+    println(hlog, "ATOI = ", m_atoi);
+    println(hlog, "ITOA = ", m_itoa, " vs ", 1/not_squared());
+    println(hlog, "PD   = ", m_pd, " vs ", factor);
+
+    ray::reset_raycaster();
+    }
+
   EX bool applicable() {
     return rotspace || among(geometry, gCell120, gECell120, gCell24, gECell24, gCell8, gECell8);
     }
 
   EX bool in() {
-    return factor && applicable();
+    return (factor || mstretch) && applicable();
     }
 
   EX transmatrix translate(hyperpoint h) {
@@ -2308,9 +2395,9 @@ EX namespace stretch {
     return translate(h);
     }
 
-  hyperpoint mulz(const hyperpoint at, const hyperpoint velocity, ld factor) {
+  hyperpoint mulz(const hyperpoint at, const hyperpoint velocity, ld zf) {
     auto vel = itranslate(at) * velocity;
-    vel[2] *= factor;
+    vel[2] *= zf;
     return translate(at) * vel;
     }
   
@@ -2323,10 +2410,16 @@ EX namespace stretch {
     }
   
   EX hyperpoint isometric_to_actual(const hyperpoint at, const hyperpoint velocity) {
+    if(mstretch)
+      return translate(at) * m_itoa * itranslate(at) * velocity;
+    else
       return mulz(at, velocity, 1/not_squared());
     }
  
   EX hyperpoint actual_to_isometric(const hyperpoint at, const hyperpoint velocity) {
+    if(mstretch)
+      return translate(at) * m_atoi * itranslate(at) * velocity;
+    else
       return mulz(at, velocity, not_squared());
     }
   
@@ -2337,14 +2430,21 @@ EX namespace stretch {
     
     hyperpoint c;
     
+    if(mstretch) {    
+      c = Hypc;
+      for(int i=0; i<3; i++)
+      for(int j=0; j<3; j++)
+      for(int k=0; k<3; k++)
+        c[i] += vel[j] * tra[k] * ms_christoffel[i][j][k];
+      }
+    
+    else {    
       auto K = factor;    
-    
+      c[0] = (sphere ? -K : K+2) * (vel[1] * tra[2] + vel[2] * tra[1]);
-    if(!sphere) K = -2 - K;
+      c[1] = (sphere ? K : -(K+2)) * (vel[0] * tra[2] + vel[2] * tra[0]);
-    
-    c[0] = -K * (vel[1] * tra[2] + vel[2] * tra[1]);
-    c[1] =  K * (vel[0] * tra[2] + vel[2] * tra[0]);
       c[2] = 0;
       c[3] = 0;
+      }
 
     return translate(at) * c;
     }  
@@ -2591,7 +2691,14 @@ EX namespace nisot {
         auto fix = [&] (hyperpoint& h, ld& m) {
           h = stretch::itranslate(at) * h;
           h[3] = 0;
-          ld m1 = h[0] * h[0] + h[1] * h[1] + h[2] * h[2] * stretch::squared();
+          ld m1;
+          if(stretch::mstretch) {
+            m1 = 0;
+            for(int i=0; i<3; i++) for(int j=0; j<3; j++)
+              m1 += h[i] * stretch::m_pd[i][j] * h[j];
+            }
+          else
+            m1 = h[0] * h[0] + h[1] * h[1] + h[2] * h[2] * stretch::squared();
           h /= sqrt(m1/m);
           h = stretch::translate(at) * h;
           };
@@ -2746,6 +2853,29 @@ EX namespace nisot {
       shift_arg_formula(stretch::factor, ray::reset_raycaster);
       return 0;
       }
+    else if(argis("-mstretch")) {
+      PHASEFROM(2);
+      auto& M = stretch::m_atoi;
+      M = Id;
+      while(true) {
+        shift();
+        string s = args();
+        if(isize(s) == 2 && among(s[0], 'a', 'b','c') && among(s[1], 'a', 'b', 'c'))
+          shift_arg_formula(M[s[0]-'a'][s[1]-'a'], stretch::enable_mstretch);
+        else break;
+        }
+      // shift_arg_formula(stretch::yfactor, ray::reset_raycaster);
+      return 0;
+      }
+    else if(argis("-mstretch1")) {
+      PHASEFROM(2);
+      auto& M = stretch::m_atoi;
+      M = Id;
+      M[2][2] = stretch::not_squared();
+      stretch::enable_mstretch();
+      // shift_arg_formula(stretch::yfactor, ray::reset_raycaster);
+      return 0;
+      }
     else if(argis("-prodturn")) {
       PHASEFROM(2);
       if(prod) stop_game();

raycaster.cpp

@@ -108,6 +108,7 @@ struct raycaster : glhr::GLprogram {
   GLint uBLevel;
   GLint uPosX, uPosY;
   GLint uWallOffset, uSides;
+  GLint uITOA, uATOI;
   
   raycaster(string vsh, string fsh);
   };
@@ -151,6 +152,9 @@ raycaster::raycaster(string vsh, string fsh) : GLprogram(vsh, fsh) {
     
     uPosX = glGetUniformLocation(_program, "uPosX");
     uPosY = glGetUniformLocation(_program, "uPosY");
+    
+    uITOA = glGetUniformLocation(_program, "uITOA");
+    uATOI = glGetUniformLocation(_program, "uATOI");
     }
 
 shared_ptr<raycaster> our_raycaster;
@@ -388,8 +392,17 @@ void enable_raycaster() {
       
     if(stretch::in()) {
       fmain += 
-        "tangent = s_itranslate(position) * tangent;\n"
+        "tangent = s_itranslate(position) * tangent;\n";
-        "tangent[2] /= " + to_glsl(stretch::not_squared()) + ";\n"
+      if(stretch::mstretch) {
+        fsh += "mediump uniform mat4 uITOA;\n";
+        fsh += "mediump uniform mat4 uATOI;\n";
+        fmain += 
+          "tangent = uITOA * tangent;\n";
+        }
+      else 
+        fmain +=
+          "tangent[2] /= " + to_glsl(stretch::not_squared()) + ";\n";
+      fmain +=
         "tangent = s_translate(position) * tangent;\n";
         ;
       }
@@ -576,22 +589,39 @@ void enable_raycaster() {
         "  }\n";
         use_christoffel = false;
         }
-      else if(sl2) {
+      else if(sl2 || stretch::in()) {
+        if(sl2) {
           fsh += "mediump mat4 s_translate(vec4 h) {\n"
             "return mat4(h.w,h.z,h.y,h.x,-h.z,h.w,-h.x,h.y,h.y,-h.x,h.w,-h.z,h.x,h.y,h.z,h.w);\n"
             "}\n";
-        fsh += "mediump mat4 s_itranslate(vec4 h) {\n"
-          "h.xyz = -h.xyz; return s_translate(h);\n"
-          "}\n";
-        use_christoffel = false;
           }
-      else if(stretch::in()) {
+        else {
           fsh += "mediump mat4 s_translate(vec4 h) {\n"
             "return mat4(h.w,h.z,-h.y,-h.x,-h.z,h.w,h.x,-h.y,h.y,-h.x,h.w,-h.z,h.x,h.y,h.z,h.w);\n"
             "}\n";
+          }
         fsh += "mediump mat4 s_itranslate(vec4 h) {\n"
           "h.xyz = -h.xyz; return s_translate(h);\n"
           "}\n";
+        if(stretch::mstretch) {
+          fsh += "mediump vec4 christoffel(mediump vec4 pos, mediump vec4 vel, mediump vec4 tra) {\n"
+            "vel = s_itranslate(pos) * vel;\n"
+            "tra = s_itranslate(pos) * tra;\n"
+            "return s_translate(pos) * vec4(\n";
+          
+          for(int i=0; i<3; i++) {
+            auto &c = stretch::ms_christoffel;
+            fsh += "  0.";
+            for(int j=0; j<3; j++) 
+            for(int k=0; k<3; k++) 
+              if(c[i][j][k])
+                fsh += "  + vel["+its(j)+"]*tra["+its(k)+"]*" + to_glsl(c[i][j][k]);
+            fsh += "  ,\n";
+            }
+          fsh += "  0);\n"
+            "}\n";
+          }
+        else
           use_christoffel = false;
         }
       else use_christoffel = false;
@@ -605,7 +635,7 @@ void enable_raycaster() {
       fmain += 
         "  dist = next < minstep ? 2.*next : next;\n";
 
-      if(nil) fsh += 
+      if(nil && !use_christoffel) fsh += 
         "mediump vec4 translate(mediump vec4 a, mediump vec4 b) {\n"
           "return vec4(a[0] + b[0], a[1] + b[1], a[2] + b[2] + a[0] * b[1], b[3]);\n"
           "}\n"
@@ -619,7 +649,7 @@ void enable_raycaster() {
           "return vec4(t[0], t[1], t[2] - a[0] * t[1], 0.);\n"
           "}\n";
                 
-      if(nil) fmain += "tangent = translate(position, itranslate(position, tangent));\n";
+      // if(nil) fmain += "tangent = translate(position, itranslate(position, tangent));\n";
       
       if(use_christoffel) fmain +=
         "mediump vec4 vel = tangent * dist;\n"
@@ -629,14 +659,23 @@ void enable_raycaster() {
         "mediump vec4 acc4 = get_acc(position + vel + acc2/2., vel + acc3/2.);\n"
         "mediump vec4 nposition = position + vel + (acc1+acc2+acc3)/6.;\n";
 
-      if(sl2 || stretch::in()) {
+      if((sl2 || stretch::in()) && use_christoffel) {
+        if(sl2) fmain += 
+          "nposition = nposition / sqrt(dot(position.zw, position.zw) - dot(nposition.xy, nposition.xy));\n";
+
+        else if(stretch::in()) fmain += 
+          "nposition = nposition / sqrt(dot(nposition, nposition));\n";
+        }
+      
+      if((sl2 || stretch::in()) && !use_christoffel) {
         ld SV = stretch::not_squared();
         ld mul = (sphere?1:-1)-1/SV/SV;
         fmain += 
           "vec4 vel = s_itranslate(position) * tangent * dist;\n"
-          "mediump float vlen = length(vel.xyz);\n"
+          "vec4 vel1 = vel; vel1.z *= " + to_glsl(stretch::not_squared()) + ";\n"
+          "mediump float vlen = length(vel1.xyz);\n"
           "if(vel.z<0.) vlen=-vlen;\n"
-          "float z_part = vel.z/vlen;\n"
+          "float z_part = vel1.z/vlen;\n"
           "float x_part = sqrt(1.-z_part*z_part);\n"
           "const float SV = " + to_glsl(SV) + ";\n"
           "float rparam = x_part / z_part / SV;\n"
@@ -718,13 +757,12 @@ void enable_raycaster() {
             "}\n";          
         
         fmain +=    
-          "ntangent = s_itranslate(nposition) * ntangent / dist;\n"
+          "ntangent = ntangent / dist;\n"
-          "ntangent.z *= SV;\n"
+          "ntangent = s_translate(position) * ntangent;\n"
-          "nposition = s_translate(position) * nposition;\n"
+          "nposition = s_translate(position) * nposition;\n";
-          "ntangent = s_translate(nposition) * ntangent;\n";
         }
       
-      if(nil) {
+      if(nil && !use_christoffel) {
         fmain +=
           "mediump vec4 xp, xt;\n"
           "mediump vec4 back = itranslatev(position, tangent);\n"
@@ -876,15 +914,31 @@ void enable_raycaster() {
           "next = maxstep;\n"
           "}\n";
       
-      fmain +=
-        "position = nposition;\n";
-      
       if(use_christoffel) fmain +=
         "tangent = tangent + (acc1+2.*acc2+2.*acc3+acc4)/(6.*dist);\n";
       else if(nil) fmain +=
         "tangent = translatev(position, xt);\n";
       else fmain +=
         "tangent = ntangent;\n";
+
+      fmain +=
+        "position = nposition;\n";
+      
+      if((stretch::in() || sl2) && use_christoffel) {
+        fmain += 
+          "tangent = s_itranslate(position) * tangent;\n"
+          "tangent[3] = 0.;\n";
+        if(stretch::mstretch)
+          fmain +=
+            "float nvelsquared = dot(tangent.xyz, (uATOI * tangent).xyz);\n";
+        else
+          fmain +=
+            "float nvelsquared = tangent.x * tangent.x + tangent.y * tangent.y + "
+              + to_glsl(stretch::squared()) + " * tangent.z * tangent.z;\n";
+        fmain +=      
+          "tangent /= sqrt(nvelsquared);\n"
+          "tangent = s_translate(position) * tangent;\n";
+        }
       }
     else fmain += 
       "position = position + tangent * dist;\n";
@@ -1270,6 +1324,11 @@ EX void cast() {
   glUniform1f(o->uIPD, vid.ipd);
   GLERR("uniform mediump IPD");
   
+  if(o->uITOA != -1) {
+    glUniformMatrix4fv(o->uITOA, 1, 0, glhr::tmtogl_transpose3(stretch::m_itoa).as_array());   
+    glUniformMatrix4fv(o->uATOI, 1, 0, glhr::tmtogl_transpose3(stretch::m_atoi).as_array());   
+    }
+  
   if(o->uWallOffset != -1) {
     glUniform1i(o->uWallOffset, wall_offset(cs));
     glUniform1i(o->uSides, cs->type + (WDIM == 2 ? 2 : 0));