From 442d14ce0c7fb818cc981e21e23ca9276ef6c2d8 Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Thu, 31 May 2018 00:17:32 +0200 Subject: [PATCH] hyperbolic tiling --- hyperbolic.c | 265 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 265 insertions(+) create mode 100644 hyperbolic.c diff --git a/hyperbolic.c b/hyperbolic.c new file mode 100644 index 0000000..53c2450 --- /dev/null +++ b/hyperbolic.c @@ -0,0 +1,265 @@ +#include +#include +#include + +#include "triangle.h" +#include "linalg.h" + +#define POINCARE 1 +#define LOOP(i) for(int i = 0; i < 3; i++) + +void cartan_matrix(gsl_matrix *cartan, double a1, double a2, double a3, double s) +{ + gsl_matrix_set(cartan, 0, 0, -2); + gsl_matrix_set(cartan, 0, 1, 2*s*cos(a3)); + gsl_matrix_set(cartan, 0, 2, 2/s*cos(a2)); + + gsl_matrix_set(cartan, 1, 0, 2/s*cos(a3)); + gsl_matrix_set(cartan, 1, 1, -2); + gsl_matrix_set(cartan, 1, 2, 2*s*cos(a1)); + + gsl_matrix_set(cartan, 2, 0, 2*s*cos(a2)); + gsl_matrix_set(cartan, 2, 1, 2/s*cos(a1)); + gsl_matrix_set(cartan, 2, 2, -2); +} + +void initialize_triangle_generators(gsl_matrix **gen, gsl_matrix *cartan) +{ + LOOP(i) { + gsl_matrix_set_identity(gen[i]); + LOOP(j) *gsl_matrix_ptr(gen[i], i, j) += gsl_matrix_get(cartan, i, j); + } +} + +typedef struct { double x[3]; } point; + +point apply(gsl_matrix *g, point x) +{ + point out; + LOOP(i) out.x[i] = 0.0; + LOOP(i) LOOP(j) out.x[i] += gsl_matrix_get(g, i, j) * x.x[j]; + return out; +} + +point incidence(point a, point b) +{ + point c; + LOOP(i) c.x[i] = a.x[(i+1)%3]*b.x[(i+2)%3] - a.x[(i+2)%3]*b.x[(i+1)%3]; + return c; +} + +double applyForm(point x, point y, gsl_matrix *form) +{ + double out = 0.0; + LOOP(i) LOOP(j) out += x.x[i] * gsl_matrix_get(form, i, j) * y.x[j]; + return out; +} + +void boundary(point x, point y, point *a, point *b, gsl_matrix *form) +{ + double formX = applyForm(x, x, form); + double formY = applyForm(y, y, form); + double formXY = applyForm(x, y, form); + double t1 = (- formXY + sqrt(formXY * formXY - formX * formY)) / formX; + double t2 = (- formXY - sqrt(formXY * formXY - formX * formY)) / formX; + LOOP(i) a->x[i] = t1 * x.x[i] + y.x[i]; + LOOP(i) b->x[i] = t2 * x.x[i] + y.x[i]; +} + +point column(gsl_matrix *m, int j) +{ + point out; + LOOP(i) out.x[i] = gsl_matrix_get(m, i, j); + return out; +} + +point row(gsl_matrix *m, int i) +{ + point out; + LOOP(j) out.x[j] = gsl_matrix_get(m, i, j); + return out; +} + +double coord(point p, int k, gsl_matrix *frame) +{ + double tmp[3] = {0, 0, 0}; + double factor = 1.0; + + LOOP(i) LOOP(j) tmp[i] += gsl_matrix_get(frame, i, j)*p.x[j]; + +#ifdef POINCARE + double norm2 = (tmp[0]*tmp[0] + tmp[1]*tmp[1]) / (tmp[2]*tmp[2]); + if(fabs(norm2) < 0.5) // just to avoid dividing by 0 + factor = 1.0 / (1.0 + sqrt(1.0 - norm2)); + else if(fabs(norm2) < 1.0) + factor = (1.0 - sqrt(1.0 - norm2)) / norm2; + else + factor = 1.0; +#endif + + return factor*tmp[k]/tmp[2]; +} + +void print_tex_header() { + char *header = + "\\documentclass{standalone}\n" + "\\usepackage[utf8]{inputenc}\n" + "\\usepackage{tikz}\n" + "\\begin{document}\n" + "\\begin{tikzpicture}[scale=5]\n"; + + printf("%s", header); +} + +void print_tex_footer() { + char *footer = + "\\draw[thick] (0,0) circle (1.0);\n" + "\\end{tikzpicture}\n" + "\\end{document}\n"; + + printf("%s", footer); +} + +char *poincare_arc(double x1, double x2, double y1, double y2, char *buffer) +{ + double t = (1 - x1*y1 - x2*y2)/(x1*y2 - x2*y1)/2; + double m1 = x1/2 + y1/2 + t*(y2 - x2); // center of circle + double m2 = x2/2 + y2/2 + t*(x1 - y1); + + double r = sqrt((x1-m1)*(x1-m1) + (x2-m2)*(x2-m2)); + double phix = atan2(x2-m2,x1-m1); + double phiy = atan2(y2-m2,y1-m1); + + if(phix - phiy > M_PI) + phiy += 2*M_PI; + else if(phiy - phix > M_PI) + phix += 2*M_PI; + + sprintf(buffer, "%f:%f:%f", phix/M_PI*180, phiy/M_PI*180, r); + return buffer; +} + +void draw_triangle(point *p, gsl_matrix *frame, const char *arguments) +{ +#ifdef POINCARE + char buffer1[100], buffer2[100], buffer3[100]; + + double x1 = coord(p[0], 0, frame); + double y1 = coord(p[0], 1, frame); + double x2 = coord(p[1], 0, frame); + double y2 = coord(p[1], 1, frame); + double x3 = coord(p[2], 0, frame); + double y3 = coord(p[2], 1, frame); + + printf("\\draw[%s] (%f, %f) arc (%s) arc (%s) arc (%s);\n", arguments, x1, y1, + poincare_arc(x1, y1, x2, y2, buffer1), + poincare_arc(x2, y2, x3, y3, buffer2), + poincare_arc(x3, y3, x1, y1, buffer3)); + +#else + printf("\\draw[%s] (%f, %f) -- (%f, %f) -- (%f, %f) -- cycle;\n", arguments, + coord(p[0], 0, frame), coord(p[0], 1, frame), + coord(p[1], 0, frame), coord(p[1], 1, frame), + coord(p[2], 0, frame), coord(p[2], 1, frame)); +#endif +} + +void draw_line(point p1, point p2, gsl_matrix *frame, const char *arguments) +{ + char buffer[100]; + + double x1 = coord(p1, 0, frame); + double y1 = coord(p1, 1, frame); + double x2 = coord(p2, 0, frame); + double y2 = coord(p2, 1, frame); + +#ifdef POINCARE + printf("\\draw[%s] (%f, %f) arc (%s);\n", arguments, x1, y1, + poincare_arc(x1, y1, x2, y2, buffer)); +#else + printf("\\draw[%s] (%f, %f) -- (%f, %f);\n", arguments, x1, y1, x2, y2); +#endif +} + +int main() +{ + groupelement_t *group; + gsl_matrix **matrices; + gsl_matrix *cartan; + gsl_matrix *gen[3]; + gsl_matrix *coxeter[3]; + gsl_matrix *coxeter_eigenvectors[3]; + gsl_matrix *frame; + workspace_t *ws; + int elements = 5000; + int p = 5, q = 5, r = 5; + + group = malloc(elements*sizeof(groupelement_t)); + matrices = malloc(elements*sizeof(gsl_matrix*)); + for(int i = 0; i < elements; i++) + matrices[i] = gsl_matrix_alloc(3, 3); + cartan = gsl_matrix_alloc(3, 3); + frame = gsl_matrix_alloc(3, 3); + LOOP(i) gen[i] = gsl_matrix_alloc(3, 3); + LOOP(i) coxeter[i] = gsl_matrix_alloc(3, 3); + LOOP(i) coxeter_eigenvectors[i] = gsl_matrix_alloc(3, 3); + ws = workspace_alloc(3); + + generate_triangle_group(group, elements, p, q, r); + cartan_matrix(cartan, M_PI/p, M_PI/q, M_PI/r, 1.0); + initialize_triangle_generators(gen, cartan); + int pos = diagonalize_symmetric_form(cartan, frame, ws); // choose frame of reference which diagonalizes the form + + gsl_matrix_set_identity(matrices[0]); + for(int i = 1; i < elements; i++) + multiply(matrices[group[i].parent->id], gen[group[i].letter], matrices[i]); + + LOOP(i) multiply_many(ws, coxeter[i], 3, gen[i%3], gen[(i+1)%3], gen[(i+2)%3]); + LOOP(i) eigenvectors(coxeter[i], coxeter_eigenvectors[i], ws); + + /* + for(int i = 0; i < elements; i++) { + printf("%4d: ", i); + for(groupelement_t *cur = &group[i]; cur->parent; cur = cur->parent) + fputc(cur->letter+'a', stdout); + fputc('\n', stdout); + } + + return 0;*/ + + point coxeter_attracting[3]; + point coxeter_repelling[3]; + point reflection_lines[3]; + point triangle_points[3]; + point transformed[3]; + point transformed2[3]; + + LOOP(i) coxeter_attracting[i] = column(coxeter_eigenvectors[i], 0); + LOOP(i) coxeter_repelling[i] = column(coxeter_eigenvectors[i], 2); + LOOP(i) reflection_lines[i] = row(cartan, i); + LOOP(i) triangle_points[i] = incidence(reflection_lines[(i+1)%3], reflection_lines[(i+2)%3]); + + print_tex_header(); + +// int indices[10] = {0, 1, 6, 10, 30, 46, 124, 185, 484, 717}; +// int indices[10] = {0, 1, 4, 10, 22}; + + for(int k = 0; k < elements; k++) { + if(group[k].length % 2) + continue; + + LOOP(i) transformed[i] = apply(matrices[k], triangle_points[i]); + draw_triangle(transformed, frame, "black,fill=black!10,line width=0pt"); + } + + for(int k = 0; k < elements; k++) { + if(group[k].length % 2) + continue; + + LOOP(i) transformed[i] = apply(matrices[k], coxeter_attracting[i]); + LOOP(i) transformed2[i] = apply(matrices[k], coxeter_repelling[i]); + LOOP(i) draw_line(transformed[i], transformed2[i], frame, "red"); + } + + print_tex_footer(); +}