movie command line arguments

draw.c

@@ -46,21 +46,22 @@ int fixedPoints(DrawingContext *ctx, const char *word, vector_t *out)
 {
 	gsl_matrix *tmp = getTempMatrix(ctx->ws);
 	gsl_matrix *ev = getTempMatrix(ctx->ws);
-	gsl_matrix **gen = getTempMatrices(ctx->ws, 3);
+	gsl_matrix **gen = getTempMatrices(ctx->ws, 6);
 
-	initializeTriangleGenerators(gen, ctx->cartan);
+	initializeTriangleGeneratorsCurrent(gen, ctx);
 
 	gsl_matrix_set_identity(tmp);
 	for(int i = 0; i < strlen(word); i++) {
-		if(word[i] == ' ')
+		if(word[i] >= 'a' && word[i] <= 'c')
-			continue;
 			multiply_right(tmp, gen[word[i]-'a'], ctx->ws);
+		else if(word[i] >= 'A' && word[i] <= 'C')
+			multiply_right(tmp, gen[word[i]-'A'+3], ctx->ws);
 	}
 	int count = real_eigenvectors(tmp, ev, ctx->ws);
 
 	LOOP(i) LOOP(j) out[i].x[j] = gsl_matrix_get(ev, j, i);
 
-	releaseTempMatrices(ctx->ws, 5);
+	releaseTempMatrices(ctx->ws, 8);
 
 	return count;
 }
@@ -69,9 +70,9 @@ int wordEigenvalues(DrawingContext *ctx, const char *word, double *out)
 {
 	gsl_matrix *tmp = getTempMatrix(ctx->ws);
 	gsl_vector *ev = getTempVector(ctx->ws);
-	gsl_matrix **gen = getTempMatrices(ctx->ws, 3);
+	gsl_matrix **gen = getTempMatrices(ctx->ws, 6);
 
-	initializeTriangleGenerators(gen, ctx->cartan);
+	initializeTriangleGeneratorsCurrent(gen, ctx);
 
 	gsl_matrix_set_identity(tmp);
 	for(int i = 0; i < strlen(word); i++) {
@@ -83,7 +84,7 @@ int wordEigenvalues(DrawingContext *ctx, const char *word, double *out)
 
 	LOOP(i) out[i] = gsl_vector_get(ev, i);
 
-	releaseTempMatrices(ctx->ws, 4);
+	releaseTempMatrices(ctx->ws, 7);
 	releaseTempVectors(ctx->ws, 1);
 
 	return count;
@@ -298,8 +299,6 @@ void drawBoxLines(DrawingContext *ctx, const char *word1, const char *word2)
 	drawPolygon(ctx, 0, 4, p[0][0], i[0], p[1][0], i[1]);
 }
 
-
-
 void drawBoxStd(DrawingContext *ctx, const char *word, char base)
 {
 	char word1[100];
@@ -329,13 +328,13 @@ void drawRotationOrbitFrame(DrawingContext *ctx, gsl_matrix *frame, vector_t sta
 	vector_t v[3], w;
 	point_t p;
 	double parameter, startangle;
-	int iterations = 200;
+	int iterations = 2000;
 	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);
+//	computeRotationMatrixFrame(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]);
@@ -371,7 +370,7 @@ void drawRotationOrbit(DrawingContext *ctx, const char *word, vector_t start)
 {
 	gsl_matrix *frame = getTempMatrix(ctx->ws);
 
-	computeRotationMatrix(ctx, frame, word);
+	computeRotationMatrixFrame(ctx, frame, word);
 	drawRotationOrbitFrame(ctx, frame, start);
 
 	releaseTempMatrices(ctx->ws, 1);
@@ -389,7 +388,7 @@ void drawDualRotationOrbit(DrawingContext *ctx, const char *word, vector_t start
 	gsl_vector *start_in_frame = getTempVector(ctx->ws);
 	cairo_t *C = ctx->cairo;
 
-	computeRotationMatrix(ctx, frame, word);
+	computeRotationMatrixFrame(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]);
@@ -434,7 +433,7 @@ void drawArcWithOutput(DrawingContext *ctx, const char *word, vector_t start, ve
 	gsl_vector *vector_in_frame = getTempVector(ctx->ws);
 	cairo_t *C = ctx->cairo;
 
-	computeRotationMatrix(ctx, frame, word);
+	computeRotationMatrixFrame(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]);
@@ -576,9 +575,9 @@ void drawAttractors(DrawingContext *ctx)
 	vector_t p[6][3];
 	vector_t l[6][3];
 
-	fixedPoints(ctx, "abc", p[0]);
+	fixedPoints(ctx, "cba", p[0]);
-	fixedPoints(ctx, "bca", p[1]);
+	fixedPoints(ctx, "bac", p[1]);
-	fixedPoints(ctx, "cab", p[2]);
+	fixedPoints(ctx, "acb", p[2]);
 	fixedPoints(ctx, "a cab a", p[3]);
 	fixedPoints(ctx, "b abc b", p[4]);
 	fixedPoints(ctx, "c bca c", p[5]);
@@ -705,200 +704,180 @@ void drawCurvedBox(DrawingContext *ctx, int base, const char *conj, int style)
 	}
 }
 
+groupelement_t *left(const char *word, groupelement_t *g)
+{
+	int n = strlen(word);
+
+	for(int i = n-1; i >= 0; i--) {
+		if(word[i] == ' ')
+			continue;
+
+		g = g->adj[word[i]-'a'];
+
+		if(!g)
+			break;
+	}
+
+	return g;
+}
 
 void drawBoxes(DrawingContext *ctx)
 {
-	gsl_matrix *rot = getTempMatrix(ctx->ws);
+	gsl_matrix *tmp = getTempMatrix(ctx->ws);
-	gsl_matrix **gen = getTempMatrices(ctx->ws, 3);
+	gsl_matrix **gen = getTempMatrices(ctx->ws, 6);
 	gsl_matrix *frame = getTempMatrix(ctx->ws);
+	gsl_matrix *frame2 = getTempMatrix(ctx->ws);
+	gsl_vector *startpoint_drawbasis = getTempVector(ctx->ws);
+	gsl_vector *startpoint_globalbasis = getTempVector(ctx->ws);
+	gsl_matrix **elements = getTempMatrices(ctx->ws, ctx->n_group_elements);
+
 	cairo_t *C = ctx->cairo;
 	cairo_save(C);
-
-	vector_t p[22][3];
-	vector_t l[22][3];
-	vector_t alpha[6];
-	vector_t ptmp[3];
-	char word[100], word2[100];
-
-	fixedPoints(ctx, "abc", p[0]);
-	fixedPoints(ctx, "bca", p[1]);
-	fixedPoints(ctx, "cab", p[2]);
-	fixedPoints(ctx, "bacabab", p[3]);
-	fixedPoints(ctx, "bcacabacb", p[4]);
-
 	cairo_set_line_width(C, 2.0/ctx->dim->scalefactor);
 	cairo_set_source_rgb(C, 0.6, 0.6, 0.6);
 
-	drawRotationOrbit(ctx, "ab", p[0][0]);
+	vector_t p[22][3];
-	drawRotationOrbit(ctx, "bc", p[0][0]);
+	vector_t fp[3];
-	drawRotationOrbit(ctx, "ca", p[0][0]);
+	vector_t l[22][3];
+	vector_t alpha[6];
+	vector_t ptmp[3];
+	vector_t start;
+	vector_t start2;
+	char word[100], word2[100];
 
-	if(ctx->mode >= 2) {
+	fixedPoints(ctx, "cba", p[0]);
-		cairo_set_source_rgb(C, 0.6, 0.6, 1);
+	fixedPoints(ctx, "acb", p[1]);
-		drawRotationOrbit(ctx, "bcabcb", p[1][0]); // bcC
+	fixedPoints(ctx, "bac", p[2]);
-		drawRotationOrbit(ctx, "abcabcba", p[0][0]); // abcC
+
-//		drawRotationOrbit(ctx, "bcabcabacb", p[1][0]); // bcabC''
+	initializeTriangleGeneratorsCurrent(gen, ctx);
-//		drawRotationOrbit(ctx, "bacabcacab", p[3][0]); // bacaC'
+	gsl_matrix_set_identity(elements[0]);
-//		drawRotationOrbit(ctx, "bcacabcacacb", p[4][0]); // bcacaC' bcacabacb
+	for(int i = 1; i < ctx->n_group_elements; i++) {
+		if(ctx->group[i].length % 2)
+			continue;
+		int letter = ROTATION_LETTER(ctx->group[i].letter, ctx->group[i].parent->letter);
+		multiply(gen[letter], elements[ctx->group[i].parent->parent->id], elements[i]);
 	}
 
-	cairo_set_source_rgb(C, 1, 0, 1);
+	gsl_vector_set(startpoint_drawbasis, 0, ctx->marking.x);
-//	drawRotationOrbit(ctx, "bacabcacab", p[3][0]); // ababcba
+	gsl_vector_set(startpoint_drawbasis, 1, ctx->marking.y);
-	cairo_set_source_rgb(C, 0, 0, 0);
+	gsl_vector_set(startpoint_drawbasis, 2, 1);
-	fixedPoints(ctx, "abababcbaba", p[3]);
+	solve(ctx->cob, startpoint_drawbasis, startpoint_globalbasis, ctx->ws);
-//	drawRotationOrbit(ctx, "abababcabcbababa", p[3][0]); //  bab abc bab
+	LOOP(i) start.x[i] = gsl_vector_get(startpoint_globalbasis, i);
-
-	cairo_set_source_rgb(C, 0, 0, 1);
-//	fixedPoints(ctx, "cab", p[3]);
-//	drawRotationOrbit(ctx, "cabc", p[3][0]);
-//	fixedPoints(ctx, "bca", p[3]);
-//	drawRotationOrbit(ctx, "bcabcb", p[3][0]);
-//	fixedPoints(ctx, "bc abc cb", p[3]);
-//	drawRotationOrbit(ctx, "bcabcb", p[3][0]);
-
-	fixedPoints(ctx, "bc bca cb", p[3]);
-//	drawRotationOrbit(ctx, "bcbcacbc", p[3][0]);
 
 	/*
-	cairo_set_source_rgb(C, 0, 0, 0);
+	fixedPoints(ctx, "ABAB", fp);
-	strncpy(word,"abc",100);
+	drawCovector(ctx, cross(fp[0], fp[2]));
-	for(int i = 0; i < 9; i++) {
+	drawVector(ctx, fp[1]);
-		conjugate_word(word, 0, "ab", word2);
+	computeRotationMatrixFrame(ctx, frame, "c");
-		strncpy(word, word2, 100);
+	drawRotationOrbit(ctx, "c", start); // ba cb ab ac = C A c B
-		fixedPoints(ctx, word, ptmp);
+	drawVector(ctx, start);
-		drawVector(ctx, ptmp[0]);
+
-//		drawVector(ctx, ptmp[2]);
+	computeMatrix(ctx, tmp, "ABABc");
+	start2 = apply(tmp, start);
+	computeRotationMatrixFrame(ctx, frame, "c");
+	multiply(tmp, frame, frame2);
+	drawVector(ctx, start2);
+	drawRotationOrbitFrame(ctx, frame2, apply(tmp, start));
+
+	computeMatrix(ctx, tmp, "ABABcABABC");
+	start2 = apply(tmp, start);
+	computeRotationMatrixFrame(ctx, frame, "c");
+	multiply(tmp, frame, frame2);
+	drawVector(ctx, start2);
+	drawRotationOrbitFrame(ctx, frame2, apply(tmp, start));
+	*/
+
+	for(groupelement_t *cur = &ctx->group[180]; cur->parent; cur = cur->parent)
+		fputc('a'+cur->letter, stdout);
+	fputc('\n', stdout);
+
+	queue_t queue;
+	queue_init(&queue);
+	queue_put(&queue, 0);
+	int current;
+	groupelement_t *cur, *next;
+
+	for(int i = 0; i < ctx->n_group_elements_combinatorial; i++)
+		ctx->group[i].visited = 0;
+
+	cur = &ctx->group[0];
+
+	computeRotationMatrixFrame(ctx, frame, "c");
+
+	for(int i = 0; i < 1000; i++) {
+		if(ctx->group[i].length % 2 != 0)
+			continue;
+
+		multiply(elements[i], frame, frame2);
+		drawRotationOrbitFrame(ctx, frame2, apply(elements[i], start));
 	}
 
-	strncpy(word,"bca",100);
+	/*
-	for(int i = 0; i < 9; i++) {
+	while((current = queue_get(&queue)) != -1) {
-		conjugate_word(word, 0, "ab", word2);
+		cur = &ctx->group[current];
-		strncpy(word, word2, 100);
+
-		fixedPoints(ctx, word, ptmp);
+		if(cur->visited > 4)
-		drawVector(ctx, ptmp[0]);
+			continue;
-//		drawVector(ctx, ptmp[2]);
+
+		if(cur->id < ctx->n_group_elements) {
+			multiply(elements[cur->id], frame, frame2);
+			drawRotationOrbitFrame(ctx, frame2, apply(elements[cur->id], start));
 		}
 
-
+	next = left("ab ab", cur);
-	strncpy(word,"abc",100);
+	if(next && next->visited == 0) {
-	for(int i = 0; i < 9; i++) {
+		queue_put(&queue, next->id);
-		conjugate_word(word, 0, "bc", word2);
+		next->visited = cur->visited+1;
-		strncpy(word, word2, 100);
-		fixedPoints(ctx, word, ptmp);
-		drawVector(ctx, ptmp[0]);
 	}
-
+	next = left("cbac cbac", cur);
-	strncpy(word,"cab",100);
+	if(next && next->visited == 0) {
-	for(int i = 0; i < 9; i++) {
+		queue_put(&queue, next->id);
-		conjugate_word(word, 0, "bc", word2);
+		next->visited = cur->visited+1;
-		strncpy(word, word2, 100);
-		fixedPoints(ctx, word, ptmp);
-		drawVector(ctx, ptmp[0]);
 	}
-
+	next = left("cacbca cacbca", cur);
-	strncpy(word,"cab",100);
+	if(next && next->visited == 0) {
-	for(int i = 0; i < 9; i++) {
+		queue_put(&queue, next->id);
-		conjugate_word(word, 0, "ca", word2);
+		next->visited = cur->visited+1;
-		strncpy(word, word2, 100);
+	}
-		fixedPoints(ctx, word, ptmp);
+	next = left("cabcacbcac cabcacbcac", cur);
-		drawVector(ctx, ptmp[0]);
+	if(next && next->visited == 0) {
+		queue_put(&queue, next->id);
+		next->visited = cur->visited+1;
+	}
+	next = left("acbcacba acbcacba", cur);
+	if(next && next->visited == 0) {
+		queue_put(&queue, next->id);
+		next->visited = cur->visited+1;
+	}
+	next = left("bcacbc bcacbc", cur);
+	if(next && next->visited == 0) {
+		queue_put(&queue, next->id);
+		next->visited = cur->visited+1;
 	}
-
-	strncpy(word,"abc",100);
-	for(int i = 0; i < 9; i++) {
-		conjugate_word(word, 0, "ca", word2);
-		strncpy(word, word2, 100);
-		fixedPoints(ctx, word, ptmp);
-		drawVector(ctx, ptmp[0]);
 	}
 	*/
 
-	/*
+//	drawRotationOrbit(ctx, "a", p[1][0]);
-	cairo_set_source_rgb(C, 1, 0, 0);
+//	drawRotationOrbit(ctx, "b", p[2][0]);
-	drawVector(ctx, p[0][0]);
-	cairo_set_source_rgb(C, 0, 0.6, 0);
-	drawVector(ctx, p[1][0]);
-	cairo_set_source_rgb(C, 0, 0, 1);
-	drawVector(ctx, p[2][0]);
-	*/
-
-	/*
-	fixedPoints(ctx, "ab abc ba", p[4]);
-	fixedPoints(ctx, "abab abc baba", p[5]);
-	fixedPoints(ctx, "ababab abc bababa", p[6]);
-	fixedPoints(ctx, "abababab abc babababa", p[7]);
-	fixedPoints(ctx, "babababa abc abababab", p[8]);
-	fixedPoints(ctx, "bababa abc ababab", p[9]);
-	fixedPoints(ctx, "baba abc abab", p[10]);
-	fixedPoints(ctx, "ba abc ab", p[11]);
-	fixedPoints(ctx, "bca", p[12]);
-	fixedPoints(ctx, "b abc b", p[13]);
-	fixedPoints(ctx, "bab abc bab", p[14]);
-	fixedPoints(ctx, "babab abc babab", p[15]);
-	fixedPoints(ctx, "bababab abc bababab", p[16]);
-	fixedPoints(ctx, "abababab bca babababa", p[17]);
-	fixedPoints(ctx, "ababab bca bababa", p[18]);
-	fixedPoints(ctx, "abab bca baba", p[19]);
-	fixedPoints(ctx, "ab bca ba", p[20]);
-	*/
-
-//	initializeTriangleGenerators(gen, ctx->cartan);
-
-//	for(int i = 0; i < 22; 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);
-
-	/*
-	gsl_matrix_set(frame, 0, 0, 2.0);
-	gsl_matrix_set(frame, 0, 1, 0.0);
-	gsl_matrix_set(frame, 0, 2, 1.0);
-	gsl_matrix_set(frame, 1, 0, -1.0);
-	gsl_matrix_set(frame, 1, 1, sqrt(3));
-	gsl_matrix_set(frame, 1, 2, 1.0);
-	gsl_matrix_set(frame, 2, 0, -1.0);
-	gsl_matrix_set(frame, 2, 1, -sqrt(3));
-	gsl_matrix_set(frame, 2, 2, 1.0);*/
-//	drawRotationOrbitFrame(ctx, frame, p[0][0]);
-
-//	drawRotationOrbit(ctx, "bc", p[0][0]);
-
-//	drawRotationOrbit(ctx, "ca", p[0][0]);
-
-/*
-	for(int i = 0; i < 18; i++) {
-		if(i == 0)
-			cairo_set_source_rgb(C, 1, 0, 0);
-		else if(i == 8)
-			cairo_set_source_rgb(C, 0, 0, 1);
-		else if(i == 9)
-			cairo_set_source_rgb(C, 0, 0.6, 0);
-		else
-			cairo_set_source_rgb(C, 0, 0, 0);
-		drawVector(ctx, p[3+i][0]);
-//		drawCovector(ctx, l[3+i][0]);
-}
-*/
-
-//	drawRotationOrbit(ctx, "ab", cross(l[0][0], l[2][1]));
-
-//	drawRotationOrbit(ctx, "abca", p[0][0]);
 
 	cairo_restore(C);
-	releaseTempMatrices(ctx->ws, 5);
+	releaseTempMatrices(ctx->ws, 9 + ctx->n_group_elements);
+	releaseTempVectors(ctx->ws, 2);
 }
 
 void drawBoxes2(DrawingContext *ctx)
 {
 	gsl_matrix *rot = getTempMatrix(ctx->ws);
-	gsl_matrix **gen = getTempMatrices(ctx->ws, 3);
+	gsl_matrix **gen = getTempMatrices(ctx->ws, 6);
 	cairo_t *C = ctx->cairo;
 	cairo_save(C);
-	initializeTriangleGenerators(gen, ctx->cartan);
+	initializeTriangleGeneratorsCurrent(gen, ctx);
 
 	vector_t p[4][3];
 
-	fixedPoints(ctx, "abc", p[0]);
+	fixedPoints(ctx, "cba", p[0]);
-	fixedPoints(ctx, "bca", p[1]);
+	fixedPoints(ctx, "acb", p[1]);
-	fixedPoints(ctx, "cab", p[2]);
+	fixedPoints(ctx, "bac", p[2]);
 
 	cairo_set_line_width(C, 2.5/ctx->dim->scalefactor);
 
@@ -1262,13 +1241,13 @@ void drawBoxes2(DrawingContext *ctx)
 */
 
 	cairo_restore(C);
-	releaseTempMatrices(ctx->ws, 4);
+	releaseTempMatrices(ctx->ws, 7);
 }
 
 void drawRotatedReflectors(DrawingContext *ctx)
 {
 	gsl_matrix *rot = getTempMatrix(ctx->ws);
-	gsl_matrix **gen = getTempMatrices(ctx->ws, 3);
+	gsl_matrix **gen = getTempMatrices(ctx->ws, 6);
 	cairo_t *C = ctx->cairo;
 	vector_t fp[3], fp2[3];
 	vector_t w;
@@ -1276,7 +1255,7 @@ void drawRotatedReflectors(DrawingContext *ctx)
 
 	cairo_set_source_rgb(C, 0.7, 0.7, 0.7);
 
-	initializeTriangleGenerators(gen, ctx->cartan);
+	initializeTriangleGeneratorsCurrent(gen, ctx);
 
 	LOOP(i) LOOP(j) v[i].x[j] = gsl_matrix_get(ctx->cartan, i, j);
 	multiply(gen[0], gen[1], rot);
@@ -1297,7 +1276,7 @@ void drawRotatedReflectors(DrawingContext *ctx)
 	fixedPoints(ctx, "cacabac", fp2);
 	drawRotationOrbit(ctx, "ac", fp[0]);
 
-	releaseTempMatrices(ctx->ws, 4);
+	releaseTempMatrices(ctx->ws, 7);
 }
 
 void drawDualLimitCurve(DrawingContext *ctx)
@@ -1307,12 +1286,36 @@ void drawDualLimitCurve(DrawingContext *ctx)
 	cairo_save(C);
 	cairo_set_source_rgb(C, 0.5, 0.5, 1);
 
-	int n = 18;
+	gsl_matrix **gen = getTempMatrices(ctx->ws, 6);
-	vector_t p[n][3];
+	gsl_matrix **elements = getTempMatrices(ctx->ws, ctx->n_group_elements);
-	vector_t l[n][3];
-	vector_t ptmp[3], ltmp[3];
 
-	fixedPoints(ctx, "abc", p[0]);
+//	wordEigenvalues(ctx, "abc", ev);
+//	LOOP(i) LOOP(j) gsl_matrix_set(coxeter_fixedpoints, j, i, cox[0][i].x[j]);
+
+	initializeTriangleGeneratorsCurrent(gen, ctx);
+	gsl_matrix_set_identity(elements[0]);
+	for(int i = 1; i < ctx->n_group_elements; i++)
+		multiply(gen[ctx->group[i].letter], elements[ctx->group[i].parent->id], elements[i]);
+
+	vector_t p[3], l[3], v;
+
+	fixedPoints(ctx, "cba", p);
+	drawVector(ctx, p[0]);
+	drawVector(ctx, p[1]);
+	drawVector(ctx, p[2]);
+
+	LOOP(i) l[i] = cross(p[(i+1)%3], p[(i+2)%3]);
+
+	for(int i = 0; i < ctx->n_group_elements; i++) {
+		v = apply_transpose(elements[i], l[0]);
+		drawCovector(ctx, v);
+	}
+
+	releaseTempMatrices(ctx->ws, 3 + ctx->n_group_elements);
+//	releaseTempVectors(ctx->ws, 4);
+
+
+/*
 	fixedPoints(ctx, "ab abc ba", p[1]);
 	fixedPoints(ctx, "abab abc baba", p[2]);
 	fixedPoints(ctx, "ababab abc bababa", p[3]);
@@ -1331,6 +1334,7 @@ void drawDualLimitCurve(DrawingContext *ctx)
 	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]);
@@ -1368,6 +1372,7 @@ void drawDualLimitCurve(DrawingContext *ctx)
 //		drawCovector(ctx, l[i][2]);
 }*/
 
+	/*
 	fixedPoints(ctx, "abc", ptmp);
 	drawCovector(ctx, cross(ptmp[0], ptmp[1]));
 	fixedPoints(ctx, "bca", ptmp);
@@ -1380,6 +1385,7 @@ void drawDualLimitCurve(DrawingContext *ctx)
 	drawCovector(ctx, cross(ptmp[0], ptmp[1]));
 	fixedPoints(ctx, "acaba", ptmp);
 	drawCovector(ctx, cross(ptmp[0], ptmp[1]));
+	*/
 
 	cairo_restore(C);
 }
@@ -1435,7 +1441,7 @@ void drawLimitCurve(DrawingContext *ctx)
 void drawCoxeterOrbit(DrawingContext *ctx)
 {
 	gsl_matrix *rot = getTempMatrix(ctx->ws);
-	gsl_matrix **gen = getTempMatrices(ctx->ws, 3);
+	gsl_matrix **gen = getTempMatrices(ctx->ws, 6);
 	gsl_vector *eval = getTempVector(ctx->ws);
 	gsl_matrix *coxeter_fixedpoints = getTempMatrix(ctx->ws);
 	gsl_vector *startpoint_coxeterbasis = getTempVector(ctx->ws);
@@ -1452,7 +1458,7 @@ void drawCoxeterOrbit(DrawingContext *ctx)
 	int first = 1;
 
 	cairo_save(C);
-	initializeTriangleGenerators(gen, ctx->cartan);
+	initializeTriangleGeneratorsCurrent(gen, ctx);
 
 	cairo_set_source_rgb(C, 0, 0, 1);
 
@@ -1463,7 +1469,7 @@ void drawCoxeterOrbit(DrawingContext *ctx)
 	wordEigenvalues(ctx, "abc", ev);
 	LOOP(i) LOOP(j) gsl_matrix_set(coxeter_fixedpoints, j, i, cox[0][i].x[j]);
 
-	initializeTriangleGenerators(gen, ctx->cartan);
+	initializeTriangleGeneratorsCurrent(gen, ctx);
 	gsl_matrix_set_identity(elements[0]);
 	for(int i = 1; i < ctx->n_group_elements; i++)
 		multiply(gen[ctx->group[i].letter], elements[ctx->group[i].parent->id], elements[i]);
@@ -1477,7 +1483,7 @@ void drawCoxeterOrbit(DrawingContext *ctx)
 	gsl_vector_set(startpoint_drawbasis, 2, 1);
 	solve(ctx->cob, startpoint_drawbasis, startpoint_globalbasis, ctx->ws);
 
-	solve(coxeter_fixedpoints, startpoint_globalbasis, startpoint_coxeterbasis, ctx->ws);
+//	solve(coxeter_fixedpoints, startpoint_globalbasis, startpoint_coxeterbasis, ctx->ws);
 //	LOOP(i) start.x[i] = gsl_vector_get(startpoint_coxeterbasis, i);
 
 	LOOP(i) start.x[i] = gsl_vector_get(startpoint_globalbasis, i);
@@ -1517,7 +1523,7 @@ void drawCoxeterOrbit(DrawingContext *ctx)
 //	LOOP(i) drawVector(ctx, abcb[i]);
 
 	cairo_restore(C);
-	releaseTempMatrices(ctx->ws, 5 + ctx->n_group_elements);
+	releaseTempMatrices(ctx->ws, 8 + ctx->n_group_elements);
 	releaseTempVectors(ctx->ws, 4);
 }
 
@@ -1526,7 +1532,7 @@ void drawText(DrawingContext *ctx)
 	cairo_move_to(ctx->cairo, 15, 30);
 	cairo_set_source_rgb(ctx->cairo, 0, 0, 0);
 	char buf[100];
-	sprintf(buf, "t = exp(%.8f) = %.8f, marking = (%.5f, %.5f)", log(ctx->parameter), ctx->parameter, ctx->marking.x, ctx->marking.y);
+	sprintf(buf, "t = exp(%.8f) = %.8f, s = exp(%.8f) = %.8f, marking = (%.5f, %.5f)", log(ctx->parameter), ctx->parameter, log(ctx->parameter2), ctx->parameter2, ctx->marking.x, ctx->marking.y);
 	cairo_show_text(ctx->cairo, buf);
 }
 

limit_set.c

@@ -5,6 +5,7 @@ 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);
@@ -20,11 +21,58 @@ 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, gsl_matrix *cartan)
+void initializeTriangleGenerators(gsl_matrix **gen, double a1, double a2, double a3, double s, double t, workspace_t *ws)
 {
-	LOOP(i) gsl_matrix_set_identity(gen[i]);
+	gsl_matrix *reflection_gen[3];
-	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);
+	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);
 }
 
 int computeLimitCurve(DrawingContext *ctx)
@@ -38,7 +86,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, 3);
+	gsl_matrix **gen = getTempMatrices(ctx->ws, 6);
 	gsl_matrix **elements = getTempMatrices(ctx->ws, ctx->n_group_elements);
 	groupelement_t *group = ctx->group;
 	int success = 0;
@@ -50,31 +98,43 @@ 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, cartan_pos);
+	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);
 	gsl_matrix_set_identity(elements[0]);
-	for(int i = 1; i < ctx->n_group_elements; i++)
+	for(int i = 1; i < ctx->n_group_elements; i++) {
-		multiply(gen[group[i].letter], elements[group[i].parent->id], 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]);
+	}
 	diagonalize_symmetric_form(cartan_pos, cob_pos, ws);
-	multiply_many(ws, coxeter_pos, 3, gen[0], gen[1], gen[2]);
+	multiply_many(ws, coxeter_pos, 3, gen[2], gen[1], gen[0]);
 	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*i+2] = atan2(
+		ctx->limit_curve[3*n+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
-	initializeTriangleGenerators(gen, ctx->cartan);
+	initializeTriangleGeneratorsCurrent(gen, ctx);
 	gsl_matrix_set_identity(elements[0]);
-	for(int i = 1; i < ctx->n_group_elements; i++)
+	for(int i = 1; i < ctx->n_group_elements; i++) {
-		multiply(gen[group[i].letter], elements[group[i].parent->id], 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]);
+	}
 
-	multiply_many(ws, coxeter, 3, gen[0], gen[1], gen[2]);
+	multiply_many(ws, coxeter, 3, gen[2], gen[1], gen[0]);
 	int ev_count = real_eigenvectors(coxeter, coxeter_fixedpoints, ws);
 
 	if(ev_count == 1)
@@ -82,11 +142,17 @@ 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*i  ] = gsl_matrix_get(fixedpoints, 0, column)/gsl_matrix_get(fixedpoints, 2, column);
+		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*i+1] = gsl_matrix_get(fixedpoints, 1, 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++;
 
 		if((x - ctx->marking.x)*(x - ctx->marking.x) + (y - ctx->marking.y)*(y - ctx->marking.y) < 25e-10)
 		{
@@ -95,19 +161,16 @@ int computeLimitCurve(DrawingContext *ctx)
 				fputc('a' + cur->letter, stdout); // bcbcbca, bacbcacab, bc bca cb
 			fputc('\n',stdout);
 		}
-
-		// bca abc acb = abc
-
 	}
 
-	qsort(ctx->limit_curve, ctx->n_group_elements, 3*sizeof(double), compareAngle);
+	qsort(ctx->limit_curve, ctx->limit_curve_count, 3*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, 11+ctx->n_group_elements);
+	releaseTempMatrices(ctx->ws, 14+ctx->n_group_elements);
 
 	return success;
 }

main.c

@@ -11,6 +11,7 @@
 #include "linalg.h"
 
 #define TOGGLE(a) do { (a) = !(a); } while(0)
+#define SIGN(x) ((x) > 0 ? 1.0 : -1.0)
 
 DrawingContext *screen_context;
 
@@ -18,6 +19,7 @@ DrawingContext *screen_context;
 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]);
@@ -28,21 +30,33 @@ void setupContext(DrawingContext *ctx, int argc, char *argv[])
 	ctx->parameter = atof(argv[7]);
     else
 	ctx->parameter = 1.0;
-//	ctx->parameter = 2.77;
+    if(argc > 8)
-//	ctx->parameter = 0.1;
+	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_limit = 0;
-	ctx->show_dual_limit= 0;
+    ctx->show_dual_limit = 0;
     ctx->show_text = 1;
     ctx->mode = 0;
-	ctx->use_rotation_basis = 0;
+    ctx->use_rotation_basis = 1;
-	ctx->limit_with_lines = 1;
+    ctx->limit_with_lines = 0;
     ctx->use_repelling = 0;
-	ctx->show_marking = 1;
+    ctx->show_marking = 0;
     ctx->marking.x = -0.73679;
     ctx->marking.y = -0.01873;
     ctx->show_coxeter_orbit = 0;
@@ -50,8 +64,8 @@ void setupContext(DrawingContext *ctx, int argc, char *argv[])
     ctx->limit_curve = malloc(3*ctx->n_group_elements*sizeof(double));
     ctx->limit_curve_count = -1;
 
-	ctx->group = malloc(ctx->n_group_elements*sizeof(groupelement_t));
+    ctx->group = malloc(ctx->n_group_elements_combinatorial*sizeof(groupelement_t));
-	generate_triangle_group(ctx->group, ctx->n_group_elements, ctx->p[0], ctx->p[1], ctx->p[2]);
+    generate_triangle_group(ctx->group, ctx->n_group_elements_combinatorial, ctx->p[0], ctx->p[1], ctx->p[2]);
 
     // the temporary stuff
     ctx->cartan = gsl_matrix_alloc(3, 3);
@@ -70,21 +84,30 @@ void destroyContext(DrawingContext *ctx)
     workspace_free(ctx->ws);
 }
 
-void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char *type)
+void computeMatrix(DrawingContext *ctx, gsl_matrix *result, const char *type)
+{
+    gsl_matrix **gen = getTempMatrices(ctx->ws, 6);
+
+    //	ERROR(strlen(type) != 2, "Invalid call of computeRotationMatrix()\n");
+
+    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);
+    }
+
+    releaseTempMatrices(ctx->ws, 6);
+}
+
+void computeRotationMatrixFrame(DrawingContext *ctx, gsl_matrix *result, const char *type)
 {
     gsl_matrix *tmp = getTempMatrix(ctx->ws);
-	gsl_matrix **gen = getTempMatrices(ctx->ws, 3);
+    computeMatrix(ctx, tmp, type);
-
-//	ERROR(strlen(type) != 2, "Invalid call of computeRotationMatrix()\n");
-
-	initializeTriangleGenerators(gen, ctx->cartan);
-	gsl_matrix_set_identity(tmp);
-	for(int i = 0; i < strlen(type); i++)
-		multiply_right(tmp, gen[type[i]-'a'], ctx->ws);
-
     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)
@@ -144,8 +167,9 @@ void updateMatrices(DrawingContext *ctx)
     cartanMatrix(ctx->cartan, angle[0], angle[1], angle[2], ctx->parameter);
 
     gsl_matrix *tmp = getTempMatrix(ctx->ws);
+    int nmodes = 5;
 
-	if(ctx->use_rotation_basis % 5 == 0) {
+    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);
@@ -156,15 +180,25 @@ void updateMatrices(DrawingContext *ctx)
 	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) {
+    } 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?
-	} else if(ctx->use_rotation_basis % 5 == 2) {
+	multiply_left(tmp, ctx->cob, ctx->ws);
-		computeRotationMatrix(ctx, tmp, "ba");
+    } else if(ctx->use_rotation_basis % nmodes == 2) {
+	computeRotationMatrixFrame(ctx, tmp, "C");
 	invert(tmp, ctx->cob, ctx->ws);
-	} else if(ctx->use_rotation_basis % 5 == 3) {
+    } else if(ctx->use_rotation_basis % nmodes == 3) {
-		computeBoxTransform(ctx, "bca", "abc", ctx->cob);
+	computeBoxTransform(ctx, "acb", "cba", ctx->cob);
-//		computeBoxTransform(ctx, "cab", "bca", ctx->cob);
+	//		computeBoxTransform(ctx, "cab", "bca", ctx->cob);
-//		computeBoxTransform(ctx, "acb", "cba", 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);
@@ -224,7 +258,7 @@ int processEvent(GraphicsInfo *info, XEvent *ev)
     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:
@@ -248,22 +282,34 @@ int processEvent(GraphicsInfo *info, XEvent *ev)
 
 	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:
-			screen_context->parameter /= exp(0.00002);
+	    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:
-			screen_context->parameter *= exp(0.00002);
+	    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;
@@ -283,7 +329,7 @@ int processEvent(GraphicsInfo *info, XEvent *ev)
 	    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);
@@ -332,37 +378,18 @@ int processEvent(GraphicsInfo *info, XEvent *ev)
 	    print(screen_context);
 	    break;
 	case 'M':
-			/*
 	    screen_context->limit_with_lines = 0;
 	    double parameter_start = screen_context->parameter;
-			for(int i = 0; i <= 1300; i++) {
+	    double parameter2_start = screen_context->parameter2;
-				if(i < 400)
+	    for(int i = 0; i <= screen_context->movie_n_frames; i++) {
-					screen_context->parameter = exp(log(parameter_start)+0.002*i);
+		screen_context->parameter = SIGN(parameter_start)*exp(log(fabs(parameter_start)) +
-				else if(i < 500)
+						i*screen_context->movie_parameter_duration/screen_context->movie_n_frames);
-					screen_context->parameter = exp(log(parameter_start)+0.002*400);
+		screen_context->parameter2 = SIGN(parameter2_start)*exp(log(fabs(parameter2_start)) +
-				else
+						 i*screen_context->movie_parameter2_duration/screen_context->movie_n_frames);
-					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);
+		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);
-			}
-			*/
-			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);
 	    }
@@ -401,14 +428,14 @@ int main(int argc, char *argv[])
     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;

main.h

@@ -11,7 +11,12 @@
 #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 50000
+#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)
 
 typedef struct {
 	double x[3];
@@ -30,7 +35,13 @@ 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;
@@ -66,7 +77,8 @@ 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, gsl_matrix *cartan);
+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);
 int computeLimitCurve(DrawingContext *ctx);
 
 // implemented in draw.c
@@ -96,7 +108,8 @@ void setupContext(DrawingContext *ctx, int argc, char *argv[]);
 void destroyContext(DrawingContext *ctx);
 void print(DrawingContext *screen);
 int processEvent(GraphicsInfo *info, XEvent *ev);
-void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char *type);
+void computeMatrix(DrawingContext *ctx, gsl_matrix *result, const char *type);
+void computeRotationMatrixFrame(DrawingContext *ctx, gsl_matrix *result, const char *type);
 void updateMatrices(DrawingContext *ctx);
 
 static vector_t vectorFromGsl(gsl_vector *v)

triangle.h

@@ -14,6 +14,7 @@ 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);