From 4daaf1ed7ca80f6fc95e76683edf12b2ae587478 Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Sun, 24 Feb 2019 08:43:52 +0100 Subject: [PATCH] start drawing rotation orbits --- draw.c | 128 +++++++++++++++++++++---------------------------------- linalg.c | 64 ++++++++++++++++++++++++++++ linalg.h | 3 +- main.c | 85 ++++++++++++++++++++++++++++-------- main.h | 14 ++++++ 5 files changed, 196 insertions(+), 98 deletions(-) diff --git a/draw.c b/draw.c index 564aa83..bf784d0 100644 --- a/draw.c +++ b/draw.c @@ -2,6 +2,12 @@ // level 0: helper functions +int isInsideBB(DrawingContext *ctx, point_t p) +{ + cairo_user_to_device(ctx->cairo, &x, &y); + return -p.x < ctx->dim->width && p.x < 3*ctx->dim->width && -p.y < ctx->dim->height && p.y < 3*ctx->dim->height; +} + vector_t cross(vector_t a, vector_t b) { vector_t result; @@ -47,26 +53,6 @@ int fixedPoints(DrawingContext *ctx, const char *word, vector_t *out) return count; } -void transformFrameStd(DrawingContext *ctx, vector_t *x, gsl_matrix *out) -{ - gsl_matrix *tmp = getTempMatrix(ctx->ws); - gsl_vector *fourth = getTempVector(ctx->ws); - gsl_vector *lambda = getTempVector(ctx->ws); - int s; - - LOOP(i) LOOP(j) gsl_matrix_set(out, j, i, x[i].x[j]); - gsl_matrix_memcpy(tmp, out); - gsl_linalg_LU_decomp(tmp, ctx->ws->permutation, &s); - gsl_linalg_LU_solve(tmp, ctx->ws->permutation, fourth, lambda); - - LOOP(i) LOOP(j) *gsl_matrix_ptr(out, i, j) *= gsl_vector_get(lambda, j); - - gsl_matrix_fprintf(stdout, out, "%f"); - - releaseTempMatrices(ctx->ws, 1); - releaseTempVectors(ctx->ws, 2); -} - // level 1: the elementary drawing functions, drawPoint, drawSegment2d void drawPoint(DrawingContext *ctx, point_t p) @@ -257,8 +243,8 @@ void drawAttractors(DrawingContext *ctx) vector_t l[3][3]; fixedPoints(ctx, "abc", p[0]); - fixedPoints(ctx, "bca", p[1]); - fixedPoints(ctx, "cab", p[2]); + fixedPoints(ctx, "bca", p[1]); // (bca, cab) -> (cab, cacabac) -> (cacabac, acacbacac) -> + fixedPoints(ctx, "cab", p[2]); // ac cab ca = bca, ac bac ca = acb double color[3][3] = {{1,0,0},{0,0.7,0},{0,0,1}}; @@ -341,71 +327,55 @@ void drawBoxes(DrawingContext *ctx) void drawBoxes2(DrawingContext *ctx) { - /* - cairo_set_source_rgb(ctx->cairo, 1, 0.5, 0); - drawBoxStd(ctx, "a", 'A'); - drawBoxStd(ctx, "ca", 'A'); - drawBoxStd(ctx, "aca", 'A'); - drawBoxStd(ctx, "caca", 'A'); - drawBoxStd(ctx, "acaca", 'A'); - cairo_set_source_rgb(ctx->cairo, 0, 0.5, 1); - drawBoxStd(ctx, "acac", 'A'); - drawBoxStd(ctx, "cac", 'A'); - drawBoxStd(ctx, "ac", 'A'); - drawBoxStd(ctx, "c", 'A'); - drawBoxStd(ctx, "", 'A'); - */ - - /* - cairo_set_source_rgb(ctx->cairo, 0, 0.5, 1); - drawBox(ctx, "abc", "cab"); -// drawBox(ctx, "ba abc ab", "ba cab ab"); - drawBox(ctx, "baba abc abab", "baba cab abab"); - drawBox(ctx, "abab abc baba", "abab cab baba"); - drawBox(ctx, "ab abc ba", "ab cab ba"); + gsl_matrix *tmp = getTempMatrix(ctx->ws); cairo_set_source_rgb(ctx->cairo, 1, 0.5, 0); - 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, "aba abc aba", "aba cab aba"); -// drawBox(ctx, "a abc a", "a cab a"); + drawBox(ctx, "ac cab ca", "ac bca ca"); + drawBox(ctx, "acac cab caca", "acac bca caca"); + drawBox(ctx, "caca cab acac", "caca bca acac"); + drawBox(ctx, "ca cab ac", "ca bca ac"); // cacabac, cab - cairo_set_source_rgb(ctx->cairo, 0, 0.5, 1); - drawBox(ctx, "bca", "abc"); -// drawBox(ctx, "cb bca bc", "cb abc bc"); - drawBox(ctx, "cbcb bca bcbc", "cbcb abc bcbc"); - drawBox(ctx, "bcbc bca cbcb", "bcbc abc cbcb"); - drawBox(ctx, "bc bca cb", "bc abc cb"); - - cairo_set_source_rgb(ctx->cairo, 1, 0.5, 0); - drawBox(ctx, "c bca c", "c abc c"); - drawBox(ctx, "cbc bca cbc", "cbc abc cbc"); - drawBox(ctx, "cbcbc bca cbcbc", "cbcbc abc cbcbc"); - drawBox(ctx, "bcb bca bcb", "bcb abc bcb"); -// drawBox(ctx, "b bca b", "b abc b"); -*/ - -// cairo_set_source_rgb(ctx->cairo, 1, 0, 0); -// drawBoxLines(ctx, "bca", "abc"); - - cairo_set_source_rgb(ctx->cairo, 1, 0.5, 0); - drawBoxLines(ctx, "ac cab ca", "ac bca ca"); - drawBoxLines(ctx, "acac cab caca", "acac bca caca"); - drawBoxLines(ctx, "caca cab acac", "caca bca acac"); - drawBoxLines(ctx, "ca cab ac", "ca bca ac"); + cairo_set_source_rgb(ctx->cairo, 1, 0, 0.5); drawBoxLines(ctx, "cab", "bca"); cairo_set_source_rgb(ctx->cairo, 0, 0, 0); - drawBoxLines(ctx, "abc", "cab"); + drawBox(ctx, "abc", "cab"); - cairo_set_source_rgb(ctx->cairo, 1, 0.5, 0); -// drawBox(ctx, "a cab a", "a bca a"); -// drawBox(ctx, "aca cab aca", "aca bca aca"); -// drawBox(ctx, "acaca cab acaca", "acaca bca acaca"); -// drawBox(ctx, "cac cab cac", "cac bca cac"); -// drawBox(ctx, "c cab c", "c bca c"); + vector_t w = {0, 1, 0}; + drawCovector(ctx, w); + vector_t v[3]; + point_t p; + cairo_set_source_rgb(ctx->cairo, 1, 0, 0); + computeRotationMatrix(ctx, tmp, "ac"); + LOOP(i) LOOP(j) v[i].x[j] = gsl_matrix_get(tmp, j, i); + + for(int k = 0; k < 50; k++) { + LOOP(i) w.x[i] = 1.8 * v[2].x[i] + cos(2*k*M_PI/50) * v[0].x[i] + sin(2*k*M_PI/50) * v[1].x[i]; + p = vectorToPoint(ctx, w); + + if(k == 0) + cairo_move_to(ctx->cairo, p.x, p.y); + else + cairo_line_to(ctx->cairo, p.x, p.y); + } + + cairo_set_source_rgb(ctx->cairo, 1, 0, 0); + cairo_stroke(ctx->cairo); + + /* vector_t p[3]; */ + /* fixedPoints(ctx, "cab", p); */ + /* drawVector(ctx, p[0]); */ + /* fixedPoints(ctx, "ca cab ac", p); */ + /* drawVector(ctx, p[0]); */ + /* fixedPoints(ctx, "caca cab acac", p); */ + /* drawVector(ctx, p[0]); */ + /* fixedPoints(ctx, "acac cab caca", p); */ + /* drawVector(ctx, p[0]); */ + /* fixedPoints(ctx, "ac cab ca", p); */ + /* drawVector(ctx, p[0]); */ + + releaseTempMatrices(ctx->ws, 1); } void drawLimitCurve(DrawingContext *ctx) diff --git a/linalg.c b/linalg.c index 704ca7c..1f0c577 100644 --- a/linalg.c +++ b/linalg.c @@ -4,7 +4,9 @@ #include #include #include +#include #include +#include #include "linalg.h" @@ -293,3 +295,65 @@ int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws) releaseTempMatrices(ws, 1); return positive; } + +// computes a matrix in SL(3, R) which projectively transforms (e1, e2, e3, e1+e2+e3) to the 4 given vectors +void projective_frame(gsl_vector **vertices, gsl_matrix *result, workspace_t *ws) +{ + gsl_matrix *tmp = getTempMatrix(ws); + gsl_vector *coeff = getTempVector(ws); + int s; + double det, scale; + + for(int i = 0; i < 3; i++) + for(int j = 0; j < 3; j++) + gsl_matrix_set(tmp, i, j, gsl_vector_get(vertices[j], i)); + + gsl_linalg_LU_decomp(tmp, ws->permutation, &s); + gsl_linalg_LU_solve(tmp, ws->permutation, vertices[3], coeff); + det = gsl_linalg_LU_det(tmp, s); + + for(int i = 0; i < 3; i++) + det *= gsl_vector_get(coeff, i); + scale = 1/cbrt(det); + + for(int i = 0; i < 3; i++) + for(int j = 0; j < 3; j++) + gsl_matrix_set(result, i, j, scale*gsl_vector_get(vertices[j], i)*gsl_vector_get(coeff, j)); + + releaseTempMatrices(ws, 1); + releaseTempVectors(ws, 1); +} + +void rotation_frame(gsl_matrix *rotation, gsl_matrix *result, workspace_t *ws) +{ + gsl_matrix *tmp = getTempMatrix(ws); + gsl_matrix *rot_basis = getTempMatrix(ws); + + gsl_matrix_memcpy(tmp, rotation); + gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv); + int r = gsl_eigen_nonsymmv(tmp, ws->eval_complex, ws->evec_complex, ws->work_nonsymmv); + ERROR(r, "gsl_eigen_nonsymmv failed!\n"); + + double arg, minarg = 5; // greater than pi + int minidx; + for(int i = 0; i < 3; i++) { + arg = gsl_complex_arg(gsl_vector_complex_get(ws->eval_complex, i)); + if(abs(arg) < minarg) + { + minidx = i; + minarg = abs(arg); + } + } + ERROR(FCMP(minarg, 0.0) != 0, "rotation_frame() failed! No eigenvalue was 1.\n"); + + for(int i = 0; i < 3; i++) { + gsl_complex x = gsl_matrix_complex_get(ws->evec_complex, i, (minidx+1)%3); + gsl_complex y = gsl_matrix_complex_get(ws->evec_complex, i, (minidx+2)%3); + gsl_complex z = gsl_matrix_complex_get(ws->evec_complex, i, minidx); + gsl_matrix_set(result, i, 0, GSL_REAL(x)+GSL_REAL(y)); + gsl_matrix_set(result, i, 1, GSL_IMAG(x)-GSL_IMAG(y)); + gsl_matrix_set(result, i, 2, GSL_REAL(z)); + } + + releaseTempMatrices(ws, 2); +} diff --git a/linalg.h b/linalg.h index 59931fa..c322898 100644 --- a/linalg.h +++ b/linalg.h @@ -51,7 +51,8 @@ int eigenvectors(gsl_matrix *g, gsl_matrix *evec, workspace_t *ws); int real_eigenvectors(gsl_matrix *g, gsl_matrix *evec, workspace_t *ws); void eigenvectors_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws); int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws); - +void projective_frame(gsl_vector **vertices, gsl_matrix *result, workspace_t *ws); +void rotation_frame(gsl_matrix *rotation, gsl_matrix *result, workspace_t *ws); // matrix allocation stuff diff --git a/main.c b/main.c index e26dde1..c5bbf5a 100644 --- a/main.c +++ b/main.c @@ -60,7 +60,6 @@ void destroyContext(DrawingContext *ctx) void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char *type) { gsl_matrix *tmp = getTempMatrix(ctx->ws); - gsl_matrix *rot_basis = getTempMatrix(ctx->ws); gsl_matrix **gen = getTempMatrices(ctx->ws, 3); ERROR(strlen(type) != 2, "Invalid call of computeRotationMatrix()\n"); @@ -68,24 +67,61 @@ void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char * initializeTriangleGenerators(gen, ctx->cartan); gsl_matrix_set_identity(tmp); multiply_right(tmp, gen[type[0]-'a'], ctx->ws); - multiply_right(tmp, gen[type[1]-'b'], ctx->ws); + multiply_right(tmp, gen[type[1]-'a'], ctx->ws); - gsl_eigen_nonsymmv_params(0, ctx->ws->work_nonsymmv); - int r = gsl_eigen_nonsymmv(tmp, ctx->ws->eval_complex, ctx->ws->evec_complex, ctx->ws->work_nonsymmv); - ERROR(r, "gsl_eigen_nonsymmv failed!\n"); + rotation_frame(tmp, result, ctx->ws); + releaseTempMatrices(ctx->ws, 4); +} + +void computeBoxTransform(DrawingContext *ctx, char *word1, char *word2, gsl_matrix *result) +{ + vector_t p[2][3],i[2]; + vector_t std[4] = { + {-1, -1, 1}, + {-1, 1, 1}, + {1, 1, 1}, + {1, -1, 1} + }; + + gsl_vector **vertices = getTempVectors(ctx->ws, 4); + gsl_vector **std_vertices = getTempVectors(ctx->ws, 4); + gsl_matrix *tmp = getTempMatrix(ctx->ws); + gsl_matrix *to_frame = getTempMatrix(ctx->ws); + gsl_matrix *to_std_frame = getTempMatrix(ctx->ws); + + fixedPoints(ctx, word1, p[0]); + fixedPoints(ctx, word2, p[1]); + + // intersect attracting line with neutral line of the other element + for(int j = 0; j < 2; j++) + i[j] = cross(cross(p[j%2][0],p[j%2][1]),cross(p[(j+1)%2][0],p[(j+1)%2][2])); + + // box consists of p[0][0], i[0], p[1][0], i[1] + + for(int i = 0; i < 4; i++) + vectorToGsl(std[i], std_vertices[i]); + + vectorToGsl(p[0][0], vertices[0]); + vectorToGsl(i[0], vertices[1]); + vectorToGsl(p[1][0], vertices[2]); + vectorToGsl(i[1], vertices[3]); + + projective_frame(std_vertices, to_std_frame, ctx->ws); + projective_frame(vertices, to_frame, ctx->ws); + invert(to_frame, tmp, ctx->ws); + multiply(to_std_frame, tmp, result); + + /* LOOP(i) { - gsl_complex x = gsl_matrix_complex_get(ctx->ws->evec_complex, i, 0); - gsl_complex y = gsl_matrix_complex_get(ctx->ws->evec_complex, i, 1); - gsl_complex z = gsl_matrix_complex_get(ctx->ws->evec_complex, i, 2); - gsl_matrix_set(tmp, i, 0, GSL_REAL(x)+GSL_REAL(y)); - gsl_matrix_set(tmp, i, 1, GSL_IMAG(x)-GSL_IMAG(y)); - gsl_matrix_set(tmp, i, 2, GSL_REAL(z)); - } + LOOP(j) { + printf("%.4f ", gsl_matrix_get(result, i, j)); + } + printf("\n"); + }*/ - invert(tmp, result, ctx->ws); - - releaseTempMatrices(ctx->ws, 5); + releaseTempVectors(ctx->ws, 8); + releaseTempMatrices(ctx->ws, 3); } void updateMatrices(DrawingContext *ctx) @@ -94,10 +130,13 @@ void updateMatrices(DrawingContext *ctx) LOOP(i) angle[i] = M_PI*ctx->k[i]/ctx->p[i]; cartanMatrix(ctx->cartan, angle[0], angle[1], angle[2], ctx->parameter); - if(ctx->use_rotation_basis) - computeRotationMatrix(ctx, ctx->cob, "ac"); - else + if(ctx->use_rotation_basis) { +// computeRotationMatrix(ctx, ctx->cob, "ac"); // and don't forget to invert! + computeBoxTransform(ctx, "abc", "cab", ctx->cob); + } else { gsl_matrix_memcpy(ctx->cob, ctx->cartan); // is this a good choice of basis for any reason? + } + } void print(DrawingContext *screen) @@ -241,12 +280,22 @@ int main() if(!info) 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; + updateDimensions(info->dim); screen_context->dim = info->dim; diff --git a/main.h b/main.h index ff86f35..bd9143d 100644 --- a/main.h +++ b/main.h @@ -1,6 +1,8 @@ #ifndef TRIANGLE_GROUP_MAIN_H #define TRIANGLE_GROUP_MAIN_H +#include + #include "triangle.h" #include "linalg.h" #include "initcairo.h" @@ -84,4 +86,16 @@ int processEvent(GraphicsInfo *info, XEvent *ev); void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char *type); void updateMatrices(DrawingContext *ctx); +static vector_t vectorFromGsl(gsl_vector *v) +{ + vector_t result; + LOOP(i) result.x[i] = gsl_vector_get(v, i); + return result; +} + +static void vectorToGsl(vector_t v, gsl_vector *out) +{ + LOOP(i) gsl_vector_set(out, i, v.x[i]); +} + #endif