99 lines
3.4 KiB
C
99 lines
3.4 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], i, j) += gsl_matrix_get(cartan, i, j);
|
|
}
|
|
|
|
int computeLimitCurve(DrawingContext *ctx)
|
|
{
|
|
workspace_t *ws = ctx->ws;
|
|
gsl_matrix *cartan_pos = getTempMatrix(ctx);
|
|
gsl_matrix *cob_pos = getTempMatrix(ctx);
|
|
gsl_matrix *coxeter_pos = getTempMatrix(ctx);
|
|
gsl_matrix *coxeter_fixedpoints_pos = getTempMatrix(ctx);
|
|
gsl_matrix *fixedpoints_pos = getTempMatrix(ctx);
|
|
gsl_matrix *coxeter = getTempMatrix(ctx);
|
|
gsl_matrix *coxeter_fixedpoints = getTempMatrix(ctx);
|
|
gsl_matrix *fixedpoints = getTempMatrix(ctx);
|
|
gsl_matrix **gen = getTempMatrices(ctx, 3);
|
|
gsl_matrix **elements = getTempMatrices(ctx, ctx->n_group_elements);
|
|
groupelement_t *group = ctx->group;
|
|
int success = 0;
|
|
|
|
int column = ctx->use_repelling ? 2 : 0;
|
|
ctx->limit_curve_valid = 0;
|
|
|
|
// 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(elements[group[i].parent->id], gen[group[i].letter], 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, 0, column)/gsl_matrix_get(fixedpoints_pos, 2, column),
|
|
gsl_matrix_get(fixedpoints_pos, 1, column)/gsl_matrix_get(fixedpoints_pos, 2, 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(elements[group[i].parent->id], gen[group[i].letter], 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);
|
|
|
|
ctx->limit_curve[3*i ] = gsl_matrix_get(fixedpoints, 0, column)/gsl_matrix_get(fixedpoints, 2, column);
|
|
ctx->limit_curve[3*i+1] = gsl_matrix_get(fixedpoints, 1, column)/gsl_matrix_get(fixedpoints, 2, column);
|
|
}
|
|
|
|
qsort(ctx->limit_curve, ctx->n_group_elements, 3*sizeof(double), compareAngle);
|
|
|
|
ctx->limit_curve_valid = 1;
|
|
|
|
success = 1;
|
|
|
|
error_out:
|
|
releaseTempMatrices(ctx, 11+ctx->n_group_elements);
|
|
|
|
return success;
|
|
}
|