draw x^alpha curves; no want to refactor limit curve

Makefile

@@ -1,4 +1,4 @@
-HEADERS=triangle.h linalg.h queue.h initcairo.h main.h
+HEADERS=triangle.h linalg.h queue.h initcairo.h main.h exp_equation.h
 
 SPECIAL_OPTIONS=-O0 -g -D_DEBUG
 #SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline
@@ -11,8 +11,8 @@ OPTIONS=$(GENERAL_OPTIONS) $(CAIRO_OPTIONS) $(SPECIAL_OPTIONS)
 
 all: limit_set
 
-limit_set: limit_set.o linalg.o triangle.o initcairo.o draw.o main.o
-	gcc $(OPTIONS) -o limit_set limit_set.o linalg.o triangle.o initcairo.o draw.o main.o -lm -lgsl -lcblas -lcairo -lX11
+limit_set: limit_set.o linalg.o triangle.o initcairo.o draw.o main.o exp_equation.o
+	gcc $(OPTIONS) -o limit_set limit_set.o linalg.o triangle.o initcairo.o draw.o main.o exp_equation.o -lm -lgsl -lcblas -lcairo -lX11
 
 linalg.o: linalg.c $(HEADERS)
 	gcc $(OPTIONS) -c linalg.c
@@ -32,5 +32,8 @@ draw.o: draw.c $(HEADERS)
 main.o: main.c $(HEADERS)
 	gcc $(OPTIONS) -c main.c
 
+exp_equation.o: exp_equation.c $(HEADERS)
+	gcc $(OPTIONS) -c exp_equation.c
+
 clean:
-	rm -f limit_set linalg.o triangle.o limit_set.o draw.o main.o
+	rm -f limit_set linalg.o triangle.o limit_set.o draw.o main.o exp_equation.o

exp_equation.c

@@ -0,0 +1,183 @@
+#include "main.h"
+#include "exp_equation.h"
+
+#include <stdio.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_roots.h>
+
+#define EPSILON 1e-9
+#define LOOP(i) for(int i = 0; i < 3; i++)
+
+struct solve_exp_plus_exp_params
+{
+	double alpha;
+	double beta;
+	double x;
+};
+
+static double solve_exp_plus_exp_f(double t, void *_params)
+{
+	struct solve_exp_plus_exp_params *params = (struct solve_exp_plus_exp_params*)_params;
+
+//	return exp(params->alpha * t) + exp(params->beta * t) - params->x;
+	return log(exp(params->alpha * t) + exp(params->beta * t)) - log(params->x);
+}
+
+static double solve_exp_plus_exp_df(double t, void *_params)
+{
+	struct solve_exp_plus_exp_params *params = (struct solve_exp_plus_exp_params*)_params;
+
+//	return params->alpha * exp(params->alpha * t) + params->beta * exp(params->beta * t);
+	return params->alpha + (params->beta - params->alpha) / (1 + exp((params->alpha-params->beta)*t));
+}
+
+static void solve_exp_plus_exp_fdf(double t, void *params, double *f, double *df)
+{
+	*f = solve_exp_plus_exp_f(t, params);
+	*df = solve_exp_plus_exp_df(t, params);
+}
+
+// solve the equation exp(alpha t) + exp(beta t) = x for t
+int solve_exp_plus_exp(double alpha, double beta, double x, double *t)
+{
+	if(alpha <= 0 && beta >= 0 || alpha >= 0 && beta <= 0) {
+		double critical =
+			pow(-beta/alpha, alpha/(alpha-beta)) +
+			pow(-beta/alpha, beta/(alpha-beta));
+
+		if(x < critical)
+			return 0;
+		else if (x < critical + EPSILON) {
+			t[0] = log(-beta/alpha)/(alpha-beta);
+			return 1;
+		}
+
+		// Newton this
+		gsl_root_fdfsolver *solver = gsl_root_fdfsolver_alloc(gsl_root_fdfsolver_newton);
+		struct solve_exp_plus_exp_params params;
+		params.alpha = alpha;
+		params.beta = beta;
+		params.x = x;
+		gsl_function_fdf FDF;
+		FDF.f = &solve_exp_plus_exp_f;
+		FDF.df = &solve_exp_plus_exp_df;
+		FDF.fdf = &solve_exp_plus_exp_fdf;
+		FDF.params = (void *)&params;
+		int status;
+		double root, lastroot;
+
+		for(int r = 0; r < 2; r++) {
+			root = r == 0 ? log(x)/beta : log(x)/alpha;
+			gsl_root_fdfsolver_set(solver, &FDF, root);
+			for(int i = 0; i < 100; i++) {
+				gsl_root_fdfsolver_iterate(solver);
+				lastroot = root;
+				root = gsl_root_fdfsolver_root(solver);
+				status = gsl_root_test_delta(root, lastroot, 0, 1e-9);
+
+				if(status == GSL_SUCCESS)
+					break;
+
+//				printf("iteration %d, root %f\n", i, root);
+			}
+			t[r] = root;
+		}
+
+		gsl_root_fdfsolver_free(solver);
+
+		return 2;
+	} else {
+		// Newton with start value 0
+		gsl_root_fdfsolver *solver = gsl_root_fdfsolver_alloc(gsl_root_fdfsolver_newton);
+		struct solve_exp_plus_exp_params params;
+		params.alpha = alpha;
+		params.beta = beta;
+		params.x = x;
+		gsl_function_fdf FDF;
+		FDF.f = &solve_exp_plus_exp_f;
+		FDF.df = &solve_exp_plus_exp_df;
+		FDF.fdf = &solve_exp_plus_exp_fdf;
+		FDF.params = (void *)&params;
+		int status;
+		double root, lastroot;
+
+		root = 0;
+		gsl_root_fdfsolver_set(solver, &FDF, root);
+		for(int i = 0; i < 100; i++) {
+			gsl_root_fdfsolver_iterate(solver);
+			lastroot = root;
+			root = gsl_root_fdfsolver_root(solver);
+			status = gsl_root_test_delta(root, lastroot, 0, 1e-9);
+
+//			printf("iteration %d, root %f\n", i, root);
+
+			if(status == GSL_SUCCESS)
+				break;
+		}
+		t[0] = root;
+
+		gsl_root_fdfsolver_free(solver);
+
+		return 1;
+	}
+}
+
+// solve the equation x1 exp(a1 t) + x2 exp(a2 t) + x3 exp(a3 t) = 0
+int solve_linear_exp(vector_t a, vector_t x, double *t)
+{
+	if(x.x[0] > 0 && x.x[1] > 0 && x.x[2] > 0 ||
+	   x.x[0] < 0 && x.x[1] < 0 && x.x[2] < 0)
+		return 0;
+
+	// ensure that y[0] < 0 and y[1], y[2] > 0
+	int j;
+	vector_t y, b;
+	for(j = 0; j < 3; j++)
+		if(x.x[(j+1)%3] * x.x[(j+2)%3] > 0)
+			break;
+	LOOP(i) y.x[i] = x.x[(i+j)%3];
+	if(y.x[0] > 0)
+		LOOP(i) y.x[i] *= -1;
+	LOOP(i) b.x[i] = a.x[(i+j)%3];
+
+	double T = (log(y.x[1]) - log(y.x[2])) / (b.x[2] - b.x[1]);
+	double rhs = - y.x[0] *
+		pow(y.x[1], (b.x[0]-b.x[2])/(b.x[2]-b.x[1])) *
+		pow(y.x[2], (b.x[0]-b.x[1])/(b.x[1]-b.x[2]));
+
+	int n = solve_exp_plus_exp(b.x[1] - b.x[0], b.x[2] - b.x[0], rhs, t);
+
+	for(int i = 0; i < n; i++)
+		t[i] += T;
+
+	return n;
+}
+
+/*
+int main(int argc, char *argv[])
+{
+
+//	int n = solve_exp_plus_exp(atof(argv[1]), atof(argv[2]), atof(argv[3]), result);
+
+	double result[2];
+	vector_t a, x;
+	a.x[0] = atof(argv[1]);
+	a.x[1] = atof(argv[2]);
+	a.x[2] = atof(argv[3]);
+	x.x[0] = atof(argv[4]);
+	x.x[1] = atof(argv[5]);
+	x.x[2] = atof(argv[6]);
+	int n = solve_linear_exp(a, x, result);
+
+	if(n == 0)
+		printf("0 results found\n");
+	else if(n == 1)
+		printf("1 result found: %.9f\n", result[0]);
+	else if(n == 2)
+		printf("2 results found: %.9f and %.9f\n", result[0], result[1]);
+	else
+		printf("%d results found\n", n);
+	return 0;
+}
+*/

exp_equation.h

@@ -0,0 +1,14 @@
+#ifndef EXP_EQUATION_H
+#define EXP_EQUATION_H
+
+#include "main.h"
+
+#include <stdio.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_roots.h>
+
+int solve_exp_plus_exp(double alpha, double beta, double x, double *t);
+int solve_linear_exp(vector_t a, vector_t x, double *t);
+
+#endif

limit_set.c

@@ -5,7 +5,6 @@ static int compareAngle(const void *x, const void *y)
 	return ((double*)x)[2] > ((double*)y)[2] ? 1 : -1;
 }
 
-// might need a rewrite
 void cartanMatrix(gsl_matrix *cartan, double a1, double a2, double a3, double s)
 {
 	gsl_matrix_set(cartan, 0, 0, 2);
@@ -21,58 +20,11 @@ void cartanMatrix(gsl_matrix *cartan, double a1, double a2, double a3, double s)
 	gsl_matrix_set(cartan, 2, 2, 2);
 }
 
-void initializeTriangleGenerators(gsl_matrix **gen, double a1, double a2, double a3, double s, double t, workspace_t *ws)
+void initializeTriangleGenerators(gsl_matrix **gen, gsl_matrix *cartan)
 {
-	gsl_matrix *reflection_gen[3];
-
-	LOOP(i) {
-		reflection_gen[i] = gsl_matrix_alloc(3, 3);
-		gsl_matrix_set_identity(reflection_gen[i]);
-	}
-
-	double rho[3];
-	rho[0] = sqrt(s*s + 2*s*cos(a1) + 1);
-	rho[1] = sqrt(s*s + 2*s*cos(a2) + 1);
-	rho[2] = sqrt(s*s + 2*s*cos(a3) + 1);
-
-	gsl_matrix_set(reflection_gen[0], 0, 0, -1.0);
-	gsl_matrix_set(reflection_gen[0], 0, 1, rho[2]*t);
-	gsl_matrix_set(reflection_gen[0], 0, 2, rho[1]/t);
-
-	gsl_matrix_set(reflection_gen[1], 1, 0, rho[2]/t);
-	gsl_matrix_set(reflection_gen[1], 1, 1, -1.0);
-	gsl_matrix_set(reflection_gen[1], 1, 2, rho[0]*t);
-
-	gsl_matrix_set(reflection_gen[2], 2, 0, rho[1]*t);
-	gsl_matrix_set(reflection_gen[2], 2, 1, rho[0]/t);
-	gsl_matrix_set(reflection_gen[2], 2, 2, -1.0);
-
-	LOOP(i) {
-		gsl_matrix_set_identity(gen[i]);
-		gsl_matrix_set(gen[i], (i+1)%3, (i+1)%3, s);
-		gsl_matrix_set(gen[i], (i+2)%3, (i+2)%3, 1/s);
-
-		gsl_matrix_set_identity(gen[i+3]);
-		gsl_matrix_set(gen[i+3], (i+1)%3, (i+1)%3, 1/s);
-		gsl_matrix_set(gen[i+3], (i+2)%3, (i+2)%3, s);
-	}
-
-	LOOP(i) {
-		multiply_left(reflection_gen[i], gen[(i+2)%3], ws);
-		multiply_right(gen[(i+2)%3], reflection_gen[(i+1)%3], ws);
-
-		multiply_left(reflection_gen[(i+1)%3], gen[(i+2)%3+3], ws);
-		multiply_right(gen[(i+2)%3+3], reflection_gen[i], ws);
-	}
-
-	LOOP(i) gsl_matrix_free(reflection_gen[i]);
-}
-
-void initializeTriangleGeneratorsCurrent(gsl_matrix **gen, DrawingContext *ctx)
-{
-	double angle[3];
-	LOOP(i) angle[i] = 2*M_PI*ctx->k[i]/ctx->p[i];
-	initializeTriangleGenerators(gen, angle[0], angle[1], angle[2], ctx->parameter2, ctx->parameter, ctx->ws);
+	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);
 }
 
 int computeLimitCurve(DrawingContext *ctx)
@@ -86,7 +38,7 @@ int computeLimitCurve(DrawingContext *ctx)
 	gsl_matrix *coxeter = getTempMatrix(ctx->ws);
 	gsl_matrix *coxeter_fixedpoints = getTempMatrix(ctx->ws);
 	gsl_matrix *fixedpoints = getTempMatrix(ctx->ws);
-	gsl_matrix **gen = getTempMatrices(ctx->ws, 6);
+	gsl_matrix **gen = getTempMatrices(ctx->ws, 3);
 	gsl_matrix **elements = getTempMatrices(ctx->ws, ctx->n_group_elements);
 	groupelement_t *group = ctx->group;
 	int success = 0;
@@ -98,43 +50,31 @@ int computeLimitCurve(DrawingContext *ctx)
 
 	// 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);
-	initializeTriangleGenerators(gen, 2*M_PI/ctx->p[0], 2*M_PI/ctx->p[1], 2*M_PI/ctx->p[2], 1.0, 1.0, ctx->ws);
+	initializeTriangleGenerators(gen, cartan_pos);
 	gsl_matrix_set_identity(elements[0]);
-	for(int i = 1; i < ctx->n_group_elements; i++) {
-		if(group[i].length % 2)
-			continue;
-		int letter = ROTATION_LETTER(group[i].letter, group[i].parent->letter);
-		multiply(gen[letter], elements[group[i].parent->parent->id], elements[i]);
-	}
+	for(int i = 1; i < ctx->n_group_elements; i++)
+		multiply(gen[group[i].letter], elements[group[i].parent->id], elements[i]);
 	diagonalize_symmetric_form(cartan_pos, cob_pos, ws);
-	multiply_many(ws, coxeter_pos, 3, gen[2], gen[1], gen[0]);
+	multiply_many(ws, coxeter_pos, 3, gen[0], gen[1], gen[2]);
 	int ev_count_pos = real_eigenvectors(coxeter_pos, coxeter_fixedpoints_pos, ws);
 
 	if(ev_count_pos != 3)
 		goto error_out;
 
-	int n = 0;
 	for(int i = 0; i < ctx->n_group_elements; i++) {
-		if(group[i].length % 2)
-			continue;
 		multiply_many(ws, fixedpoints_pos, 3, cob_pos, elements[i], coxeter_fixedpoints_pos);
-		ctx->limit_curve[3*n+2] = atan2(
+		ctx->limit_curve[12*i+2] = atan2(
 			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));
-		n++;
 	}
 
 	// now do it again to calculate x and y coordinates
-	initializeTriangleGeneratorsCurrent(gen, ctx);
+	initializeTriangleGenerators(gen, ctx->cartan);
 	gsl_matrix_set_identity(elements[0]);
-	for(int i = 1; i < ctx->n_group_elements; i++) {
-		if(group[i].length % 2)
-			continue;
-		int letter = ROTATION_LETTER(group[i].letter, group[i].parent->letter);
-		multiply(gen[letter], elements[group[i].parent->parent->id], elements[i]);
-	}
+	for(int i = 1; i < ctx->n_group_elements; i++)
+		multiply(gen[group[i].letter], elements[group[i].parent->id], elements[i]);
 
-	multiply_many(ws, coxeter, 3, gen[2], gen[1], gen[0]);
+	multiply_many(ws, coxeter, 3, gen[0], gen[1], gen[2]);
 	int ev_count = real_eigenvectors(coxeter, coxeter_fixedpoints, ws);
 
 	if(ev_count == 1)
@@ -142,17 +82,11 @@ int computeLimitCurve(DrawingContext *ctx)
 	if(ev_count == 0)
 		goto error_out;
 
-	ctx->limit_curve_count = 0;
 	for(int i = 0; i < ctx->n_group_elements; i++) {
-		if(group[i].length % 2)
-			continue;
-
 		multiply_many(ws, fixedpoints, 3, ctx->cob, elements[i], coxeter_fixedpoints);
 
-		x = ctx->limit_curve[3*ctx->limit_curve_count  ] = gsl_matrix_get(fixedpoints, 0, column)/gsl_matrix_get(fixedpoints, 2, column);
-		y = ctx->limit_curve[3*ctx->limit_curve_count+1] = gsl_matrix_get(fixedpoints, 1, column)/gsl_matrix_get(fixedpoints, 2, column);
-
-		ctx->limit_curve_count++;
+		x = ctx->limit_curve[12*i  ] = gsl_matrix_get(fixedpoints, 0, column)/gsl_matrix_get(fixedpoints, 2, column);
+		y = ctx->limit_curve[12*i+1] = gsl_matrix_get(fixedpoints, 1, column)/gsl_matrix_get(fixedpoints, 2, column);
 
 		if((x - ctx->marking.x)*(x - ctx->marking.x) + (y - ctx->marking.y)*(y - ctx->marking.y) < 25e-10)
 		{
@@ -161,16 +95,22 @@ int computeLimitCurve(DrawingContext *ctx)
 				fputc('a' + cur->letter, stdout); // bcbcbca, bacbcacab, bc bca cb
 			fputc('\n',stdout);
 		}
+
+		multiply_many(ws, fixedpoints, 2, elements[i], coxeter_fixedpoints);
+		LOOP(j) LOOP(k) ctx->limit_curve[12*i+3+3*j+k] = gsl_matrix_get(fixedpoints, k, j);
+
+		// bca abc acb = abc
+
 	}
 
-	qsort(ctx->limit_curve, ctx->limit_curve_count, 3*sizeof(double), compareAngle);
+	qsort(ctx->limit_curve, ctx->n_group_elements, 12*sizeof(double), compareAngle);
 
-//	ctx->limit_curve_count = ctx->n_group_elements;
+	ctx->limit_curve_count = ctx->n_group_elements;
 
 	success = 1;
 
 error_out:
-	releaseTempMatrices(ctx->ws, 14+ctx->n_group_elements);
+	releaseTempMatrices(ctx->ws, 11+ctx->n_group_elements);
 
 	return success;
 }

linalg.c

@@ -315,7 +315,6 @@ int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws)
 	ERROR(r, "gsl_eigen_symmv failed!\n");
 
 	gsl_eigen_symmv_sort(ws->eval_real, cob, GSL_EIGEN_SORT_VAL_ASC);
-
 	gsl_matrix_transpose(cob);
 
 	int positive = 0;
@@ -332,6 +331,32 @@ int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws)
 	return positive;
 }
 
+int diagonalize_symmetric_matrix(gsl_matrix *A, gsl_matrix *cob, double *eval, workspace_t *ws)
+{
+	gsl_matrix *A_ = getTempMatrix(ws);
+
+	gsl_matrix_memcpy(A_, A);
+	int r = gsl_eigen_symmv (A_, ws->eval_real, cob, ws->work_symmv);
+	ERROR(r, "gsl_eigen_symmv failed!\n");
+
+	gsl_eigen_symmv_sort(ws->eval_real, cob, GSL_EIGEN_SORT_VAL_ASC);
+
+	gsl_matrix_transpose(cob);
+
+	int positive = 0;
+
+	for(int i = 0; i < ws->n; i++) {
+		if(eval)
+			eval[i] = gsl_vector_get(ws->eval_real, i);
+
+		if(gsl_vector_get(ws->eval_real, i) > 0)
+			positive++;
+	}
+
+	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)
 {

linalg.h

@@ -53,6 +53,7 @@ int real_eigenvectors(gsl_matrix *g, gsl_matrix *evec, workspace_t *ws);
 int real_eigenvalues(gsl_matrix *g, gsl_vector *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);
+int diagonalize_symmetric_matrix(gsl_matrix *A, gsl_matrix *cob, double *eval, 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);
 

main.c

@@ -11,474 +11,481 @@
 #include "linalg.h"
 
 #define TOGGLE(a) do { (a) = !(a); } while(0)
-#define SIGN(x) ((x) > 0 ? 1.0 : -1.0)
 
 DrawingContext *screen_context;
 
 // setup everything except cairo and dim, which will be provided by the graphics system
 void setupContext(DrawingContext *ctx, int argc, char *argv[])
 {
-    ctx->n_group_elements = NUM_GROUP_ELEMENTS;
-    ctx->n_group_elements_combinatorial = NUM_GROUP_ELEMENTS_COMBINATORIAL;
-    ctx->p[0] = atoi(argv[1]);
-    ctx->p[1] = atoi(argv[2]);
-    ctx->p[2] = atoi(argv[3]);
-    ctx->k[0] = atoi(argv[4]);
-    ctx->k[1] = atoi(argv[5]);
-    ctx->k[2] = atoi(argv[6]);
-    if(argc > 7)
-	ctx->parameter = atof(argv[7]);
-    else
-	ctx->parameter = 1.0;
-    if(argc > 8)
-	ctx->parameter2 = atof(argv[8]);
-    else
-	ctx->parameter2 = 1.0;
-    if(argc > 12) {
-	ctx->movie_filename = argv[9];
-	ctx->movie_parameter_duration = atof(argv[10]);
-	ctx->movie_parameter2_duration = atof(argv[11]);
-	ctx->movie_n_frames = atoi(argv[12]);
-    } else {
-	ctx->movie_n_frames = 0;
-    }
-    //	ctx->parameter = 2.77;
-    //	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 = 0;
-    ctx->show_dual_limit = 0;
-    ctx->show_text = 1;
-    ctx->mode = 0;
-    ctx->use_rotation_basis = 1;
-    ctx->limit_with_lines = 0;
-    ctx->use_repelling = 0;
-    ctx->show_marking = 0;
-    ctx->marking.x = -0.73679;
-    ctx->marking.y = -0.01873;
-    ctx->show_coxeter_orbit = 0;
+	ctx->n_group_elements = NUM_GROUP_ELEMENTS;
+	ctx->p[0] = atoi(argv[1]);
+	ctx->p[1] = atoi(argv[2]);
+	ctx->p[2] = atoi(argv[3]);
+	ctx->k[0] = atoi(argv[4]);
+	ctx->k[1] = atoi(argv[5]);
+	ctx->k[2] = atoi(argv[6]);
+	if(argc > 7)
+		ctx->parameter = atof(argv[7]);
+	else
+		ctx->parameter = 1.0;
+//	ctx->parameter = 2.77;
+//	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= 0;
+	ctx->show_dual_limit= 0;
+	ctx->show_text = 1;
+	ctx->mode = 0;
+	ctx->use_rotation_basis = 0;
+	ctx->limit_with_lines = 1;
+	ctx->use_repelling = 0;
+	ctx->show_marking = 1;
+	ctx->marking.x = -0.73679;
+	ctx->marking.y = -0.01873;
+	ctx->marking2.x = -0.73679;
+	ctx->marking2.y = -0.11873;
+	ctx->marking3.x = -0.73679;
+	ctx->marking3.y = -0.21873;
+	ctx->distance_parameter1 = 0.45;
+	ctx->distance_parameter2 = 0.2;
+	ctx->show_coxeter_orbit = 0;
+	ctx->extra_text = malloc(1000*sizeof(char));
+	memset(ctx->extra_text, 0, 1000*sizeof(char));
 
-    ctx->limit_curve = malloc(3*ctx->n_group_elements*sizeof(double));
-    ctx->limit_curve_count = -1;
+	ctx->limit_curve = malloc(12*ctx->n_group_elements*sizeof(double));
+	ctx->limit_curve_count = -1;
 
-    ctx->group = malloc(ctx->n_group_elements_combinatorial*sizeof(groupelement_t));
-    generate_triangle_group(ctx->group, ctx->n_group_elements_combinatorial, ctx->p[0], ctx->p[1], ctx->p[2]);
+	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]);
 
-    // the temporary stuff
-    ctx->cartan = gsl_matrix_alloc(3, 3);
-    ctx->cob = gsl_matrix_alloc(3, 3);
-    ctx->ws = workspace_alloc(3);
+	// the temporary stuff
+	ctx->cartan = gsl_matrix_alloc(3, 3);
+	ctx->cob = gsl_matrix_alloc(3, 3);
+	ctx->ws = workspace_alloc(3);
 }
 
 void destroyContext(DrawingContext *ctx)
 {
-    free(ctx->limit_curve);
-    free(ctx->group);
+	free(ctx->limit_curve);
+	free(ctx->group);
+	free(ctx->extra_text);
 
-    gsl_matrix_free(ctx->cartan);
-    gsl_matrix_free(ctx->cob);
+	gsl_matrix_free(ctx->cartan);
+	gsl_matrix_free(ctx->cob);
 
-    workspace_free(ctx->ws);
+	workspace_free(ctx->ws);
 }
 
-void computeMatrix(DrawingContext *ctx, gsl_matrix *result, const char *type)
+void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char *word)
 {
-    gsl_matrix **gen = getTempMatrices(ctx->ws, 6);
+	gsl_matrix *tmp = getTempMatrix(ctx->ws);
+	gsl_matrix **gen = getTempMatrices(ctx->ws, 3);
 
-    //	ERROR(strlen(type) != 2, "Invalid call of computeRotationMatrix()\n");
+	initializeTriangleGenerators(gen, ctx->cartan);
+	gsl_matrix_set_identity(tmp);
+	for(int i = 0; i < strlen(word); i++)
+		multiply_right(tmp, gen[word[i]-'a'], ctx->ws);
 
-    initializeTriangleGeneratorsCurrent(gen, ctx);
-    gsl_matrix_set_identity(result);
-    for(int i = 0; i < strlen(type); i++) {
-	if(type[i] >= 'a' && type[i] <= 'c')
-	    multiply_right(result, gen[type[i]-'a'], ctx->ws);
-	else if(type[i] >= 'A' && type[i] <= 'C')
-	    multiply_right(result, gen[type[i]-'A'+3], ctx->ws);
-    }
+	rotation_frame(tmp, result, ctx->ws);
 
-    releaseTempMatrices(ctx->ws, 6);
-}
-
-void computeRotationMatrixFrame(DrawingContext *ctx, gsl_matrix *result, const char *type)
-{
-    gsl_matrix *tmp = getTempMatrix(ctx->ws);
-    computeMatrix(ctx, tmp, type);
-    rotation_frame(tmp, result, ctx->ws);
-    releaseTempMatrices(ctx->ws, 1);
+	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}
-    };
+	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);
+	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]);
+	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]));
+	// 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]
+	// 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]);
+	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]);
+	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);
+	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) {
-      LOOP(j) {
-      printf("%.4f ", gsl_matrix_get(result, i, j));
-      }
-      printf("\n");
-      }*/
+	/*
+	LOOP(i) {
+		LOOP(j) {
+			printf("%.4f ", gsl_matrix_get(result, i, j));
+		}
+		printf("\n");
+		}*/
 
-    releaseTempVectors(ctx->ws, 8);
-    releaseTempMatrices(ctx->ws, 3);
+	releaseTempVectors(ctx->ws, 8);
+	releaseTempMatrices(ctx->ws, 3);
 }
 
 void updateMatrices(DrawingContext *ctx)
 {
-    double angle[3];
-    LOOP(i) angle[i] = M_PI*ctx->k[i]/ctx->p[i];
-    cartanMatrix(ctx->cartan, angle[0], angle[1], angle[2], ctx->parameter);
+	double angle[3];
+	LOOP(i) angle[i] = M_PI*ctx->k[i]/ctx->p[i];
+	cartanMatrix(ctx->cartan, angle[0], angle[1], angle[2], ctx->parameter);
 
-    gsl_matrix *tmp = getTempMatrix(ctx->ws);
-    int nmodes = 5;
+	gsl_matrix *tmp = getTempMatrix(ctx->ws);
 
-    if(ctx->use_rotation_basis % nmodes == 0) {
-	gsl_matrix_set(tmp, 0, 0, 0.0);
-	gsl_matrix_set(tmp, 0, 1, sqrt(3.0)/2.0);
-	gsl_matrix_set(tmp, 0, 2, -sqrt(3.0)/2.0);
-	gsl_matrix_set(tmp, 1, 0, 1.0);
-	gsl_matrix_set(tmp, 1, 1, -0.5);
-	gsl_matrix_set(tmp, 1, 2, -0.5);
-	gsl_matrix_set(tmp, 2, 0, 1.0);
-	gsl_matrix_set(tmp, 2, 1, 1.0);
-	gsl_matrix_set(tmp, 2, 2, 1.0);
-	gsl_matrix_memcpy(ctx->cob, tmp);
-    } else if(ctx->use_rotation_basis % nmodes == 1) {
-	gsl_matrix_set(tmp, 0, 0, 1.0);
-	gsl_matrix_set(tmp, 0, 1, -1.0);
-	gsl_matrix_set(tmp, 0, 2, 0.0);
-	gsl_matrix_set(tmp, 1, 0, 1.0);
-	gsl_matrix_set(tmp, 1, 1, 1.0);
-	gsl_matrix_set(tmp, 1, 2, 0.0);
-	gsl_matrix_set(tmp, 2, 0, 0.0);
-	gsl_matrix_set(tmp, 2, 1, 0.0);
-	gsl_matrix_set(tmp, 2, 2, 1.0);
-	gsl_matrix_memcpy(ctx->cob, ctx->cartan); // is this a good choice of basis for any reason?
-	multiply_left(tmp, ctx->cob, ctx->ws);
-    } else if(ctx->use_rotation_basis % nmodes == 2) {
-	computeRotationMatrixFrame(ctx, tmp, "C");
-	invert(tmp, ctx->cob, ctx->ws);
-    } else if(ctx->use_rotation_basis % nmodes == 3) {
-	computeBoxTransform(ctx, "acb", "cba", ctx->cob);
-	//		computeBoxTransform(ctx, "cab", "bca", ctx->cob);
-	//		computeBoxTransform(ctx, "acb", "cba", ctx->cob);
-    } else {
-	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);
-    }
+	if(ctx->use_rotation_basis % 5 == 0) {
+		gsl_matrix_set(tmp, 0, 0, 0.0);
+		gsl_matrix_set(tmp, 0, 1, sqrt(3.0)/2.0);
+		gsl_matrix_set(tmp, 0, 2, -sqrt(3.0)/2.0);
+		gsl_matrix_set(tmp, 1, 0, 1.0);
+		gsl_matrix_set(tmp, 1, 1, -0.5);
+		gsl_matrix_set(tmp, 1, 2, -0.5);
+		gsl_matrix_set(tmp, 2, 0, 1.0);
+		gsl_matrix_set(tmp, 2, 1, 1.0);
+		gsl_matrix_set(tmp, 2, 2, 1.0);
+		gsl_matrix_memcpy(ctx->cob, tmp);
+	} else if(ctx->use_rotation_basis % 5 == 1) {
+		gsl_matrix_memcpy(ctx->cob, ctx->cartan); // is this a good choice of basis for any reason?
+	} else if(ctx->use_rotation_basis % 5 == 2) {
+		computeRotationMatrix(ctx, tmp, "ba");
+		invert(tmp, ctx->cob, ctx->ws);
+	} else if(ctx->use_rotation_basis % 5 == 3) {
+		computeBoxTransform(ctx, "bca", "abc", ctx->cob);
+//		computeBoxTransform(ctx, "cab", "bca", ctx->cob);
+//		computeBoxTransform(ctx, "acb", "cba", ctx->cob);
+	} else {
+		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);
+	releaseTempMatrices(ctx->ws, 1);
 }
 
 void output_info(DrawingContext *ctx)
 {
-    vector_t p[4][3];
-    point_t pt;
+	vector_t p[4][3];
+	point_t pt;
 
-    fixedPoints(ctx, "abc", p[0]);
-    fixedPoints(ctx, "bca", p[1]);
-    fixedPoints(ctx, "cab", p[2]);
+	fixedPoints(ctx, "abc", p[0]);
+	fixedPoints(ctx, "bca", p[1]);
+	fixedPoints(ctx, "cab", p[2]);
 
-    pt = vectorToPoint(ctx, p[0][0]);
-    printf("(abc)-+ = (%f %f)\n", pt.x, pt.y);
-    pt = vectorToPoint(ctx, p[1][0]);
-    printf("(bca)-+ = (%f %f)\n", pt.x, pt.y);
+	pt = vectorToPoint(ctx, p[0][0]);
+	printf("(abc)-+ = (%f %f)\n", pt.x, pt.y);
+	pt = vectorToPoint(ctx, p[1][0]);
+	printf("(bca)-+ = (%f %f)\n", pt.x, pt.y);
 }
 
 void print(DrawingContext *screen)
 {
-    DrawingContext file;
-    DimensionsInfo dim;
-    cairo_surface_t *surface;
+	DrawingContext file;
+	DimensionsInfo dim;
+	cairo_surface_t *surface;
 
-    char filename[100];
-    time_t t = time(NULL);
-    strftime(filename, sizeof(filename), "screenshot_%Y%m%d_%H%M%S.pdf", localtime(&t));
+	char filename[100];
+	time_t t = time(NULL);
+	strftime(filename, sizeof(filename), "screenshot_%Y%m%d_%H%M%S.pdf", localtime(&t));
 
-    memcpy(&file, screen, sizeof(file));
+	memcpy(&file, screen, sizeof(file));
 
-    dim.width = screen->dim->width;
-    dim.height = screen->dim->width / sqrt(2.0);
-    dim.matrix = screen->dim->matrix;
-    dim.matrix.y0 += ((double)dim.height - (double)screen->dim->height) / 2.0; // recenter vertically
-    updateDimensions(&dim);
-    file.dim = &dim;
+	dim.width = screen->dim->width;
+	dim.height = screen->dim->width / sqrt(2.0);
+	dim.matrix = screen->dim->matrix;
+	dim.matrix.y0 += ((double)dim.height - (double)screen->dim->height) / 2.0; // recenter vertically
+	updateDimensions(&dim);
+	file.dim = &dim;
 
-    surface = cairo_pdf_surface_create(filename, (double)dim.width, (double)dim.height);
+	surface = cairo_pdf_surface_create(filename, (double)dim.width, (double)dim.height);
     file.cairo = cairo_create(surface);
 
-    draw(&file);
+	draw(&file);
 
-    cairo_destroy(file.cairo);
-    cairo_surface_destroy(surface);
+	cairo_destroy(file.cairo);
+	cairo_surface_destroy(surface);
 
-    printf("Wrote sceenshot to file: %s\n", filename);
+	printf("Wrote sceenshot to file: %s\n", filename);
 }
 
 int processEvent(GraphicsInfo *info, XEvent *ev)
 {
-    int state;
-    unsigned long key;
-    char filename[100];
+	int state;
+	unsigned long key;
+	char filename[100];
 
-    //	fprintf(stderr, "Event: %d\n", ev->type);
+//	fprintf(stderr, "Event: %d\n", ev->type);
 
-    switch(ev->type) {
-    case ButtonPress:
-	state = ev->xbutton.state & (ShiftMask | LockMask | ControlMask);
+	switch(ev->type) {
+	case ButtonPress:
+		state = ev->xbutton.state & (ShiftMask | LockMask | ControlMask);
 
-	if(ev->xbutton.button == 1 && state & ShiftMask) {
-	    screen_context->marking.x = (double)ev->xbutton.x;
-	    screen_context->marking.y = (double)ev->xbutton.y;
-	    printf("mouse button pressed: %f, %f\n", screen_context->marking.x, screen_context->marking.y);
-	    cairo_set_matrix(screen_context->cairo, &screen_context->dim->matrix);
-	    cairo_device_to_user(screen_context->cairo, &screen_context->marking.x, &screen_context->marking.y);
-	    printf("mouse button pressed transformed: %f, %f\n", screen_context->marking.x, screen_context->marking.y);
-	    return STATUS_REDRAW;
-	}
-	break;
+		if(ev->xbutton.button == 1 && state & ShiftMask && state & ControlMask) {
+			screen_context->marking.x = (double)ev->xbutton.x;
+			screen_context->marking.y = (double)ev->xbutton.y;
+			printf("mouse button pressed: %f, %f\n", screen_context->marking.x, screen_context->marking.y);
+			cairo_set_matrix(screen_context->cairo, &screen_context->dim->matrix);
+			cairo_device_to_user(screen_context->cairo, &screen_context->marking.x, &screen_context->marking.y);
+			printf("mouse button pressed transformed: %f, %f\n", screen_context->marking.x, screen_context->marking.y);
+			return STATUS_REDRAW;
+		} else if(ev->xbutton.button == 1 && state & ControlMask) {
+			screen_context->marking2.x = (double)ev->xbutton.x;
+			screen_context->marking2.y = (double)ev->xbutton.y;
+			printf("mouse button pressed: %f, %f\n", screen_context->marking2.x, screen_context->marking2.y);
+			cairo_set_matrix(screen_context->cairo, &screen_context->dim->matrix);
+			cairo_device_to_user(screen_context->cairo, &screen_context->marking2.x, &screen_context->marking2.y);
+			printf("mouse button pressed transformed: %f, %f\n", screen_context->marking2.x, screen_context->marking2.y);
+			return STATUS_REDRAW;
+		} else if(ev->xbutton.button == 1 && state & ShiftMask) {
+			screen_context->marking3.x = (double)ev->xbutton.x;
+			screen_context->marking3.y = (double)ev->xbutton.y;
+			printf("mouse button pressed: %f, %f\n", screen_context->marking3.x, screen_context->marking3.y);
+			cairo_set_matrix(screen_context->cairo, &screen_context->dim->matrix);
+			cairo_device_to_user(screen_context->cairo, &screen_context->marking3.x, &screen_context->marking3.y);
+			printf("mouse button pressed transformed: %f, %f\n", screen_context->marking3.x, screen_context->marking3.y);
+			return STATUS_REDRAW;
+		}
 
-    case KeyPress:
-	state = ev->xkey.state & (ShiftMask | LockMask | ControlMask);
-	key = XkbKeycodeToKeysym(ev->xkey.display, ev->xkey.keycode, 0, !!(state & ShiftMask));
-	printf("Key pressed: %ld\n", key);
+		break;
 
-	switch(key) {
-	case XK_Down:
-	    if(ev->xkey.state & ShiftMask)
-		screen_context->parameter /= exp(0.00005);
-	    else
-		screen_context->parameter /= exp(0.002);
-	    updateMatrices(screen_context);
-	    computeLimitCurve(screen_context);
-	    break;
-	case XK_Up:
-	    if(ev->xkey.state & ShiftMask)
-		screen_context->parameter *= exp(0.00005);
-	    else
-		screen_context->parameter *= exp(0.002);
-	    updateMatrices(screen_context);
-	    computeLimitCurve(screen_context);
-	    break;
-	case XK_Left:
-	    if(ev->xkey.state & ShiftMask)
-		screen_context->parameter2 /= exp(0.00005);
-	    else
-		screen_context->parameter2 /= exp(0.002);
-	    updateMatrices(screen_context);
-	    computeLimitCurve(screen_context);
-	    break;
-	case XK_Right:
-	    if(ev->xkey.state & ShiftMask)
-		screen_context->parameter2 *= exp(0.00005);
-	    else
-		screen_context->parameter2 *= exp(0.002);
-	    updateMatrices(screen_context);
-	    computeLimitCurve(screen_context);
-	    break;
-	case XK_Page_Down:
-	    screen_context->parameter /= exp(0.02);
-	    updateMatrices(screen_context);
-	    computeLimitCurve(screen_context);
-	    break;
-	case XK_Page_Up:
-	    screen_context->parameter *= exp(0.02);
-	    updateMatrices(screen_context);
-	    computeLimitCurve(screen_context);
-	    break;
-	case ' ':
-	    screen_context->parameter = 5.57959706;
-	    updateMatrices(screen_context);
-	    computeLimitCurve(screen_context);
-	    break;
-	case XK_Return:
-	    //			screen_context->parameter = 2.76375163;
-	    screen_context->parameter = 5.29063366;
-	    updateMatrices(screen_context);
-	    computeLimitCurve(screen_context);
-	    break;
-	case 'm':
-	    printf("matrix.xx = %f;\n", info->dim->matrix.xx);
-	    printf("matrix.xy = %f;\n", info->dim->matrix.xy);
-	    printf("matrix.x0 = %f;\n", info->dim->matrix.x0);
-	    printf("matrix.yx = %f;\n", info->dim->matrix.yx);
-	    printf("matrix.yy = %f;\n", info->dim->matrix.yy);
-	    printf("matrix.y0 = %f;\n", info->dim->matrix.y0);
-	    break;
-	case 'i':
-	    output_info(screen_context);
-	    break;
-	case 'b':
-	    TOGGLE(screen_context->show_boxes);
-	    break;
-	case 'B':
-	    TOGGLE(screen_context->show_boxes2);
-	    break;
-	case 'a':
-	    TOGGLE(screen_context->show_attractors);
-	    break;
-	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);
-	    computeLimitCurve(screen_context);
-	    break;
-	case 'p':
-	    print(screen_context);
-	    break;
-	case 'M':
-	    screen_context->limit_with_lines = 0;
-	    double parameter_start = screen_context->parameter;
-	    double parameter2_start = screen_context->parameter2;
-	    for(int i = 0; i <= screen_context->movie_n_frames; i++) {
-		screen_context->parameter = SIGN(parameter_start)*exp(log(fabs(parameter_start)) +
-						i*screen_context->movie_parameter_duration/screen_context->movie_n_frames);
-		screen_context->parameter2 = SIGN(parameter2_start)*exp(log(fabs(parameter2_start)) +
-						 i*screen_context->movie_parameter2_duration/screen_context->movie_n_frames);
-		updateMatrices(screen_context);
-		computeLimitCurve(screen_context);
-		draw(screen_context);
-		sprintf(filename, "output/%s%03d.png", screen_context->movie_filename, i);
-		cairo_surface_write_to_png(info->buffer_surface, filename);
-		printf("Finished drawing %s\n", filename);
-	    }
+	case KeyPress:
+		state = ev->xkey.state & (ShiftMask | LockMask | ControlMask);
+		key = XkbKeycodeToKeysym(ev->xkey.display, ev->xkey.keycode, 0, !!(state & ShiftMask));
+//		printf("Key pressed: %ld\n", key);
 
-	case 'f':
-	    TOGGLE(screen_context->use_repelling);
-	    computeLimitCurve(screen_context);
-	    break;
-	case 't':
-	    TOGGLE(screen_context->show_text);
-	    break;
-	case 'c':
-	    TOGGLE(screen_context->show_coxeter_orbit);
-	    break;
-	case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9': case '0':
-	    screen_context->mode = key - '0';
-	    break;
+		switch(key) {
+		case XK_Down:
+			screen_context->parameter /= exp(0.002);
+			updateMatrices(screen_context);
+			computeLimitCurve(screen_context);
+			break;
+		case XK_Up:
+			screen_context->parameter *= exp(0.002);
+			updateMatrices(screen_context);
+			computeLimitCurve(screen_context);
+			break;
+		case XK_Left:
+			screen_context->parameter /= exp(0.00002);
+			updateMatrices(screen_context);
+			computeLimitCurve(screen_context);
+			break;
+		case XK_Right:
+			screen_context->parameter *= exp(0.00002);
+			updateMatrices(screen_context);
+			computeLimitCurve(screen_context);
+			break;
+		case XK_Page_Down:
+			screen_context->parameter /= exp(0.02);
+			updateMatrices(screen_context);
+			computeLimitCurve(screen_context);
+			break;
+		case XK_Page_Up:
+			screen_context->parameter *= exp(0.02);
+			updateMatrices(screen_context);
+			computeLimitCurve(screen_context);
+			break;
+		case ' ':
+			screen_context->parameter = 5.57959706;
+			updateMatrices(screen_context);
+			computeLimitCurve(screen_context);
+			break;
+		case XK_Return:
+//			screen_context->parameter = 2.76375163;
+			screen_context->parameter = 5.29063366;
+			updateMatrices(screen_context);
+			computeLimitCurve(screen_context);
+			break;
+		case 'm':
+			printf("matrix.xx = %f;\n", info->dim->matrix.xx);
+			printf("matrix.xy = %f;\n", info->dim->matrix.xy);
+			printf("matrix.x0 = %f;\n", info->dim->matrix.x0);
+			printf("matrix.yx = %f;\n", info->dim->matrix.yx);
+			printf("matrix.yy = %f;\n", info->dim->matrix.yy);
+			printf("matrix.y0 = %f;\n", info->dim->matrix.y0);
+			break;
+		case 'i':
+			output_info(screen_context);
+			break;
+		case 'b':
+			TOGGLE(screen_context->show_boxes);
+			break;
+		case 'B':
+			TOGGLE(screen_context->show_boxes2);
+			break;
+		case 'a':
+			TOGGLE(screen_context->show_attractors);
+			break;
+		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);
+			computeLimitCurve(screen_context);
+			break;
+		case 'p':
+			print(screen_context);
+			break;
+		case 'P':
+			for(int i = 0; i <= 100; i++) {
+				for(int j = 0; j <= 100; j++) {
+					screen_context->distance_parameter1 = 0.02*i;
+					screen_context->distance_parameter2 = 0.02*j;
+					draw(screen_context);
+				}
+				fflush(stdout);
+			}
+			break;
+		case 'M':
+			/*
+			screen_context->limit_with_lines = 0;
+			double parameter_start = screen_context->parameter;
+			for(int i = 0; i <= 1300; i++) {
+				if(i < 400)
+					screen_context->parameter = exp(log(parameter_start)+0.002*i);
+				else if(i < 500)
+					screen_context->parameter = exp(log(parameter_start)+0.002*400);
+				else
+					screen_context->parameter = exp(log(parameter_start)+0.002*(900-i));
+				updateMatrices(screen_context);
+				computeLimitCurve(screen_context);
+				draw(screen_context);
+				sprintf(filename, "movie3/test%03d.png", i);
+				cairo_surface_write_to_png(info->buffer_surface, filename);
+				printf("Finished drawing %s\n", filename);
+			}
+			*/
+			screen_context->limit_with_lines = 0;
+			double parameter_start = screen_context->parameter;
+			for(int i = 0; i <= 1300; i++) {
+				if(i < 400)
+					screen_context->parameter = exp(0.003*i);
+				else if(i < 500)
+					screen_context->parameter = exp(0.003*400);
+				else
+					screen_context->parameter = exp(0.003*(900-i));
+				updateMatrices(screen_context);
+				computeLimitCurve(screen_context);
+				draw(screen_context);
+				sprintf(filename, "movie5/test%03d.png", i);
+				cairo_surface_write_to_png(info->buffer_surface, filename);
+				printf("Finished drawing %s\n", filename);
+			}
+
+		case 'f':
+			TOGGLE(screen_context->use_repelling);
+			computeLimitCurve(screen_context);
+			break;
+		case 't':
+			TOGGLE(screen_context->show_text);
+			break;
+		case 'c':
+			TOGGLE(screen_context->show_coxeter_orbit);
+			break;
+		case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9': case '0':
+			screen_context->mode = key - '0';
+			break;
+		}
+
+		return STATUS_REDRAW;
 	}
 
-	return STATUS_REDRAW;
-    }
-
-    return STATUS_NOTHING;
+	return STATUS_NOTHING;
 }
 
 int main(int argc, char *argv[])
 {
-    GraphicsInfo *info;
+	GraphicsInfo *info;
 
-    screen_context = malloc(sizeof(DrawingContext));
-    setupContext(screen_context, argc, argv);
-    updateMatrices(screen_context);
-    computeLimitCurve(screen_context);
+	screen_context = malloc(sizeof(DrawingContext));
+	setupContext(screen_context, argc, argv);
+	updateMatrices(screen_context);
+	computeLimitCurve(screen_context);
 
-    info = initCairo(0, KeyPressMask, 200, 200, "Triangle group");
-    if(!info)
-	return 1;
+	info = initCairo(0, KeyPressMask, 200, 200, "Triangle group");
+	if(!info)
+		return 1;
 
-    /*
-      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;
-    */
+/*
+	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;
+*/
 
-    info->dim->matrix.xx = 274.573171;
-    info->dim->matrix.xy = 0.000000;
-    info->dim->matrix.x0 = 910.073462;
-    info->dim->matrix.yx = 0.000000;
-    info->dim->matrix.yy = 274.573171;
-    info->dim->matrix.y0 = 509.225293;
+	info->dim->matrix.xx = 274.573171;
+	info->dim->matrix.xy = 0.000000;
+	info->dim->matrix.x0 = 910.073462;
+	info->dim->matrix.yx = 0.000000;
+	info->dim->matrix.yy = 274.573171;
+	info->dim->matrix.y0 = 509.225293;
 
-    updateDimensions(info->dim);
+	updateDimensions(info->dim);
 
-    screen_context->dim = info->dim;
-    screen_context->cairo = info->buffer_context;
+	screen_context->dim = info->dim;
+	screen_context->cairo = info->buffer_context;
 
-    startTimer(info);
+	startTimer(info);
 
-    while(1) {
-	int result = checkEvents(info, processEvent, NULL);
-	if(result == STATUS_QUIT)
-	    return 0;
-	else if(result == STATUS_REDRAW) {
-	    struct timeval current_time;
-	    double start_time, intermediate_time, end_time;
-	    gettimeofday(&current_time, 0);
-	    start_time = current_time.tv_sec + current_time.tv_usec*1e-6;
+	while(1) {
+		int result = checkEvents(info, processEvent, NULL);
+		if(result == STATUS_QUIT)
+			return 0;
+		else if(result == STATUS_REDRAW) {
+			struct timeval current_time;
+			double start_time, intermediate_time, end_time;
+			gettimeofday(&current_time, 0);
+			start_time = current_time.tv_sec + current_time.tv_usec*1e-6;
 
-	    draw(screen_context);
+			draw(screen_context);
 
-	    gettimeofday(&current_time, 0);
-	    intermediate_time = current_time.tv_sec + current_time.tv_usec*1e-6;
+			gettimeofday(&current_time, 0);
+			intermediate_time = current_time.tv_sec + current_time.tv_usec*1e-6;
 
-	    cairo_set_source_surface(info->front_context, info->buffer_surface, 0, 0);
-	    cairo_paint(info->front_context);
+			cairo_set_source_surface(info->front_context, info->buffer_surface, 0, 0);
+			cairo_paint(info->front_context);
 
-	    gettimeofday(&current_time, 0);
-	    end_time = current_time.tv_sec + current_time.tv_usec*1e-6;
-	    printf("drawing finished in %.2f milliseconds, of which %.2f milliseconds were buffer switching\n", (end_time - start_time) * 1000, (end_time - intermediate_time) * 1000);
+			gettimeofday(&current_time, 0);
+			end_time = current_time.tv_sec + current_time.tv_usec*1e-6;
+//			printf("drawing finished in %.2f milliseconds, of which %.2f milliseconds were buffer switching\n", (end_time - start_time) * 1000, (end_time - intermediate_time) * 1000);
+		}
+		waitUpdateTimer(info);
 	}
-	waitUpdateTimer(info);
-    }
 
-    free(screen_context);
-    destroyCairo(info);
-    destroyContext(screen_context);
+	free(screen_context);
+	destroyCairo(info);
+	destroyContext(screen_context);
 
-    return 0;
+	return 0;
 }

main.h

@@ -11,12 +11,7 @@
 #define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);}
 #define LOOP(i) for(int i = 0; i < 3; i++)
 
-#define NUM_GROUP_ELEMENTS 10000
-#define NUM_GROUP_ELEMENTS_COMBINATORIAL 100000
-
-// (0,1) -> 2, (1,2) -> 0, (2,0) -> 1
-// (1,0) -> 5, (2,1) -> 3, (0,2) -> 4
-#define ROTATION_LETTER(x,y) (((y)-(x)+3)%3 == 1 ? ((y)+1)%3 : ((x)+1)%3+3)
+#define NUM_GROUP_ELEMENTS 50000
 
 typedef struct {
 	double x[3];
@@ -27,6 +22,14 @@ typedef struct {
 	double y;
 } point_t;
 
+typdef struct {
+	double x;
+	double y;
+	double angle;
+	vector_t fp[3];
+	groupelement_t *g;
+} limit_curve_t;
+
 typedef struct {
 	// infos about the screen to draw on
 	cairo_t *cairo;
@@ -35,13 +38,7 @@ typedef struct {
 	// a priori parameter
 	int p[3],k[3];
 	double parameter;
-	double parameter2;
-	char *movie_filename;
-	double movie_parameter_duration;
-	double movie_parameter2_duration;
-	int movie_n_frames;
 	int n_group_elements;
-	int n_group_elements_combinatorial;
 	int show_boxes;
 	int show_boxes2;
 	int show_attractors;
@@ -55,9 +52,14 @@ typedef struct {
 	int limit_with_lines;
 	int show_marking;
 	point_t marking;
+	point_t marking2;
+	point_t marking3;
+	double distance_parameter1;
+	double distance_parameter2;
 	int show_coxeter_orbit;
 	int use_repelling;
 	gsl_matrix *cartan, *cob;
+	char *extra_text;
 
     // precomputed and constant
 	groupelement_t *group;
@@ -77,8 +79,7 @@ typedef enum {
 
 // implemented in limit_set.c
 void cartanMatrix(gsl_matrix *cartan, double a1, double a2, double a3, double s);
-void initializeTriangleGenerators(gsl_matrix **gen, double a1, double a2, double a3, double s, double t, workspace_t *ws);
-void initializeTriangleGeneratorsCurrent(gsl_matrix **gen, DrawingContext *ctx);
+void initializeTriangleGenerators(gsl_matrix **gen, gsl_matrix *cartan);
 int computeLimitCurve(DrawingContext *ctx);
 
 // implemented in draw.c
@@ -108,8 +109,7 @@ void setupContext(DrawingContext *ctx, int argc, char *argv[]);
 void destroyContext(DrawingContext *ctx);
 void print(DrawingContext *screen);
 int processEvent(GraphicsInfo *info, XEvent *ev);
-void computeMatrix(DrawingContext *ctx, gsl_matrix *result, const char *type);
-void computeRotationMatrixFrame(DrawingContext *ctx, gsl_matrix *result, const char *type);
+void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char *type);
 void updateMatrices(DrawingContext *ctx);
 
 static vector_t vectorFromGsl(gsl_vector *v)

triangle.h

@@ -14,7 +14,6 @@ typedef struct _groupelement {
 	struct _groupelement *parent;
 	struct _groupelement *inverse;
 	int letter;
-	int visited;
 } groupelement_t;
 
 int generate_triangle_group(groupelement_t *group, int size, int k1, int k2, int k3);