From 0cde265d1ede849e37d0ec4dd7bc3d2ece0d9430 Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Mon, 23 Dec 2019 12:29:50 +0100 Subject: [PATCH] draw arcs --- draw.c | 619 ++++++++++++++++++++++++++++++++++++++++++---------- limit_set.c | 29 +-- linalg.h | 2 +- main.c | 82 +++++-- main.h | 10 +- queue.h | 2 +- 6 files changed, 592 insertions(+), 152 deletions(-) diff --git a/draw.c b/draw.c index d0ad184..1a873a0 100644 --- a/draw.c +++ b/draw.c @@ -1,5 +1,10 @@ #include "main.h" +#define FMOD(x,y) (fmod(x,y) < 0 ? fmod(x,y) + y : fmod(x,y)) +#define ANGLE_DIFF(x,y) (FMOD((x)-(y), 2*M_PI)) +#define ANGLE_IN_INTERVAL(a,b,x) (ANGLE_DIFF(x,a) < ANGLE_DIFF(b,a)) +#define FLIP(x,y) do {double tmp = x; x = y; y = tmp;} while(0) + // level 0: helper functions int isInsideBB(DrawingContext *ctx, point_t p) @@ -27,6 +32,16 @@ vector_t apply(gsl_matrix *m, vector_t x) return out; } +vector_t apply_transpose(gsl_matrix *m, vector_t x) +{ + vector_t out; + + LOOP(i) out.x[i] = 0.0; + LOOP(i) LOOP(j) out.x[i] += gsl_matrix_get(m, j, i) * x.x[j]; + + return out; +} + int fixedPoints(DrawingContext *ctx, const char *word, vector_t *out) { gsl_matrix *tmp = getTempMatrix(ctx->ws); @@ -56,6 +71,13 @@ void drawPoint(DrawingContext *ctx, point_t p) { cairo_t *C = ctx->cairo; + cairo_save(C); + cairo_arc(C, p.x, p.y, 5.0/ctx->dim->scalefactor, 0, 2*M_PI); + cairo_close_path(C); + cairo_fill(C); + cairo_restore(C); + + /* cairo_save(C); cairo_move_to(C, p.x, p.y); cairo_close_path(C); @@ -63,6 +85,7 @@ void drawPoint(DrawingContext *ctx, point_t p) cairo_set_line_width(C, 10.0/ctx->dim->scalefactor); cairo_stroke(C); cairo_restore(C); + */ } void drawSegment2d(DrawingContext *ctx, point_t a, point_t b) @@ -223,7 +246,7 @@ void drawRotationOrbit(DrawingContext *ctx, const char *word, vector_t start) { vector_t v[3], w; point_t p; - double parameter; + double parameter, startangle; int iterations = 200; gsl_matrix *frame = getTempMatrix(ctx->ws); gsl_matrix *inverse = getTempMatrix(ctx->ws); @@ -239,6 +262,7 @@ void drawRotationOrbit(DrawingContext *ctx, const char *word, vector_t start) parameter = gsl_vector_get(start_in_frame, 2); parameter /= sqrt(gsl_vector_get(start_in_frame, 0)*gsl_vector_get(start_in_frame, 0) + gsl_vector_get(start_in_frame, 1)*gsl_vector_get(start_in_frame, 1)); + startangle = atan2(gsl_vector_get(start_in_frame, 1), gsl_vector_get(start_in_frame, 0)); int previous_inside = 0; for(int k = 0; k <= iterations; k++) { @@ -262,6 +286,167 @@ void drawRotationOrbit(DrawingContext *ctx, const char *word, vector_t start) releaseTempVectors(ctx->ws, 2); } +void drawDualRotationOrbit(DrawingContext *ctx, const char *word, vector_t start) +{ + vector_t v[3], w; + point_t p; + double parameter, startangle; + int iterations = 200; + gsl_matrix *frame = getTempMatrix(ctx->ws); + gsl_matrix *inverse = getTempMatrix(ctx->ws); + gsl_vector *start_v = getTempVector(ctx->ws); + gsl_vector *start_in_frame = getTempVector(ctx->ws); + cairo_t *C = ctx->cairo; + + computeRotationMatrix(ctx, frame, word); + LOOP(i) LOOP(j) v[i].x[j] = gsl_matrix_get(frame, j, i); + + LOOP(i) gsl_vector_set(start_v, i, start.x[i]); +// solve(frame, start_v, start_in_frame, ctx->ws); + gsl_blas_dgemv(CblasTrans, 1.0, frame, start_v, 0.0, start_in_frame); + parameter = sqrt(gsl_vector_get(start_in_frame, 0)*gsl_vector_get(start_in_frame, 0) + + gsl_vector_get(start_in_frame, 1)*gsl_vector_get(start_in_frame, 1)); + parameter /= gsl_vector_get(start_in_frame, 2); + startangle = atan2(gsl_vector_get(start_in_frame, 1), gsl_vector_get(start_in_frame, 0)); + + int previous_inside = 0; + for(int k = 0; k <= iterations; k++) { + LOOP(i) w.x[i] = parameter * v[2].x[i] + cos(2*k*M_PI/iterations) * v[0].x[i] + sin(2*k*M_PI/iterations) * v[1].x[i]; + p = vectorToPoint(ctx, w); + + if(isInsideBB(ctx, p)) { + if(!previous_inside) + cairo_move_to(C, p.x, p.y); + else + cairo_line_to(C, p.x, p.y); + previous_inside = 1; + } else { + previous_inside = 0; + } + } + + cairo_stroke(C); + + releaseTempMatrices(ctx->ws, 2); + releaseTempVectors(ctx->ws, 2); +} + +void drawArc(DrawingContext *ctx, const char *word, vector_t start, vector_type_t starttype, vector_t end, vector_type_t endtype, vector_t third, int contain) +{ + vector_t v[3], w; + point_t p; + double radius, angle_start, angle_end, angle_third, angle, angle_end_delta, sign, angle_start_final, angle_end_final, angle_end_other; + int iterations = 200; + gsl_matrix *frame = getTempMatrix(ctx->ws); + gsl_matrix *inverse = getTempMatrix(ctx->ws); + gsl_vector *vector = getTempVector(ctx->ws); + gsl_vector *vector_in_frame = getTempVector(ctx->ws); + cairo_t *C = ctx->cairo; + + computeRotationMatrix(ctx, frame, word); + LOOP(i) LOOP(j) v[i].x[j] = gsl_matrix_get(frame, j, i); + + LOOP(i) gsl_vector_set(vector, i, start.x[i]); + if(starttype == VT_POINT) { + solve(frame, vector, vector_in_frame, ctx->ws); + radius = sqrt(gsl_vector_get(vector_in_frame, 0)*gsl_vector_get(vector_in_frame, 0) + + gsl_vector_get(vector_in_frame, 1)*gsl_vector_get(vector_in_frame, 1)); + radius /= fabs(gsl_vector_get(vector_in_frame, 2)); + angle_start = atan2(gsl_vector_get(vector_in_frame, 1)/gsl_vector_get(vector_in_frame, 2), + gsl_vector_get(vector_in_frame, 0)/gsl_vector_get(vector_in_frame, 2)); + } else { + gsl_blas_dgemv(CblasTrans, 1.0, frame, vector, 0.0, vector_in_frame); + radius = fabs(gsl_vector_get(vector_in_frame, 2)); + radius /= sqrt(gsl_vector_get(vector_in_frame, 0)*gsl_vector_get(vector_in_frame, 0) + + gsl_vector_get(vector_in_frame, 1)*gsl_vector_get(vector_in_frame, 1)); + + angle_start = atan2(gsl_vector_get(vector_in_frame, 1)/gsl_vector_get(vector_in_frame, 2), + gsl_vector_get(vector_in_frame, 0)/gsl_vector_get(vector_in_frame, 2)); + } + + LOOP(i) gsl_vector_set(vector, i, third.x[i]); + solve(frame, vector, vector_in_frame, ctx->ws); + angle_third = atan2(gsl_vector_get(vector_in_frame, 1)/gsl_vector_get(vector_in_frame, 2), + gsl_vector_get(vector_in_frame, 0)/gsl_vector_get(vector_in_frame, 2)); + + LOOP(i) gsl_vector_set(vector, i, end.x[i]); + if(endtype == VT_POINT) { + solve(frame, vector, vector_in_frame, ctx->ws); + angle_end = atan2(gsl_vector_get(vector_in_frame, 1)/gsl_vector_get(vector_in_frame, 2), + gsl_vector_get(vector_in_frame, 0)/gsl_vector_get(vector_in_frame, 2)); + } else { + gsl_blas_dgemv(CblasTrans, 1.0, frame, vector, 0.0, vector_in_frame); + + // this is only the average angle + angle_end = atan2(gsl_vector_get(vector_in_frame, 1)/gsl_vector_get(vector_in_frame, 2), + gsl_vector_get(vector_in_frame, 0)/gsl_vector_get(vector_in_frame, 2)); + + angle_end_delta = acos(-fabs(gsl_vector_get(vector_in_frame, 2))/radius/ + sqrt(gsl_vector_get(vector_in_frame, 0)*gsl_vector_get(vector_in_frame, 0) + + gsl_vector_get(vector_in_frame, 1)*gsl_vector_get(vector_in_frame, 1))); + + printf("angle_end_delta: %f\n", angle_end_delta); + } + + int previous_inside = 0; + + for(int i = 0; i < 4; i++) { + angle_start_final = angle_start; + + if(endtype == VT_POINT) { + angle_end_final = angle_end; + } else { + if(i >= 2) { + angle_end_final = angle_end - angle_end_delta; + angle_end_other = angle_end + angle_end_delta; + } else { + angle_end_final = angle_end + angle_end_delta; + angle_end_other = angle_end - angle_end_delta; + } + } + + if(i%2) + FLIP(angle_start_final, angle_end_final); + + printf("start: %f, end: %f, other: %f, third: %f\n", angle_start_final, angle_end_final, angle_end_other, angle_third); + printf("interval1: %d, interval2: %d, diff: %f %f %f\n", ANGLE_IN_INTERVAL(angle_start, angle_end_final, angle_end_other), ANGLE_IN_INTERVAL(angle_start_final, angle_end_final, angle_third), ANGLE_DIFF(angle_end_final, angle_start_final), ANGLE_DIFF(angle_end_other, angle_start_final), ANGLE_DIFF(angle_third, angle_start_final)); + + if(endtype == VT_LINE && ANGLE_IN_INTERVAL(angle_start_final, angle_end_final, angle_end_other)) + continue; + if(contain && !ANGLE_IN_INTERVAL(angle_start_final, angle_end_final, angle_third)) + continue; + if(!contain && ANGLE_IN_INTERVAL(angle_start_final, angle_end_final, angle_third)) + continue; + + break; + } + + +// printf("angle_start: %f, angle_end: %f, angle_third: %f\n", angle_start, angle_end, angle_third); + + for(int k = 0; k <= iterations; k++) { + angle = angle_start_final + (double)k/(double)iterations * ANGLE_DIFF(angle_end_final, angle_start_final); + + LOOP(i) w.x[i] = v[2].x[i] / radius + cos(angle) * v[0].x[i] + sin(angle) * v[1].x[i]; + p = vectorToPoint(ctx, w); + + if(isInsideBB(ctx, p)) { + if(!previous_inside) + cairo_move_to(C, p.x, p.y); + else + cairo_line_to(C, p.x, p.y); + previous_inside = 1; + } else { + previous_inside = 0; + } + } + + cairo_stroke(C); + + releaseTempMatrices(ctx->ws, 2); + releaseTempVectors(ctx->ws, 2); +} + // level 4: draw the actual image components void drawReflectors(DrawingContext *ctx) @@ -279,90 +464,133 @@ void drawReflectors(DrawingContext *ctx) void drawAttractors(DrawingContext *ctx) { - vector_t p[3][3]; - vector_t l[3][3]; + int n = 3; + vector_t p[6][3]; + vector_t l[6][3]; fixedPoints(ctx, "abc", p[0]); - fixedPoints(ctx, "bca", p[1]); // (bca, cab) -> (cab, cacabac) -> (cacabac, acacbacac) -> - fixedPoints(ctx, "cab", p[2]); // ac cab ca = bca, ac bac ca = acb + fixedPoints(ctx, "bca", p[1]); + fixedPoints(ctx, "cab", p[2]); + fixedPoints(ctx, "a cab a", p[3]); + fixedPoints(ctx, "b abc b", p[4]); + fixedPoints(ctx, "c bca c", p[5]); - double color[3][3] = {{1,0,0},{0,0.7,0},{0,0,1}}; + double color[6][3] = {{1,0,0},{0,0.7,0},{0,0,1},{0,1,1},{0,1,1},{0,1,1}}; - LOOP(i) LOOP(j) l[i][j] = cross(p[i][(3-j)%3], p[i][(4-j)%3]); + for(int i = 0; i < n; i++) LOOP(j) l[i][j] = cross(p[i][(3-j)%3], p[i][(4-j)%3]); - LOOP(i) LOOP(j) if(i != 1) { + for(int i = 0; i < n; i++) LOOP(j) { cairo_set_source_rgb(ctx->cairo, color[i][0], color[i][1], color[i][2]); drawVector(ctx, p[i][j]); } - LOOP(i) LOOP(j) if(i != 1) { + for(int i = 0; i < n; i++) LOOP(j) { cairo_set_source_rgb(ctx->cairo, color[i][0], color[i][1], color[i][2]); drawCovector(ctx, l[i][j]); } } +void drawCurvedBox(DrawingContext *ctx) +{ + vector_t p[6][3]; + vector_t l[2][3]; + + fixedPoints(ctx, "abc", p[0]); + fixedPoints(ctx, "bca", p[1]); + fixedPoints(ctx, "babcb", p[2]); + fixedPoints(ctx, "ababcba", p[3]); + fixedPoints(ctx, "bacac abc cacab", p[4]); + fixedPoints(ctx, "abacac abc cacaba", p[5]); + + LOOP(j) l[0][j] = cross(p[0][(3-j)%3], p[0][(4-j)%3]); + LOOP(j) l[1][j] = cross(p[1][(3-j)%3], p[1][(4-j)%3]); + + drawArc(ctx, "ab", p[0][0], VT_POINT, p[2][0], VT_POINT, p[1][0], 0); + drawArc(ctx, "bcab", p[2][0], VT_POINT, l[1][0], VT_LINE, p[4][0], 1); + drawArc(ctx, "ab", p[1][0], VT_POINT, p[3][0], VT_POINT, p[0][0], 0); + drawArc(ctx, "abcaba", p[3][0], VT_POINT, l[0][0], VT_LINE, p[5][0], 1); + + drawCovector(ctx, l[0][0]); + drawCovector(ctx, l[1][0]); +} + + void drawBoxes(DrawingContext *ctx) { + gsl_matrix *rot = getTempMatrix(ctx->ws); + gsl_matrix **gen = getTempMatrices(ctx->ws, 3); cairo_t *C = ctx->cairo; + vector_t p[6][3]; + vector_t l[6][3]; + vector_t alpha[6]; + + fixedPoints(ctx, "abc", p[0]); + fixedPoints(ctx, "bca", p[1]); + fixedPoints(ctx, "cab", p[2]); + + initializeTriangleGenerators(gen, ctx->cartan); + + for(int i = 0; i < 6; i++) LOOP(j) l[i][j] = cross(p[i][(3-j)%3], p[i][(4-j)%3]); + LOOP(i) LOOP(j) alpha[i].x[j] = gsl_matrix_get(ctx->cartan, i, j); + + cairo_set_source_rgb(C, 0.5, 0.5, 0.5); + + drawCurvedBox(ctx); + +/* + drawRotationOrbit(ctx, "ca", p[0][0]); + drawRotationOrbit(ctx, "ab", p[1][0]); + drawRotationOrbit(ctx, "bc", p[2][0]); + drawRotationOrbit(ctx, "cbac", p[4][0]); + drawRotationOrbit(ctx, "abca", p[5][0]); +*/ + /* - cairo_set_source_rgb(C, 1, 0, 0); - drawTriangle(ctx, "abc"); - cairo_set_source_rgb(C, 0, 0, 1); - drawTriangle(ctx, "aca abc aca"); - drawTriangle(ctx, "acac abc caca"); - drawTriangle(ctx, "acaca abc acaca"); - cairo_set_source_rgb(C, 0, 0.8, 0); - drawTriangle(ctx, "cac abc cac"); - drawTriangle(ctx, "caca abc acac"); - drawTriangle(ctx, "cacac abc cacac"); + cairo_set_source_rgb(C, 0.0, 0.0, 1.0); + drawDualRotationOrbit(ctx, "ca", l[0][2]); + drawDualRotationOrbit(ctx, "ab", l[1][2]); + drawDualRotationOrbit(ctx, "bc", l[2][2]); + drawDualRotationOrbit(ctx, "bacb", l[3][2]); + drawDualRotationOrbit(ctx, "cbac", l[4][2]); + drawDualRotationOrbit(ctx, "abca", l[5][2]); */ - cairo_set_source_rgb(C, 1, 0, 0); - drawBoxStd(ctx, "c", 'C'); - drawBoxStd(ctx, "", 'B'); - drawBoxStd(ctx, "a", 'A'); - drawBoxStd(ctx, "", 'C'); - drawBoxStd(ctx, "b", 'B'); + /* + cairo_set_source_rgb(C, 0, 0, 1); + fixedPoints(ctx, "abc", p[0]); // + fixedPoints(ctx, "ababcba", p[1]); + drawRotationOrbit(ctx, "abcaba", p[0][0]); // abcaba a rar- r-ar a rar- a + drawRotationOrbit(ctx, "ab", p[0][0]); // ab a rar- a rar- rar- a + drawRotationOrbit(ctx, "abcabcba", p[0][0]); // abc abc ba + drawCovector(ctx, cross(p[0][0],p[0][1])); + drawCovector(ctx, cross(p[1][0],p[1][1])); - cairo_set_source_rgb(C, 0, 0, 0); - drawBoxStd(ctx, "ca", 'A'); - drawBoxStd(ctx, "cac", 'C'); - drawBoxStd(ctx, "caca", 'A'); - drawBoxStd(ctx, "acac", 'C'); - drawBoxStd(ctx, "aca", 'A'); - drawBoxStd(ctx, "ac", 'C'); + cairo_set_source_rgb(C, 0, 0.7, 0); + fixedPoints(ctx, "ab abc ba", p[0]); // + fixedPoints(ctx, "ab ababcba ba", p[1]); + drawRotationOrbit(ctx, "ababcababa", p[0][0]); // abcaba a rar- r-ar a rar- a + drawRotationOrbit(ctx, "ab", p[0][0]); // ab a rar- a rar- rar- a + drawRotationOrbit(ctx, "ababcabcbaba", p[0][0]); // abc abc ba + drawCovector(ctx, cross(p[0][0],p[0][1])); + drawCovector(ctx, cross(p[1][0],p[1][1])); + */ - drawBoxStd(ctx, "aca cb", 'B'); - drawBoxStd(ctx, "aca cbc", 'C'); - drawBoxStd(ctx, "aca cbcb", 'B'); - drawBoxStd(ctx, "aca bcbc", 'C'); - drawBoxStd(ctx, "aca bcb", 'B'); - drawBoxStd(ctx, "aca bc", 'C'); + cairo_set_source_rgb(C, 0, 0, 1.0); +// drawRotationOrbit(ctx, "bcab", p[3][0]); +// drawRotationOrbit(ctx, "abca", p[5][0]); +// drawRotationOrbit(ctx, "bacb", p[3][0]); - drawBoxStd(ctx, "caca cb", 'B'); - drawBoxStd(ctx, "caca cbc", 'C'); - drawBoxStd(ctx, "caca cbcb", 'B'); - drawBoxStd(ctx, "caca bcbc", 'C'); - drawBoxStd(ctx, "caca bcb", 'B'); - drawBoxStd(ctx, "caca bc", 'C'); +// LOOP(i) drawVector(ctx, p[i+3][0]); - cairo_set_source_rgb(C, 1, 0, 1); - drawBoxStd(ctx, "ca bc", 'C'); - drawBoxStd(ctx, "ca bcb", 'B'); - drawBoxStd(ctx, "ca bcbc", 'C'); - drawBoxStd(ctx, "ca cbcb", 'B'); - drawBoxStd(ctx, "ca cbc", 'C'); - drawBoxStd(ctx, "ca cb", 'B'); - - cairo_set_source_rgb(C, 0, 1, 0); -// drawBoxStd(ctx, "ca bc", 'C'); - drawBoxStd(ctx, "cabc ba", 'A'); - drawBoxStd(ctx, "cabc bab", 'B'); - drawBoxStd(ctx, "cabc baba", 'A'); - drawBoxStd(ctx, "cabc abab", 'B'); - drawBoxStd(ctx, "cabc aba", 'A'); - drawBoxStd(ctx, "cabc ab", 'B'); + cairo_set_source_rgb(C, 0.5, 0.5, 0.5); +// drawCovector(ctx, alpha[1]); +// drawCovector(ctx, alpha[2]); + /* + alpha[3] = apply_transpose(gen[1], alpha[0]); + drawCovector(ctx, alpha[3]); + */ + releaseTempMatrices(ctx->ws, 4); } void drawBoxes2(DrawingContext *ctx) @@ -372,77 +600,236 @@ void drawBoxes2(DrawingContext *ctx) cairo_t *C = ctx->cairo; initializeTriangleGenerators(gen, ctx->cartan); + // abc, ababcba, abababcbaba, ..., cab + // bca, acaba, abacababa, ..., babcb + + vector_t v[4][3]; + vector_t i[4]; + + fixedPoints(ctx, "abc", v[0]); + fixedPoints(ctx, "bca", v[1]); + fixedPoints(ctx, "cab", v[2]); + fixedPoints(ctx, "acaba", v[3]); + i[0] = cross(cross(v[0][1],v[0][2]),cross(v[2][0],v[2][2])); + i[1] = cross(cross(v[1][1],v[1][2]),cross(v[3][0],v[3][2])); + i[2] = cross(cross(v[2][0],v[2][2]),cross(v[3][0],v[3][2])); + i[3] = cross(cross(v[0][0],v[0][2]),cross(v[1][0],v[1][2])); + + cairo_set_source_rgb(C, 0, 0, 1); +// drawPolygon(ctx, 1, 6, v[2][2], v[1][1], v[0][1], v[3][2], v[3][1], v[2][1]); +// drawPolygon(ctx, 1, 6, v[1][2], i[1], i[2], i[0], v[0][2], i[3]); + +/* + cairo_set_source_rgb(C, 1, 0, 0); + drawBox(ctx, "cab", "bca"); cairo_set_source_rgb(C, 1, 0.5, 0); -// drawBox(ctx, "caca cab acac", "caca bca acac"); -// drawBox(ctx, "ca cab ac", "ca bca ac"); // ac abc ca -// drawBox(ctx, "aca cab aca", "aca bca aca"); -// drawBox(ctx, "a cab a", "a bca a"); + drawBox(ctx, "bc bca cb", "bc abc cb"); + drawBox(ctx, "bcbc bca cbcb", "bcbc abc cbcb"); + drawBox(ctx, "bcbcbc bca cbcbcb", "bcbcbc abc cbcbcb"); + drawBox(ctx, "bcbcbcbc bca cbcbcbcb", "bcbcbcbc abc cbcbcbcb"); + drawBox(ctx, "bcbcbcbcbc bca cbcbcbcbcb", "bcbcbcbcbc abc cbcbcbcbcb"); +*/ - vector_t fp[3], fp2[3]; - vector_t w; + cairo_set_source_rgb(C, 1, 0, 0); + drawBox(ctx, "bca", "abc"); + cairo_set_source_rgb(C, 1, 0.5, 0); + drawBox(ctx, "ab abc ba", "ab cab ba"); + drawBox(ctx, "abab abc baba", "abab cab baba"); + drawBox(ctx, "ababab abc bababa", "ababab cab bababa"); + drawBox(ctx, "abababab abc babababa", "abababab cab babababa"); + drawBox(ctx, "ababababab abc bababababa", "ababababab cab bababababa"); -// fixedPoints(ctx, "ac abc ca", fp); // -// drawVector(ctx, fp[0]); - fixedPoints(ctx, "cab", fp); - fixedPoints(ctx, "abc", fp2); - w = cross(cross(fp[0], fp[2]), cross(fp2[0], fp2[1])); - drawRotationOrbit(ctx, "ac", w); -// drawVector(ctx, cross(cross(fp[0], fp[1]), cross(fp2[0], fp2[2]))); -// drawVector(ctx, fp[0]); + cairo_set_source_rgb(C, 1, 0, 0); + drawBox(ctx, "ab ca cab ac ba", "ab ca bca ac ba"); + drawBox(ctx, "ab caca cab acac ba", "ab caca bca acac ba"); + drawBox(ctx, "ab cacaca cab acacac ba", "ab cacaca bca acacac ba"); + drawBox(ctx, "ab cacacaca cab acacacac ba", "ab cacacaca bca acacacac ba"); + drawBox(ctx, "ab cacacacaca cab acacacacac ba", "ab cacacacaca bca acacacacac ba"); + /* + drawBox(ctx, "ab abc ba", "ab cab ba"); + drawBox(ctx, "abab abc baba", "abab cab baba"); + drawBox(ctx, "ababab abc bababa", "ababab cab bababa"); + drawBox(ctx, "abababab abc babababa", "abababab cab babababa"); + drawBox(ctx, "ababababab abc bababababa", "ababababab cab bababababa"); + */ - cairo_set_source_rgb(C, 1, 0, 0.5); -// drawBoxLines(ctx, "cab", "bca"); + // abc -> arar abc r-ar-a = (ar)³ + // bca -> (ar)² (ra)³ (ar)-² -> a rarr a rar- r-ar rar- a = ababcba + // a rar- r-ar a r-ar a rar- a r-ar ra-r a = abcacabacba - cairo_set_source_rgb(C, 0, 0, 0); -// drawBox(ctx, "abc", "cab"); +// drawBox(ctx, "ab ca cab ac ba", "ab ca bca ac ba"); +// drawBox(ctx, "ab ca bc bca cb ac ba", "ab ca bc abc cb ac ba"); +// drawBox(ctx, "ab ca bc ab abc ba cb ac ba", "ab ca bc ab cab ba cb ac ba"); -// fixedPoints(ctx, "abc", fp); - cairo_set_source_rgb(C, 0.7, 0.7, 0.7); - drawRotationOrbit(ctx, "ac", fp2[0]); -// drawRotationOrbit(ctx, "ac", fp[1]); -// drawRotationOrbit(ctx, "ac", fp[2]); -// drawRotationOrbit(ctx, "ab", fp[0]); -// drawRotationOrbit(ctx, "bc", fp[0]); + // these are the right boxes + /* + cairo_set_source_rgb(C, 1, 0.5, 0); + drawBox(ctx, "abababab abc babababa", "abababab cab babababa"); + drawBox(ctx, "ababab abc bababa", "ababab cab bababa"); + drawBox(ctx, "abab abc baba", "abab cab baba"); + drawBox(ctx, "ab abc ba", "ab cab ba"); + drawBox(ctx, "b abc b", "b cab b"); + drawBox(ctx, "bab abc bab", "bab cab bab"); + drawBox(ctx, "babab abc babab", "babab cab babab"); + drawBox(ctx, "bababab abc bababab", "bababab cab bababab"); + */ + + +// drawBox(ctx, "ab a cab a ba", "ab a bca a ba"); +// drawBox(ctx, "ab ca cab ac ba", "ab ca bca ac ba"); + +// drawBox(ctx, "cab", "bca"); +// drawBox(ctx, "ca cab ac", "ca bca ac"); releaseTempMatrices(ctx->ws, 4); } +void drawRotatedReflectors(DrawingContext *ctx) +{ + gsl_matrix *rot = getTempMatrix(ctx->ws); + gsl_matrix **gen = getTempMatrices(ctx->ws, 3); + cairo_t *C = ctx->cairo; + vector_t fp[3], fp2[3]; + vector_t w; + vector_t v[3]; + + cairo_set_source_rgb(C, 0.7, 0.7, 0.7); + + initializeTriangleGenerators(gen, ctx->cartan); + + LOOP(i) LOOP(j) v[i].x[j] = gsl_matrix_get(ctx->cartan, i, j); + multiply(gen[0], gen[1], rot); + + for(int j = 0; j < ctx->p[2]; j++) { + drawCovector(ctx, v[0]); + v[0] = apply_transpose(rot, v[0]); + } + + LOOP(i) LOOP(j) { v[i].x[j] = (i==j) ? 1.0 : 0.0; } + + for(int j = 0; j < ctx->p[2]; j++) { + drawVector(ctx, v[0]); + v[0] = apply(rot, v[0]); + } + + fixedPoints(ctx, "cab", fp); + fixedPoints(ctx, "cacabac", fp2); + drawRotationOrbit(ctx, "ac", fp[0]); + + releaseTempMatrices(ctx->ws, 4); +} + +void drawDualLimitCurve(DrawingContext *ctx) +{ + cairo_t *C = ctx->cairo; + + cairo_save(C); + cairo_set_source_rgb(C, 0, 0, 0); + + int n = 18; + vector_t p[18][3]; + vector_t l[n][3]; + + /* + fixedPoints(ctx, "abc", p[0]); + fixedPoints(ctx, "ab abc ba", p[1]); + fixedPoints(ctx, "abab abc baba", p[2]); + fixedPoints(ctx, "ababab abc bababa", p[3]); + fixedPoints(ctx, "abababab abc babababa", p[4]); + fixedPoints(ctx, "babababa abc abababab", p[5]); + fixedPoints(ctx, "bababa abc ababab", p[6]); + fixedPoints(ctx, "baba abc abab", p[7]); + fixedPoints(ctx, "ba abc ab", p[8]); + + fixedPoints(ctx, "bca", p[9]); + fixedPoints(ctx, "ab bca ba", p[10]); + fixedPoints(ctx, "abab bca baba", p[11]); + fixedPoints(ctx, "ababab bca bababa", p[12]); + fixedPoints(ctx, "abababab bca babababa", p[13]); + fixedPoints(ctx, "babababa bca abababab", p[14]); + fixedPoints(ctx, "bababa bca ababab", p[15]); + fixedPoints(ctx, "baba bca abab", p[16]); + fixedPoints(ctx, "ba bca ab", p[17]); + */ + + fixedPoints(ctx, "abc", p[0]); + fixedPoints(ctx, "ac abc ca", p[1]); + fixedPoints(ctx, "acac abc caca", p[2]); + fixedPoints(ctx, "acacac abc cacaca", p[3]); + fixedPoints(ctx, "acacacac abc cacacaca", p[4]); + fixedPoints(ctx, "cacacaca abc acacacac", p[5]); + fixedPoints(ctx, "cacaca abc acacac", p[6]); + fixedPoints(ctx, "caca abc acac", p[7]); + fixedPoints(ctx, "ca abc ac", p[8]); + + fixedPoints(ctx, "bca", p[9]); + fixedPoints(ctx, "ac bca ca", p[10]); + fixedPoints(ctx, "acac bca caca", p[11]); + fixedPoints(ctx, "acacac bca cacaca", p[12]); + fixedPoints(ctx, "acacacac bca cacacaca", p[13]); + fixedPoints(ctx, "cacacaca bca acacacac", p[14]); + fixedPoints(ctx, "cacaca bca acacac", p[15]); + fixedPoints(ctx, "caca bca acac", p[16]); + fixedPoints(ctx, "ca bca ac", p[17]); + + + /* + fixedPoints(ctx, "cab", p[2]); + fixedPoints(ctx, "b abc b", p[3]); + fixedPoints(ctx, "c bca c", p[4]); + fixedPoints(ctx, "a cab a", p[5]); + */ + + for(int i = 0; i < n; i++) { + LOOP(j) l[i][j] = cross(p[i][(3-j)%3], p[i][(4-j)%3]); +// drawCovector(ctx, l[i][0]); + drawCovector(ctx, l[i][2]); + } + + cairo_restore(C); +} + void drawLimitCurve(DrawingContext *ctx) { cairo_t *C = ctx->cairo; cairo_save(C); - int previous_inside = 0; - for(int i = 0; i < ctx->n_group_elements; i++) { - point_t p; - p.x = ctx->limit_curve[3*i]; - p.y = ctx->limit_curve[3*i+1]; + cairo_set_source_rgb(C, 0, 0, 0); - if(isInsideBB(ctx, p)) { - if(ctx->limit_with_lines) { + if(ctx->limit_with_lines) { + int previous_inside = 0; + for(int i = 0; i < ctx->limit_curve_count; i++) { + point_t p; + p.x = ctx->limit_curve[3*i]; + p.y = ctx->limit_curve[3*i+1]; + + if(isInsideBB(ctx, p)) { if(!previous_inside) cairo_move_to(C, p.x, p.y); else cairo_line_to(C, p.x, p.y); + previous_inside = 1; } else { - cairo_move_to(C, p.x, p.y); - cairo_close_path(C); + previous_inside = 0; + } + } + + cairo_stroke(C); + } else { + for(int i = 0; i < ctx->limit_curve_count; i++) { + point_t p; + p.x = ctx->limit_curve[3*i]; + p.y = ctx->limit_curve[3*i+1]; + + if(isInsideBB(ctx, p)) { + cairo_arc(C, p.x, p.y, 2.0/ctx->dim->scalefactor, 0, 2*M_PI); + cairo_close_path(C); + cairo_fill(C); } - previous_inside = 1; - } else { - previous_inside = 0; } } - if(!ctx->limit_with_lines) { // draw dots instead of lines - cairo_set_line_cap(C, CAIRO_LINE_CAP_ROUND); - cairo_set_line_width(C, 3.0/ctx->dim->scalefactor); - } - - cairo_set_source_rgb(C, 0, 0, 0); - cairo_stroke(C); - cairo_restore(C); } @@ -472,26 +859,34 @@ void draw(DrawingContext *ctx) cairo_set_line_join(C, CAIRO_LINE_JOIN_BEVEL); cairo_set_line_cap(C, CAIRO_LINE_CAP_ROUND); - if(ctx->limit_curve_valid) { + if(ctx->limit_curve_count >= 0) { if(ctx->show_limit) drawLimitCurve(ctx); + if(ctx->show_dual_limit) + drawDualLimitCurve(ctx); + + if(ctx->show_attractors) + drawAttractors(ctx); + + if(ctx->show_rotated_reflectors) + drawRotatedReflectors(ctx); + + if(ctx->show_reflectors) + drawReflectors(ctx); + if(ctx->show_boxes) drawBoxes(ctx); if(ctx->show_boxes2) drawBoxes2(ctx); - if(ctx->show_attractors) - drawAttractors(ctx); - - if(ctx->show_reflectors) - drawReflectors(ctx); } cairo_identity_matrix(C); // text is in screen coordinates - drawText(ctx); + if(ctx->show_text) + drawText(ctx); cairo_surface_flush(cairo_get_target(C)); } diff --git a/limit_set.c b/limit_set.c index 270098c..bbe648a 100644 --- a/limit_set.c +++ b/limit_set.c @@ -7,22 +7,23 @@ static int compareAngle(const void *x, const void *y) void cartanMatrix(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, 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, 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); + 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 initializeTriangleGenerators(gsl_matrix **gen, gsl_matrix *cartan) { LOOP(i) gsl_matrix_set_identity(gen[i]); + LOOP(i) LOOP(j) *gsl_matrix_ptr(gen[i], j, j) = -1.0; LOOP(i) LOOP(j) *gsl_matrix_ptr(gen[i], i, j) += gsl_matrix_get(cartan, i, j); } @@ -43,7 +44,8 @@ int computeLimitCurve(DrawingContext *ctx) int success = 0; int column = ctx->use_repelling ? 2 : 0; - ctx->limit_curve_valid = 0; +// int column = 1; + ctx->limit_curve_count = -1; // do first in the Fuchsian positive case to get the angles cartanMatrix(cartan_pos, M_PI/ctx->p[0], M_PI/ctx->p[1], M_PI/ctx->p[2], 1.0); @@ -61,8 +63,8 @@ int computeLimitCurve(DrawingContext *ctx) for(int i = 0; i < ctx->n_group_elements; i++) { multiply_many(ws, fixedpoints_pos, 3, cob_pos, elements[i], coxeter_fixedpoints_pos); ctx->limit_curve[3*i+2] = atan2( - gsl_matrix_get(fixedpoints_pos, 0, column)/gsl_matrix_get(fixedpoints_pos, 2, column), - gsl_matrix_get(fixedpoints_pos, 1, column)/gsl_matrix_get(fixedpoints_pos, 2, column)); + gsl_matrix_get(fixedpoints_pos, 2, column)/gsl_matrix_get(fixedpoints_pos, 0, column), + gsl_matrix_get(fixedpoints_pos, 1, column)/gsl_matrix_get(fixedpoints_pos, 0, column)); } // now do it again to calculate x and y coordinates @@ -70,6 +72,7 @@ int computeLimitCurve(DrawingContext *ctx) gsl_matrix_set_identity(elements[0]); for(int i = 1; i < ctx->n_group_elements; i++) multiply(elements[group[i].parent->id], gen[group[i].letter], elements[i]); + multiply_many(ws, coxeter, 3, gen[0], gen[1], gen[2]); int ev_count = real_eigenvectors(coxeter, coxeter_fixedpoints, ws); @@ -87,7 +90,7 @@ int computeLimitCurve(DrawingContext *ctx) qsort(ctx->limit_curve, ctx->n_group_elements, 3*sizeof(double), compareAngle); - ctx->limit_curve_valid = 1; + ctx->limit_curve_count = ctx->n_group_elements; success = 1; diff --git a/linalg.h b/linalg.h index 38477bb..c8c75fe 100644 --- a/linalg.h +++ b/linalg.h @@ -13,7 +13,7 @@ #define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);} #define FCMP(x, y) gsl_fcmp(x, y, 1e-10) -#define MAX_TEMP_MATRICES 60000 +#define MAX_TEMP_MATRICES 600000 #define MAX_TEMP_VECTORS 100 typedef struct _workspace { diff --git a/main.c b/main.c index 6cc250d..057ad25 100644 --- a/main.c +++ b/main.c @@ -18,24 +18,28 @@ DrawingContext *screen_context; void setupContext(DrawingContext *ctx) { ctx->n_group_elements = NUM_GROUP_ELEMENTS; - ctx->p[0] = 5; - ctx->p[1] = 5; - ctx->p[2] = 5; - ctx->k[0] = 2; - ctx->k[1] = 2; - ctx->k[2] = 2; - ctx->parameter = 3.0; + ctx->p[0] = 9; + ctx->p[1] = 9; + ctx->p[2] = 9; + ctx->k[0] = 4; + ctx->k[1] = 4; + ctx->k[2] = 4; + ctx->parameter = 5.35; +// ctx->parameter = 0.1; ctx->show_boxes = 0; ctx->show_boxes2 = 0; ctx->show_attractors = 0; ctx->show_reflectors = 0; + ctx->show_rotated_reflectors = 0; ctx->show_limit= 1; + ctx->show_dual_limit= 0; + ctx->show_text = 1; ctx->use_rotation_basis = 0; ctx->limit_with_lines = 1; ctx->use_repelling = 0; ctx->limit_curve = malloc(3*ctx->n_group_elements*sizeof(double)); - ctx->limit_curve_valid = 0; + ctx->limit_curve_count = -1; ctx->group = malloc(ctx->n_group_elements*sizeof(groupelement_t)); generate_triangle_group(ctx->group, ctx->n_group_elements, ctx->p[0], ctx->p[1], ctx->p[2]); @@ -132,13 +136,19 @@ void updateMatrices(DrawingContext *ctx) gsl_matrix *tmp = getTempMatrix(ctx->ws); - if(ctx->use_rotation_basis % 3 == 0) { + if(ctx->use_rotation_basis % 4 == 0) { gsl_matrix_memcpy(ctx->cob, ctx->cartan); // is this a good choice of basis for any reason? - } else if(ctx->use_rotation_basis % 3 == 1) { - computeRotationMatrix(ctx, tmp, "ac"); +// gsl_matrix_set_identity(ctx->cob); // is this a good choice of basis for any reason? + } else if(ctx->use_rotation_basis % 4 == 1) { + computeRotationMatrix(ctx, tmp, "ba"); invert(tmp, ctx->cob, ctx->ws); + } else if(ctx->use_rotation_basis % 4 == 2) { + computeBoxTransform(ctx, "bca", "abc", ctx->cob); +// computeBoxTransform(ctx, "cab", "bca", ctx->cob); +// computeBoxTransform(ctx, "acb", "cba", ctx->cob); } else { - computeBoxTransform(ctx, "abc", "cab", ctx->cob); + cartanMatrix(tmp, M_PI/ctx->p[0], M_PI/ctx->p[1], M_PI/ctx->p[2], 1.0); + diagonalize_symmetric_form(tmp, ctx->cob, ctx->ws); } releaseTempMatrices(ctx->ws, 1); @@ -178,6 +188,7 @@ int processEvent(GraphicsInfo *info, XEvent *ev) { int state; unsigned long key; + char filename[100]; switch(ev->type) { @@ -198,12 +209,12 @@ int processEvent(GraphicsInfo *info, XEvent *ev) computeLimitCurve(screen_context); break; case XK_Left: - screen_context->parameter /= exp(0.0001); + screen_context->parameter /= exp(0.00002); updateMatrices(screen_context); computeLimitCurve(screen_context); break; case XK_Right: - screen_context->parameter *= exp(0.0001); + screen_context->parameter *= exp(0.00002); updateMatrices(screen_context); computeLimitCurve(screen_context); break; @@ -218,12 +229,13 @@ int processEvent(GraphicsInfo *info, XEvent *ev) computeLimitCurve(screen_context); break; case ' ': - screen_context->parameter = 2.890053638; + screen_context->parameter = 5.57959706; updateMatrices(screen_context); computeLimitCurve(screen_context); break; case XK_Return: - screen_context->parameter = 2.76375163; +// screen_context->parameter = 2.76375163; + screen_context->parameter = 5.29063366; updateMatrices(screen_context); computeLimitCurve(screen_context); break; @@ -247,12 +259,18 @@ int processEvent(GraphicsInfo *info, XEvent *ev) case 'r': TOGGLE(screen_context->show_reflectors); break; + case 'x': + TOGGLE(screen_context->show_rotated_reflectors); + break; case 'L': TOGGLE(screen_context->limit_with_lines); break; case 'l': TOGGLE(screen_context->show_limit); break; + case 'd': + TOGGLE(screen_context->show_dual_limit); + break; case 'R': screen_context->use_rotation_basis++; updateMatrices(screen_context); @@ -260,10 +278,26 @@ int processEvent(GraphicsInfo *info, XEvent *ev) break; case 'p': print(screen_context); + /* + screen_context->limit_with_lines = 0; + for(int i = 0; i <= 800; i++) { + screen_context->parameter = exp(0.005*i-2); + updateMatrices(screen_context); + computeLimitCurve(screen_context); + draw(screen_context); + sprintf(filename, "test%03d.png", i); + cairo_surface_write_to_png(info->buffer_surface, filename); + printf("Finished drawing %s\n", filename); + } + */ break; case 'f': TOGGLE(screen_context->use_repelling); computeLimitCurve(screen_context); + break; + case 't': + TOGGLE(screen_context->show_text); + break; } return STATUS_REDRAW; @@ -286,20 +320,20 @@ int main() return 1; /* - info->dim->matrix.xx = 837.930824; - info->dim->matrix.xy = -712.651341; - info->dim->matrix.x0 = 180.427716; - info->dim->matrix.yx = 712.651341; - info->dim->matrix.yy = 837.930824; - info->dim->matrix.y0 = 1412.553240; - */ - info->dim->matrix.xx = 112.465171; info->dim->matrix.xy = 0.000000; info->dim->matrix.x0 = 891.180490; info->dim->matrix.yx = 0.000000; info->dim->matrix.yy = 112.465171; info->dim->matrix.y0 = 506.676280; + */ + + info->dim->matrix.xx = 274.573171; + info->dim->matrix.xy = 0.000000; + info->dim->matrix.x0 = 583.073462; + info->dim->matrix.yx = 0.000000; + info->dim->matrix.yy = 274.573171; + info->dim->matrix.y0 = 777.225293; updateDimensions(info->dim); diff --git a/main.h b/main.h index bd9143d..ba6a8dd 100644 --- a/main.h +++ b/main.h @@ -35,7 +35,10 @@ typedef struct { int show_boxes2; int show_attractors; int show_reflectors; + int show_rotated_reflectors; int show_limit; + int show_dual_limit; + int show_text; int use_rotation_basis; int limit_with_lines; int use_repelling; @@ -46,12 +49,17 @@ typedef struct { // computed stuff double *limit_curve; // x, y, angle triples - int limit_curve_valid; + int limit_curve_count; // temporary; matrices can only be freed from the top, but that's enough for us workspace_t *ws; } DrawingContext; +typedef enum { + VT_POINT, + VT_LINE +} vector_type_t; + // implemented in limit_set.c void cartanMatrix(gsl_matrix *cartan, double a1, double a2, double a3, double s); void initializeTriangleGenerators(gsl_matrix **gen, gsl_matrix *cartan); diff --git a/queue.h b/queue.h index 88ba8d0..ea5e9e5 100644 --- a/queue.h +++ b/queue.h @@ -4,7 +4,7 @@ #include #include -#define QUEUE_SIZE 100000 +#define QUEUE_SIZE 1000000 #define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);} #ifdef _DEBUG