1 // Hyperbolic Rogue -- binary tilings
2 // Copyright (C) 2011-2019 Zeno Rogue, see 'hyper.cpp' for details
3
4 /** \file binary-tiling.cpp
5 * \brief Binary tilings in 2D and 3D
6 */
7
8 #include "hyper.h"
9 namespace hr {
10
11 EX namespace bt {
12
13 /** note: nihsolv and kd3 tilings return bt::in(). They are defined elsewhere, although some of bt:: functions are used for them */
in()14 EX bool in() {
15 #if CAP_BT
16 return cgflags & qBINARY;
17 #else
18 return false;
19 #endif
20 }
21
22 #if !CAP_BT
updir()23 EX int updir() { return 0; }
24 #endif
25
26 #if CAP_BT
27 #if HDR
28 enum bindir {
29 bd_right = 0,
30 bd_up_right = 1,
31 bd_up = 2,
32 bd_up_left = 3,
33 bd_left = 4,
34 bd_down = 5, /* for cells of degree 6 */
35 bd_down_left = 5, /* for cells of degree 7 */
36 bd_down_right = 6 /* for cells of degree 7 */
37 };
38 #endif
39
type_of(heptagon * h)40 EX int type_of(heptagon *h) {
41 return h->c7->type;
42 }
43
44 // 0 - central, -1 - left, +1 - right
mapside(heptagon * h)45 EX int mapside(heptagon *h) {
46 return h->zebraval;
47 }
48
49 #if DEBUG_BINARY_TILING
50 map<heptagon*, long long> xcode;
51 map<long long, heptagon*> rxcode;
52
expected_xcode(heptagon * h,int d)53 long long expected_xcode(heptagon *h, int d) {
54 auto r =xcode[h];
55 if(d == 0) return r + 1;
56 if(d == 1) return 2*r + 1;
57 if(d == 2) return 2*r;
58 if(d == 3) return 2*r - 1;
59 if(d == 4) return r-1;
60 if(d == 5 && type_of(h) == 6) return r / 2;
61 if(d == 5 && type_of(h) == 7) return (r-1) / 2;
62 if(d == 6 && type_of(h) == 7) return (r+1) / 2;
63 breakhere();
64 }
65 #endif
66
path(heptagon * h,int d,int d1,std::initializer_list<int> p)67 EX heptagon *path(heptagon *h, int d, int d1, std::initializer_list<int> p) {
68 static int rec = 0;
69 rec++; if(rec>100) exit(1);
70 // printf("{generating path from %p (%d/%d) dir %d:", h, type_of(h), mapside(h), d);
71 heptagon *h1 = h;
72 for(int d0: p) {
73 // printf(" [%d]", d0);
74 h1 = currentmap->may_create_step(h1, d0);
75 // printf(" %p", h1);
76 }
77
78 #if DEBUG_BINARY_TILING
79 if(xcode[h1] != expected_xcode(h, d)) {
80 printf("expected_xcode mismatch\n");
81 breakhere();
82 }
83 #endif
84 // printf("}\n");
85 if(h->move(d) && h->move(d) != h1) {
86 printf("already connected to something else (1)\n");
87 breakhere();
88 }
89 if(h1->move(d1) && h1->move(d1) != h) {
90 printf("already connected to something else (2)\n");
91 breakhere();
92 }
93 h->c.connect(d, h1, d1, false);
94 rec--;
95 return h1;
96 }
97
pathc(heptagon * h,int d,int d1,std::vector<std::initializer_list<int>> p)98 EX heptagon *pathc(heptagon *h, int d, int d1, std::vector<std::initializer_list<int>> p) {
99 h->cmove(S7-1);
100 int z = h->c.spin(S7-1);
101 return path(h, d, d1, p[z]);
102 }
103
104 EX ld hororec_scale = 0.25;
105 EX ld horohex_scale = 0.6;
106
make_binary_lands(heptagon * parent,heptagon * h)107 EX void make_binary_lands(heptagon *parent, heptagon *h) {
108 if(!parent->emeraldval) parent->emeraldval = currentmap->gamestart()->land;
109 eLand z = eLand(parent->emeraldval);
110 int chance = 0;
111 if(ls::no_walls() || parent->emeraldval == laCrossroads4) {
112 eLand x = parent->c7->land;
113 parent->c7->land = z;
114 chance = wallchance(parent->c7, deep_ocean_at(parent->c7, parent->c7));
115 parent->c7->land = x;
116 }
117 if(ls::std_chaos()) chance = 1000;
118 if(chance && hrand(40000) < chance)
119 h->emeraldval = getNewLand(z);
120 else
121 h->emeraldval = z;
122 }
123
build(heptagon * parent,int d,int d1,int t,int side,int delta)124 EX heptagon *build(heptagon *parent, int d, int d1, int t, int side, int delta) {
125 auto h = buildHeptagon1(init_heptagon(t), parent, d, hsA, d1);
126 h->distance = parent->distance + delta;
127 h->dm4 = parent->dm4 + delta;
128 h->c7 = NULL;
129 if(parent->c7) h->c7 = newCell(t, h);
130 h->zebraval = side;
131 switch(geometry) {
132 case gBinary4:
133 if(d < 2)
134 h->emeraldval = gmod(parent->emeraldval * 2 + d, 15015);
135 else
136 h->emeraldval = gmod((parent->emeraldval - d1) * 7508, 15015);
137 break;
138 case gTernary:
139 if(d < 2)
140 h->emeraldval = gmod(parent->emeraldval * 3 + d, 10010);
141 else
142 h->emeraldval = gmod((parent->emeraldval - d1) * 3337, 10010);
143 break;
144 case gHoroRec: {
145 int x = parent->fieldval & 4095;
146 int y = (parent->fieldval >> 12) & 4095;
147 if(d < 2) tie(x, y) = make_pair(y, gmod(x * 2 + d, 1155));
148 else tie(x,y) = make_pair(gmod((y-d1)*578, 1155), x);
149 h->fieldval = x + (y << 12);
150 break;
151 }
152 case gBinary3: {
153 int x = parent->fieldval & 4095;
154 int y = (parent->fieldval >> 12) & 4095;
155 if(d < 4) x = gmod(x * 2 + (d&1), 1155), y = gmod(y * 2 + (d>>1), 1155);
156 else x = gmod((x-(d1&1))*578, 1155), y = gmod((y-(d1>>1))*578, 1155);
157 h->fieldval = x + (y << 12);
158 break;
159 }
160 default:
161 break;
162 }
163 if(WDIM == 3 && h->c7) make_binary_lands(parent, h);
164 #if DEBUG_BINARY_TILING
165 xcode[h] = expected_xcode(parent, d);
166 if(rxcode.count(xcode[h])) {
167 printf("xcode clash\n");
168 breakhere();
169 }
170 rxcode[xcode[h]] = h;
171 #endif
172 return h;
173 }
174
175 #if MAXMDIM==4
build3(heptagon * parent,int d,int d1,int delta)176 EX heptagon *build3(heptagon *parent, int d, int d1, int delta) {
177 int side = 0;
178 if(geometry == gBinary3) {
179 if(d < 4) side = (parent->zebraval * 2 + d) % 5;
180 if(d == S7-1) side = ((5+parent->zebraval-d1) * 3) % 5;
181 }
182 if(geometry == gHoroHex) {
183 if(d < 3) side = (parent->zebraval + d) % 3;
184 if(d == S7-1) side = (parent->zebraval + 3 - d1) % 3;
185 }
186 return build(parent, d, d1, S7, side, delta);
187 }
188 #endif
189
190 struct hrmap_binary : hrmap {
191
192 heptagon *origin;
193 std::mt19937 directions_generator;
194
hrmap_binaryhr::bt::hrmap_binary195 hrmap_binary(heptagon *o) : origin(o) { set_seed(); }
196
set_seedhr::bt::hrmap_binary197 void set_seed() { directions_generator.seed(137137137); }
198
nextdirhr::bt::hrmap_binary199 int nextdir(int choices) { return directions_generator() % choices; }
200
getOriginhr::bt::hrmap_binary201 heptagon *getOrigin() override { return origin; }
202
hrmap_binaryhr::bt::hrmap_binary203 hrmap_binary() {
204 set_seed();
205 origin = hyperbolic_origin();
206 #if DEBUG_BINARY_TILING
207 bt::xcode.clear();
208 bt::rxcode.clear();
209 bt::xcode[&h] = (1 << 16);
210 bt::rxcode[1<<16] = &h;
211 #endif
212 origin->zebraval = 0;
213 origin->emeraldval = 0;
214 }
215
create_stephr::bt::hrmap_binary216 heptagon *create_step(heptagon *parent, int d) override {
217 auto h = parent;
218 switch(geometry) {
219 case gBinaryTiling: {
220 switch(d) {
221 case bd_right: {
222 if(mapside(h) > 0 && type_of(h) == 7)
223 return path(h, d, bd_left, {bd_left, bd_down, bd_right, bd_up});
224 else if(mapside(h) >= 0)
225 return build(parent, bd_right, bd_left, type_of(parent) ^ 1, 1, 0);
226 else if(type_of(h) == 6)
227 return path(h, d, bd_left, {bd_down, bd_right, bd_up, bd_left});
228 else
229 return path(h, d, bd_left, {bd_down_right, bd_up});
230 }
231 case bd_left: {
232 if(mapside(h) < 0 && type_of(h) == 7)
233 return path(h, d, bd_right, {bd_right, bd_down, bd_left, bd_up});
234 else if(mapside(h) <= 0)
235 return build(parent, bd_left, bd_right, type_of(parent) ^ 1, -1, 0);
236 else if(type_of(h) == 6)
237 return path(h, d, bd_right, {bd_down, bd_left, bd_up, bd_right});
238 else
239 return path(h, d, bd_right, {bd_down_left, bd_up});
240 }
241 case bd_up_right: {
242 return path(h, d, bd_down_left, {bd_up, bd_right});
243 }
244 case bd_up_left: {
245 return path(h, d, bd_down_right, {bd_up, bd_left});
246 }
247 case bd_up:
248 return build(parent, bd_up, bd_down, 6, mapside(parent), 1);
249 default:
250 /* bd_down */
251 if(type_of(h) == 6) {
252 if(mapside(h) == 0)
253 return build(parent, bd_down, bd_up, 6, 0, -1);
254 else if(mapside(h) == 1)
255 return path(h, d, bd_up, {bd_left, bd_left, bd_down, bd_right});
256 else if(mapside(h) == -1)
257 return path(h, d, bd_up, {bd_right, bd_right, bd_down, bd_left});
258 }
259 /* bd_down_left */
260 else if(d == bd_down_left) {
261 return path(h, d, bd_up_right, {bd_left, bd_down});
262 }
263 else if(d == bd_down_right) {
264 return path(h, d, bd_up_left, {bd_right, bd_down});
265 }
266 }
267 throw hr_exception("wrong dir");
268 }
269 case gBinary4: {
270 switch(d) {
271 case 0: case 1:
272 return build(parent, d, 3, 5, d, 1);
273 case 3:
274 return build(parent, 3, parent->zebraval, 5, nextdir(2), -1);
275 case 2:
276 if(parent->zebraval == 0)
277 return path(h, 2, 4, {3, 1});
278 else
279 return path(h, 2, 4, {3, 2, 0});
280 case 4:
281 if(parent->zebraval == 1)
282 return path(h, 4, 2, {3, 0});
283 else
284 return path(h, 4, 2, {3, 4, 1});
285 default:
286 throw hr_exception("wrong dir");
287 }
288 }
289 case gTernary: {
290 switch(d) {
291 case 0: case 1: case 2:
292 return build(parent, d, 4, 6, d, 1);
293 case 4:
294 return build(parent, 4, parent->zebraval, 6, nextdir(3), -1);
295 case 3:
296 if(parent->zebraval < 2)
297 return path(h, 3, 5, {4, parent->zebraval + 1});
298 else
299 return path(h, 3, 5, {4, 3, 0});
300 case 5:
301 if(parent->zebraval > 0)
302 return path(h, 5, 3, {4, parent->zebraval - 1});
303 else
304 return path(h, 5, 3, {4, 5, 2});
305 default:
306 throw hr_exception("wrong dir");
307 }
308 }
309 #if MAXMDIM >= 4
310 case gBinary3: {
311 switch(d) {
312 case 0: case 1:
313 case 2: case 3:
314 return build3(parent, d, 8, 1);
315 case 8:
316 return build3(parent, 8, nextdir(4), -1);
317 case 4:
318 parent->cmove(8);
319 if(parent->c.spin(8) & 1)
320 return path(h, 4, 5, {8, parent->c.spin(8) ^ 1});
321 else
322 return path(h, 4, 5, {8, 4, parent->c.spin(8) ^ 1});
323 case 5:
324 parent->cmove(8);
325 if(!(parent->c.spin(8) & 1))
326 return path(h, 5, 4, {8, parent->c.spin(8) ^ 1});
327 else
328 return path(h, 5, 4, {8, 5, parent->c.spin(8) ^ 1});
329 case 6:
330 parent->cmove(8);
331 if(parent->c.spin(8) & 2)
332 return path(h, 6, 7, {8, parent->c.spin(8) ^ 2});
333 else
334 return path(h, 6, 7, {8, 6, parent->c.spin(8) ^ 2});
335 case 7:
336 parent->cmove(8);
337 if(!(parent->c.spin(8) & 2))
338 return path(h, 7, 6, {8, parent->c.spin(8) ^ 2});
339 else
340 return path(h, 7, 6, {8, 7, parent->c.spin(8) ^ 2});
341 default:
342 throw hr_exception("wrong dir");
343 }
344 }
345 case gHoroRec: {
346 switch(d) {
347 case 0: case 1:
348 return build3(parent, d, 6, 1);
349 case 6:
350 return build3(parent, 6, nextdir(2), -1);
351 case 2:
352 parent->cmove(6);
353 if(parent->c.spin(6) == 0)
354 return path(h, 2, 4, {6, 1});
355 else
356 return path(h, 2, 4, {6, 3, 0});
357 case 4:
358 parent->cmove(6);
359 if(parent->c.spin(6) == 0)
360 return path(h, 4, 2, {6, 5, 1});
361 else
362 return path(h, 4, 2, {6, 0});
363 case 3:
364 parent->cmove(6);
365 return path(h, 3, 5, {6, 4, parent->c.spin(6)});
366 case 5:
367 parent->cmove(6);
368 return path(h, 5, 3, {6, 2, parent->c.spin(6)});
369 default:
370 throw hr_exception("wrong dir");
371 }
372 }
373 case gHoroTris: {
374 switch(d) {
375 case 0: case 1: case 2: case 3:
376 return build3(parent, d, 7, 1);
377 case 7:
378 return build3(parent, 7, nextdir(3), -1);
379 case 4: case 5: case 6: {
380 parent->cmove(7);
381 int s = parent->c.spin(7);
382 if(s == 0) return path(h, d, d, {7, d-3});
383 else if(s == d-3) return path(h, d, d, {7, 0});
384 else return path(h, d, d, {7, d, 9-d-s});
385 }
386 default:
387 throw hr_exception("wrong dir");
388 }
389 }
390 case gHoroHex: {
391 // the comment is a picture...
392 // generated with the help of hexb.cpp
393 switch(d) {
394 case 0: case 1: case 2:
395 return build3(parent, d, 13, 1);
396 case 13:
397 return build3(parent, 13, nextdir(3), -1);
398 case 3:
399 return pathc(h, 3, 12, {{13,4,2}, {13,5,2}, {13,3,2}});
400 case 4:
401 return pathc(h, 4, 12, {{13,6,2,0}, {13,7,0,0}, {13,8,1,0}});
402 case 5:
403 return pathc(h, 5, 12, {{13,1,1}, {13,2,1}, {13,0,1}});
404 case 6:
405 return pathc(h, 6, 10, {{13,5}, {13,3}, {13,4}});
406 case 7:
407 return pathc(h, 7, 11, {{13,2}, {13,0}, {13,1}});
408 case 8:
409 return pathc(h, 8, 9, {{13,6,0}, {13,7,1}, {13,8,2}});
410 case 9:
411 return pathc(h, 9, 8, {{13,4}, {13,5}, {13,3}});
412 case 10:
413 return pathc(h, 10, 6, {{13,6,2}, {13,7,0}, {13,8,1}});
414 case 11:
415 return pathc(h, 11, 7, {{13,1}, {13,2}, {13,0}});
416 case 12: {
417 h->cmove(13);
418 int z = h->c.spin(13);
419 return path(h, 12, (z+1)%3+3, {13, z+6});
420 }
421 default:
422 throw hr_exception("wrong dir");
423 }
424 }
425 #endif
426 default:
427 throw hr_exception("wrong geometry");
428 }
429 }
430
shvidhr::bt::hrmap_binary431 int shvid(cell *c) override {
432 if(geometry == gBinaryTiling)
433 return c->type-6;
434 else if(geometry == gBinary4 || geometry == gTernary)
435 return c->master->zebraval;
436 else
437 return 0;
438 }
439
get_cornerhr::bt::hrmap_binary440 hyperpoint get_corner(cell *c, int cid, ld cf) override {
441 if(WDIM == 3) {
442 println(hlog, "get_corner_position called");
443 return C0;
444 }
445 return mid_at_actual(bt::get_horopoint(bt::get_corner_horo_coordinates(c, cid)), 3/cf);
446 }
447
updir_athr::bt::hrmap_binary448 int updir_at(heptagon *h) {
449 if(geometry != gBinaryTiling) return updir();
450 else if(type_of(h) == 6) return bd_down;
451 else if(mapside(h) == 1) return bd_left;
452 else if(mapside(h) == -1) return bd_right;
453 else throw hr_exception("wrong dir");
454 }
455
relative_matrixhhr::bt::hrmap_binary456 transmatrix relative_matrixh(heptagon *h2, heptagon *h1, const hyperpoint& hint) override {
457 if(gmatrix0.count(h2->c7) && gmatrix0.count(h1->c7))
458 return inverse_shift(gmatrix0[h1->c7], gmatrix0[h2->c7]);
459 transmatrix gm = Id, where = Id;
460 while(h1 != h2) {
461 if(h1->distance <= h2->distance) {
462 int d = updir_at(h2);
463 where = iadj(h2, d) * where;
464 h2 = may_create_step(h2, d);
465 }
466 else {
467 int d = updir_at(h1);
468 gm = gm * adj(h1, d);
469 h1 = may_create_step(h1, d);
470 }
471 }
472 return gm * where;
473 }
474
spin_anglehr::bt::hrmap_binary475 ld spin_angle(cell *c, int d) override {
476 if(WDIM == 3 || geometry == gBinary4 || geometry == gTernary) {
477 return hrmap::spin_angle(c, d);
478 }
479 if(d == NODIR) return 0;
480 if(d == c->type-1) d++;
481 return -(d+2)*M_PI/4;
482 }
483
adjhr::bt::hrmap_binary484 transmatrix adj(heptagon *h, int dir) override {
485 if(geometry == gBinaryTiling) switch(dir) {
486 case bd_up: return xpush(-log(2));
487 case bd_left: return parabolic(-1);
488 case bd_right: return parabolic(+1);
489 case bd_down:
490 if(h->type == 6) return xpush(log(2));
491 /* case bd_down_left: */
492 return parabolic(-1) * xpush(log(2));
493 case bd_down_right:
494 return parabolic(+1) * xpush(log(2));
495 case bd_up_left:
496 return xpush(-log(2)) * parabolic(-1);
497 case bd_up_right:
498 return xpush(-log(2)) * parabolic(1);
499 default:
500 throw hr_exception("unknown direction");
501 }
502 else if(use_direct_for(dir))
503 return cgi.direct_tmatrix[dir];
504 else {
505 h->cmove(dir);
506 return cgi.inverse_tmatrix[h->c.spin(dir)];
507 }
508 }
509
iadjhr::bt::hrmap_binary510 const transmatrix iadj(heptagon *h, int dir) { heptagon *h1 = h->cmove(dir); return adj(h1, h->c.spin(dir)); }
511
virtualRebasehr::bt::hrmap_binary512 void virtualRebase(heptagon*& base, transmatrix& at) override {
513
514 while(true) {
515
516 double currz = at[LDIM][LDIM];
517
518 heptagon *h = base;
519
520 heptagon *newbase = NULL;
521
522 transmatrix bestV;
523
524 for(int d=0; d<S7; d++) {
525 transmatrix V2 = iadj(h, d) * at;
526 double newz = V2[LDIM][LDIM];
527 if(newz < currz) {
528 currz = newz;
529 bestV = V2;
530 newbase = h->cmove(d);
531 }
532 }
533
534 if(newbase) {
535 base = newbase;
536 at = bestV;
537 continue;
538 }
539
540 return;
541 }
542 }
543
~hrmap_binaryhr::bt::hrmap_binary544 ~hrmap_binary() { if(origin) clearfrom(origin); }
545 };
546
new_map()547 EX hrmap *new_map() { return new hrmap_binary; }
548
549 struct hrmap_alternate_binary : hrmap_binary {
550 heptagon *origin;
hrmap_alternate_binaryhr::bt::hrmap_alternate_binary551 hrmap_alternate_binary(heptagon *o) { origin = o; }
~hrmap_alternate_binaryhr::bt::hrmap_alternate_binary552 ~hrmap_alternate_binary() { clearfrom(origin); }
553 };
554
new_alt_map(heptagon * o)555 EX hrmap *new_alt_map(heptagon *o) { return new hrmap_binary(o); }
556
557 /** \brief return if ew should use direct_tmatrix[dir] to get the adjacent cell the given direction
558 *
559 * Otherwise, this is the 'up' direction and thus we should use inverse_tmatrix for the inverse direction
560 */
use_direct_for(int dir)561 EX bool use_direct_for(int dir) {
562 return (cgi.use_direct >> dir) & 1;
563 }
564
565 /** \brief which coordinate is expanding */
expansion_coordinate()566 EX int expansion_coordinate() {
567 if(WDIM == 2) return 0;
568 return 2;
569 }
570
571 /** \brief by what factor does the area expand after moving one level in hr::bt::expansion_coordinate() */
area_expansion_rate()572 EX ld area_expansion_rate() {
573 switch(geometry) {
574 case gBinaryTiling: case gBinary4:
575 return 2;
576 case gTernary:
577 return 3;
578 case gBinary3: case gHoroTris:
579 return 4;
580 case gHoroRec:
581 return 2;
582 case gHoroHex:
583 return 3;
584 case gNil:
585 return 1;
586 case gEuclidSquare:
587 return 1;
588 case gKiteDart3:
589 return pow(golden_phi, 2);
590 case gSol:
591 return 1;
592 case gNIH:
593 return 6;
594 case gSolN:
595 return 3/2.;
596 case gArnoldCat:
597 return 1;
598
599 default:
600 return 0;
601 }
602 }
603
604 /** \brief by what factor do the lengths expand after moving one level in hr::bt::expansion_coordinate() */
expansion()605 EX ld expansion() {
606 if(WDIM == 2) return area_expansion_rate();
607 else return sqrt(area_expansion_rate());
608 }
609
610 /** \brief Get a point in the current cell, normalized to [-1,1]^WDIM
611 *
612 * This function returns the matrix moving point (0,0,0) to the given point in a parallelogram-like box
613 * Dimensions of the box are normalized to [-1,1], and directions are the same as usual (i.e., expansion_coordinate() is the correct one)
614 *
615 * This should works for all geometries which actually have boxes.
616 *
617 * For binary-based tessellations which are not based on square sections (e.g. gKiteDart3), 'x' and 'y' coordinates are not given in [-1,1], but take binary_width into account
618 *
619 * Otherwise: just return h
620 *
621 * See also: in devmods/tests.cpp, -bt-test tests whether this works correctly
622 *
623 */
624
normalized_at(hyperpoint h)625 EX transmatrix normalized_at(hyperpoint h) {
626 ld z2 = -log(2) / 2;
627 ld z3 = -log(3) / 2;
628 ld bwhn = vid.binary_width / 2;
629 ld bwh = vid.binary_width * z2;
630 ignore(bwh); ignore(bwhn);
631 ld r2 = sqrt(2);
632 const ld hs = hororec_scale;
633 auto &x = h[0], &y = h[1], &z = h[2];
634 switch(geometry) {
635 case gBinaryTiling: case gBinary4:
636 return bt::parabolic(y/2) * xpush(x*z2);
637 case gTernary:
638 return bt::parabolic(y/2) * xpush(x*z3);
639 #if CAP_SOLV
640 case gSol:
641 return xpush(bwh*x) * ypush(bwh*y) * zpush(z2*z);
642 case gSolN: case gNIH:
643 return xpush(bwhn*x) * ypush(bwhn*y) * zpush(-z*.5);
644 case gArnoldCat:
645 return rgpushxto0(asonov::tx*x/2 + asonov::ty*y/2 + asonov::tz*z/2);
646 #endif
647 case gNil:
648 return rgpushxto0(point31(x/2, y/2, z/2));
649 case gEuclidSquare:
650 return rgpushxto0(hpxy(x, y));
651 case gBinary3:
652 return parabolic3(x,y) * xpush(z*z2);
653 case gHoroRec:
654 return parabolic3(r2*hs*x, 2*hs*y) * xpush(z*z2/2);
655 case gHoroTris:
656 return parabolic3(x,y) * xpush(z*z2);
657 case gHoroHex:
658 return parabolic3(x,y) * xpush(z*z3/2);
659 case gKiteDart3:
660 return parabolic3(x,y) * xpush(-z*log_golden_phi/2);
661 default:
662 return rgpushxto0(h);
663 }
664 }
665
666 EX transmatrix normalized_at(ld x, ld y, ld z IS(0)) {
667 return normalized_at(point3(x, y, z));
668 }
669
updir()670 EX int updir() {
671 if(geometry == gBinary4) return 3;
672 if(geometry == gTernary) return 4;
673 if(geometry == gBinaryTiling) return 5;
674 if(kite::in()) return 0;
675 if(!bt::in()) return 0;
676 return S7-1;
677 }
678
dirs_outer()679 EX int dirs_outer() {
680 switch(geometry) {
681 case gBinary3: return 4;
682 case gHoroTris: return 4;
683 case gHoroRec: return 2;
684 case gHoroHex: return 6;
685 default: return -1;
686 }
687 }
688
dirs_inner()689 EX int dirs_inner() {
690 if(among(geometry, gBinaryTiling, gHoroHex)) return 2;
691 return 1;
692 }
693
build_tmatrix()694 EX void build_tmatrix() {
695 if(among(geometry, gBinaryTiling, gSol, gArnoldCat)) return; // unused
696 auto& direct_tmatrix = cgi.direct_tmatrix;
697 auto& inverse_tmatrix = cgi.inverse_tmatrix;
698 auto& use_direct = cgi.use_direct;
699 use_direct = (1 << (S7-1)) - 1;
700 if(geometry == gBinary4) {
701 use_direct = 3;
702 direct_tmatrix[0] = xpush(-log(2)) * parabolic(-0.5);
703 direct_tmatrix[1] = xpush(-log(2)) * parabolic(+0.5);
704 direct_tmatrix[2] = parabolic(1);
705 direct_tmatrix[4] = parabolic(-1);
706 use_direct = 1+2+4+16;
707 }
708 if(geometry == gTernary) {
709 direct_tmatrix[0] = xpush(-log(3)) * parabolic(-1);
710 direct_tmatrix[1] = xpush(-log(3));
711 direct_tmatrix[2] = xpush(-log(3)) * parabolic(+1);
712 direct_tmatrix[3] = parabolic(1);
713 direct_tmatrix[5] = parabolic(-1);
714 use_direct = 1+2+4+8+32;
715 }
716 if(geometry == gBinary3) {
717 direct_tmatrix[0] = xpush(-log(2)) * parabolic3(-1, -1);
718 direct_tmatrix[1] = xpush(-log(2)) * parabolic3(1, -1);
719 direct_tmatrix[2] = xpush(-log(2)) * parabolic3(-1, 1);
720 direct_tmatrix[3] = xpush(-log(2)) * parabolic3(1, 1);
721 direct_tmatrix[4] = parabolic3(-2, 0);
722 direct_tmatrix[5] = parabolic3(+2, 0);
723 direct_tmatrix[6] = parabolic3(0, -2);
724 direct_tmatrix[7] = parabolic3(0, +2);
725 }
726 if(geometry == gHoroTris) {
727 ld r3 = sqrt(3);
728 direct_tmatrix[0] = xpush(-log(2)) * cspin(1,2, M_PI);
729 direct_tmatrix[1] = parabolic3(0, +r3/3) * xpush(-log(2));
730 direct_tmatrix[2] = parabolic3(-0.5, -r3/6) * xpush(-log(2));
731 direct_tmatrix[3] = parabolic3(+0.5, -r3/6) * xpush(-log(2));
732 direct_tmatrix[4] = parabolic3(0, -r3*2/3) * cspin(1,2, M_PI);
733 direct_tmatrix[5] = parabolic3(1, r3/3) * cspin(1,2,M_PI);
734 direct_tmatrix[6] = parabolic3(-1, r3/3) * cspin(1,2,M_PI);
735 }
736 if(geometry == gHoroRec) {
737 ld r2 = sqrt(2);
738 ld l = -log(2)/2;
739 ld z = hororec_scale;
740 direct_tmatrix[0] = parabolic3(0, -z) * xpush(l) * cspin(2,1,M_PI/2);
741 direct_tmatrix[1] = parabolic3(0, +z) * xpush(l) * cspin(2,1,M_PI/2);
742 direct_tmatrix[2] = parabolic3(+2*r2*z, 0);
743 direct_tmatrix[3] = parabolic3(0, +4*z);
744 direct_tmatrix[4] = parabolic3(-2*r2*z, 0);
745 direct_tmatrix[5] = parabolic3(0, -4*z);
746 }
747 if(geometry == gHoroHex) {
748 // also generated with the help of hexb.cpp
749 ld l = log(3)/2;
750 auto& t = direct_tmatrix;
751 t[0] = parabolic3(horohex_scale, 0) * xpush(-l) * cspin(1, 2, M_PI/2);
752 t[1] = cspin(1, 2, 2*M_PI/3) * t[0];
753 t[2] = cspin(1, 2, 4*M_PI/3) * t[0];
754 auto it = iso_inverse(t[0]);
755
756 t[5] = it * t[1] * t[1];
757 t[6] = it * t[5];
758 t[4] = it * t[6] * t[2] * t[0];
759 t[3] = it * t[4] * t[2];
760
761 t[7] = it * t[2];
762 t[8] = it * t[6] * t[0];
763 t[9] = it * t[4];
764 t[10] = it * t[6] * t[2];
765 t[11] = it * t[1];
766
767 if(debugflags & DF_GEOM)
768 for(int a=0; a<12; a++)
769 println(hlog, t[a]);
770
771 use_direct >>= 1;
772 }
773 for(int i=0; i<S7; i++) if(use_direct_for(i))
774 inverse_tmatrix[i] = iso_inverse(direct_tmatrix[i]);
775 }
776
777 #if MAXMDIM == 4
778
queuecube(const shiftmatrix & V,ld size,color_t linecolor,color_t facecolor)779 EX void queuecube(const shiftmatrix& V, ld size, color_t linecolor, color_t facecolor) {
780 ld yy = log(2) / 2;
781 const int STEP=3;
782 const ld MUL = 1. / STEP;
783 auto at = [&] (ld x, ld y, ld z) { curvepoint(parabolic3(size*x, size*y) * xpush0(size*yy*z)); };
784 for(int a:{-1,1}) {
785 for(ld t=-STEP; t<STEP; t++) at(a, 1,t*MUL);
786 for(ld t=-STEP; t<STEP; t++) at(a, -t*MUL,1);
787 for(ld t=-STEP; t<STEP; t++) at(a, -1,-t*MUL);
788 for(ld t=-STEP; t<STEP; t++) at(a, t*MUL,-1);
789 at(a, 1,-1);
790 queuecurve(V, linecolor, facecolor, PPR::LINE);
791
792 for(ld t=-STEP; t<STEP; t++) at(1,t*MUL,a);
793 for(ld t=-STEP; t<STEP; t++) at(-t*MUL,1,a);
794 for(ld t=-STEP; t<STEP; t++) at(-1,-t*MUL,a);
795 for(ld t=-STEP; t<STEP; t++) at(t*MUL,-1,a);
796 at(1,-1,a);
797 queuecurve(V, linecolor, facecolor, PPR::LINE);
798
799 for(ld t=-STEP; t<STEP; t++) at(1,a,t*MUL);
800 for(ld t=-STEP; t<STEP; t++) at(-t*MUL,a,1);
801 for(ld t=-STEP; t<STEP; t++) at(-1,a,-t*MUL);
802 for(ld t=-STEP; t<STEP; t++) at(t*MUL,a,-1);
803 at(1,a,-1);
804 queuecurve(V, linecolor, facecolor, PPR::LINE);
805 }
806 /*for(int a:{-1,1}) for(int b:{-1,1}) for(int c:{-1,1}) {
807 at(0,0,0); at(a,b,c); queuecurve(linecolor, facecolor, PPR::LINE);
808 }*/
809 }
810 #endif
811
parabolic(ld u)812 EX transmatrix parabolic(ld u) {
813 return parabolic1(u * vid.binary_width / log(2) / 2);
814 }
815
parabolic3(ld y,ld z)816 EX transmatrix parabolic3(ld y, ld z) {
817 ld co = vid.binary_width / log(2) / 4;
818 return hr::parabolic13(y * co, z * co);
819 }
820
821 // on which horocycle are we
horo_level(hyperpoint h)822 EX ld horo_level(hyperpoint h) {
823 h /= (1 + h[LDIM]);
824 h[0] -= 1;
825 h /= sqhypot_d(GDIM, h);
826 h[0] += .5;
827 return log(2) + log(-h[0]);
828 }
829
deparabolic3(hyperpoint h)830 EX hyperpoint deparabolic3(hyperpoint h) {
831 h /= (1 + h[3]);
832 hyperpoint one = point3(1,0,0);
833 h -= one;
834 h /= sqhypot_d(3, h);
835 h[0] += .5;
836 ld co = vid.binary_width / log(2) / 8;
837 return point3(log(2) + log(-h[0]), h[1] / co, h[2] / co);
838 }
839
840 #if CAP_COMMANDLINE
__anona40cbd760202null841 auto bt_config = arg::add2("-btwidth", [] {arg::shift_arg_formula(vid.binary_width); });
842 #endif
843
pseudohept(cell * c)844 EX bool pseudohept(cell *c) {
845 if(WDIM == 2)
846 return c->type & c->master->distance & 1;
847 else if(geometry == gHoroRec)
848 return c->c.spin(S7-1) == 0 && (c->master->distance & 1) && c->cmove(S7-1)->c.spin(S7-1) == 0;
849 else if(geometry == gHoroTris)
850 return c->c.spin(S7-1) == 0 && (c->master->distance & 1);
851 else
852 return (c->master->zebraval == 1) && (c->master->distance & 1);
853 }
854
gpvalue(heptagon * h)855 EX pair<gp::loc, gp::loc> gpvalue(heptagon *h) {
856 int d = h->c.spin(S7-1);
857 if(d == 0) return make_pair(gp::loc(0,0), gp::loc(-1,0));
858 else return make_pair(gp::eudir((d-1)*2), gp::loc(1,0));
859 }
860
861 // distance in a triangular grid
tridist(gp::loc v)862 EX int tridist(gp::loc v) {
863 using namespace gp;
864 int d = v.first - v.second;
865 int d0 = d % 3;
866 if(d0 == 1 || d0 == -2) return 1 + min(tridist(v - eudir(0)), min(tridist(v - eudir(2)), tridist(v - eudir(4))));
867 if(d0 == 2 || d0 == -1) return 1 + min(tridist(v + eudir(0)), min(tridist(v + eudir(2)), tridist(v + eudir(4))));
868 return length(v * loc(1,1)) * 2 / 3;
869 }
870
equalize(heptagon * & c1,heptagon * & c2)871 EX int equalize(heptagon*& c1, heptagon*& c2) {
872 int steps = 0;
873 int d1 = c1->distance;
874 int d2 = c2->distance;
875 while(d1 > d2) c1 = c1->cmove(S7-1), steps++, d1--;
876 while(d2 > d1) c2 = c2->cmove(S7-1), steps++, d2--;
877 return steps;
878 }
879
celldistance3_tri(heptagon * c1,heptagon * c2)880 EX int celldistance3_tri(heptagon *c1, heptagon *c2) {
881 using namespace gp;
882 int steps = equalize(c1, c2);
883 vector<pair<loc, loc> > m1, m2;
884 while(c1 != c2) {
885 m2.push_back(gpvalue(c2));
886 m1.push_back(gpvalue(c1));
887 c1 = c1->cmove(S7-1);
888 c2 = c2->cmove(S7-1);
889 steps += 2;
890 }
891 loc T1(0,0), T2(0,0), inv1(1,0), inv2(1,0);
892 int xsteps = steps;
893 while(isize(m1)) {
894 xsteps -= 2;
895 inv1 = inv1 * m1.back().second;
896 inv2 = inv2 * m2.back().second;
897 T1 = T1 + T1 + m1.back().first * inv1;
898 T2 = T2 + T2 + m2.back().first * inv2;
899 m1.pop_back(); m2.pop_back();
900 loc T0 = T2 - T1;
901 if(T0.first > 3 || T0.second > 3 || T0.first < -3 || T0.second < -3) break;
902 steps = min(steps, xsteps + tridist(T0));
903 }
904 return steps;
905 }
906
celldistance3_rec(heptagon * c1,heptagon * c2)907 EX int celldistance3_rec(heptagon *c1, heptagon *c2) {
908 int steps = equalize(c1, c2);
909 vector<int> dx;
910 while(c1 != c2) {
911 dx.push_back(c1->c.spin(S7-1) - c2->c.spin(S7-1));
912 c1 = c1->cmove(S7-1);
913 c2 = c2->cmove(S7-1);
914 steps += 2;
915 }
916 int xsteps = steps, sx = 0, sy = 0;
917 while(isize(dx)) {
918 xsteps -= 2;
919 tie(sx, sy) = make_pair(-sy, 2 * sx + dx.back());
920 dx.pop_back();
921 int ysteps = xsteps + abs(sx) + abs(sy);
922 if(ysteps < steps) steps = ysteps;
923 if(sx >= 8 || sx <= -8 || sy >= 8 || sy <= -8) break;
924 }
925 return steps;
926 }
927
celldistance3_square(heptagon * c1,heptagon * c2)928 EX int celldistance3_square(heptagon *c1, heptagon *c2) {
929 int steps = equalize(c1, c2);
930 vector<int> dx, dy;
931 while(c1 != c2) {
932 dx.push_back((c1->c.spin(S7-1) & 1) - (c2->c.spin(S7-1) & 1));
933 dy.push_back((c1->c.spin(S7-1) >> 1) - (c2->c.spin(S7-1) >> 1));
934 c1 = c1->cmove(S7-1);
935 c2 = c2->cmove(S7-1);
936 steps += 2;
937 }
938 int xsteps = steps, sx = 0, sy = 0;
939 while(isize(dx)) {
940 xsteps -= 2;
941 sx *= 2;
942 sy *= 2;
943 sx += dx.back(); sy += dy.back();
944 dx.pop_back(); dy.pop_back();
945 int ysteps = xsteps + abs(sx) + abs(sy);
946 if(ysteps < steps) steps = ysteps;
947 if(sx >= 8 || sx <= -8 || sy >= 8 || sy <= -8) break;
948 }
949 return steps;
950 }
951
952 // this algorithm is wrong: it never considers the "narrow gap" moves
celldistance3_hex(heptagon * c1,heptagon * c2)953 EX int celldistance3_hex(heptagon *c1, heptagon *c2) {
954 int steps = equalize(c1, c2);
955 vector<int> d1, d2;
956 while(c1 != c2) {
957 d1.push_back(c1->c.spin(S7-1));
958 d2.push_back(c2->c.spin(S7-1));
959 c1 = c1->cmove(S7-1);
960 c2 = c2->cmove(S7-1);
961 steps += 2;
962 }
963 int xsteps = steps;
964 dynamicval<eGeometry> g(geometry, gEuclid);
965 transmatrix T = Id;
966 while(isize(d1)) {
967 xsteps -= 2;
968
969 T = euscalezoom(hpxy(0,sqrt(3))) * eupush(1,0) * spin(-d2.back() * 2 * M_PI/3) * T * spin(d1.back() * 2 * M_PI/3) * eupush(-1,0) * euscalezoom(hpxy(0,-1/sqrt(3)));
970
971 d1.pop_back(); d2.pop_back();
972
973 hyperpoint h = tC0(T);
974 int sx = int(floor(h[0] - h[1] / sqrt(3) + .5)) / 3;
975 int sy = int(floor(h[1] * 2 / sqrt(3) + .5)) / 3;
976
977 int ysteps = xsteps + euc::dist(sx, sy);
978 if(ysteps < steps) steps = ysteps;
979 if(sx >= 8 || sx <= -8 || sy >= 8 || sy <= -8) break;
980 }
981 return steps;
982 }
983
celldistance3_approx(heptagon * c1,heptagon * c2)984 EX int celldistance3_approx(heptagon *c1, heptagon *c2) {
985 int d = 0;
986 while(true) {
987 if(d > 1000000) return d; /* sanity check */
988 if(c1 == c2) return d;
989 for(int i=0; i<c1->type; i++)
990 if(c1->move(i) == c2) return d + 1;
991 for(int i=0; i<c1->type; i++) {
992 heptagon *c3 = c1->move(i);
993 for(int j=0; j<c3->type; j++)
994 if(c3->move(j) == c2) return d+2;
995 }
996 if(c1->distance > c2->distance) c1=c1->cmove(updir()), d++;
997 else c2=c2->cmove(updir()), d++;
998 }
999 }
1000
celldistance3(heptagon * c1,heptagon * c2)1001 EX int celldistance3(heptagon *c1, heptagon *c2) {
1002 switch(geometry) {
1003 case gBinary3: return celldistance3_square(c1, c2);
1004 case gHoroTris: return celldistance3_tri(c1, c2);
1005 case gHoroRec: return celldistance3_rec(c1, c2);
1006 case gHoroHex: return celldistance3_hex(c1, c2);
1007 default:
1008 if(sol || !bt::in()) {
1009 println(hlog, "called celldistance3 for wrong geometry"); return 0;
1010 }
1011 return celldistance3_approx(c1, c2);
1012 }
1013 }
1014
celldistance3(cell * c1,cell * c2)1015 EX int celldistance3(cell *c1, cell *c2) { return celldistance3(c1->master, c2->master); }
1016
get_horopoint(ld y,ld x)1017 EX hyperpoint get_horopoint(ld y, ld x) {
1018 return xpush(-y) * bt::parabolic(x) * C0;
1019 }
1020
get_horopoint(hyperpoint h)1021 EX hyperpoint get_horopoint(hyperpoint h) {
1022 return get_horopoint(h[0], h[1]);
1023 }
1024
get_corner_horo_coordinates(cell * c,int i)1025 EX hyperpoint get_corner_horo_coordinates(cell *c, int i) {
1026 ld yx = log(2) / 2;
1027 ld yy = yx;
1028 ld xx = 1 / sqrt(2)/2;
1029 switch(geometry) {
1030 case gBinaryTiling:
1031 switch(gmod(i, c->type)) {
1032 case 0: return point2(-yy, xx);
1033 case 1: return point2(yy, 2*xx);
1034 case 2: return point2(yy, xx);
1035 case 3: return point2(yy, -xx);
1036 case 4: return point2(yy, -2*xx);
1037 case 5: return point2(-yy, -xx);
1038 case 6: return point2(-yy, 0);
1039 default: return point2(0, 0);
1040 }
1041
1042 case gBinary4:
1043 switch(gmod(i, c->type)) {
1044 case 0: return point2(yy, -2*xx);
1045 case 1: return point2(yy, +0*xx);
1046 case 2: return point2(yy, +2*xx);
1047 case 3: return point2(-yy, xx);
1048 case 4: return point2(-yy, -xx);
1049 default: return point2(0, 0);
1050 }
1051
1052 case gTernary:
1053 yy = log(3) / 2;
1054 xx = 1 / sqrt(3) / 2;
1055 switch(gmod(i, c->type)) {
1056 case 0: return point2(yy, -3*xx);
1057 case 1: return point2(yy, -1*xx);
1058 case 2: return point2(yy, +1*xx);
1059 case 3: return point2(yy, +3*xx);
1060 case 4: return point2(-yy, xx);
1061 case 5: return point2(-yy, -xx);
1062 default: return point2(0, 0);
1063 }
1064
1065 default:
1066 return point2(0, 0);
1067 }
1068 return point2(0, 0);
1069 }
1070
make4(hyperpoint a,hyperpoint b,hyperpoint c)1071 vector<hyperpoint> make4(hyperpoint a, hyperpoint b, hyperpoint c) {
1072 return {a, b, b+c-a, c};
1073 }
1074
make5(hyperpoint a,hyperpoint b,hyperpoint c)1075 vector<hyperpoint> make5(hyperpoint a, hyperpoint b, hyperpoint c) {
1076 return {a, (a+b)/2, b, b+c-a, c};
1077 }
1078
create_faces()1079 EX void create_faces() {
1080
1081 if(geometry == gBinary3) {
1082 hyperpoint h00 = point3(-1,-1,-1);
1083 hyperpoint h01 = point3(-1,0,-1);
1084 hyperpoint h02 = point3(-1,+1,-1);
1085 hyperpoint h10 = point3(0,-1,-1);
1086 hyperpoint h11 = point3(0,0,-1);
1087 hyperpoint h12 = point3(0,+1,-1);
1088 hyperpoint h20 = point3(+1,-1,-1);
1089 hyperpoint h21 = point3(+1,0,-1);
1090 hyperpoint h22 = point3(+1,+1,-1);
1091 hyperpoint down = point3(0,0,2);
1092
1093 add_wall(0, make4(h11, h01, h10));
1094 add_wall(1, make4(h11, h21, h10));
1095 add_wall(2, make4(h11, h01, h12));
1096 add_wall(3, make4(h11, h21, h12));
1097 add_wall(4, make5(h00, h02, h00+down));
1098 add_wall(5, make5(h20, h22, h20+down));
1099 add_wall(6, make5(h00, h20, h00+down));
1100 add_wall(7, make5(h02, h22, h02+down));
1101 add_wall(8, make4(h22+down, h02+down, h20+down));
1102 }
1103
1104 if(GDIM == 3 && bt::in() && geometry == gHoroTris) {
1105 ld r = sqrt(3)/6;
1106 ld r1 = r;
1107 ld r2 = r * 2;
1108
1109 hyperpoint t0 = point3(0,-r2,-1);
1110 hyperpoint t1 = point3(+.5,r1,-1);
1111 hyperpoint t2 = point3(-.5,r1,-1);
1112 hyperpoint shift = point3(0,0,-3);
1113 hyperpoint down = point3(0,0,2);
1114 hyperpoint d0 = -2 * t0 + shift;
1115 hyperpoint d1 = -2 * t1 + shift;
1116 hyperpoint d2 = -2 * t2 + shift;
1117
1118 add_wall(0, {t0, t1, t2});
1119 add_wall(1, {d0, t1, t2});
1120 add_wall(2, {t0, d1, t2});
1121 add_wall(3, {t0, t1, d2});
1122 add_wall(4, make5(d2, d1, d2 + down));
1123 add_wall(5, make5(d0, d2, d0 + down));
1124 add_wall(6, make5(d1, d0, d1 + down));
1125 add_wall(7, {d0+down, d1+down, d2+down});
1126 }
1127
1128 if(geometry == gHoroRec) {
1129 ld r2 = sqrt(2);
1130 ld z = bt::hororec_scale;
1131
1132 hyperpoint a00 = point3(-r2*z,-2*z,-.5);
1133 hyperpoint a01 = point3(+r2*z,-2*z,-.5);
1134 hyperpoint a10 = point3(-r2*z, 0*z,-.5);
1135 hyperpoint a11 = point3(+r2*z, 0*z,-.5);
1136 hyperpoint a20 = point3(-r2*z,+2*z,-.5);
1137 hyperpoint a21 = point3(+r2*z,+2*z,-.5);
1138
1139 hyperpoint down = point3(0,0,1);
1140
1141 add_wall(0, make4(a00, a01, a10));
1142 add_wall(1, make4(a10, a11, a20));
1143 add_wall(2, make5(a01, a21, a01+down));
1144 add_wall(3, make4(a21, a20, a21+down));
1145 add_wall(4, make5(a20, a00, a20+down));
1146 add_wall(5, make4(a00, a01, a00+down));
1147 add_wall(6, make4(a00+down, a01+down, a20+down));
1148 }
1149
1150 if(geometry == gHoroHex) {
1151 ld z = log(3) / log(2) / 2;
1152 ld r3 = sqrt(3) / 2 * bt::horohex_scale;
1153 ld h = bt::horohex_scale / 2;
1154 hyperpoint down = point3(0,0,2*z);
1155
1156 for(int j=0; j<4; j++) for(int i=0; i<3; i++) {
1157 transmatrix T = cspin(0, 1, 2*M_PI*i/3);
1158
1159 hyperpoint hcenter = point3(0,0,-z);
1160 hyperpoint hu0 = T*point3(+h, +r3,-z);
1161 hyperpoint hu1 = T*point3(+h*3,+r3,-z);
1162 hyperpoint hd0 = T*point3(+h, -r3,-z);
1163 hyperpoint hd1 = T*point3(+h*3,-r3,-z);
1164 hyperpoint hcn = T*point3(-h*2,0, -z);
1165 hyperpoint hun = T*point3(-h*3,+r3,-z);
1166 hyperpoint hdn = T*point3(-h*3,-r3,-z);
1167 if(j == 0) add_wall(i, {hcenter, hu0, hu1, hd1, hd0});
1168 if(j == 1) add_wall(i+3, {hcn, hun, hdn});
1169 if(j == 2) add_wall(i+6, make4(hd1, hu1, hd1+down));
1170 if(j == 3) add_wall(i+9, make4(hun, hdn, hun+down));
1171 }
1172
1173 add_wall(12, {point3(3*h,r3,z), point3(0,2*r3,z), point3(-3*h,r3,z)});
1174 add_wall(13, {point3(3*h,r3,z), point3(3*h,-r3,z), point3(0,-2*r3,z), point3(-3*h,-r3,z), point3(-3*h,r3,z)});
1175 }
1176
1177 if(kite::in()) {
1178 auto kv = kite::make_walls();
1179 for(auto& v: kv.first) for(auto& h: v) {
1180 h = bt::deparabolic3(h);
1181 h = point3(h[1], h[2], h[0] / (log(2)/2));
1182 }
1183 for(int i=0; i<isize(kv.first); i++) {
1184 add_wall(i, kv.first[i]);
1185 }
1186 get_hsh().weights = kv.second;
1187 }
1188
1189 get_hsh().compute_hept();
1190 }
1191
__anona40cbd760302null1192 auto hooksw = addHook(hooks_swapdim, 100, [] {
1193 if(bt::in()) build_tmatrix();
1194 });
1195 #endif
1196
1197 }
1198
1199 }
1200