transform into rotation basis

draw.c

@@ -151,7 +151,7 @@ void drawSegment(DrawingContext *ctx, vector_t a, vector_t b)
 
 // level 3: boxes and polygons
 
-void drawPolygon(DrawingContext *ctx, int sides, ...)
+void drawPolygon(DrawingContext *ctx, int segments, int sides, ...)
 {
 	va_list args;
 	vector_t first, prev, current;
@@ -163,9 +163,15 @@ void drawPolygon(DrawingContext *ctx, int sides, ...)
 	for(int i = 0; i < sides-1; i++) {
 		prev = current;
 		current = va_arg(args, vector_t);
+		if(segments)
 			drawSegment(ctx, prev, current);
+		else
+			drawCovector(ctx, cross(prev, current));
 	}
+	if(segments)
 		drawSegment(ctx, current, first);
+	else
+		drawCovector(ctx, cross(current, first));
 
 	va_end(args);
 }
@@ -175,7 +181,7 @@ void drawTriangle(DrawingContext *ctx, const char *word)
 	vector_t p[3];
 
 	fixedPoints(ctx, word, p);
-	drawPolygon(ctx, 3, p[0], p[1], p[2]);
+	drawPolygon(ctx, 1, 3, p[0], p[1], p[2]);
 }
 
 void drawBox(DrawingContext *ctx, const char *word1, const char *word2)
@@ -189,7 +195,21 @@ void drawBox(DrawingContext *ctx, const char *word1, const char *word2)
 	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]));
 
-	drawPolygon(ctx, 4, p[0][0], i[0], p[1][0], i[1]);
+	drawPolygon(ctx, 1, 4, p[0][0], i[0], p[1][0], i[1]);
+}
+
+void drawBoxLines(DrawingContext *ctx, const char *word1, const char *word2)
+{
+	vector_t p[2][3],i[2];
+
+	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]));
+
+	drawPolygon(ctx, 0, 4, p[0][0], i[0], p[1][0], i[1]);
 }
 
 void drawBoxStd(DrawingContext *ctx, const char *word, char base)
@@ -366,19 +386,22 @@ void drawBoxes2(DrawingContext *ctx)
 //	drawBox(ctx, "b bca b", "b abc b");
 */
 
-	cairo_set_source_rgb(ctx->cairo, 0, 0, 0);
+//	cairo_set_source_rgb(ctx->cairo, 1, 0, 0);
-	drawBox(ctx, "abc", "cab");
+//	drawBoxLines(ctx, "bca", "abc");
-
-	cairo_set_source_rgb(ctx->cairo, 0, 0.5, 1);
-//	drawBox(ctx, "cab", "bca");
-//	drawBox(ctx, "ac cab ca", "ac bca ca");
-//	drawBox(ctx, "acac cab caca", "acac bca caca");
-	drawBox(ctx, "caca cab acac", "caca bca acac");
-	drawBox(ctx, "ca cab ac", "ca bca ac");
 
 	cairo_set_source_rgb(ctx->cairo, 1, 0.5, 0);
-	drawBox(ctx, "a cab a", "a bca a");
+	drawBoxLines(ctx, "ac cab ca", "ac bca ca");
-	drawBox(ctx, "aca cab aca", "aca bca aca");
+	drawBoxLines(ctx, "acac cab caca", "acac bca caca");
+	drawBoxLines(ctx, "caca cab acac", "caca bca acac");
+	drawBoxLines(ctx, "ca cab ac", "ca bca ac");
+	drawBoxLines(ctx, "cab", "bca");
+
+	cairo_set_source_rgb(ctx->cairo, 0, 0, 0);
+	drawBoxLines(ctx, "abc", "cab");
+
+	cairo_set_source_rgb(ctx->cairo, 1, 0.5, 0);
+//	drawBox(ctx, "a cab a", "a bca a");
+//	drawBox(ctx, "aca cab aca", "aca bca aca");
 //	drawBox(ctx, "acaca cab acaca", "acaca bca acaca");
 //	drawBox(ctx, "cac cab cac", "cac bca cac");
 //	drawBox(ctx, "c cab c", "c bca c");

main.c

@@ -30,6 +30,7 @@ void setupContext(DrawingContext *ctx)
 	ctx->show_attractors = 0;
 	ctx->show_reflectors = 0;
 	ctx->show_limit= 1;
+	ctx->use_rotation_basis = 0;
 	ctx->limit_with_lines = 1;
 	ctx->use_repelling = 0;
 
@@ -61,7 +62,36 @@ 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);
+
+	if(!ctx->use_rotation_basis)
 		gsl_matrix_memcpy(ctx->cob, ctx->cartan); // is this a good choice of basis
+	else {
+		gsl_matrix *tmp = getTempMatrix(ctx->ws);
+		gsl_matrix *rot_basis = getTempMatrix(ctx->ws);
+		gsl_matrix **gen = getTempMatrices(ctx->ws, 3);
+
+		initializeTriangleGenerators(gen, ctx->cartan);
+		gsl_matrix_set_identity(tmp);
+		multiply_right(tmp, gen[0], ctx->ws);
+		multiply_right(tmp, gen[2], ctx->ws);
+
+		gsl_eigen_nonsymmv_params(0, ctx->ws->work_nonsymmv);
+		int r = gsl_eigen_nonsymmv(tmp, ctx->ws->eval_complex, ctx->ws->evec_complex, ctx->ws->work_nonsymmv);
+		ERROR(r, "gsl_eigen_nonsymmv failed!\n");
+
+		LOOP(i) {
+			gsl_complex x = gsl_matrix_complex_get(ctx->ws->evec_complex, i, 0);
+			gsl_complex y = gsl_matrix_complex_get(ctx->ws->evec_complex, i, 1);
+			gsl_complex z = gsl_matrix_complex_get(ctx->ws->evec_complex, i, 2);
+			gsl_matrix_set(tmp, i, 0, GSL_REAL(x)+GSL_REAL(y));
+			gsl_matrix_set(tmp, i, 1, GSL_IMAG(x)-GSL_IMAG(y));
+			gsl_matrix_set(tmp, i, 2, GSL_REAL(z));
+		}
+
+		invert(tmp, ctx->cob, ctx->ws);
+
+		releaseTempMatrices(ctx->ws, 5);
+	}
 }
 
 void print(DrawingContext *screen)
@@ -173,6 +203,11 @@ int processEvent(GraphicsInfo *info, XEvent *ev)
 		case 'l':
 			TOGGLE(screen_context->show_limit);
 			break;
+		case 'R':
+			TOGGLE(screen_context->use_rotation_basis);
+			updateMatrices(screen_context);
+			computeLimitCurve(screen_context);
+			break;
 		case 'p':
 			print(screen_context);
 			break;

main.h

@@ -34,6 +34,7 @@ typedef struct {
 	int show_attractors;
 	int show_reflectors;
 	int show_limit;
+	int use_rotation_basis;
 	int limit_with_lines;
 	int use_repelling;
 	gsl_matrix *cartan, *cob;
@@ -62,7 +63,7 @@ void drawSegment2d(DrawingContext *ctx, point_t a, point_t b);
 void drawVector(DrawingContext *ctx, vector_t v);
 void drawCovector(DrawingContext *ctx, vector_t v);
 void drawSegment(DrawingContext *ctx, vector_t a, vector_t b);
-void drawPolygon(DrawingContext *ctx, int sides, ...);
+void drawPolygon(DrawingContext *ctx, int segments, int sides, ...);
 void drawTriangle(DrawingContext *ctx, const char *word);
 void drawBox(DrawingContext *ctx, const char *word1, const char *word2);
 void drawBoxStd(DrawingContext *ctx, const char *word, char base);