401 lines
12 KiB
C
401 lines
12 KiB
C
#include <complex.h>
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
|
|
#include "triangle.h"
|
|
#include "linalg.h"
|
|
|
|
#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);
|
|
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 initialize_triangle_generators(gsl_matrix **gen, gsl_matrix *cartan)
|
|
{
|
|
LOOP(i) {
|
|
gsl_matrix_set_identity(gen[i]);
|
|
LOOP(j) *gsl_matrix_ptr(gen[i], i, j) += gsl_matrix_get(cartan, i, j);
|
|
}
|
|
}
|
|
|
|
typedef struct { double x[3]; } point;
|
|
|
|
point apply(gsl_matrix *g, point x)
|
|
{
|
|
point out;
|
|
LOOP(i) out.x[i] = 0.0;
|
|
LOOP(i) LOOP(j) out.x[i] += gsl_matrix_get(g, i, j) * x.x[j];
|
|
return out;
|
|
}
|
|
|
|
point incidence(point a, point b)
|
|
{
|
|
point c;
|
|
LOOP(i) c.x[i] = a.x[(i+1)%3]*b.x[(i+2)%3] - a.x[(i+2)%3]*b.x[(i+1)%3];
|
|
return c;
|
|
}
|
|
|
|
double applyForm(point x, point y, gsl_matrix *form)
|
|
{
|
|
double out = 0.0;
|
|
LOOP(i) LOOP(j) out += x.x[i] * gsl_matrix_get(form, i, j) * y.x[j];
|
|
return out;
|
|
}
|
|
|
|
void boundary(point x, point y, point *a, point *b, gsl_matrix *form)
|
|
{
|
|
double formX = applyForm(x, x, form);
|
|
double formY = applyForm(y, y, form);
|
|
double formXY = applyForm(x, y, form);
|
|
double t1 = (- formXY + sqrt(formXY * formXY - formX * formY)) / formX;
|
|
double t2 = (- formXY - sqrt(formXY * formXY - formX * formY)) / formX;
|
|
LOOP(i) a->x[i] = t1 * x.x[i] + y.x[i];
|
|
LOOP(i) b->x[i] = t2 * x.x[i] + y.x[i];
|
|
}
|
|
|
|
point column(gsl_matrix *m, int j)
|
|
{
|
|
point out;
|
|
LOOP(i) out.x[i] = gsl_matrix_get(m, i, j);
|
|
return out;
|
|
}
|
|
|
|
point row(gsl_matrix *m, int i)
|
|
{
|
|
point out;
|
|
LOOP(j) out.x[j] = gsl_matrix_get(m, i, j);
|
|
return out;
|
|
}
|
|
|
|
double coord(point p, int k, gsl_matrix *frame)
|
|
{
|
|
double tmp[3] = {0, 0, 0};
|
|
double factor = 1.0;
|
|
|
|
LOOP(i) LOOP(j) tmp[i] += gsl_matrix_get(frame, i, j)*p.x[j];
|
|
|
|
#ifdef POINCARE
|
|
double norm2 = (tmp[0]*tmp[0] + tmp[1]*tmp[1]) / (tmp[2]*tmp[2]);
|
|
if(fabs(norm2) < 0.5) // just to avoid dividing by 0
|
|
factor = 1.0 / (1.0 + sqrt(1.0 - norm2));
|
|
else if(fabs(norm2) < 1.0)
|
|
factor = (1.0 - sqrt(1.0 - norm2)) / norm2;
|
|
else
|
|
factor = 1.0;
|
|
#endif
|
|
|
|
return factor*tmp[k]/tmp[2];
|
|
}
|
|
|
|
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";
|
|
|
|
printf("%s", header);
|
|
}
|
|
|
|
void print_tex_footer() {
|
|
char *footer =
|
|
"\\draw[thick] (0,0) circle (1.0);\n"
|
|
"\\end{tikzpicture}\n"
|
|
"\\end{document}\n";
|
|
|
|
printf("%s", footer);
|
|
}
|
|
|
|
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)
|
|
{
|
|
// 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);
|
|
|
|
// 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 0 0 %d %f %f", DCONV(r), DCONV(r), t>0, CONV(bx), CONV(by));
|
|
return buffer;
|
|
}
|
|
|
|
void draw_triangle(point *p, gsl_matrix *frame, const char *arguments)
|
|
{
|
|
#ifdef POINCARE
|
|
char buffer1[100], buffer2[100], buffer3[100];
|
|
|
|
double x1 = coord(p[0], 0, frame);
|
|
double y1 = coord(p[0], 1, frame);
|
|
double x2 = coord(p[1], 0, frame);
|
|
double y2 = coord(p[1], 1, frame);
|
|
double x3 = coord(p[2], 0, frame);
|
|
double y3 = coord(p[2], 1, frame);
|
|
|
|
// 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)); */
|
|
|
|
// 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,
|
|
coord(p[0], 0, frame), coord(p[0], 1, frame),
|
|
coord(p[1], 0, frame), coord(p[1], 1, frame),
|
|
coord(p[2], 0, frame), coord(p[2], 1, frame));
|
|
#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=\"6\" 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];
|
|
|
|
double x1 = coord(p1, 0, frame);
|
|
double y1 = coord(p1, 1, frame);
|
|
double x2 = coord(p2, 0, frame);
|
|
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("<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
|
|
}
|
|
|
|
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[])
|
|
{
|
|
groupelement_t *group;
|
|
gsl_matrix **matrices;
|
|
gsl_matrix *cartan;
|
|
gsl_matrix *gen[3];
|
|
gsl_matrix **special;
|
|
gsl_matrix **special_eigenvectors;
|
|
point *special_attracting, *special_repelling, *special_rotation;
|
|
gsl_matrix *frame;
|
|
workspace_t *ws;
|
|
int p,q,r;
|
|
int elements, nspecial, nspecial_hyp, nspecial_rot;
|
|
|
|
if(argc < 5) {
|
|
fprintf(stderr, "Usage: %s <p> <q> <r> <n_elements> <word1> <word2> ...\n", argv[0]);
|
|
exit(1);
|
|
}
|
|
|
|
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*));
|
|
for(int i = 0; i < elements; i++)
|
|
matrices[i] = gsl_matrix_alloc(3, 3);
|
|
cartan = gsl_matrix_alloc(3, 3);
|
|
frame = gsl_matrix_alloc(3, 3);
|
|
LOOP(i) gen[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));
|
|
special_rotation = 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);
|
|
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]);
|
|
|
|
nspecial_hyp = nspecial_rot = 0;
|
|
for(int i = 0; i < nspecial; i++) {
|
|
int j = 0;
|
|
int nreal;
|
|
|
|
compute_word(ws, special[3*i+j], gen, argv[i+5], j, 0);
|
|
nreal = eigenvectors(special[3*i+j], special_eigenvectors[3*i+j], ws);
|
|
if(nreal == 3) {
|
|
special_attracting[nspecial_hyp] = column(special_eigenvectors[3*i+j], 0);
|
|
|
|
// 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[nspecial_hyp] = column(special_eigenvectors[3*i+j], 0);
|
|
nspecial_hyp++;
|
|
} else {
|
|
special_rotation[nspecial_rot] = column(special_eigenvectors[3*i+j], 0);
|
|
nspecial_rot++;
|
|
}
|
|
}
|
|
|
|
fprintf(stderr, "%d special elements, %d rotations, %d hyperbolic\n", nspecial, nspecial_rot, nspecial_hyp);
|
|
|
|
/*
|
|
// 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);
|
|
*/
|
|
|
|
// 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, "fill:#cfcfcf;");
|
|
}
|
|
|
|
char stylestring[100];
|
|
char colors[3][20] = {"red", "blue", "green"};
|
|
|
|
// draw special elements
|
|
for(int k = 0; k < elements; k++) {
|
|
// if(group[k].length % 2)
|
|
// continue;
|
|
|
|
for(int i = 0; i < nspecial_hyp; i+=3) {
|
|
// draw_dot(apply(matrices[k], special_repelling[i]), frame, "fill:red;stroke-width:1;");
|
|
// draw_dot(apply(matrices[k], special_attracting[i]), frame, "fill:blue;stroke-width:1;");
|
|
snprintf(stylestring, sizeof(stylestring), "fill:none;stroke:%s;stroke-width:1;", colors[i%3]);
|
|
draw_line(apply(matrices[k], special_repelling[i]),
|
|
apply(matrices[k], special_attracting[i]),
|
|
frame, stylestring);
|
|
}
|
|
|
|
for(int i = 0; i < nspecial_rot; i+=3) {
|
|
snprintf(stylestring, sizeof(stylestring), "fill:%s;stroke:%s;stroke-width:1;", colors[i%3], colors[i%3]);
|
|
// fprintf(stderr, "%f %f\n", special_rotation[i].x[0], special_rotation[i].x[1]);
|
|
draw_dot(apply(matrices[k], special_rotation[i]), frame, stylestring);
|
|
}
|
|
}
|
|
|
|
|
|
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);
|
|
free(special_rotation);
|
|
}
|