start drawing rotation orbits

draw.c

@@ -2,6 +2,12 @@
 
 // level 0: helper functions
 
+int isInsideBB(DrawingContext *ctx, point_t p)
+{
+	cairo_user_to_device(ctx->cairo, &x, &y);
+	return -p.x < ctx->dim->width && p.x < 3*ctx->dim->width && -p.y < ctx->dim->height && p.y < 3*ctx->dim->height;
+}
+
 vector_t cross(vector_t a, vector_t b)
 {
 	vector_t result;
@@ -47,26 +53,6 @@ int fixedPoints(DrawingContext *ctx, const char *word, vector_t *out)
 	return count;
 }
 
-void transformFrameStd(DrawingContext *ctx, vector_t *x, gsl_matrix *out)
-{
-	gsl_matrix *tmp = getTempMatrix(ctx->ws);
-	gsl_vector *fourth = getTempVector(ctx->ws);
-	gsl_vector *lambda = getTempVector(ctx->ws);
-	int s;
-
-	LOOP(i) LOOP(j) gsl_matrix_set(out, j, i, x[i].x[j]);
-	gsl_matrix_memcpy(tmp, out);
-	gsl_linalg_LU_decomp(tmp, ctx->ws->permutation, &s);
-	gsl_linalg_LU_solve(tmp, ctx->ws->permutation, fourth, lambda);
-
-	LOOP(i) LOOP(j) *gsl_matrix_ptr(out, i, j) *= gsl_vector_get(lambda, j);
-
-	gsl_matrix_fprintf(stdout, out, "%f");
-
-	releaseTempMatrices(ctx->ws, 1);
-	releaseTempVectors(ctx->ws, 2);
-}
-
 // level 1: the elementary drawing functions, drawPoint, drawSegment2d
 
 void drawPoint(DrawingContext *ctx, point_t p)
@@ -257,8 +243,8 @@ void drawAttractors(DrawingContext *ctx)
 	vector_t l[3][3];
 
 	fixedPoints(ctx, "abc", p[0]);
-	fixedPoints(ctx, "bca", p[1]);
+	fixedPoints(ctx, "bca", p[1]); // (bca, cab) -> (cab, cacabac) -> (cacabac, acacbacac) ->
-	fixedPoints(ctx, "cab", p[2]);
+	fixedPoints(ctx, "cab", p[2]); // ac cab ca = bca, ac bac ca = acb
 
 	double color[3][3] = {{1,0,0},{0,0.7,0},{0,0,1}};
 
@@ -341,71 +327,55 @@ void drawBoxes(DrawingContext *ctx)
 
 void drawBoxes2(DrawingContext *ctx)
 {
-	/*
+	gsl_matrix *tmp = getTempMatrix(ctx->ws);
-	cairo_set_source_rgb(ctx->cairo, 1, 0.5, 0);
-	drawBoxStd(ctx, "a", 'A');
-	drawBoxStd(ctx, "ca", 'A');
-	drawBoxStd(ctx, "aca", 'A');
-	drawBoxStd(ctx, "caca", 'A');
-	drawBoxStd(ctx, "acaca", 'A');
-	cairo_set_source_rgb(ctx->cairo, 0, 0.5, 1);
-	drawBoxStd(ctx, "acac", 'A');
-	drawBoxStd(ctx, "cac", 'A');
-	drawBoxStd(ctx, "ac", 'A');
-	drawBoxStd(ctx, "c", 'A');
-	drawBoxStd(ctx, "", 'A');
-	*/
-
-	/*
-	cairo_set_source_rgb(ctx->cairo, 0, 0.5, 1);
-	drawBox(ctx, "abc", "cab");
-//	drawBox(ctx, "ba abc ab", "ba cab ab");
-	drawBox(ctx, "baba abc abab", "baba cab abab");
-	drawBox(ctx, "abab abc baba", "abab cab baba");
-	drawBox(ctx, "ab abc ba", "ab cab ba");
 
 	cairo_set_source_rgb(ctx->cairo, 1, 0.5, 0);
-	drawBox(ctx, "b abc b", "b cab b");
+	drawBox(ctx, "ac cab ca", "ac bca ca");
-	drawBox(ctx, "bab abc bab", "bab cab bab");
+	drawBox(ctx, "acac cab caca", "acac bca caca");
-	drawBox(ctx, "babab abc babab", "babab cab babab");
+	drawBox(ctx, "caca cab acac", "caca bca acac");
-	drawBox(ctx, "aba abc aba", "aba cab aba");
+	drawBox(ctx, "ca cab ac", "ca bca ac"); // cacabac, cab
-//	drawBox(ctx, "a abc a", "a cab a");
 
-	cairo_set_source_rgb(ctx->cairo, 0, 0.5, 1);
+	cairo_set_source_rgb(ctx->cairo, 1, 0, 0.5);
-	drawBox(ctx, "bca", "abc");
-//	drawBox(ctx, "cb bca bc", "cb abc bc");
-	drawBox(ctx, "cbcb bca bcbc", "cbcb abc bcbc");
-	drawBox(ctx, "bcbc bca cbcb", "bcbc abc cbcb");
-	drawBox(ctx, "bc bca cb", "bc abc cb");
-
-	cairo_set_source_rgb(ctx->cairo, 1, 0.5, 0);
-	drawBox(ctx, "c bca c", "c abc c");
-	drawBox(ctx, "cbc bca cbc", "cbc abc cbc");
-	drawBox(ctx, "cbcbc bca cbcbc", "cbcbc abc cbcbc");
-	drawBox(ctx, "bcb bca bcb", "bcb abc bcb");
-//	drawBox(ctx, "b bca b", "b abc b");
-*/
-
-//	cairo_set_source_rgb(ctx->cairo, 1, 0, 0);
-//	drawBoxLines(ctx, "bca", "abc");
-
-	cairo_set_source_rgb(ctx->cairo, 1, 0.5, 0);
-	drawBoxLines(ctx, "ac cab ca", "ac bca ca");
-	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");
+	drawBox(ctx, "abc", "cab");
 
-	cairo_set_source_rgb(ctx->cairo, 1, 0.5, 0);
+	vector_t w = {0, 1, 0};
-//	drawBox(ctx, "a cab a", "a bca a");
+	drawCovector(ctx, w);
-//	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");
 
+	vector_t v[3];
+	point_t p;
+	cairo_set_source_rgb(ctx->cairo, 1, 0, 0);
+	computeRotationMatrix(ctx, tmp, "ac");
+	LOOP(i) LOOP(j) v[i].x[j] = gsl_matrix_get(tmp, j, i);
+
+	for(int k = 0; k < 50; k++) {
+		LOOP(i) w.x[i] = 1.8 * v[2].x[i] + cos(2*k*M_PI/50) * v[0].x[i] + sin(2*k*M_PI/50) * v[1].x[i];
+		p = vectorToPoint(ctx, w);
+
+		if(k == 0)
+			cairo_move_to(ctx->cairo, p.x, p.y);
+		else
+			cairo_line_to(ctx->cairo, p.x, p.y);
+	}
+
+	cairo_set_source_rgb(ctx->cairo, 1, 0, 0);
+	cairo_stroke(ctx->cairo);
+
+	/* vector_t p[3]; */
+	/* fixedPoints(ctx, "cab", p); */
+	/* drawVector(ctx, p[0]); */
+	/* fixedPoints(ctx, "ca cab ac", p); */
+	/* drawVector(ctx, p[0]); */
+	/* fixedPoints(ctx, "caca cab acac", p); */
+	/* drawVector(ctx, p[0]); */
+	/* fixedPoints(ctx, "acac cab caca", p); */
+	/* drawVector(ctx, p[0]); */
+	/* fixedPoints(ctx, "ac cab ca", p); */
+	/* drawVector(ctx, p[0]); */
+
+	releaseTempMatrices(ctx->ws, 1);
 }
 
 void drawLimitCurve(DrawingContext *ctx)

linalg.c

@@ -4,7 +4,9 @@
 #include <gsl/gsl_eigen.h>
 #include <gsl/gsl_blas.h>
 #include <gsl/gsl_linalg.h>
+#include <gsl/gsl_complex_math.h>
 #include <memory.h>
+#include <math.h>
 
 #include "linalg.h"
 
@@ -293,3 +295,65 @@ int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws)
 	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)
+{
+	gsl_matrix *tmp = getTempMatrix(ws);
+	gsl_vector *coeff = getTempVector(ws);
+	int s;
+	double det, scale;
+
+	for(int i = 0; i < 3; i++)
+		for(int j = 0; j < 3; j++)
+			gsl_matrix_set(tmp, i, j, gsl_vector_get(vertices[j], i));
+
+	gsl_linalg_LU_decomp(tmp, ws->permutation, &s);
+	gsl_linalg_LU_solve(tmp, ws->permutation, vertices[3], coeff);
+	det = gsl_linalg_LU_det(tmp, s);
+
+	for(int i = 0; i < 3; i++)
+		det *= gsl_vector_get(coeff, i);
+	scale = 1/cbrt(det);
+
+	for(int i = 0; i < 3; i++)
+		for(int j = 0; j < 3; j++)
+			gsl_matrix_set(result, i, j, scale*gsl_vector_get(vertices[j], i)*gsl_vector_get(coeff, j));
+
+	releaseTempMatrices(ws, 1);
+	releaseTempVectors(ws, 1);
+}
+
+void rotation_frame(gsl_matrix *rotation, gsl_matrix *result, workspace_t *ws)
+{
+	gsl_matrix *tmp = getTempMatrix(ws);
+	gsl_matrix *rot_basis = getTempMatrix(ws);
+
+	gsl_matrix_memcpy(tmp, rotation);
+	gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv);
+	int r = gsl_eigen_nonsymmv(tmp, ws->eval_complex, ws->evec_complex, ws->work_nonsymmv);
+	ERROR(r, "gsl_eigen_nonsymmv failed!\n");
+
+	double arg, minarg = 5; // greater than pi
+	int minidx;
+	for(int i = 0; i < 3; i++) {
+		arg = gsl_complex_arg(gsl_vector_complex_get(ws->eval_complex, i));
+		if(abs(arg) < minarg)
+		{
+			minidx = i;
+			minarg = abs(arg);
+		}
+	}
+	ERROR(FCMP(minarg, 0.0) != 0, "rotation_frame() failed! No eigenvalue was 1.\n");
+
+	for(int i = 0; i < 3; i++) {
+		gsl_complex x = gsl_matrix_complex_get(ws->evec_complex, i, (minidx+1)%3);
+		gsl_complex y = gsl_matrix_complex_get(ws->evec_complex, i, (minidx+2)%3);
+		gsl_complex z = gsl_matrix_complex_get(ws->evec_complex, i, minidx);
+		gsl_matrix_set(result, i, 0, GSL_REAL(x)+GSL_REAL(y));
+		gsl_matrix_set(result, i, 1, GSL_IMAG(x)-GSL_IMAG(y));
+		gsl_matrix_set(result, i, 2, GSL_REAL(z));
+	}
+
+	releaseTempMatrices(ws, 2);
+}

linalg.h

@@ -51,7 +51,8 @@ int eigenvectors(gsl_matrix *g, gsl_matrix *evec, workspace_t *ws);
 int real_eigenvectors(gsl_matrix *g, gsl_matrix *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);
-
+void projective_frame(gsl_vector **vertices, gsl_matrix *result, workspace_t *ws);
+void rotation_frame(gsl_matrix *rotation, gsl_matrix *result, workspace_t *ws);
 
 // matrix allocation stuff
 

main.c

@@ -60,7 +60,6 @@ void destroyContext(DrawingContext *ctx)
 void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char *type)
 {
 	gsl_matrix *tmp = getTempMatrix(ctx->ws);
-	gsl_matrix *rot_basis = getTempMatrix(ctx->ws);
 	gsl_matrix **gen = getTempMatrices(ctx->ws, 3);
 
 	ERROR(strlen(type) != 2, "Invalid call of computeRotationMatrix()\n");
@@ -68,24 +67,61 @@ void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char *
 	initializeTriangleGenerators(gen, ctx->cartan);
 	gsl_matrix_set_identity(tmp);
 	multiply_right(tmp, gen[type[0]-'a'], ctx->ws);
-	multiply_right(tmp, gen[type[1]-'b'], ctx->ws);
+	multiply_right(tmp, gen[type[1]-'a'], ctx->ws);
 
-	gsl_eigen_nonsymmv_params(0, ctx->ws->work_nonsymmv);
+	rotation_frame(tmp, result, ctx->ws);
-	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) {
+	releaseTempMatrices(ctx->ws, 4);
-		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, result, ctx->ws);
+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}
+	};
 
-	releaseTempMatrices(ctx->ws, 5);
+	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]);
+
+	// 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]
+
+	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]);
+
+	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");
+		}*/
+
+	releaseTempVectors(ctx->ws, 8);
+	releaseTempMatrices(ctx->ws, 3);
 }
 
 void updateMatrices(DrawingContext *ctx)
@@ -94,12 +130,15 @@ void updateMatrices(DrawingContext *ctx)
 	LOOP(i) angle[i] = M_PI*ctx->k[i]/ctx->p[i];
 	cartanMatrix(ctx->cartan, angle[0], angle[1], angle[2], ctx->parameter);
 
-	if(ctx->use_rotation_basis)
+	if(ctx->use_rotation_basis) {
-		computeRotationMatrix(ctx, ctx->cob, "ac");
+//		computeRotationMatrix(ctx, ctx->cob, "ac"); // and don't forget to invert!
-	else
+		computeBoxTransform(ctx, "abc", "cab", ctx->cob);
+	} else {
 		gsl_matrix_memcpy(ctx->cob, ctx->cartan); // is this a good choice of basis for any reason?
 	}
 
+}
+
 void print(DrawingContext *screen)
 {
 	DrawingContext file;
@@ -241,12 +280,22 @@ int main()
 	if(!info)
 		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;
+
 	updateDimensions(info->dim);
 
 	screen_context->dim = info->dim;

main.h

@@ -1,6 +1,8 @@
 #ifndef TRIANGLE_GROUP_MAIN_H
 #define TRIANGLE_GROUP_MAIN_H
 
+#include <gsl/gsl_linalg.h>
+
 #include "triangle.h"
 #include "linalg.h"
 #include "initcairo.h"
@@ -84,4 +86,16 @@ int processEvent(GraphicsInfo *info, XEvent *ev);
 void computeRotationMatrix(DrawingContext *ctx, gsl_matrix *result, const char *type);
 void updateMatrices(DrawingContext *ctx);
 
+static vector_t vectorFromGsl(gsl_vector *v)
+{
+	vector_t result;
+	LOOP(i) result.x[i] = gsl_vector_get(v, i);
+	return result;
+}
+
+static void vectorToGsl(vector_t v, gsl_vector *out)
+{
+	LOOP(i) gsl_vector_set(out, i, v.x[i]);
+}
+
 #endif