From 598603c9370061dfd19266ffafcc0825d7d36ce1 Mon Sep 17 00:00:00 2001 From: Zeno Rogue Date: Sun, 28 Jul 2019 11:07:21 +0200 Subject: [PATCH] geodesics in Sol --- basegraph.cpp | 56 ++++++- binary-tiling.cpp | 136 ----------------- classes.h | 6 +- compileunits.h | 1 + config.cpp | 4 +- conformal.cpp | 12 +- drawing.cpp | 3 +- geom-exp.cpp | 4 +- graph.cpp | 15 +- hyper.h | 14 +- hyperpoint.cpp | 4 +- hypgraph.cpp | 14 +- shaders.cpp | 25 ++- sol.cpp | 376 ++++++++++++++++++++++++++++++++++++++++++++++ system.cpp | 11 +- 15 files changed, 519 insertions(+), 162 deletions(-) create mode 100644 sol.cpp diff --git a/basegraph.cpp b/basegraph.cpp index 3fa5c5d3..3aa4d21d 100644 --- a/basegraph.cpp +++ b/basegraph.cpp @@ -233,10 +233,10 @@ void display_data::set_projection(int ed) { shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::halfplane3; if(DIM == 3 && hyperbolic && apply_models && pmodel == mdPerspective) shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardH3, pers3 = true; - if(DIM == 3 && euclid && apply_models && pmodel == mdPerspective) - shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardR3, pers3 = true; - if(DIM == 3 && sol && apply_models && pmodel == mdPerspective) + if(DIM == 3 && (euclid || sol) && apply_models && pmodel == mdPerspective) shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardR3, pers3 = true; + if(DIM == 3 && apply_models && pmodel == mdSolPerspective) + shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardSolv, pers3 = true; if(DIM == 3 && sphere && apply_models && pmodel == mdPerspective) { shaderside_projection = true; pers3 = true; int sp = spherephase & 3; @@ -250,6 +250,54 @@ void display_data::set_projection(int ed) { start_projection(ed, shaderside_projection); if(pmodel == mdRug) return; + + if(glhr::new_shader_projection == glhr::shader_projection::standardSolv) { + + static bool toload = true; + + static GLuint invexpid = 0; + + if(toload) { // if(!has_table.count(_program)) { + solv::load_table(); + if(!solv::table_loaded) { pmodel = mdPerspective; set_projection(ed); return; } + + println(hlog, "installing table"); + using namespace solv; + toload = false; + + if(invexpid == 0) glGenTextures(1, &invexpid); + + glBindTexture( GL_TEXTURE_3D, invexpid); + + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); + + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE); + // unsigned int loc = glGetUniformLocation(_program, "inv_exp_table"); + // glUniform4fv(loc, PRECX*PRECX*PRECZ, &(xbuffer[0])); + + auto xbuffer = new glvertex[PRECZ][PRECX][PRECX]; + + // for(int y=0; yuILP, 1, 0, glhr::tmtogl_transpose(solv::ilocal_perspective).as_array()); + glUniform1i(glhr::current->tInvExpTable, glhr::INVERSE_EXP_BINDING); + glActiveTexture(GL_TEXTURE0 + glhr::INVERSE_EXP_BINDING); + glBindTexture(GL_TEXTURE_3D, invexpid); + + glActiveTexture(GL_TEXTURE0 + 0); + } auto cd = current_display; @@ -284,6 +332,8 @@ void display_data::set_projection(int ed) { if(pers3) { glhr::projection_multiply(glhr::frustum(current_display->tanfov, current_display->tanfov * cd->ysize / cd->xsize)); glhr::projection_multiply(glhr::scale(1, -1, -1)); + if(glhr::new_shader_projection == glhr::shader_projection::standardSolv) + glhr::projection_multiply(glhr::tmtogl_transpose(solv::local_perspective)); } else if(DIM == 3) { glhr::glmatrix M = glhr::ortho(cd->xsize/current_display->radius/2, -cd->ysize/current_display->radius/2, 1); diff --git a/binary-tiling.cpp b/binary-tiling.cpp index acd76e7f..19954270 100644 --- a/binary-tiling.cpp +++ b/binary-tiling.cpp @@ -892,140 +892,4 @@ auto hooksw = addHook(hooks_swapdim, 100, [] { } -namespace solv { - - struct hrmap_sol : hrmap { - hrmap *binary_map; - unordered_map, heptagon*> at; - unordered_map> coords; - - heptagon *origin; - - heptagon *getOrigin() override { return origin; } - - heptagon *get_at(heptagon *x, heptagon *y) { - auto& h = at[make_pair(x, y)]; - if(h) return h; - h = tailored_alloc (S7); - h->c7 = newCell(S7, h); - coords[h] = make_pair(x, y); - h->distance = x->distance; - h->dm4 = 0; - h->zebraval = x->emeraldval; - h->emeraldval = y->emeraldval; - h->fieldval = 0; - h->cdata = NULL; - return h; - } - - hrmap_sol() { - - heptagon *alt; - - if(true) { - dynamicval g(geometry, gBinary4); - alt = tailored_alloc (S7); - alt->s = hsOrigin; - alt->alt = alt; - alt->cdata = NULL; - alt->c7 = NULL; - alt->zebraval = 0; - alt->distance = 0; - alt->emeraldval = 0; - binary_map = binary::new_alt_map(alt); - } - - origin = get_at(alt, alt); - } - - heptagon *altstep(heptagon *h, int d) { - dynamicval g(geometry, gBinary4); - dynamicval cm(currentmap, binary_map); - return h->cmove(d); - } - - heptagon *create_step(heptagon *parent, int d) override { - auto p = coords[parent]; - auto pf = p.first, ps = p.second; - auto rule = [&] (heptagon *c1, heptagon *c2, int d1) { - auto g = get_at(c1, c2); - parent->c.connect(d, g, d1, false); - return g; - }; - - switch(d) { - case 0: // right - return rule(altstep(pf, 2), ps, 4); - case 1: // up - return rule(pf, altstep(ps, 2), 5); - case 2: // front left - return rule(altstep(pf, 0), altstep(ps, 3), ps->zebraval ? 7 : 6); - case 3: // front right - return rule(altstep(pf, 1), altstep(ps, 3), ps->zebraval ? 7 : 6); - case 4: // left - return rule(altstep(pf, 4), ps, 0); - case 5: // down - return rule(pf, altstep(ps, 4), 1); - case 6: // back down - return rule(altstep(pf, 3), altstep(ps, 0), pf->zebraval ? 3 : 2); - case 7: // back up - return rule(altstep(pf, 3), altstep(ps, 1), pf->zebraval ? 3 : 2); - default: - return NULL; - } - } - - ~hrmap_sol() { - delete binary_map; - } - - transmatrix adjmatrix(int i, int j) { - ld z = log(2); - ld bw = vid.binary_width * z; - ld bwh = bw / 4; - switch(i) { - case 0: return xpush(+bw); - case 1: return ypush(+bw); - case 2: return xpush(-bwh) * zpush(+z) * ypush(j == 6 ? +bwh : -bwh); - case 3: return xpush(+bwh) * zpush(+z) * ypush(j == 6 ? +bwh : -bwh); - case 4: return xpush(-bw); - case 5: return ypush(-bw); - case 6: return ypush(-bwh) * zpush(-z) * xpush(j == 2 ? +bwh : -bwh); - case 7: return ypush(+bwh) * zpush(-z) * xpush(j == 2 ? +bwh : -bwh); - default:return Id; - } - } - - virtual transmatrix relative_matrix(heptagon *h2, heptagon *h1) override { - for(int i=0; itype; i++) if(h1->move(i) == h2) return adjmatrix(i, h1->c.spin(i)); - return Id; - } - - void draw() override { - dq::visited.clear(); - dq::enqueue(viewctr.at, cview()); - - while(!dq::drawqueue.empty()) { - auto& p = dq::drawqueue.front(); - heptagon *h = get<0>(p); - transmatrix V = get<1>(p); - dq::drawqueue.pop(); - - cell *c = h->c7; - if(!do_draw(c, V)) continue; - drawcell(c, V, 0, false); - - for(int i=0; icmove(i); - dq::enqueue(h1, V * adjmatrix(i, h->c.spin(i))); - } - } - } - - }; - - hrmap *new_map() { return new hrmap_sol; } - } - } diff --git a/classes.h b/classes.h index bf91b86f..e1727e8e 100644 --- a/classes.h +++ b/classes.h @@ -267,9 +267,9 @@ enum eModel { mdFisheye, mdJoukowsky, mdJoukowskyInverted, // 17..19 mdRotatedHyperboles, mdSpiral, mdPerspective, - // 20..22 - mdEquivolume, mdCentralInversion, mdSimulatedPerspective, mdTwoHybrid, - // 24.. + // 20..24 + mdEquivolume, mdCentralInversion, mdSimulatedPerspective, mdTwoHybrid, mdSolPerspective, + // 25.. mdGUARD, mdUnchanged, mdHyperboloidFlat, mdPolynomial, mdRug, mdFlatten }; diff --git a/compileunits.h b/compileunits.h index 31618384..7ec83b57 100644 --- a/compileunits.h +++ b/compileunits.h @@ -30,6 +30,7 @@ #include "fieldpattern.cpp" #include "heptagon.cpp" #include "binary-tiling.cpp" +#include "sol.cpp" #include "penrose.cpp" #include "archimedean.cpp" #include "euclid.cpp" diff --git a/config.cpp b/config.cpp index 993584c5..6782077b 100644 --- a/config.cpp +++ b/config.cpp @@ -464,6 +464,8 @@ void initConfig() { addsaver(bounded_mine_percentage, "bounded_mine_percentage"); + addsaver(solv::geodesic_movement, "solv_geodesic_movement", true); + #if CAP_SHMUP shmup::initConfig(); #endif @@ -1305,7 +1307,7 @@ void show3D() { if(GDIM == 2) dialog::addSelItem(XLAT("Projection at the ground level"), fts(vid.alpha), 'p'); - else if(pmodel != mdPerspective) + else if(!in_perspective()) dialog::addSelItem(XLAT("Projection distance"), fts(vid.alpha), 'p'); dialog::addBreak(50); diff --git a/conformal.cpp b/conformal.cpp index 3761c343..f2ebef65 100644 --- a/conformal.cpp +++ b/conformal.cpp @@ -630,6 +630,8 @@ namespace conformal { #endif bool model_available(eModel pm) { + if(sol) return among(pm, mdDisk, mdPerspective, mdSolPerspective); + if(pm == mdSolPerspective && !sol) return false; if(sphere && (pm == mdHalfplane || pm == mdBall)) return false; if(DIM == 2 && pm == mdPerspective) return false; @@ -652,6 +654,11 @@ namespace conformal { string get_model_name(eModel m) { if(m == mdDisk && DIM == 3 && hyperbolic) return XLAT("ball model/Gans"); + if(sol) { + if(m == mdDisk) return XLAT("simple model: projection"); + if(m == mdPerspective) return XLAT("simple model: perspective"); + if(m == mdSolPerspective) return XLAT("native perspective"); + } if(m == mdDisk && DIM == 3) return XLAT("perspective in 4D"); if(m == mdHalfplane && DIM == 3 && hyperbolic) return XLAT("half-space"); if(sphere) @@ -770,6 +777,9 @@ namespace conformal { } dialog::addBreak(100); + + if(sol) + dialog::addBoolItem_action(XLAT("geodesic movement in Sol"), solv::geodesic_movement, 'G'); dialog::addBoolItem(XLAT("rotation"), do_rotate == 2, 'r'); if(do_rotate == 0) dialog::lastItem().value = XLAT("NEVER"); @@ -780,7 +790,7 @@ namespace conformal { dialog::add_action([] { edit_rotation(conformal::rotation); }); // if(pmodel == mdBand && sphere) - if(pmodel != mdPerspective) + if(!in_perspective()) dialog::addSelItem(XLAT("scale factor"), fts(vid.scale), 'z'); if(abs(vid.alpha-1) > 1e-3 && pmodel != mdBall && pmodel != mdHyperboloid && pmodel != mdHemisphere && pmodel != mdDisk) { diff --git a/drawing.cpp b/drawing.cpp index 7ebf0fd2..1a907b3e 100644 --- a/drawing.cpp +++ b/drawing.cpp @@ -382,7 +382,8 @@ void drawTexturedTriangle(SDL_Surface *s, int *px, int *py, glvertex *tv, color_ void glapplymatrix(const transmatrix& V) { GLfloat mat[16]; int id = 0; - if(pmodel == mdPerspective && DIM == 3) { + + if(in_perspective() && DIM == 3) { if(spherephase & 4) { for(int y=0; y<4; y++) { for(int x=0; x<4; x++) mat[id++] = -V[x][y]; diff --git a/geom-exp.cpp b/geom-exp.cpp index 1ca0bf77..b8ee5e4e 100644 --- a/geom-exp.cpp +++ b/geom-exp.cpp @@ -448,7 +448,7 @@ void ge_select_tiling(const vector& lst) { } string current_proj_name() { - if(pmodel != mdDisk) + if(pmodel != mdDisk || sol) return conformal::get_model_name(pmodel); else if(hyperbolic && vid.alpha == 1) return XLAT("Poincaré model"); @@ -724,6 +724,8 @@ void showEuclideanMenu() { dialog::add_action_push(show3D); } dialog::addSelItem(XLAT("projection"), current_proj_name(), '1'); + if(sol) + dialog::addBoolItem_action(XLAT("geodesic movement in Sol"), solv::geodesic_movement, 'G'); dialog::add_action_push(conformal::model_menu); #if CAP_CRYSTAL && MAXMDIM >= 4 crystal::add_crystal_transform('x'); diff --git a/graph.cpp b/graph.cpp index db3d88d4..620acf0b 100644 --- a/graph.cpp +++ b/graph.cpp @@ -18,8 +18,12 @@ int detaillevel = 0; bool first_cell_to_draw = true; +bool in_perspective() { + return among(pmodel, mdPerspective, mdSolPerspective); + } + bool hide_player() { - return DIM == 3 && playermoved && vid.yshift == 0 && vid.sspeed > -5 && pmodel == mdPerspective && (first_cell_to_draw || elliptic) && (WDIM == 3 || vid.camera == 0) && !inmirrorcount + return DIM == 3 && playermoved && vid.yshift == 0 && vid.sspeed > -5 && in_perspective() && (first_cell_to_draw || elliptic) && (WDIM == 3 || vid.camera == 0) && !inmirrorcount #if CAP_RACING && !(racing::on && !racing::standard_centering && !racing::player_relative) #endif @@ -809,6 +813,7 @@ color_t kind_outline(eItem it) { transmatrix face_the_player(const transmatrix V) { if(DIM == 2) return V; + if(sol) return V * cspin(0, 2, ptick(618, 0)); return rgpushxto0(tC0(V)); } @@ -863,7 +868,7 @@ bool drawItemType(eItem it, cell *c, const transmatrix& V, color_t icol, int pti if(DIM == 3 && mapeditor::drawUserShape(V, mapeditor::sgItem, it, darkena(icol, 0, 0xFF), c)) return false; - if(WDIM == 3 && c == viewctr.at->c7 && pmodel == mdPerspective && hdist0(tC0(V)) < cgi.orbsize * 0.25) return false; + if(WDIM == 3 && c == viewctr.at->c7 && in_perspective() && hdist0(tC0(V)) < cgi.orbsize * 0.25) return false; transmatrix Vit = V; if(GDIM == 3 && WDIM == 2 && c && it != itBabyTortoise) Vit = mscale(V, cgi.STUFF); @@ -4586,7 +4591,7 @@ void draw_grid_at(cell *c, const transmatrix& V) { void queue_transparent_wall(const transmatrix& V, hpcshape& sh, color_t color) { auto& poly = queuepolyat(V, sh, color, PPR::TRANSPARENT_WALL); hyperpoint h = V * sh.intester; - if(pmodel == mdPerspective) + if(in_perspective()) poly.subprio = int(hdist0(h) * 100000); else { hyperpoint h2; @@ -7406,7 +7411,7 @@ void calcparam() { current_display->sidescreen = false; - if(vid.xres < vid.yres - 2 * vid.fsize && !inHighQual && pmodel != mdPerspective) { + if(vid.xres < vid.yres - 2 * vid.fsize && !inHighQual && !in_perspective()) { cd->ycenter = vid.yres - cd->scrsize - vid.fsize; } else { @@ -7421,7 +7426,7 @@ void calcparam() { } cd->radius = vid.scale * cd->scrsize; - if(DIM == 3 && pmodel == mdPerspective) cd->radius = cd->scrsize; + if(DIM == 3 && in_perspective()) cd->radius = cd->scrsize; realradius = min(realradius, cd->radius); if(dronemode) { cd->ycenter -= cd->radius; cd->ycenter += vid.fsize/2; cd->ycenter += vid.fsize/2; cd->radius *= 2; } diff --git a/hyper.h b/hyper.h index cf32aacf..df76f879 100644 --- a/hyper.h +++ b/hyper.h @@ -4023,7 +4023,7 @@ namespace glhr { enum class shader_projection { standard, band, halfplane, standardH3, standardR3, standardS30, standardS31, standardS32, standardS33, - ball, halfplane3, band3, flatten, + ball, halfplane3, band3, flatten, standardSolv, MAX }; @@ -5625,4 +5625,16 @@ namespace kite { hyperpoint get_corner(cell *c, int d, ld cf); } #endif + +namespace solv { + extern transmatrix local_perspective, ilocal_perspective; + hrmap *new_map(); + bool inverse_exp(hyperpoint h, hyperpoint& res); + + transmatrix get_solmul(const transmatrix T, const transmatrix V); + extern string solshader; + } + +bool in_perspective(); + } \ No newline at end of file diff --git a/hyperpoint.cpp b/hyperpoint.cpp index 0b380099..e5f158b4 100644 --- a/hyperpoint.cpp +++ b/hyperpoint.cpp @@ -791,8 +791,8 @@ bool asign(ld y1, ld y2) { return signum(y1) != signum(y2); } ld xcross(ld x1, ld y1, ld x2, ld y2) { return x1 + (x2 - x1) * y1 / (y1 - y2); } -transmatrix solmul(transmatrix T, transmatrix V) { - if(sol) return V * eupush(inverse(V) * T * V * C0); +transmatrix solmul(const transmatrix T, const transmatrix V) { + if(sol) return solv::get_solmul(T, V); else return T * V; } diff --git a/hypgraph.cpp b/hypgraph.cpp index 079609af..6f75552d 100644 --- a/hypgraph.cpp +++ b/hypgraph.cpp @@ -318,6 +318,15 @@ void applymodel(hyperpoint H, hyperpoint& ret) { ret[2] = 1; return; } + + case mdSolPerspective: { + auto S = solv::inverse_exp(H); + ld ratio = vid.xres / current_display->tanfov / current_display->radius / 2; + ret[0] = S[0]/S[2] * ratio; + ret[1] = S[1]/S[2] * ratio; + ret[2] = 1; + return; + } case mdUnchanged: case mdFlatten: @@ -897,6 +906,7 @@ bool invalid_point(const transmatrix T) { bool in_smart_range(const transmatrix& T) { if(invalid_point(T)) return false; + if(pmodel == mdSolPerspective) return solv::in_table_range(solv::ilocal_perspective * tC0(T)); hyperpoint h1; applymodel(tC0(T), h1); for(int i=0; i vid.cells_drawn_limit) return false; - if(vid.use_smart_range) { + if(sol) + return solv::in_table_range(solv::ilocal_perspective * tC0(T)); + else if(vid.use_smart_range) { if(cells_drawn >= 50 && !in_smart_range(T)) return false; if(!limited_generation(c)) return false; } diff --git a/shaders.cpp b/shaders.cpp index 7affa097..7794bb88 100644 --- a/shaders.cpp +++ b/shaders.cpp @@ -107,6 +107,14 @@ glmatrix tmtogl(const transmatrix& T) { return tmp; } +glmatrix tmtogl_transpose(const transmatrix& T) { + glmatrix tmp; + for(int i=0; i<4; i++) + for(int j=0; j<4; j++) + tmp[i][j] = T[j][i]; + return tmp; + } + glmatrix ortho(ld x, ld y, ld z) { return scale(1/x, 1/y, 1/z); } @@ -204,11 +212,13 @@ static const int aPosition = 0; static const int aColor = 3; static const int aTexture = 8; +const int INVERSE_EXP_BINDING = 2; + struct GLprogram { GLuint _program; GLuint vertShader, fragShader; - GLint uMVP, uFog, uFogColor, uColor, tTexture, uMV, uProjection, uAlpha, uFogBase; + GLint uMVP, uFog, uFogColor, uColor, tTexture, tInvExpTable, uMV, uProjection, uAlpha, uFogBase, uILP; GLprogram(string vsh, string fsh) { _program = glCreateProgram(); @@ -263,8 +273,10 @@ struct GLprogram { uFogBase = glGetUniformLocation(_program, "uFogBase"); uAlpha = glGetUniformLocation(_program, "uAlpha"); uColor = glGetUniformLocation(_program, "uColor"); + uILP = glGetUniformLocation(_program, "uILP"); tTexture = glGetUniformLocation(_program, "tTexture"); - + tInvExpTable = glGetUniformLocation(_program, "tInvExpTable"); + #if DEBUG_GL printf("uniforms: %d %d %d %d\n", uMVP, uFog, uColor, tTexture); #endif @@ -546,6 +558,7 @@ void init() { bool band = among(sp, shader_projection::band, shader_projection::band3); bool hp = among(sp, shader_projection::halfplane, shader_projection::halfplane3); bool sh3 = (sp == shader_projection::standardH3); + bool ssol = (sp == shader_projection::standardSolv); bool sr3 = (sp == shader_projection::standardR3); bool ss30 = (sp == shader_projection::standardS30); bool ss31 = (sp == shader_projection::standardS31); @@ -553,7 +566,7 @@ void init() { bool ss33 = (sp == shader_projection::standardS33); bool ss3 = ss30 || ss31 || ss32 || ss33; - bool s3 = (sh3 || sr3 || ss3); + bool s3 = (sh3 || sr3 || ss3 || ssol); bool dim3 = s3 || among(sp, shader_projection::ball, shader_projection::halfplane3, shader_projection::band3); bool dim2 = !dim3; bool ball = (sp == shader_projection::ball); @@ -601,6 +614,8 @@ void init() { 1, "float zlevel(vec4 h) {", 1, " return (h[2] < 0.0 ? -1.0 : 1.0) * sqrt(h[2]*h[2] - h[0]*h[0] - h[1]*h[1]);", 1, " }", + + ssol, solv::solshader, 1, "void main() {", texture, "vTexCoord = aTexture;", @@ -635,8 +650,10 @@ void init() { hp && dim3, "t.x /= -rads; t.y /= -rads; t.z /= -rads; t[3] = 1.0;", s3, "vec4 t = uMV * aPosition;", + ssol, "t = inverse_exp(uILP * t);", sh3, "float fogs = (uFogBase - acosh(t[3]) / uFog);", sr3, "float fogs = (uFogBase - sqrt(t[0]*t[0] + t[1]*t[1] + t[2]*t[2]) / uFog);", + ssol, "float fogs = (uFogBase - sqrt(t[0]*t[0] + t[1]*t[1] + t[2]*t[2]) / uFog);", ss30, "float fogs = (uFogBase - (6.284 - acos(t[3])) / uFog); t = -t; ", ss31, "float fogs = (uFogBase - (6.284 - acos(t[3])) / uFog); t.xyz = -t.xyz; ", @@ -645,7 +662,7 @@ void init() { s3, "vColor.xyz = vColor.xyz * fogs + uFogColor.xyz * (1.0-fogs);", - sh3 || sr3 || ball,"t[3] = 1.0;", + sh3 || sr3 || ssol || ball,"t[3] = 1.0;", band || hp || s3 || ball,"gl_Position = uP * t;", dim3 && !s3, "vColor.xyz = vColor.xyz * (0.5 - gl_Position.z / 2.0) + uFogColor.xyz * (0.5 + gl_Position.z / 2.0);", diff --git a/sol.cpp b/sol.cpp new file mode 100644 index 00000000..ca148704 --- /dev/null +++ b/sol.cpp @@ -0,0 +1,376 @@ + +// implementation of the Solv space + +// Copyright (C) 2011-2019 Zeno Rogue, see 'hyper.cpp' for details + +#include "extra/solv/common.cpp" + +namespace hr { + +namespace solv { + + const int PRECX = 64; + const int PRECZ = 32; + const ld SXY = 32.; + const ld SZ = 8.; + + using pt = hyperpoint; + using ptlow = array; + + ptlow inverse_exp_table[PRECZ][PRECX][PRECX]; + + bool table_loaded; + + void load_table() { + if(table_loaded) return; + FILE *f = fopen("soltable.dat", "rb"); + if(!f) f = fopen("/usr/lib/soltable.dat", "rb"); + if(!f) { addMessage(XLAT("geodesic table missing")); pmodel = mdPerspective; return; } + fread(inverse_exp_table, sizeof(inverse_exp_table), 1, f); + fclose(f); + table_loaded = true; + } + + void step2(ld& x, ld &y, ld &z, ld &vx, ld &vy, ld &vz) { + + ld ax = -2 * vx * vz; + ld ay = 2 * vy * vz; + ld az = exp(2*z) * vx * vx - exp(-2*z) * vy * vy; + + // ld x2 = x + vx/2; + // ld y2 = y + vy/2; + ld z2 = z + vz/2; + + ld vx2 = vx + ax/2; + ld vy2 = vy + ay/2; + ld vz2 = vz + az/2; + + ld ax2 = -2 * vx2 * vz2; + ld ay2 = 2 * vy2 * vz2; + ld az2 = exp(2*z2) * vx2 * vx2 - exp(-2*z2) * vy2 * vy2; + + x += vx + ax/2; + y += vy + ay/2; + z += vz + az/2; + + vx += ax2; + vy += ay2; + vz += az2; + } + + pt direct_exp(pt v, int steps) { + pt at = hpxy(0, 0); + v[0] /= steps; + v[1] /= steps; + v[2] /= steps; + for(int i=0; i transported) { + ld x = 0, y = 0, z = 0; + v[0] /= steps; + v[1] /= steps; + v[2] /= steps; + for(int i=0; i= PRECX-1) ix = PRECX-2; + if(iy >= PRECX-1) iy = PRECX-2; + if(iz >= PRECZ-1) iz = PRECZ-2; + + int ax = ix, bx = ax+1; + int ay = iy, by = ay+1; + int az = iz, bz = az+1; + + pt res; + + // println(hlog, "inv_table", make_tuple(iz,iy,ix), " = ", make_tuple(inverse_exp_table[z][y][x][0], inverse_exp_table[z][y][x][1], inverse_exp_table[z][y][x][2])); + + #define S0(x,y,z) inverse_exp_table[z][y][x][t] + #define S1(x,y) (S0(x,y,az) * (bz-iz) + S0(x,y,bz) * (iz-az)) + #define S2(x) (S1(x,ay) * (by-iy) + S1(x,by) * (iy-ay)) + + for(int t=0; t<3; t++) + res[t] = S2(ax) * (bx-ix) + S2(bx) * (ix-ax); + if(nz) swap(res[0], res[1]), res[2] = -res[2]; + res[0] *= sx; res[1] *= sy; + res[3] = 1; + + #undef S0 + #undef S1 + #undef S2 + + // println(hlog, kz(h), " => ", kz(res), " => ", kz(direct_exp(res, 1000)), " [", ix, ",", iy, ",", iz, " | ", sx, "/", sy, "/", nz, "]"); + + return res; + } + + transmatrix local_perspective, ilocal_perspective; + + bool geodesic_movement = true; + + struct hrmap_sol : hrmap { + hrmap *binary_map; + unordered_map, heptagon*> at; + unordered_map> coords; + + heptagon *origin; + + heptagon *getOrigin() override { return origin; } + + heptagon *get_at(heptagon *x, heptagon *y) { + auto& h = at[make_pair(x, y)]; + if(h) return h; + h = tailored_alloc (S7); + h->c7 = newCell(S7, h); + coords[h] = make_pair(x, y); + h->distance = x->distance; + h->dm4 = 0; + h->zebraval = x->emeraldval; + h->emeraldval = y->emeraldval; + h->fieldval = 0; + h->cdata = NULL; + return h; + } + + hrmap_sol() { + + heptagon *alt; + + if(true) { + dynamicval g(geometry, gBinary4); + alt = tailored_alloc (S7); + alt->s = hsOrigin; + alt->alt = alt; + alt->cdata = NULL; + alt->c7 = NULL; + alt->zebraval = 0; + alt->distance = 0; + alt->emeraldval = 0; + binary_map = binary::new_alt_map(alt); + } + + origin = get_at(alt, alt); + } + + heptagon *altstep(heptagon *h, int d) { + dynamicval g(geometry, gBinary4); + dynamicval cm(currentmap, binary_map); + return h->cmove(d); + } + + heptagon *create_step(heptagon *parent, int d) override { + auto p = coords[parent]; + auto pf = p.first, ps = p.second; + auto rule = [&] (heptagon *c1, heptagon *c2, int d1) { + auto g = get_at(c1, c2); + parent->c.connect(d, g, d1, false); + return g; + }; + + switch(d) { + case 0: // right + return rule(altstep(pf, 2), ps, 4); + case 1: // up + return rule(pf, altstep(ps, 2), 5); + case 2: // front left + return rule(altstep(pf, 0), altstep(ps, 3), ps->zebraval ? 7 : 6); + case 3: // front right + return rule(altstep(pf, 1), altstep(ps, 3), ps->zebraval ? 7 : 6); + case 4: // left + return rule(altstep(pf, 4), ps, 0); + case 5: // down + return rule(pf, altstep(ps, 4), 1); + case 6: // back down + return rule(altstep(pf, 3), altstep(ps, 0), pf->zebraval ? 3 : 2); + case 7: // back up + return rule(altstep(pf, 3), altstep(ps, 1), pf->zebraval ? 3 : 2); + default: + return NULL; + } + } + + ~hrmap_sol() { + delete binary_map; + } + + transmatrix adjmatrix(int i, int j) { + ld z = log(2); + ld bw = vid.binary_width * z; + ld bwh = bw / 4; + switch(i) { + case 0: return xpush(+bw); + case 1: return ypush(+bw); + case 2: return xpush(-bwh) * zpush(+z) * ypush(j == 6 ? +bwh : -bwh); + case 3: return xpush(+bwh) * zpush(+z) * ypush(j == 6 ? +bwh : -bwh); + case 4: return xpush(-bw); + case 5: return ypush(-bw); + case 6: return ypush(-bwh) * zpush(-z) * xpush(j == 2 ? +bwh : -bwh); + case 7: return ypush(+bwh) * zpush(-z) * xpush(j == 2 ? +bwh : -bwh); + default:return Id; + } + } + + virtual transmatrix relative_matrix(heptagon *h2, heptagon *h1) override { + for(int i=0; itype; i++) if(h1->move(i) == h2) return adjmatrix(i, h1->c.spin(i)); + if(gmatrix0.count(h2->c7) && gmatrix0.count(h1->c7)) + return inverse(gmatrix0[h1->c7]) * gmatrix0[h2->c7]; + return Id; // not implemented yet + } + + void draw() override { + dq::visited.clear(); + + transmatrix T = eupush( tC0(inverse(View)) ); + local_perspective = View * T; + ilocal_perspective = inverse(local_perspective); + + dq::enqueue(viewctr.at, cview()); + + while(!dq::drawqueue.empty()) { + auto& p = dq::drawqueue.front(); + heptagon *h = get<0>(p); + transmatrix V = get<1>(p); + dq::drawqueue.pop(); + + cell *c = h->c7; + if(!do_draw(c, V)) continue; + drawcell(c, V, 0, false); + + for(int i=0; icmove(i); + dq::enqueue(h1, V * adjmatrix(i, h->c.spin(i))); + } + } + } + + }; + + hrmap *new_map() { return new hrmap_sol; } + + bool in_table_range(hyperpoint h) { + ld ix = asinh(h[0]) * SXY; + ld iy = asinh(h[1]) * SXY; + ld iz = h[2] * SZ / log(2); + return abs(ix) < PRECX && abs(iy) < PRECX && abs(iz) < PRECZ; + } + + transmatrix get_solmul(const transmatrix T, const transmatrix V) { + if(!geodesic_movement) return V * eupush(inverse(V) * T * V * C0); + + transmatrix push_back = eupush( tC0(inverse(V)) ); + transmatrix space_to_view = V * push_back; + + transmatrix view_to_space = inverse(space_to_view); + using namespace hyperpoint_vec; + hyperpoint shift = /* inverse(V) * T * V * C0; */ view_to_space * T * C0; + + transmatrix push_to = inverse(push_back); + + int steps = 100; + + hyperpoint units[3] = { point3(1,0,0), point3(0,1,0), point3(0,0,1) }; + + /* + println(hlog, "shift = ", kz(shift)); + println(hlog, "space_to_view = ", kz(space_to_view)); + + for(int i=0; i<3; i++) + println(hlog, "view_to_space[", i, "] = ", kz(view_to_space * units[i])); + */ + + for(int s=0; s<3; s++) shift[s] /= steps; + ld x = 0, y = 0, z = 0; + + + for(int i=0; i