From 7e85f074584767f8fc2666c5bfc6bbe86bc54438 Mon Sep 17 00:00:00 2001 From: Zeno Rogue Date: Wed, 22 Jul 2020 00:19:13 +0200 Subject: [PATCH] primitive rendering now implemented for stretched H2xE; removed the old implementation of SL(2,R) --- drawing.cpp | 1 - geom-exp.cpp | 3 +- glhr.cpp | 7 +- hyperpoint.cpp | 11 +- nonisotropic.cpp | 625 +++++++++++++++++++++++++++-------------------- raycaster.cpp | 2 +- shaders.cpp | 1 + 7 files changed, 381 insertions(+), 269 deletions(-) diff --git a/drawing.cpp b/drawing.cpp index 4a7b4efa..73ba645a 100644 --- a/drawing.cpp +++ b/drawing.cpp @@ -2039,7 +2039,6 @@ EX void draw_main() { if(ray::in_use && !ray::comparison_mode) { ray::cast(); reset_projection(); - if(stretch::in()) return; /*primitive not implemented */ } DEBB(DF_GRAPH, ("outcircle")); diff --git a/geom-exp.cpp b/geom-exp.cpp index 22c6b287..128fea35 100644 --- a/geom-exp.cpp +++ b/geom-exp.cpp @@ -874,8 +874,7 @@ EX void showEuclideanMenu() { dialog::editNumber(stretch::factor, -1, 9, 0.1, 0, XLAT("stretched geometry"), XLAT( "Stretch the metric along the fibers. This can currently be done in rotation spaces and in 8-cell, 24-cell and 120-cell. " - "Value of 0 means not stretched, -1 means S2xE or H2xE (works only in the limit). " - "Only the raycaster is implemented for stretched geometry, so you will see only walls. (Must be > -1)" + "Value of 0 means not stretched, -1 means S2xE or H2xE (works only in the limit). (Must be > -1)" ) ); dialog::reaction = ray::reset_raycaster; diff --git a/glhr.cpp b/glhr.cpp index 65a162cc..f050d17c 100644 --- a/glhr.cpp +++ b/glhr.cpp @@ -274,7 +274,7 @@ struct GLprogram { GLuint vertShader, fragShader; GLint uFog, uFogColor, uColor, tTexture, tInvExpTable, uMV, uProjection, uAlpha, uFogBase, uPP; - GLint uPRECX, uPRECY, uPRECZ, uIndexSL, uIterations, uLevelLines; + GLint uPRECX, uPRECY, uPRECZ, uIndexSL, uIterations, uLevelLines, uSV; flagtype shader_flags; @@ -397,6 +397,7 @@ GLprogram::GLprogram(string vsh, string fsh) { uPRECY = glGetUniformLocation(_program, "PRECY"); uPRECZ = glGetUniformLocation(_program, "PRECZ"); uIndexSL = glGetUniformLocation(_program, "uIndexSL"); + uSV = glGetUniformLocation(_program, "uSV"); uIterations = glGetUniformLocation(_program, "uIterations"); uLevelLines = glGetUniformLocation(_program, "uLevelLines"); } @@ -411,6 +412,10 @@ EX void set_index_sl(ld x) { glUniform1f(glhr::current_glprogram->uIndexSL, x); } +EX void set_sv(ld x) { + glUniform1f(glhr::current_glprogram->uSV, x); + } + EX void set_sl_iterations(int steps) { glUniform1i(glhr::current_glprogram->uIterations, steps); } diff --git a/hyperpoint.cpp b/hyperpoint.cpp index 9c9ea8dd..20eba419 100644 --- a/hyperpoint.cpp +++ b/hyperpoint.cpp @@ -211,6 +211,7 @@ EX ld sin_auto(ld x) { case gcHyperbolic: return sinh(x); case gcSphere: return sin(x); case gcProduct: return PIU(sin_auto(x)); + case gcSL2: return sinh(x); default: return x; } } @@ -221,6 +222,7 @@ EX ld asin_auto(ld x) { case gcHyperbolic: return asinh(x); case gcSphere: return asin(x); case gcProduct: return PIU(asin_auto(x)); + case gcSL2: return asinh(x); default: return x; } } @@ -230,6 +232,7 @@ EX ld acos_auto(ld x) { case gcHyperbolic: return acosh(x); case gcSphere: return acos(x); case gcProduct: return PIU(acos_auto(x)); + case gcSL2: return acosh(x); default: return x; } } @@ -266,6 +269,7 @@ EX ld asin_auto_clamp(ld x) { switch(cgclass) { case gcEuclid: return x; case gcHyperbolic: return asinh(x); + case gcSL2: return asinh(x); case gcSphere: return asin_clamp(x); default: return x; } @@ -274,6 +278,7 @@ 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 gcSL2: 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; @@ -284,6 +289,7 @@ EX ld cos_auto(ld x) { switch(cgclass) { case gcEuclid: return 1; case gcHyperbolic: return cosh(x); + case gcSL2: return cosh(x); case gcSphere: return cos(x); case gcProduct: return PIU(cos_auto(x)); default: return 1; @@ -296,6 +302,7 @@ EX ld tan_auto(ld x) { case gcHyperbolic: return tanh(x); case gcSphere: return tan(x); case gcProduct: return PIU(tan_auto(x)); + case gcSL2: return tanh(x); default: return 1; } } @@ -306,6 +313,7 @@ EX ld atan_auto(ld x) { case gcHyperbolic: return atanh(x); case gcSphere: return atan(x); case gcProduct: return PIU(atan_auto(x)); + case gcSL2: return atanh(x); default: return x; } } @@ -314,6 +322,7 @@ EX ld atan2_auto(ld y, ld x) { switch(cgclass) { case gcEuclid: return y/x; case gcHyperbolic: return atanh(y/x); + case gcSL2: return atanh(y/x); case gcSphere: return atan2(y, x); case gcProduct: return PIU(atan2_auto(y, x)); default: return y/x; @@ -1207,7 +1216,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) return slr::formula_exp(v); + if(sl2 || stretch::in()) return rots::formula_exp(v); if(prod) return product::direct_exp(v); ld d = hypot_d(GDIM, v); if(d > 0) for(int i=0; i2 and 1<->3 are swapped, - // then coordinates 2<->3 are swapped - */ - EX ld range_xy = 2; EX int steps = 15; - EX hyperpoint from_phigans(hyperpoint h) { - ld r = asinh(hypot_d(2, h)); - ld x = h[0]; - ld y = h[1]; - ld z = h[2]; - return hyperpoint(x * cos(z) + y * sin(z), y * cos(z) - x * sin(z), cosh(r) * sin(z), cosh(r) * cos(z)); - } - - EX hyperpoint to_phigans(hyperpoint h) { - ld z = atan2(h[2], h[3]); - ld x = h[0]; - ld y = h[1]; - return point31(x * cos(z) - y * sin(z), y * cos(z) + x * sin(z), z); - } - - /** in the 'phigans' model */ - hyperpoint christoffel(const hyperpoint Position, const hyperpoint Velocity, const hyperpoint Transported) { - ld x = Position[0]; - ld y = Position[1]; - ld s = x*x + y*y + 1; - ld x2 = x * x; - ld y2 = y * y; - ld x4 = x2 * x2; - ld y4 = y2 * y2; - return point3( - + Velocity[ 0 ] * Transported[ 0 ] * (x*(4*x2*y2 + 4*y4 + 9*y2 + 1)) - + Velocity[ 0 ] * Transported[ 1 ] * (-y*(4*x4 + 4*x2*y2 + 9*x2 + 2)) - + Velocity[ 0 ] * Transported[ 2 ] * (-x*y*(x2 + y2 + 2)) - + Velocity[ 1 ] * Transported[ 0 ] * (-y*(4*x4 + 4*x2*y2 + 9*x2 + 2)) - + Velocity[ 1 ] * Transported[ 1 ] * (x*(4*x4 + 4*x2*y2 + 9*x2 + 5)) - + Velocity[ 1 ] * Transported[ 2 ] * (x4 + x2*y2 + 2*x2 + 1) - + Velocity[ 2 ] * Transported[ 0 ] * (-x*y*(x2 + y2 + 2)) - + Velocity[ 2 ] * Transported[ 1 ] * (x4 + x2*y2 + 2*x2 + 1), - - + Velocity[ 0 ] * Transported[ 0 ] * (y*(4*x2*y2 + 4*y4 + 9*y2 + 5)) - + Velocity[ 0 ] * Transported[ 1 ] * (-x*(4*x2*y2 + 4*y4 + 9*y2 + 2)) - + Velocity[ 0 ] * Transported[ 2 ] * (-x2*y2 - y4 - 2*y2 - 1) - + Velocity[ 1 ] * Transported[ 0 ] * (-x*(4*x2*y2 + 4*y4 + 9*y2 + 2)) - + Velocity[ 1 ] * Transported[ 1 ] * (y*(4*x4 + 4*x2*y2 + 9*x2 + 1)) - + Velocity[ 1 ] * Transported[ 2 ] * (x*y*(x2 + y2 + 2)) - + Velocity[ 2 ] * Transported[ 0 ] * (-x2*y2 - y4 - 2*y2 - 1) - + Velocity[ 2 ] * Transported[ 1 ] * (x*y*(x2 + y2 + 2)), - - + Velocity[ 0 ] * Transported[ 0 ] * (-4*x*y) - + Velocity[ 0 ] * Transported[ 1 ] * (2*x2 - 2*y2) - + Velocity[ 0 ] * Transported[ 2 ] * x - + Velocity[ 1 ] * Transported[ 0 ] * (2*x2 - 2*y2) - + Velocity[ 1 ] * Transported[ 1 ] * 4*x*y - + Velocity[ 1 ] * Transported[ 2 ] * y - + Velocity[ 2 ] * Transported[ 0 ] * x - + Velocity[ 2 ] * Transported[ 1 ] * y - ) / s; - } - EX transmatrix translate(hyperpoint h) { return matrix4( h[3], -h[2], h[1], h[0], @@ -1662,204 +1601,324 @@ EX namespace slr { EX hyperpoint polar(ld r, ld theta, ld phi) { return hyperpoint(sinh(r) * cos(theta-phi), sinh(r) * sin(theta-phi), cosh(r) * sin(phi), cosh(r) * cos(phi)); } - + EX hyperpoint xyz_point(ld x, ld y, ld z) { ld r = hypot(x, y); ld f = r ? sinh(r) / r : 1; return hyperpoint(x * f * cos(z) + y * f * sin(z), y * f * cos(z) - x * f * sin(z), cosh(r) * sin(z), cosh(r) * cos(z)); } - ld rootsin(ld square, ld s) { - if(square > 0) return sinh(sqrt(square) * s) / sqrt(square); - else if(square < 0) return sin(sqrt(-square) * s) / sqrt(-square); - else return s; - } - - /** it==0 is standard asin, it==1 is the next solution (PI-asin) */ - ld asin_it(ld z, int it) { - auto ans = asin(z); - if(it & 1) ans = M_PI - ans; - return ans; - } - - ld arootsin(ld square, ld v, int it) { - if(square > 0) return asinh(v * sqrt(square)) / sqrt(square); - else if(square < 0) return asin_it(v * sqrt(-square), it) / sqrt(-square); - else return v; - } - - ld roottan(ld square, ld s) { - if(square > 0) return tanh(sqrt(square) * s) / sqrt(square); - else if(square < 0) return tan(sqrt(-square) * s) / sqrt(-square); - else return s; - } - - hyperpoint geodesic_polar(ld alpha, ld beta, ld s) { - auto c = cos(2*alpha); - - ld t; - if(c > 0) - t = atan(sin(alpha) * tanh(sqrt(c) * s) / sqrt(c)); - else if(c < 0) { - /* the formula in the paper is roughly atan(k*tan(s)) - * however, atan is not always to be taken in [-PI/2,PI/2]: - * if s is in [kPI-PI/2, kPI+PI/2], we should also increase the result by kPI - */ - ld x = sqrt(-c) * s; - ld steps = floor(x/M_PI + 0.5); - t = atan(sin(alpha) * tan(sqrt(-c) * s) / sqrt(-c)) + M_PI * steps; - } - else t = atan(sin(alpha) * s); - - return polar( - asinh(cos(alpha) * rootsin(c, s)), - beta - t, - 2*sin(alpha)*s - t - ); - } - - EX hyperpoint formula_exp(hyperpoint h) { - ld s = hypot_d(3, h); - ld beta = atan2(h[1], h[0]); - ld alpha = asin(h[2] / s); - return geodesic_polar(alpha, beta, s); - } - - void find_alpha(ld phi, ld r, ld theta, ld &alpha, ld &s, ld &beta) { - if(phi < 0) { find_alpha(-phi, r, -theta, alpha, s, beta); alpha = -alpha; beta = -beta; return; } - ld mina = 0, maxa = M_PI/2; - - bool next_nan = true; - ld c; - - for(int it=0; it<40; it++) { - alpha = (mina + maxa) / 2; - - c = cos(2 * alpha); - s = arootsin(c, sinh(r) / cos(alpha), 0); - if(isnan(s)) { next_nan = true, maxa = alpha; continue; } - ld got_phi = 2*sin(alpha)*s - atan(sin(alpha) * roottan(c, s)); - if(got_phi > phi) next_nan = false, maxa = alpha; - else mina = alpha; - } - - if(next_nan) { - mina = M_PI/4; - - for(int it=0; it<40; it++) { - alpha = (mina + maxa) / 2; - c = cos(2 * alpha); - s = arootsin(c, sinh(r) / cos(alpha), 1); - ld got_phi = 2*sin(alpha)*s - atan(sin(alpha) * roottan(c, s)) - M_PI; - if(got_phi < phi) maxa = alpha; - else mina = alpha; - } - beta = theta + atan(sin(alpha) * roottan(c, s)) + M_PI; - } - else beta = theta + atan(sin(alpha) * roottan(c, s)); - } - EX hyperpoint get_inverse_exp(hyperpoint h, ld index IS(0)) { - if(sqhypot_d(2, h) < 1e-12) return point3(0, 0, atan2(h[2], h[3]) + index); - ld r = asinh(hypot_d(2, h)); + ld xy = hypot_d(2, h); ld phi = atan2(h[2], h[3]) + index; - ld theta = atan2(h[1], h[0]) + phi + index; + + bool flipped = phi > 0; + if(flipped) phi = -phi, h[2] *= -1, h[0] *= -1; - ld alpha, s, beta; - find_alpha(phi, r, theta, alpha, s, beta); + ld SV = stretch::not_squared(); + ld K = -1; - return point3(s * cos(beta) * cos(alpha), s * sin(beta) * cos(alpha), s * sin(alpha)); + ld alpha = atan2(h[1], -h[0]); + + hyperpoint res; + + ld fiber_barrier = atan(1/SV); + + ld flip_barrier = atan( 1 / tanh(asinh(xy)) / SV); + + // test the side of the flip barrier + + int part = -1; + + if(1) { + ld kk = flip_barrier; + + ld x_part = cos(kk); + ld z_part = sin(kk); + + ld rparam = x_part / z_part / SV; + + ld r = atanh(rparam); + + ld cr = cosh(r); + ld sr = sinh(r); + + // sinh(r) = xy + // r = tanh(sinh(xy)) + + + ld z = cr * (K - 1/SV/SV); + + ld k = M_PI/2; + ld a = k / K; + ld zw = xy * cr / sr; + ld u = z * a; + + ld phi1 = atan2(zw, cos(k)) - u; + + if(phi < phi1) part = 2; + } + + if(part == -1) { + ld zw = xy; + + ld u = xy * (K - 1/SV/SV) / K; + ld phi1 = atan2(zw, 1) - u; + + if(phi > phi1) part = 0; else part = 1; + } + + if(part == 2) { + ld min_k = fiber_barrier; + ld max_k = flip_barrier; + + for(int it=0; it<30; it++) { + ld kk = (min_k + max_k) / 2; + + ld x_part = cos(kk); + ld z_part = sin(kk); + + ld rparam = x_part / z_part / SV; + + assert(rparam <= 1); + + ld r = atanh(rparam); + ld cr = cosh(r); + ld sr = sinh(r); + + ld z = cr * (K - 1/SV/SV); + + ld k = M_PI - asin(xy / sr); + ld a = k / K; + ld len = a * hypot(sr, cr/SV); + ld zw = xy * cr / sr; + ld u = z * a; + + ld phi1 = atan2(zw, cos(k)) - u; + + if(phi < phi1) max_k = kk; + else min_k = kk; + + ld r_angle = alpha + u; + res = point3(cos(r_angle) * x_part * len, -sin(r_angle) * x_part * len, z_part * len); + } + } + + if(part == 0) { + ld min_k = 0; + ld max_k = fiber_barrier; + + for(int it=0; it<30; it++) { + ld kk = (min_k + max_k) / 2; + + ld x_part = cos(kk); + ld z_part = sin(kk); + + ld rparam = x_part / z_part / SV; + + ld cr = 1 / sqrt(rparam*rparam - 1); + ld sr = rparam * cr; + + ld z = cr * (K - 1/SV/SV); + + ld k = asinh(xy / sr); + ld a = k / K; + ld len = a * hypot(sr, cr/SV); + + ld zw = xy * cr / sr; + + ld u = z * a; + ld phi1 = atan2(zw, cosh(k)) - u; + + if(phi > phi1) max_k = kk; else min_k = kk; + + ld r_angle = alpha + u; + res = point3(cos(r_angle) * x_part * len, -sin(r_angle) * x_part * len, z_part * len); + } + } + + if(part == 1) { + ld min_k = fiber_barrier; + ld max_k = flip_barrier; + + for(int it=0; it<30; it++) { + ld kk = (min_k + max_k) / 2; + + ld x_part = cos(kk); + ld z_part = sin(kk); + + ld rparam = x_part / z_part / SV; + + ld r = atanh(rparam); + ld cr = cosh(r); + ld sr = sinh(r); + + ld z = cr * (K - 1/SV/SV); + + ld k = asin(xy / sr); + ld a = k / K; + ld len = a * hypot(sr, cr/SV); + ld zw = xy * cr / sr; + ld u = z * a; + + ld phi1 = atan2(zw, cos(k)) - u; + + if(isnan(phi1)) max_k = kk; + else if(phi > phi1) max_k = kk; + else min_k = kk; + + ld r_angle = alpha + u; + res = point3(cos(r_angle) * x_part * len, -sin(r_angle) * x_part * len, z_part * len); + } + } + + if(flipped) res[0] *= -1, res[2] *= -1; + + return res; } +#if ISWEB +#define ITERATE " for(int it=0; it<50; it++) { if(it >= uIterations) break; " +#else +#define ITERATE " for(int it=0; it phi) {" - " float mina = 0.;" - " float maxa = PI/4.;" - " lo_gphi = 0.; lo_s = r; lo_alpha = 0.;" -#if ISWEB - " for(int it=0; it<50; it++) { if(it >= uIterations) break; " -#else - " for(int it=0; it phi) { maxa = alpha; hi_alpha = alpha; hi_s = s; hi_gphi = gphi; }" - " else { mina = alpha; lo_alpha = alpha; lo_s = s; lo_gphi = gphi; }" - " }" - " }" - "else {" - " hi_gphi = phi; hi_s = phi; hi_alpha = 9.;" - " int next_nan = 1;" - " float mina = PI/4.;" - " float maxa = PI/2.;" -#if ISWEB - " for(int it=0; it<50; it++) { if(it >= uIterations) break; " -#else - " for(int it=0; it bound * cos(alpha)) { next_nan = 1; maxa = alpha; continue; }" - " s = asin(sinh(r) * c / cos(alpha)) / c;" - " gphi = 2.*sin(alpha)*s - atan(sin(alpha) * tan(c*s) / c);" - " if(gphi > phi) { next_nan = 0; maxa = alpha; hi_gphi = gphi; hi_s = s; hi_alpha = alpha; }" - " else { mina = alpha; lx_gphi = lo_gphi; lx_s = lo_s; lx_alpha = lo_alpha; lo_gphi = gphi; lo_s = s; lo_alpha = alpha; }" - " }" - " if(next_nan != 0) {" - " mina = PI/4.; " -#if ISWEB - " for(int it=0; it<50; it++) { if(it >= uIterations) break; " -#else - " for(int it=0; itbound) { maxa = alpha; next_nan = 1; continue; }" - " float s1 = PI - asin(z);" - " s = s1 / c;" - " gphi = 2.*sin(alpha)*s - atan(sin(alpha) * tan(s1) / c) - PI;" - " if(gphi < phi) { next_nan = 0; maxa = alpha; hi_gphi = gphi; hi_s = s; hi_alpha = alpha; }" - " else { mina = alpha; lo_gphi = gphi; lo_s = s; lo_alpha = alpha; }" - " }" - " }" - " }" - "if(hi_alpha <= 9.) { hi_gphi = lx_gphi; hi_s = lx_s; hi_alpha = lx_alpha; } " - "float fr = (phi-lo_gphi) / (hi_gphi-lo_gphi);" - "alpha = lo_alpha + (hi_alpha-lo_alpha) * fr;" - "s = lo_s + (hi_s-lo_s) * fr;" - "beta = theta - phi + 2.*sin(alpha)*s;" - "alpha = alpha * sgn; beta = beta * sgn;" - "return vec4(s * cos(beta) * cos(alpha), s * sin(beta) * cos(alpha), s * sin(alpha), 1.);" + + "vec4 res = vec4(sqrt(-1.),sqrt(-1.),sqrt(-1.),sqrt(-1.));" + + "bool flipped = phi > 0.;" + "if(flipped) phi = -phi, h[2] *= -1, h[0] *= -1;" + + "float alpha = atan2(h[1], -h[0]) + uIndexSL;" + + "float fiber_barrier = atan(1./uSV);" + + "float flip_barrier = atan(1. / tanh(asinh(xy)) / uSV);" + + "int part = 0;" + + "if(true) {" + "float x_part = cos(flip_barrier);" + "float z_part = sin(flip_barrier);" + "float rparam = x_part / z_part / uSV;" + "float r = atanh(rparam);" + "float cr = cosh(r);" + "float sr = sinh(r);" + "float z = cr * (-1.-1./uSV/uSV);" + "float k = PI/2.;" + "float a = -k;" + "float zw = xy * cr / sr;" + "float u = z * a;" + "float phi1 = atan2(zw, cos(k)) - u;" + "if(phi < phi1) part = 2;" + "}\n" + + "if(part == 0) {" + "float zw = xy;" + "float u = xy * (1. + 1./uSV/uSV);" + "float phi1 = atan2(zw, 1.) - u;" + "if(phi > phi1) part = 0; else part = 1;" + "}\n" + + "if(part == 2) {" + "float min_k = fiber_barrier;" + "float max_k = flip_barrier;" + + ITERATE + "float kk = (min_k + max_k) / 2.;" + "float x_part = cos(kk);" + "float z_part = sin(kk);" + "float rparam = x_part / z_part / uSV;" + "float r = atanh(rparam);" + "float cr = cosh(r);" + "float sr = sinh(r);" + + "float z = cr * (-1. - 1./uSV/uSV);" + "float k = PI - asin(xy / sr);" + "float a = -k;" + "float len = a * length(vec2(sr, cr/uSV));" + "float zw = xy * cr / sr;" + "float u = z * a;" + "float phi1 = atan2(zw, cos(k)) - u;" + "if(phi < phi1) max_k = kk; else min_k = kk;" + "float r_angle = alpha + u;" + "res = vec4(cos(r_angle) * x_part * len, -sin(r_angle) * x_part * len, z_part * len, 1);" + "}" + "}\n" + + "if(part == 0) {" + "float min_k = 0.;" + "float max_k = fiber_barrier;" + + ITERATE + "float kk = (min_k + max_k) / 2.;" + "float x_part = cos(kk);" + "float z_part = sin(kk);" + "float rparam = x_part / z_part / uSV;" + "float cr = 1. / sqrt(rparam*rparam - 1.);" + "float sr = rparam * cr;" + "float z = cr * (-1. - 1./uSV/uSV);" + "float k = asinh(xy / sr);" + "float a = -k;" + "float len = a * length(vec2(sr, cr/uSV));" + "float zw = xy * cr / sr;" + "float u = z * a;" + "float phi1 = atan2(zw, cosh(k)) - u;" + + "if(phi > phi1) max_k = kk; else min_k = kk;" + + "float r_angle = alpha + u;" + "res = vec4(cos(r_angle) * x_part * len, -sin(r_angle) * x_part * len, z_part * len, 1);" + "}" + "}\n" + + "if(part == 1) {" + "float min_k = fiber_barrier;" + "float max_k = flip_barrier;" + + ITERATE + "float kk = (min_k + max_k) / 2.;" + + "float x_part = cos(kk);" + "float z_part = sin(kk);" + + "float rparam = x_part / z_part / uSV;" + + "float r = atanh(rparam);" + "float cr = cosh(r);" + "float sr = sinh(r);" + + "float z = cr * (-1. - 1./uSV/uSV);" + + "float k = asin(xy / sr);" + "float a = -k;" + "float len = a * length(vec2(sr, cr/uSV));" + "float zw = xy * cr / sr;" + "float u = z * a;" + + "float phi1 = atan2(zw, cos(k)) - u;" + + "if(phi > phi1) max_k = kk;" + "else min_k = kk;" + + "float r_angle = alpha + u;" + "res = vec4(cos(r_angle) * x_part * len, -sin(r_angle) * x_part * len, z_part * len, 1);" + "}" + "}\n" + + "if(flipped) res[0] *= -1., res[2] *= -1.;" + + "return res;" "}"; EX } @@ -2058,6 +2117,63 @@ EX namespace rots { reset_projection(); current_display->set_all(0); } + /** @brief exponential function for both slr and Berger sphere */ + + EX hyperpoint formula_exp(hyperpoint vel) { + bool sp = sphere; + ld K = sp ? 1 : -1; + + ld len = hypot_d(3, vel); + + if(vel[2] < 0) len = -len; + + ld z_part = vel[2]/len; + ld x_part = sqrt(1 - z_part * z_part); + + ld SV = stretch::not_squared(); + + ld rparam = x_part / z_part / SV; + + ld beta = atan2(vel[1], vel[0]); + if(len < 0) beta += M_PI; + + if(sl2 && rparam > 1) { + ld cr = 1 / sqrt(rparam*rparam - 1); // *i + ld sr = rparam * cr; // *i + + ld z = cr * (K - 1/SV/SV); // *i + + ld a = len / hypot(sr, cr/SV); // /i + + ld k = K*a; // /i + ld u = z*a; + + ld xy = sr * sinh(k); + ld zw = cr * sinh(k); + + return hyperpoint(K*xy * cos(u+beta), K*xy * sin(u+beta), zw * cos(u) - cosh(k) * sin(u), zw * sin(u) + cosh(k)*cos(u)); + } + + else { + ld r = atan_auto(rparam); + ld cr = cos_auto(r); + ld sr = sin_auto(r); + + ld z = cr * (K - 1/SV/SV); + + ld a = len / hypot(sr, cr/SV); + + ld k = K*a; + ld u = z*a; + + ld xy = sr * sin(k); + ld zw = cr * sin(k); + + return hyperpoint(K*xy * cos(u+beta), K*xy * sin(u+beta), zw * cos(u) - cos(k) * sin(u), zw * sin(u) + cos(k)*cos(u)); + } + } + + EX } /** stretched rotation space (S3 or SLR) */ @@ -2105,15 +2221,15 @@ EX namespace stretch { return sqrt(squared()); } - hyperpoint isometric_to_actual(const hyperpoint at, const hyperpoint velocity) { + EX hyperpoint isometric_to_actual(const hyperpoint at, const hyperpoint velocity) { return mulz(at, velocity, 1/not_squared()); } - hyperpoint actual_to_isometric(const hyperpoint at, const hyperpoint velocity) { + EX hyperpoint actual_to_isometric(const hyperpoint at, const hyperpoint velocity) { return mulz(at, velocity, not_squared()); } - hyperpoint christoffel(const hyperpoint at, const hyperpoint velocity, const hyperpoint transported) { + EX hyperpoint christoffel(const hyperpoint at, const hyperpoint velocity, const hyperpoint transported) { auto vel = itranslate(at) * velocity; auto tra = itranslate(at) * transported; @@ -2147,8 +2263,7 @@ EX namespace nisot { #if CAP_SOLV else if(sn::in()) return sn::christoffel(at, velocity, transported); #endif - else if(stretch::in()) return stretch::christoffel(at, velocity, transported); - else if(sl2) return slr::christoffel(at, velocity, transported); + else if(stretch::in() || sl2) return stretch::christoffel(at, velocity, transported); else return point3(0, 0, 0); } @@ -2186,22 +2301,13 @@ EX namespace nisot { EX transmatrix parallel_transport_bare(transmatrix Pos, hyperpoint h) { - bool stretch = stretch::in(); + bool stretch = stretch::in() || sl2; h[3] = 0; auto tPos = transpose(Pos); - const ld eps = 1e-4; - - if(sl2 && !stretch) { - hyperpoint p = slr::to_phigans(tPos[3]); - for(int i=0; i<3; i++) - tPos[i] = (slr::to_phigans(tPos[3] + tPos[i] * eps) - p) / eps; - tPos[3] = p; - h = transpose(tPos) * h; - } - else h = Pos * h; + h = Pos * h; int steps = rk_steps; h /= steps; @@ -2266,13 +2372,6 @@ EX namespace nisot { for(int i=0; i<3; i++) tPos[i] = stretch::actual_to_isometric(at, tPos[i]); } - else if(sl2) { - hyperpoint p = slr::from_phigans(tPos[3]); - for(int i=0; i<3; i++) - tPos[i] = (slr::from_phigans(tPos[3] + tPos[i] * eps) - p) / eps; - tPos[3] = p; - } - return transpose(tPos); } diff --git a/raycaster.cpp b/raycaster.cpp index ba5f9d67..975e7dc0 100644 --- a/raycaster.cpp +++ b/raycaster.cpp @@ -87,7 +87,7 @@ EX bool available() { /** do we want to use the raycaster? */ EX bool requested() { if(cgflags & qRAYONLY) return true; - if(stretch::in()) return true; + if(stretch::in() && sphere) return true; if(!want_use) return false; #if CAP_TEXTURE if(texture::config.tstate == texture::tsActive) return false; diff --git a/shaders.cpp b/shaders.cpp index 4cce19ef..79add50a 100644 --- a/shaders.cpp +++ b/shaders.cpp @@ -359,6 +359,7 @@ void display_data::set_projection(int ed) { if(selected->uIterations != -1) { glhr::set_index_sl(0); + glhr::set_sv(stretch::not_squared()); glhr::set_sl_iterations(slr::steps); }