allow rotations

This commit is contained in:
Florian Stecker 2023-01-08 10:22:53 -05:00
parent 040f994d6f
commit c62044d637
3 changed files with 45 additions and 38 deletions

View File

@ -204,7 +204,7 @@ void draw_dot(point p, gsl_matrix *frame, const char *arguments)
double y = coord(p, 1, frame); double y = coord(p, 1, frame);
#ifdef POINCARE #ifdef POINCARE
printf("<circle cx=\"%f\" cy=\"%f\" r=\"1\" style=\"%s\"/>\n", CONV(x), CONV(y), arguments); printf("<circle cx=\"%f\" cy=\"%f\" r=\"6\" style=\"%s\"/>\n", CONV(x), CONV(y), arguments);
#endif #endif
} }
@ -247,11 +247,11 @@ int main(int argc, const char *argv[])
gsl_matrix *gen[3]; gsl_matrix *gen[3];
gsl_matrix **special; gsl_matrix **special;
gsl_matrix **special_eigenvectors; gsl_matrix **special_eigenvectors;
point *special_attracting, *special_repelling; point *special_attracting, *special_repelling, *special_rotation;
gsl_matrix *frame; gsl_matrix *frame;
workspace_t *ws; workspace_t *ws;
int p,q,r; int p,q,r;
int elements, nspecial; int elements, nspecial, nspecial_hyp, nspecial_rot;
if(argc < 5) { if(argc < 5) {
fprintf(stderr, "Usage: %s <p> <q> <r> <n_elements> <word1> <word2> ...\n", argv[0]); fprintf(stderr, "Usage: %s <p> <q> <r> <n_elements> <word1> <word2> ...\n", argv[0]);
@ -274,6 +274,7 @@ int main(int argc, const char *argv[])
special_eigenvectors = malloc(3*nspecial*sizeof(gsl_matrix*)); special_eigenvectors = malloc(3*nspecial*sizeof(gsl_matrix*));
special_attracting = malloc(3*nspecial*sizeof(point)); special_attracting = malloc(3*nspecial*sizeof(point));
special_repelling = 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++) { for(int i = 0; i < 3*nspecial; i++) {
special[i] = gsl_matrix_alloc(3, 3); special[i] = gsl_matrix_alloc(3, 3);
special_eigenvectors[i] = gsl_matrix_alloc(3, 3); special_eigenvectors[i] = gsl_matrix_alloc(3, 3);
@ -288,17 +289,29 @@ int main(int argc, const char *argv[])
for(int i = 1; i < elements; i++) for(int i = 1; i < elements; i++)
multiply(matrices[group[i].parent->id], gen[group[i].letter], matrices[i]); multiply(matrices[group[i].parent->id], gen[group[i].letter], matrices[i]);
for(int i = 0; i < nspecial; i++) nspecial_hyp = nspecial_rot = 0;
for(int i = 0; i < nspecial; i++) {
LOOP(j) { LOOP(j) {
int nreal;
compute_word(ws, special[3*i+j], gen, argv[i+5], j, 0); compute_word(ws, special[3*i+j], gen, argv[i+5], j, 0);
eigenvectors(special[3*i+j], special_eigenvectors[3*i+j], ws); nreal = eigenvectors(special[3*i+j], special_eigenvectors[3*i+j], ws);
special_attracting[3*i+j] = column(special_eigenvectors[3*i+j], 0); if(nreal == 3) {
special_attracting[nspecial_hyp] = column(special_eigenvectors[3*i+j], 0);
// repelling = attracting of inverse // repelling = attracting of inverse
compute_word(ws, special[3*i+j], gen, argv[i+5], j, 1); compute_word(ws, special[3*i+j], gen, argv[i+5], j, 1);
eigenvectors(special[3*i+j], special_eigenvectors[3*i+j], ws); eigenvectors(special[3*i+j], special_eigenvectors[3*i+j], ws);
special_repelling[3*i+j] = column(special_eigenvectors[3*i+j], 0); 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 // let's correct the frame of reference by using hyperbolic transformations
@ -348,24 +361,21 @@ int main(int argc, const char *argv[])
// if(group[k].length % 2) // if(group[k].length % 2)
// continue; // continue;
for(int i = 0; i < nspecial; i++) { for(int i = 0; i < nspecial_hyp; i+=3) {
// draw_dot(apply(matrices[k], special_repelling[i]), frame, "fill:red;stroke-width:0;"); // 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:0;"); // 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]); snprintf(stylestring, sizeof(stylestring), "fill:none;stroke:%s;stroke-width:1;", colors[i%3]);
draw_line(apply(matrices[k], special_repelling[3*i]), draw_line(apply(matrices[k], special_repelling[i]),
apply(matrices[k], special_attracting[3*i]), apply(matrices[k], special_attracting[i]),
frame, stylestring); frame, stylestring);
} }
}
/* for(int i = 0; i < nspecial_rot; i+=3) {
int k = 0; snprintf(stylestring, sizeof(stylestring), "fill:%s;stroke:%s;stroke-width:1;", colors[i%3], colors[i%3]);
for(int i = 0; i < nspecial; i++) { // fprintf(stderr, "%f %f\n", special_rotation[i].x[0], special_rotation[i].x[1]);
draw_line(apply(matrices[k], special_repelling[3*i]), draw_dot(apply(matrices[k], special_rotation[i]), frame, stylestring);
apply(matrices[k], special_attracting[3*i]), }
frame, "fill:none;stroke:black;stroke-width:1;");
} }
*/
print_svg_footer(); print_svg_footer();
@ -387,4 +397,5 @@ int main(int argc, const char *argv[])
free(special_eigenvectors); free(special_eigenvectors);
free(special_attracting); free(special_attracting);
free(special_repelling); free(special_repelling);
free(special_rotation);
} }

View File

@ -169,7 +169,7 @@ int jordan_calc(gsl_matrix *g, double *evs, workspace_t *ws)
return real; return real;
} }
void eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws) int eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws)
{ {
gsl_matrix_memcpy(ws->stack[0], g); gsl_matrix_memcpy(ws->stack[0], g);
gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv); gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv);
@ -178,20 +178,16 @@ void eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws)
gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_ABS_DESC); gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_ABS_DESC);
int real = 1; int real = 0;
for(int j = 0; j < ws->n; j++) {
if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, j)), 0) == 0) {
for(int i = 0; i < ws->n; i++) for(int i = 0; i < ws->n; i++)
if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) != 0) gsl_matrix_set(evec_real, i, real, GSL_REAL(gsl_matrix_complex_get(ws->evec_complex, i, j)));
real = 0; real++;
}
if(!real) {
fprintf(stderr,"We have non-real eigenvalues!\n");
// exit(1);
} }
for(int i = 0; i < ws->n; i++) return real;
for(int j = 0; j < ws->n; j++)
gsl_matrix_set(evec_real, i, j, GSL_REAL(gsl_matrix_complex_get(ws->evec_complex, i, j)));
} }
void eigensystem_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws) void eigensystem_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws)

View File

@ -41,7 +41,7 @@ void rotation_matrix(gsl_matrix *g, double *vector);
int 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 trace(gsl_matrix *g);
double determinant(gsl_matrix *g, workspace_t *ws); double determinant(gsl_matrix *g, workspace_t *ws);
void eigenvectors(gsl_matrix *g, gsl_matrix *evec, workspace_t *ws); int eigenvectors(gsl_matrix *g, gsl_matrix *evec, workspace_t *ws);
void eigenvectors_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws); void eigenvectors_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws);
int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws); int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws);