draw arcs

This commit is contained in:
Florian Stecker 2019-12-23 12:29:50 +01:00
parent c43e9d89e0
commit 0cde265d1e
6 changed files with 592 additions and 152 deletions

619
draw.c
View File

@ -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); // <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]);
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));
}

View File

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

View File

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

82
main.c
View File

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

10
main.h
View File

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

View File

@ -4,7 +4,7 @@
#include <stdio.h>
#include <stdlib.h>
#define QUEUE_SIZE 100000
#define QUEUE_SIZE 1000000
#define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);}
#ifdef _DEBUG