177 lines
5.8 KiB
C
177 lines
5.8 KiB
C
#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;
|
|
}
|