diff --git a/draw.c b/draw.c index bf784d0..d0ad184 100644 --- a/draw.c +++ b/draw.c @@ -4,7 +4,7 @@ int isInsideBB(DrawingContext *ctx, point_t p) { - cairo_user_to_device(ctx->cairo, &x, &y); + cairo_user_to_device(ctx->cairo, &p.x, &p.y); return -p.x < ctx->dim->width && p.x < 3*ctx->dim->width && -p.y < ctx->dim->height && p.y < 3*ctx->dim->height; } @@ -17,17 +17,14 @@ vector_t cross(vector_t a, vector_t b) return result; } -vector_t apply(DrawingContext *ctx, gsl_matrix *m, vector_t x) +vector_t apply(gsl_matrix *m, vector_t x) { - gsl_vector *tmp = getTempVector(ctx->ws); - gsl_vector *tmp2 = getTempVector(ctx->ws); vector_t out; - LOOP(i) gsl_vector_set(tmp, i, x.x[i]); - gsl_blas_dgemv(CblasNoTrans, 1.0, m, tmp, 0.0, tmp2); - LOOP(i) out.x[i] = gsl_vector_get(tmp2, i); + LOOP(i) out.x[i] = 0.0; + LOOP(i) LOOP(j) out.x[i] += gsl_matrix_get(m, i, j) * x.x[j]; - releaseTempVectors(ctx->ws, 2); + return out; } int fixedPoints(DrawingContext *ctx, const char *word, vector_t *out) @@ -222,6 +219,49 @@ void drawBoxStd(DrawingContext *ctx, const char *word, char base) drawBox(ctx, word1, word2); } +void drawRotationOrbit(DrawingContext *ctx, const char *word, vector_t start) +{ + vector_t v[3], w; + point_t p; + double parameter; + 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); + 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)); + + 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); +} + // level 4: draw the actual image components void drawReflectors(DrawingContext *ctx) @@ -250,12 +290,12 @@ void drawAttractors(DrawingContext *ctx) LOOP(i) LOOP(j) l[i][j] = cross(p[i][(3-j)%3], p[i][(4-j)%3]); - LOOP(i) LOOP(j) { + LOOP(i) LOOP(j) if(i != 1) { cairo_set_source_rgb(ctx->cairo, color[i][0], color[i][1], color[i][2]); drawVector(ctx, p[i][j]); } - LOOP(i) LOOP(j) { + LOOP(i) LOOP(j) if(i != 1) { cairo_set_source_rgb(ctx->cairo, color[i][0], color[i][1], color[i][2]); drawCovector(ctx, l[i][j]); } @@ -328,54 +368,43 @@ void drawBoxes(DrawingContext *ctx) void drawBoxes2(DrawingContext *ctx) { gsl_matrix *tmp = getTempMatrix(ctx->ws); + gsl_matrix **gen = getTempMatrices(ctx->ws, 3); + cairo_t *C = ctx->cairo; + initializeTriangleGenerators(gen, ctx->cartan); - cairo_set_source_rgb(ctx->cairo, 1, 0.5, 0); - 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(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"); - cairo_set_source_rgb(ctx->cairo, 1, 0, 0.5); - drawBoxLines(ctx, "cab", "bca"); + vector_t fp[3], fp2[3]; + vector_t w; - cairo_set_source_rgb(ctx->cairo, 0, 0, 0); - drawBox(ctx, "abc", "cab"); +// 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]); - vector_t w = {0, 1, 0}; - drawCovector(ctx, w); + cairo_set_source_rgb(C, 1, 0, 0.5); +// drawBoxLines(ctx, "cab", "bca"); - 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); + cairo_set_source_rgb(C, 0, 0, 0); +// drawBox(ctx, "abc", "cab"); - 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); +// 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]); - 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); + releaseTempMatrices(ctx->ws, 4); } void drawLimitCurve(DrawingContext *ctx) @@ -386,19 +415,18 @@ void drawLimitCurve(DrawingContext *ctx) int previous_inside = 0; for(int i = 0; i < ctx->n_group_elements; i++) { - double x = ctx->limit_curve[3*i]; - double y = ctx->limit_curve[3*i+1]; + point_t p; + p.x = ctx->limit_curve[3*i]; + p.y = ctx->limit_curve[3*i+1]; - cairo_user_to_device(C, &x, &y); - - if(-x < ctx->dim->width && x < 3*ctx->dim->width && -y < ctx->dim->height && y < 3*ctx->dim->height) { + if(isInsideBB(ctx, p)) { if(ctx->limit_with_lines) { if(!previous_inside) - cairo_move_to(C, ctx->limit_curve[3*i], ctx->limit_curve[3*i+1]); + cairo_move_to(C, p.x, p.y); else - cairo_line_to(C, ctx->limit_curve[3*i], ctx->limit_curve[3*i+1]); + cairo_line_to(C, p.x, p.y); } else { - cairo_move_to(C, ctx->limit_curve[3*i], ctx->limit_curve[3*i+1]); + cairo_move_to(C, p.x, p.y); cairo_close_path(C); } previous_inside = 1; diff --git a/linalg.c b/linalg.c index 1f0c577..6256dca 100644 --- a/linalg.c +++ b/linalg.c @@ -72,6 +72,18 @@ void invert(gsl_matrix *in, gsl_matrix *out, workspace_t *ws) releaseTempMatrices(ws, 1); } +void solve(gsl_matrix *A, gsl_vector *b, gsl_vector *result, workspace_t *ws) +{ + int s; + gsl_matrix *tmp = getTempMatrix(ws); + + gsl_matrix_memcpy(tmp, A); + gsl_linalg_LU_decomp(tmp, ws->permutation, &s); + gsl_linalg_LU_solve(tmp, ws->permutation, b, result); + + releaseTempMatrices(ws, 1); +} + void conjugate(gsl_matrix *in, gsl_matrix *conjugator, gsl_matrix *out, workspace_t *ws) { gsl_matrix *tmp = getTempMatrix(ws); diff --git a/linalg.h b/linalg.h index c322898..38477bb 100644 --- a/linalg.h +++ b/linalg.h @@ -35,6 +35,7 @@ typedef struct _workspace { workspace_t *workspace_alloc(int n); void workspace_free(workspace_t *workspace); +void solve(gsl_matrix *A, gsl_vector *b, gsl_vector *result, workspace_t *ws); void invert(gsl_matrix *in, gsl_matrix *out, workspace_t *ws); void conjugate(gsl_matrix *in, gsl_matrix *conjugator, gsl_matrix *out, workspace_t *ws); void multiply(gsl_matrix *a, gsl_matrix *b, gsl_matrix *out); diff --git a/main.c b/main.c index c5bbf5a..6cc250d 100644 --- a/main.c +++ b/main.c @@ -62,12 +62,12 @@ void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char * gsl_matrix *tmp = getTempMatrix(ctx->ws); gsl_matrix **gen = getTempMatrices(ctx->ws, 3); - ERROR(strlen(type) != 2, "Invalid call of computeRotationMatrix()\n"); +// ERROR(strlen(type) != 2, "Invalid call of computeRotationMatrix()\n"); initializeTriangleGenerators(gen, ctx->cartan); gsl_matrix_set_identity(tmp); - multiply_right(tmp, gen[type[0]-'a'], ctx->ws); - multiply_right(tmp, gen[type[1]-'a'], ctx->ws); + for(int i = 0; i < strlen(type); i++) + multiply_right(tmp, gen[type[i]-'a'], ctx->ws); rotation_frame(tmp, result, ctx->ws); @@ -130,13 +130,18 @@ 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"); // and don't forget to invert! - computeBoxTransform(ctx, "abc", "cab", ctx->cob); - } else { + gsl_matrix *tmp = getTempMatrix(ctx->ws); + + if(ctx->use_rotation_basis % 3 == 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"); + invert(tmp, ctx->cob, ctx->ws); + } else { + computeBoxTransform(ctx, "abc", "cab", ctx->cob); } + releaseTempMatrices(ctx->ws, 1); } void print(DrawingContext *screen) @@ -249,7 +254,7 @@ int processEvent(GraphicsInfo *info, XEvent *ev) TOGGLE(screen_context->show_limit); break; case 'R': - TOGGLE(screen_context->use_rotation_basis); + screen_context->use_rotation_basis++; updateMatrices(screen_context); computeLimitCurve(screen_context); break;