florian
/
triangle_group_limit_set
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

draw rotation orbit

This commit is contained in:
Florian Stecker 2019-04-12 11:51:57 +02:00
parent 4daaf1ed7c
commit c43e9d89e0
4 changed files with 113 additions and 67 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

146
draw.c
View File

@ -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); // <acaba+, cab+>
// 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;

12
linalg.c
View File

@ -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);

1
linalg.h
View File

@ -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);

21
main.c
View File

@ -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;