geodesics in Sol

This commit is contained in:
Zeno Rogue 2019-07-28 11:07:21 +02:00
parent 9326b9594b
commit 598603c937
15 changed files with 519 additions and 162 deletions

View File

@ -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; y<PRECX*PRECX*PRECZ; y++) xbuffer[y] = glvertex({inv_exp_table[0][0][y][0]/64, inv_exp_table[0][0][y][1]/64, inv_exp_table[0][0][y][2]/64, 1});
for(int z=0; z<PRECZ; z++)
for(int y=0; y<PRECX; y++)
for(int x=0; x<PRECX; x++)
xbuffer[z][y][x] = glvertex({inverse_exp_table[z][y][x][0], inverse_exp_table[z][y][x][1], inverse_exp_table[z][y][x][2], 1});
glTexImage3D(GL_TEXTURE_3D, 0, GL_RGBA32F, PRECX, PRECX, PRECZ, 0, GL_RGBA, GL_FLOAT, xbuffer);
delete xbuffer;
}
glUniformMatrix4fv(glhr::current->uILP, 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);

View File

@ -892,140 +892,4 @@ auto hooksw = addHook(hooks_swapdim, 100, [] {
}
namespace solv {
struct hrmap_sol : hrmap {
hrmap *binary_map;
unordered_map<pair<heptagon*, heptagon*>, heptagon*> at;
unordered_map<heptagon*, pair<heptagon*, heptagon*>> 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<heptagon> (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<eGeometry> g(geometry, gBinary4);
alt = tailored_alloc<heptagon> (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<eGeometry> g(geometry, gBinary4);
dynamicval<hrmap*> 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; i<h1->type; 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; i<S7; i++) {
// note: need do cmove before c.spin
heptagon *h1 = h->cmove(i);
dq::enqueue(h1, V * adjmatrix(i, h->c.spin(i)));
}
}
}
};
hrmap *new_map() { return new hrmap_sol; }
}
}

View File

@ -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
};

View File

@ -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"

View File

@ -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);

View File

@ -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) {

View File

@ -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];

View File

@ -448,7 +448,7 @@ void ge_select_tiling(const vector<eGeometry>& 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');

View File

@ -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; }

14
hyper.h
View File

@ -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();
}

View File

@ -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;
}

View File

@ -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<DIM; i++)
@ -1863,7 +1873,9 @@ bool limited_generation(cell *c) {
bool do_draw(cell *c, const transmatrix& T) {
if(WDIM == 3) {
if(cells_drawn > 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;
}

View File

@ -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);",

376
sol.cpp Normal file
View File

@ -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<float, 3>;
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<steps; i++) step2(at[0], at[1], at[2], v[0],v[1],v[2]);
return at;
}
void parallel(ld x, ld y, ld z, ld vx, ld vy, ld vz, ld& wx, ld& wy, ld& wz, ld t) {
// ld dax = -dz * gxyz * ay;
ld dwx = -vz * wx - vx * wz;
ld dwy = vz * wy + vy * wz;
ld dwz = vx * wx * exp(2*z) - vy * wy * exp(-2*z);
wx += dwx * t;
wy += dwy * t;
wz += dwz * t;
}
pt direct_exp(pt v, int steps, vector<pt> transported) {
ld x = 0, y = 0, z = 0;
v[0] /= steps;
v[1] /= steps;
v[2] /= steps;
for(int i=0; i<steps; i++) {
for(auto t: transported) parallel(x,y,z, v[0],v[1],v[2], t[0], t[1], t[2], 1);
step2(x,y,z, v[0],v[1],v[2]);
}
return hpxy3(x, y, z);
}
pt inverse_exp(pt h) {
load_table();
ld sx = 1, sy = 1;
bool nz = false;
ld ix = asinh(h[0]) * SXY;
ld iy = asinh(h[1]) * SXY;
ld iz = h[2] * SZ / log(2);
if(ix < 0) ix = -ix, sx = -sx;
if(iy < 0) iy = -iy, sy = -sy;
if(iz < 0) iz = -iz, nz = true, swap(ix, iy);
if(ix >= 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<pair<heptagon*, heptagon*>, heptagon*> at;
unordered_map<heptagon*, pair<heptagon*, heptagon*>> 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<heptagon> (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<eGeometry> g(geometry, gBinary4);
alt = tailored_alloc<heptagon> (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<eGeometry> g(geometry, gBinary4);
dynamicval<hrmap*> 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; i<h1->type; 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; i<S7; i++) {
// note: need do cmove before c.spin
heptagon *h1 = h->cmove(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<steps; i++) {
for(int j=0; j<3; j++)
parallel(x, y, z, shift[0], shift[1], shift[2], view_to_space[0][j], view_to_space[1][j], view_to_space[2][j], -1);
step2(x, y, z, shift[0], shift[1], shift[2]);
}
space_to_view = inverse(view_to_space);
/*
println(hlog, "arrive at = ", kz(pt({x, y, z})), " with vel = ", (shift * steps));
for(int i=0; i<3; i++)
println(hlog, "view_to_space[", i, "] = ", kz(view_to_space * units[i]));
println(hlog, "space_to_view = ", kz(space_to_view));
*/
for(int i=0; i<3; i++) for(int j=0; j<3; j++) {
print(hlog, (view_to_space*units[i] | view_to_space*units[j]), " | ");
}
println(hlog);
// println(hlog, "view_to_space = ", view_to_space);
transmatrix npush = Id;
npush[0][GDIM] = x;
npush[1][GDIM] = y;
npush[2][GDIM] = z;
// println(hlog, "result = ", space_to_view * npush * push_to);
// println(hlog, "naive = ", V * eupush(inverse(V) * T * V * C0));
return space_to_view * npush * push_to;
}
string solshader =
"uniform mediump sampler3D tInvExpTable;"
"uniform mediump mat4 uILP;"
"vec4 inverse_exp(vec4 h) {"
"float ix = asinh(h[0]) * 32.;"
"float iy = asinh(h[1]) * 32.;"
"float iz = h[2] * 8. / log(2.);"
"if(ix < 0.) ix = -ix;"
"if(iy < 0.) iy = -iy;"
"if(iz < 0.) { iz = -iz; float s = ix; ix = iy; iy = s; }"
"vec4 res = texture3D(tInvExpTable, vec3((ix + 0.5) / 64., (iy + 0.5) / 64., (iz + 0.5) / 32.));"
"if(h[2] < 0.) { res.xy = res.yx; res[2] = -res[2]; }"
"if(h[0] < 0.) res[0] = -res[0];"
"if(h[1] < 0.) res[1] = -res[1];"
"return res;"
"}";
}
}

View File

@ -1152,14 +1152,18 @@ void set_geometry(eGeometry target) {
if(binarytiling || WDIM == 3 || penrose) variation = eVariation::pure;
#endif
if(DIM == 3 && old_DIM == 2 && pmodel == mdDisk) pmodel = mdPerspective;
if(DIM == 2 && pmodel == mdPerspective) pmodel = mdDisk;
if(sol && old_DIM == 2) pmodel = mdSolPerspective;
if(DIM == 2 && among(pmodel, mdPerspective, mdSolPerspective)) pmodel = mdDisk;
if(sol && old_DIM == 2 && vid.texture_step < 4) vid.texture_step = 4;
if(sol && shmup::on) shmup::on = false, racing::on = false;
}
}
void set_variation(eVariation target) {
if(variation != target) {
stop_game();
if(euclid6 || binarytiling) geometry = gNormal;
if(euclid6 || binarytiling || sol || penrose) geometry = gNormal;
auto& cd = ginf[gCrystal];
if(target == eVariation::bitruncated && geometry == gCrystal && cd.sides == 8 && cd.vertex == 4) {
cd.vertex = 3;
@ -1250,7 +1254,7 @@ void switch_game_mode(char switchWhat) {
inv::on = false;
chaosmode = false;
princess::challenge = false;
if(bounded) set_geometry(gNormal);
if(sol || bounded) set_geometry(gNormal);
dual::disable();
break;
#endif
@ -1272,6 +1276,7 @@ void switch_game_mode(char switchWhat) {
shmup::on = !shmup::on;
princess::challenge = false;
if(!shmup::on) racing::on = false;
if(sol) set_geometry(gNormal);
break;
case rg::randpattern: