#include "main.h"

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);
	gsl_matrix_set(cartan, 0, 1, -2*s*cos(a3));
	gsl_matrix_set(cartan, 0, 2, -2/s*cos(a2));

	gsl_matrix_set(cartan, 1, 0, -2/s*cos(a3));
	gsl_matrix_set(cartan, 1, 1, 2);
	gsl_matrix_set(cartan, 1, 2, -2*s*cos(a1));

	gsl_matrix_set(cartan, 2, 0, -2*s*cos(a2));
	gsl_matrix_set(cartan, 2, 1, -2/s*cos(a1));
	gsl_matrix_set(cartan, 2, 2, 2);
}

void initializeTriangleGenerators(gsl_matrix **gen, double a1, double a2, double a3, double s, double t, workspace_t *ws)
{
	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);
}

int computeLimitCurve(DrawingContext *ctx)
{
	workspace_t *ws = ctx->ws;
	gsl_matrix *cartan_pos = getTempMatrix(ctx->ws);
	gsl_matrix *cob_pos = getTempMatrix(ctx->ws);
	gsl_matrix *coxeter_pos = getTempMatrix(ctx->ws);
	gsl_matrix *coxeter_fixedpoints_pos = getTempMatrix(ctx->ws);
	gsl_matrix *fixedpoints_pos = getTempMatrix(ctx->ws);
	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 **elements = getTempMatrices(ctx->ws, ctx->n_group_elements);
	groupelement_t *group = ctx->group;
	int success = 0;

	int column = ctx->use_repelling ? 2 : 0;
	double x,y;
//	int column = 1;
	ctx->limit_curve_count = -1;

	// do first in the Fuchsian positive case to get the angles
	cartanMatrix(cartan_pos, M_PI/ctx->p[0], M_PI/ctx->p[1], M_PI/ctx->p[2], 1.0);
	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++) {
		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[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*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
	initializeTriangleGeneratorsCurrent(gen, ctx);
	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]);
	}

	multiply_many(ws, coxeter, 3, gen[2], gen[1], gen[0]);
	int ev_count = real_eigenvectors(coxeter, coxeter_fixedpoints, ws);

	if(ev_count == 1)
		column = 0;
	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++;

		if((x - ctx->marking.x)*(x - ctx->marking.x) + (y - ctx->marking.y)*(y - ctx->marking.y) < 25e-10)
		{
			printf("limit point %d is close: length %d, ", i, group[i].length);
			for(groupelement_t *cur = &group[i]; cur->parent; cur = cur->parent)
				fputc('a' + cur->letter, stdout); // bcbcbca, bacbcacab, bc bca cb
			fputc('\n',stdout);
		}
	}

	qsort(ctx->limit_curve, ctx->limit_curve_count, 3*sizeof(double), compareAngle);

//	ctx->limit_curve_count = ctx->n_group_elements;

	success = 1;

error_out:
	releaseTempMatrices(ctx->ws, 14+ctx->n_group_elements);

	return success;
}