redo hyperbolic pictures
This commit is contained in:
Florian Stecker
199
hyperbolic.c
199
hyperbolic.c
@@ -198,6 +198,16 @@ void draw_triangle(point *p, gsl_matrix *frame, const char *arguments)
|
||||
#endif
|
||||
}
|
||||
|
||||
void draw_dot(point p, gsl_matrix *frame, const char *arguments)
|
||||
{
|
||||
double x = coord(p, 0, frame);
|
||||
double y = coord(p, 1, frame);
|
||||
|
||||
#ifdef POINCARE
|
||||
printf("<circle cx=\"%f\" cy=\"%f\" r=\"1\" style=\"%s\"/>\n", CONV(x), CONV(y), arguments);
|
||||
#endif
|
||||
}
|
||||
|
||||
void draw_line(point p1, point p2, gsl_matrix *frame, const char *arguments)
|
||||
{
|
||||
char buffer[100];
|
||||
@@ -217,11 +227,16 @@ void draw_line(point p1, point p2, gsl_matrix *frame, const char *arguments)
|
||||
#endif
|
||||
}
|
||||
|
||||
void compute_word(workspace_t *ws, gsl_matrix *result, gsl_matrix **gen, const char *word, int modifier)
|
||||
void compute_word(workspace_t *ws, gsl_matrix *result, gsl_matrix **gen, const char *word, int modifier, int inverse)
|
||||
{
|
||||
gsl_matrix_set_identity(result);
|
||||
if(inverse) {
|
||||
for(int i = 0; word[i] != 0; i++)
|
||||
multiply_left(gen[(word[i]-'a'+modifier)%3], result, ws);
|
||||
} else {
|
||||
for(int i = 0; word[i] != 0; i++)
|
||||
multiply_right(result, gen[(word[i]-'a'+modifier)%3], ws);
|
||||
}
|
||||
}
|
||||
|
||||
int main(int argc, const char *argv[])
|
||||
@@ -230,20 +245,22 @@ int main(int argc, const char *argv[])
|
||||
gsl_matrix **matrices;
|
||||
gsl_matrix *cartan;
|
||||
gsl_matrix *gen[3];
|
||||
gsl_matrix *coxeter[3];
|
||||
gsl_matrix *coxeter2[3];
|
||||
gsl_matrix *coxeter_eigenvectors[3];
|
||||
gsl_matrix *coxeter_eigenvectors2[3];
|
||||
gsl_matrix **special;
|
||||
gsl_matrix **special_eigenvectors;
|
||||
point *special_attracting, *special_repelling;
|
||||
gsl_matrix *frame;
|
||||
workspace_t *ws;
|
||||
int p,q,r;
|
||||
int elements, nspecial;
|
||||
|
||||
if(argc < 5) {
|
||||
fprintf(stderr, "Usage: %s <p> <q> <r> <n_elements>\n", argv[0]);
|
||||
fprintf(stderr, "Usage: %s <p> <q> <r> <n_elements> <word1> <word2> ...\n", argv[0]);
|
||||
exit(1);
|
||||
}
|
||||
|
||||
int elements = atoi(argv[4]);
|
||||
int p = atoi(argv[1]), q = atoi(argv[2]), r = atoi(argv[3]);
|
||||
nspecial = argc - 5;
|
||||
elements = atoi(argv[4]);
|
||||
p = atoi(argv[1]), q = atoi(argv[2]), r = atoi(argv[3]);
|
||||
|
||||
group = malloc(elements*sizeof(groupelement_t));
|
||||
matrices = malloc(elements*sizeof(gsl_matrix*));
|
||||
@@ -252,79 +269,44 @@ int main(int argc, const char *argv[])
|
||||
cartan = gsl_matrix_alloc(3, 3);
|
||||
frame = gsl_matrix_alloc(3, 3);
|
||||
LOOP(i) gen[i] = gsl_matrix_alloc(3, 3);
|
||||
LOOP(i) coxeter[i] = gsl_matrix_alloc(3, 3);
|
||||
LOOP(i) coxeter_eigenvectors[i] = gsl_matrix_alloc(3, 3);
|
||||
LOOP(i) coxeter2[i] = gsl_matrix_alloc(3, 3);
|
||||
LOOP(i) coxeter_eigenvectors2[i] = gsl_matrix_alloc(3, 3);
|
||||
ws = workspace_alloc(3);
|
||||
special = malloc(3*nspecial*sizeof(gsl_matrix*));
|
||||
special_eigenvectors = malloc(3*nspecial*sizeof(gsl_matrix*));
|
||||
special_attracting = malloc(3*nspecial*sizeof(point));
|
||||
special_repelling = malloc(3*nspecial*sizeof(point));
|
||||
for(int i = 0; i < 3*nspecial; i++) {
|
||||
special[i] = gsl_matrix_alloc(3, 3);
|
||||
special_eigenvectors[i] = gsl_matrix_alloc(3, 3);
|
||||
}
|
||||
|
||||
generate_triangle_group(group, elements, p, q, r);
|
||||
cartan_matrix(cartan, M_PI/p, M_PI/q, M_PI/r, 1.0);
|
||||
initialize_triangle_generators(gen, cartan);
|
||||
int pos = diagonalize_symmetric_form(cartan, frame, ws); // choose frame of reference which diagonalizes the form
|
||||
diagonalize_symmetric_form(cartan, frame, ws); // choose frame of reference which diagonalizes the form
|
||||
|
||||
gsl_matrix_set_identity(matrices[0]);
|
||||
for(int i = 1; i < elements; i++)
|
||||
multiply(matrices[group[i].parent->id], gen[group[i].letter], matrices[i]);
|
||||
|
||||
// LOOP(i) multiply_many(ws, coxeter[i], 3, gen[i%3], gen[(i+1)%3], gen[(i+2)%3]); // coxeter
|
||||
// LOOP(i) multiply_many(ws, coxeter[i], 4, gen[i%3], gen[(i+1)%3], gen[i%3], gen[(i+2)%3]); // abcb
|
||||
LOOP(i) compute_word(ws, coxeter[i], gen, "abcb", i);
|
||||
/* LOOP(i) multiply_many(ws, coxeter[i], 10,
|
||||
gen[i%3], gen[(i+1)%3], gen[(i+2)%3],
|
||||
gen[i%3], gen[(i+1)%3], gen[(i+2)%3],
|
||||
gen[i%3], gen[(i+1)%3], gen[(i+2)%3],
|
||||
gen[(i+1)%3]); // (abc)^3 b */
|
||||
LOOP(i) eigenvectors(coxeter[i], coxeter_eigenvectors[i], ws);
|
||||
LOOP(i) eigenvectors(coxeter2[i], coxeter_eigenvectors2[i], ws);
|
||||
for(int i = 0; i < nspecial; i++)
|
||||
LOOP(j) {
|
||||
compute_word(ws, special[3*i+j], gen, argv[i+5], j, 0);
|
||||
eigenvectors(special[3*i+j], special_eigenvectors[3*i+j], ws);
|
||||
special_attracting[3*i+j] = column(special_eigenvectors[3*i+j], 0);
|
||||
|
||||
/*
|
||||
for(int i = 0; i < elements; i++) {
|
||||
printf("%4d: ", i);
|
||||
for(groupelement_t *cur = &group[i]; cur->parent; cur = cur->parent)
|
||||
fputc(cur->letter+'a', stdout);
|
||||
fputc('\n', stdout);
|
||||
// repelling = attracting of inverse
|
||||
compute_word(ws, special[3*i+j], gen, argv[i+5], j, 1);
|
||||
eigenvectors(special[3*i+j], special_eigenvectors[3*i+j], ws);
|
||||
special_repelling[3*i+j] = column(special_eigenvectors[3*i+j], 0);
|
||||
}
|
||||
|
||||
return 0;*/
|
||||
|
||||
point coxeter_attracting[3];
|
||||
point coxeter_repelling[3];
|
||||
point coxeter_axes[3];
|
||||
point coxeter2_attracting[3];
|
||||
point coxeter2_repelling[3];
|
||||
point coxeter2_axes[3];
|
||||
point edge_midpoints[3];
|
||||
point reflection_lines[3];
|
||||
point triangle_points[3];
|
||||
point transformed[3];
|
||||
point transformed2[3];
|
||||
point transformed3[3];
|
||||
point transformed4[3];
|
||||
point transformed5[3];
|
||||
point transformed6[3];
|
||||
point center;
|
||||
|
||||
LOOP(i) coxeter_attracting[i] = column(coxeter_eigenvectors[i], 0);
|
||||
LOOP(i) coxeter_repelling[i] = column(coxeter_eigenvectors[i], 2);
|
||||
LOOP(i) coxeter_axes[i] = incidence(coxeter_attracting[i], coxeter_repelling[i]);
|
||||
LOOP(i) edge_midpoints[i] = incidence(coxeter_axes[(i+1)%3], coxeter_axes[(i+2)%3]);
|
||||
LOOP(i) reflection_lines[i] = row(cartan, i);
|
||||
LOOP(i) triangle_points[i] = incidence(reflection_lines[(i+1)%3], reflection_lines[(i+2)%3]);
|
||||
LOOP(i) coxeter2_attracting[i] = column(coxeter_eigenvectors2[i], 0);
|
||||
LOOP(i) coxeter2_repelling[i] = column(coxeter_eigenvectors2[i], 2);
|
||||
LOOP(i) coxeter2_axes[i] = incidence(coxeter2_attracting[i], coxeter2_repelling[i]);
|
||||
|
||||
print_svg_header();
|
||||
|
||||
/*
|
||||
// let's correct the frame of reference by using hyperbolic transformations
|
||||
center = apply(frame, triangle_points[2]);
|
||||
double angle = atan2(center.x[1], center.x[0]);
|
||||
double boost = atanh(-sqrt(center.x[0]*center.x[0]+center.x[1]*center.x[1])/center.x[2]);
|
||||
gsl_matrix *frame_correction = gsl_matrix_alloc(3, 3);
|
||||
gsl_matrix_set_identity(frame_correction);
|
||||
|
||||
/*
|
||||
gsl_matrix_set(frame_correction, 0, 0, cos(angle-M_PI/2));
|
||||
gsl_matrix_set(frame_correction, 0, 1, sin(angle-M_PI/2));
|
||||
gsl_matrix_set(frame_correction, 1, 0, -sin(angle-M_PI/2));
|
||||
@@ -338,80 +320,57 @@ int main(int argc, const char *argv[])
|
||||
gsl_matrix_set(frame_correction, 2, 0, sinh(-boost));
|
||||
gsl_matrix_set(frame_correction, 2, 2, cosh(-boost));
|
||||
multiply_left(frame_correction, frame, ws);
|
||||
*/
|
||||
gsl_matrix_free(frame_correction);
|
||||
*/
|
||||
|
||||
// int indices[10] = {0, 1, 6, 10, 30, 46, 124, 185, 484, 717};
|
||||
// int indices[10] = {0, 1, 4, 10, 22};
|
||||
// the actual drawing
|
||||
point transformed[3];
|
||||
point reflection_lines[3];
|
||||
point triangle_points[3];
|
||||
|
||||
print_svg_header();
|
||||
|
||||
for(int k = 0; k < elements; k++) {
|
||||
if(group[k].length % 2)
|
||||
continue;
|
||||
|
||||
LOOP(i) reflection_lines[i] = row(cartan, i);
|
||||
LOOP(i) triangle_points[i] = incidence(reflection_lines[(i+1)%3], reflection_lines[(i+2)%3]);
|
||||
LOOP(i) transformed[i] = apply(matrices[k], triangle_points[i]);
|
||||
// draw_triangle(transformed, frame, "black,fill=black!10,line width=0pt");
|
||||
draw_triangle(transformed, frame, "fill:#cfcfcf;");
|
||||
// draw_triangle(transformed, frame, "fill:#000000;");
|
||||
}
|
||||
|
||||
// draw special elements
|
||||
for(int k = 0; k < elements; k++) {
|
||||
// if(group[k].length % 2)
|
||||
// continue;
|
||||
|
||||
LOOP(i) transformed[i] = apply(matrices[k], edge_midpoints[(i+2)%3]);
|
||||
LOOP(i) transformed2[i] = apply(matrices[k], coxeter_repelling[i]);
|
||||
LOOP(i) transformed3[i] = apply(matrices[k], coxeter_repelling[(i+1)%3]);
|
||||
LOOP(i) transformed4[i] = apply(matrices[k], coxeter_attracting[i%3]);
|
||||
LOOP(i) transformed5[i] = apply(matrices[k], coxeter2_repelling[i]);
|
||||
LOOP(i) transformed6[i] = apply(matrices[k], coxeter2_attracting[i%3]);
|
||||
//LOOP(i) draw_line(transformed2[i], transformed4[i], frame, "red");
|
||||
draw_line(transformed2[0], transformed4[0], frame, "fill:none;stroke:red;stroke-width:1;");
|
||||
// draw_line(transformed2[1], transformed4[1], frame, "fill:none;stroke:blue;stroke-width:1;");
|
||||
// draw_line(transformed2[2], transformed4[2], frame, "fill:none;stroke:darkgreen;stroke-width:1;");
|
||||
|
||||
// draw_line(transformed5[0], transformed6[0], frame, "fill:none;stroke:blue;stroke-width:1;");
|
||||
// draw_line(transformed5[1], transformed6[1], frame, "fill:none;stroke:blue;stroke-width:1;");
|
||||
// draw_line(transformed5[2], transformed6[2], frame, "fill:none;stroke:blue;stroke-width:1;");
|
||||
// draw_line(transformed2[1], transformed4[1], frame, "fill:none;stroke:darkgreen;stroke-width:1;");
|
||||
// draw_line(transformed2[2], transformed4[2], frame, "fill:none;stroke:blue;stroke-width:1;");
|
||||
// LOOP(i) draw_line(transformed[i], transformed3[i], frame, "red");
|
||||
|
||||
// LOOP(i) transformed[i] = apply(matrices[k], coxeter_attracting[i]);
|
||||
// draw_line(transformed[1], transformed[2], frame, "red");
|
||||
for(int i = 0; i < 3*nspecial; i++) {
|
||||
// draw_dot(apply(matrices[k], special_repelling[i]), frame, "fill:red;stroke-width:0;");
|
||||
// draw_dot(apply(matrices[k], special_attracting[i]), frame, "fill:blue;stroke-width:0;");
|
||||
draw_line(apply(matrices[k], special_repelling[i]),
|
||||
apply(matrices[k], special_attracting[i]),
|
||||
frame, "fill:none;stroke:red;stroke-width:1;");
|
||||
}
|
||||
}
|
||||
|
||||
/*
|
||||
draw_line(apply(matrices[0], coxeter_repelling[0]),
|
||||
apply(matrices[0], coxeter_attracting[0]),
|
||||
frame, "fill:none;stroke:red;stroke-width:1;");
|
||||
draw_line(apply(matrices[1], coxeter_repelling[0]),
|
||||
apply(matrices[1], coxeter_attracting[0]),
|
||||
frame, "fill:none;stroke:red;stroke-width:1;");
|
||||
draw_line(apply(matrices[2], coxeter_repelling[0]),
|
||||
apply(matrices[2], coxeter_attracting[0]),
|
||||
frame, "fill:none;stroke:red;stroke-width:1;");
|
||||
draw_line(apply(matrices[3], coxeter_repelling[0]),
|
||||
apply(matrices[3], coxeter_attracting[0]),
|
||||
frame, "fill:none;stroke:red;stroke-width:1;");
|
||||
draw_line(apply(matrices[4], coxeter_repelling[0]),
|
||||
apply(matrices[4], coxeter_attracting[0]),
|
||||
frame, "fill:none;stroke:red;stroke-width:1;");
|
||||
draw_line(apply(matrices[5], coxeter_repelling[0]),
|
||||
apply(matrices[5], coxeter_attracting[0]),
|
||||
frame, "fill:none;stroke:red;stroke-width:1;");
|
||||
draw_line(apply(matrices[14], coxeter_repelling[0]),
|
||||
apply(matrices[14], coxeter_attracting[0]),
|
||||
frame, "fill:none;stroke:red;stroke-width:1;");
|
||||
draw_line(apply(matrices[15], coxeter_repelling[0]),
|
||||
apply(matrices[15], coxeter_attracting[0]),
|
||||
frame, "fill:none;stroke:red;stroke-width:1;");
|
||||
draw_line(apply(matrices[22], coxeter_repelling[0]),
|
||||
apply(matrices[22], coxeter_attracting[0]),
|
||||
frame, "fill:none;stroke:red;stroke-width:1;");
|
||||
draw_line(apply(matrices[23], coxeter_repelling[0]),
|
||||
apply(matrices[23], coxeter_attracting[0]),
|
||||
frame, "fill:none;stroke:red;stroke-width:1;");
|
||||
*/
|
||||
|
||||
print_svg_footer();
|
||||
|
||||
// clean up
|
||||
free(group);
|
||||
for(int i = 0; i < elements; i++)
|
||||
gsl_matrix_free(matrices[i]);
|
||||
free(matrices);
|
||||
gsl_matrix_free(cartan);
|
||||
gsl_matrix_free(frame);
|
||||
LOOP(i) gsl_matrix_free(gen[i]);
|
||||
workspace_free(ws);
|
||||
for(int i = 0; i < 3*nspecial; i++) {
|
||||
gsl_matrix_free(special[i]);
|
||||
gsl_matrix_free(special_eigenvectors[i]);
|
||||
}
|
||||
free(special);
|
||||
free(special_eigenvectors);
|
||||
free(special_attracting);
|
||||
free(special_repelling);
|
||||
}
|
||||
|
||||
10
linalg.c
10
linalg.c
@@ -145,11 +145,11 @@ double determinant(gsl_matrix *g, workspace_t *ws)
|
||||
return gsl_linalg_LU_det(ws->tmp, s);
|
||||
}
|
||||
|
||||
void jordan_calc(gsl_matrix *g, double *evs, workspace_t *ws)
|
||||
int jordan_calc(gsl_matrix *g, double *evs, workspace_t *ws)
|
||||
{
|
||||
gsl_eigen_nonsymmv_params(1, ws->work_nonsymmv);
|
||||
gsl_eigen_nonsymmv(g, ws->eval_complex, ws->evec_complex, ws->work_nonsymmv);
|
||||
gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_VAL_DESC);
|
||||
gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_ABS_DESC);
|
||||
|
||||
int real = 1;
|
||||
for(int i = 0; i < ws->n; i++)
|
||||
@@ -165,6 +165,8 @@ void jordan_calc(gsl_matrix *g, double *evs, workspace_t *ws)
|
||||
evs[2*i] = GSL_REAL(gsl_vector_complex_get(ws->eval_complex, i));
|
||||
evs[2*i+1] = GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i));
|
||||
}
|
||||
|
||||
return real;
|
||||
}
|
||||
|
||||
void eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws)
|
||||
@@ -182,8 +184,8 @@ void eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws)
|
||||
real = 0;
|
||||
|
||||
if(!real) {
|
||||
printf("We have non-real eigenvalues!\n");
|
||||
exit(1);
|
||||
fprintf(stderr,"We have non-real eigenvalues!\n");
|
||||
// exit(1);
|
||||
}
|
||||
|
||||
for(int i = 0; i < ws->n; i++)
|
||||
|
||||
2
linalg.h
2
linalg.h
@@ -38,7 +38,7 @@ void multiply_many(workspace_t *ws, gsl_matrix *out, int n, ...);
|
||||
void cartan_calc(gsl_matrix *g, double *mu, workspace_t *ws);
|
||||
void initialize(gsl_matrix *g, double *data, int x, int y);
|
||||
void rotation_matrix(gsl_matrix *g, double *vector);
|
||||
void jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws);
|
||||
int jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws);
|
||||
double trace(gsl_matrix *g);
|
||||
double determinant(gsl_matrix *g, workspace_t *ws);
|
||||
void eigenvectors(gsl_matrix *g, gsl_matrix *evec, workspace_t *ws);
|
||||
|
||||
Reference in New Issue
Block a user