switch hyperbolic program to svg

This commit is contained in:
Florian Stecker 2022-02-04 09:17:46 -06:00
parent 7fb5c1f442
commit bea92930ca
1 changed files with 176 additions and 28 deletions

View File

@ -8,6 +8,10 @@
#define POINCARE 1
#define LOOP(i) for(int i = 0; i < 3; i++)
#define CUTOFF(x) ((x)>1000.0?1000.0:(x)<-1000.0?-1000.0:(x))
#define CONV(x) (CUTOFF(x)*490+500)
#define DCONV(x) (CUTOFF(x)*490)
void cartan_matrix(gsl_matrix *cartan, double a1, double a2, double a3, double s)
{
gsl_matrix_set(cartan, 0, 0, -2);
@ -104,6 +108,7 @@ void print_tex_header() {
char *header =
"\\documentclass{standalone}\n"
"\\usepackage[utf8]{inputenc}\n"
"\\usepackage[dvipsnames]{xcolor}\n"
"\\usepackage{tikz}\n"
"\\begin{document}\n"
"\\begin{tikzpicture}[scale=5]\n";
@ -120,22 +125,43 @@ void print_tex_footer() {
printf("%s", footer);
}
char *poincare_arc(double x1, double x2, double y1, double y2, char *buffer)
void print_svg_header() {
char *header =
"<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n"
"<svg width=\"1000\" height=\"1000\" xmlns=\"http://www.w3.org/2000/svg\">\n";
printf("%s", header);
}
void print_svg_footer() {
printf("<circle cx=\"%f\" cy=\"%f\" r=\"%f\" style=\"fill:none;stroke-width:3;stroke:black;\" />\n",
CONV(0.0), CONV(0.0), DCONV(1.0));
printf("</svg>\n");
}
char *poincare_arc(double ax, double ay, double bx, double by, char *buffer)
{
double t = (1 - x1*y1 - x2*y2)/(x1*y2 - x2*y1)/2;
double m1 = x1/2 + y1/2 + t*(y2 - x2); // center of circle
double m2 = x2/2 + y2/2 + t*(x1 - y1);
// find m = (a+b)/2 + tc with c orthogonal to a-b and |m|^2 = |m-a|^2 + 1
// m^2 = m^2 - 2<m,a> + a^2 + 1
// 0 = -2<m,a> + a^2 + 1 = -<a+b,a> - 2t<c,a> + a^2 + 1 = -<b,a> - 2t<c,a> + 1
// t = (1 - <a,b>)/<a,c>/2
double cx = ay - by;
double cy = bx - ax;
double t = (1 - ax*bx - ay*by)/(ax*cx + ay*cy)/2;
double mx = ax/2 + bx/2 + t*cx; // center of circle
double my = ay/2 + by/2 + t*cy;
double r = sqrt((ax-mx)*(ax-mx) + (ay-my)*(ay-my));
//double phix = atan2(x2-m2,x1-m1);
//double phiy = atan2(y2-m2,y1-m1);
double r = sqrt((x1-m1)*(x1-m1) + (x2-m2)*(x2-m2));
double phix = atan2(x2-m2,x1-m1);
double phiy = atan2(y2-m2,y1-m1);
// if(phix - phiy > M_PI)
// phiy += 2*M_PI;
// else if(phiy - phix > M_PI)
// phix += 2*M_PI;
if(phix - phiy > M_PI)
phiy += 2*M_PI;
else if(phiy - phix > M_PI)
phix += 2*M_PI;
sprintf(buffer, "%f:%f:%f", phix/M_PI*180, phiy/M_PI*180, r);
// sprintf(buffer, "%f:%f:%f", phix/M_PI*180, phiy/M_PI*180, r);
sprintf(buffer, "%f %f 0 0 %d %f %f", DCONV(r), DCONV(r), t>0, CONV(bx), CONV(by));
return buffer;
}
@ -151,10 +177,18 @@ void draw_triangle(point *p, gsl_matrix *frame, const char *arguments)
double x3 = coord(p[2], 0, frame);
double y3 = coord(p[2], 1, frame);
printf("\\draw[%s] (%f, %f) arc (%s) arc (%s) arc (%s);\n", arguments, x1, y1,
// TikZ version
/* printf("\\draw[%s] (%f, %f) arc (%s) arc (%s) arc (%s);\n", arguments, x1, y1,
poincare_arc(x1, y1, x2, y2, buffer1),
poincare_arc(x2, y2, x3, y3, buffer2),
poincare_arc(x3, y3, x1, y1, buffer3));
poincare_arc(x3, y3, x1, y1, buffer3)); */
// SVG version
printf("<path d=\"M %f %f A %s A %s A %s Z\" style=\"%s\" />\n", CONV(x1), CONV(y1),
poincare_arc(x1, y1, x2, y2, buffer1),
poincare_arc(x2, y2, x3, y3, buffer2),
poincare_arc(x3, y3, x1, y1, buffer3),
arguments);
#else
printf("\\draw[%s] (%f, %f) -- (%f, %f) -- (%f, %f) -- cycle;\n", arguments,
@ -174,25 +208,42 @@ void draw_line(point p1, point p2, gsl_matrix *frame, const char *arguments)
double y2 = coord(p2, 1, frame);
#ifdef POINCARE
printf("\\draw[%s] (%f, %f) arc (%s);\n", arguments, x1, y1,
poincare_arc(x1, y1, x2, y2, buffer));
// printf("\\draw[%s] (%f, %f) arc (%s);\n", arguments, x1, y1,
// poincare_arc(x1, y1, x2, y2, buffer));
printf("<path d=\"M %f %f A %s\" style=\"%s\"/>\n", CONV(x1), CONV(y1),
poincare_arc(x1, y1, x2, y2, buffer), arguments);
#else
printf("\\draw[%s] (%f, %f) -- (%f, %f);\n", arguments, x1, y1, x2, y2);
#endif
}
int main()
void compute_word(workspace_t *ws, gsl_matrix *result, gsl_matrix **gen, const char *word, int modifier)
{
gsl_matrix_set_identity(result);
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[])
{
groupelement_t *group;
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 *frame;
workspace_t *ws;
int elements = 2000;
int p = 4, q = 4, r = 4;
if(argc < 5) {
fprintf(stderr, "Usage: %s <p> <q> <r> <n_elements>\n", argv[0]);
exit(1);
}
int elements = atoi(argv[4]);
int p = atoi(argv[1]), q = atoi(argv[2]), r = atoi(argv[3]);
group = malloc(elements*sizeof(groupelement_t));
matrices = malloc(elements*sizeof(gsl_matrix*));
@ -203,6 +254,8 @@ int main()
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);
generate_triangle_group(group, elements, p, q, r);
@ -214,8 +267,16 @@ int main()
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]);
// 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 < elements; i++) {
@ -230,12 +291,19 @@ int main()
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);
@ -243,8 +311,35 @@ int main()
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_tex_header();
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));
gsl_matrix_set(frame_correction, 1, 1, cos(angle-M_PI/2));
gsl_matrix_set(frame_correction, 2, 2, 1);
// multiply_left(frame_correction, frame, ws);
gsl_matrix_set_identity(frame_correction);
gsl_matrix_set(frame_correction, 0, 0, cosh(-boost));
gsl_matrix_set(frame_correction, 0, 2, sinh(-boost));
gsl_matrix_set(frame_correction, 1, 1, 1);
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};
@ -254,16 +349,69 @@ int main()
continue;
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, "black,fill=black!10,line width=0pt");
draw_triangle(transformed, frame, "fill:#cfcfcf;");
// draw_triangle(transformed, frame, "fill:#000000;");
}
for(int k = 0; k < 2000; k++) {
if(group[k].length % 2)
continue;
for(int k = 0; k < elements; k++) {
// if(group[k].length % 2)
// continue;
LOOP(i) transformed[i] = apply(matrices[k], coxeter_attracting[i]);
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");
}
print_tex_footer();
/*
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();
}