diff --git a/hyperbolic.c b/hyperbolic.c index 7cb7c74..ce2b132 100644 --- a/hyperbolic.c +++ b/hyperbolic.c @@ -204,7 +204,7 @@ void draw_dot(point p, gsl_matrix *frame, const char *arguments) double y = coord(p, 1, frame); #ifdef POINCARE - printf("\n", CONV(x), CONV(y), arguments); + printf("\n", CONV(x), CONV(y), arguments); #endif } @@ -247,11 +247,11 @@ int main(int argc, const char *argv[]) gsl_matrix *gen[3]; gsl_matrix **special; gsl_matrix **special_eigenvectors; - point *special_attracting, *special_repelling; + point *special_attracting, *special_repelling, *special_rotation; gsl_matrix *frame; workspace_t *ws; int p,q,r; - int elements, nspecial; + int elements, nspecial, nspecial_hyp, nspecial_rot; if(argc < 5) { fprintf(stderr, "Usage: %s

...\n", argv[0]); @@ -274,6 +274,7 @@ int main(int argc, const char *argv[]) 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); @@ -288,17 +289,29 @@ int main(int argc, const char *argv[]) for(int i = 1; i < elements; 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) { - 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); + int nreal; - // 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); + 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 @@ -348,24 +361,21 @@ int main(int argc, const char *argv[]) // if(group[k].length % 2) // continue; - for(int i = 0; i < 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;"); + 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[3*i]), - apply(matrices[k], special_attracting[3*i]), + draw_line(apply(matrices[k], special_repelling[i]), + apply(matrices[k], special_attracting[i]), frame, stylestring); } - } - /* - int k = 0; - for(int i = 0; i < nspecial; i++) { - draw_line(apply(matrices[k], special_repelling[3*i]), - apply(matrices[k], special_attracting[3*i]), - frame, "fill:none;stroke:black;stroke-width:1;"); + 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(); @@ -387,4 +397,5 @@ int main(int argc, const char *argv[]) free(special_eigenvectors); free(special_attracting); free(special_repelling); + free(special_rotation); } diff --git a/linalg.c b/linalg.c index 7d06f43..3f3a275 100644 --- a/linalg.c +++ b/linalg.c @@ -169,7 +169,7 @@ int jordan_calc(gsl_matrix *g, double *evs, workspace_t *ws) 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_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); - int real = 1; - for(int i = 0; i < ws->n; i++) - if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) != 0) - real = 0; - - if(!real) { - fprintf(stderr,"We have non-real eigenvalues!\n"); -// exit(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++) + gsl_matrix_set(evec_real, i, real, GSL_REAL(gsl_matrix_complex_get(ws->evec_complex, i, j))); + real++; + } } - for(int i = 0; i < ws->n; i++) - 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))); - + return real; } void eigensystem_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws) diff --git a/linalg.h b/linalg.h index 1077090..880d160 100644 --- a/linalg.h +++ b/linalg.h @@ -41,7 +41,7 @@ void rotation_matrix(gsl_matrix *g, double *vector); 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); +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); int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws);