114 lines
3.9 KiB
C
114 lines
3.9 KiB
C
#include "main.h"
|
|
|
|
static int compareAngle(const void *x, const void *y)
|
|
{
|
|
return ((double*)x)[2] > ((double*)y)[2] ? 1 : -1;
|
|
}
|
|
|
|
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, gsl_matrix *cartan)
|
|
{
|
|
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)
|
|
{
|
|
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, 3);
|
|
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, cartan_pos);
|
|
gsl_matrix_set_identity(elements[0]);
|
|
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[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;
|
|
|
|
for(int i = 0; i < ctx->n_group_elements; i++) {
|
|
multiply_many(ws, fixedpoints_pos, 3, cob_pos, elements[i], coxeter_fixedpoints_pos);
|
|
ctx->limit_curve[3*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));
|
|
}
|
|
|
|
// now do it again to calculate x and y coordinates
|
|
initializeTriangleGenerators(gen, ctx->cartan);
|
|
gsl_matrix_set_identity(elements[0]);
|
|
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[0], gen[1], gen[2]);
|
|
int ev_count = real_eigenvectors(coxeter, coxeter_fixedpoints, ws);
|
|
|
|
if(ev_count == 1)
|
|
column = 0;
|
|
if(ev_count == 0)
|
|
goto error_out;
|
|
|
|
for(int i = 0; i < ctx->n_group_elements; i++) {
|
|
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);
|
|
y = ctx->limit_curve[3*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)
|
|
{
|
|
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);
|
|
}
|
|
|
|
// bca abc acb = abc
|
|
|
|
}
|
|
|
|
qsort(ctx->limit_curve, ctx->n_group_elements, 3*sizeof(double), compareAngle);
|
|
|
|
ctx->limit_curve_count = ctx->n_group_elements;
|
|
|
|
success = 1;
|
|
|
|
error_out:
|
|
releaseTempMatrices(ctx->ws, 11+ctx->n_group_elements);
|
|
|
|
return success;
|
|
}
|