| commit | author | age | ||
| 7f6076 | 1 | /* vim: set ts=4 sts=4 sw=4 noet : */ |
| 37d14a | 2 | #include <stdio.h> |
| SP | 3 | #include <time.h> |
| 4 | //#include <gsl/gsl_math.h> | |
| 5 | //#include <gsl/gsl_blas.h> | |
| 6 | #include <gsl/gsl_vector.h> | |
| 7 | #include <gsl/gsl_matrix.h> | |
| 8 | #include <gsl/gsl_eigen.h> | |
| 9 | #include <gsl/gsl_sort_vector.h> | |
| 10 | #include "general.h" | |
| 11 | ||
| 12 | void gyration_eigen(ts_vesicle *vesicle, ts_double *l1, ts_double *l2, ts_double *l3){ | |
| 13 | gsl_matrix *matrix = gsl_matrix_alloc(3,3); | |
| 14 | gsl_vector *eig = gsl_vector_alloc(3); | |
| 15 | gsl_eigen_symm_workspace *workspace = gsl_eigen_symm_alloc(3); | |
| 16 | ts_uint i,j,k; | |
| 17 | ts_double mat[3][3]; | |
| 18 | ts_double vec[3]; | |
| 611412 | 19 | |
| 37d14a | 20 | for(i = 0; i < 3; i++) |
| SP | 21 | for(j = 0; j < 3; j++) |
| 22 | mat[i][j]=0; | |
| 23 | ||
| 24 | for(k=0;k<vesicle->vlist->n;k++){ | |
| 25 | vec[0]=vesicle->vlist->vtx[k]->x; | |
| 26 | vec[1]=vesicle->vlist->vtx[k]->y; | |
| 27 | vec[2]=vesicle->vlist->vtx[k]->z; | |
| 28 | for(i = 0; i < 3; i++) | |
| 29 | for(j = 0; j <= i; j++) | |
| 611412 | 30 | mat[i][j]+=vec[i]*vec[j]; |
| 37d14a | 31 | } |
| SP | 32 | |
| 611412 | 33 | // Normalize gyration tensor: |
| M | 34 | for(i = 0; i < 3; i++) |
| 35 | for(j = 0; j <= i; j++) | |
| 36 | mat[i][j]=mat[i][j]/(ts_double)vesicle->vlist->n; | |
| 37 | ||
| 38 | ||
| 37d14a | 39 | // diagonal elements are copied twice! |
| SP | 40 | for(i = 0; i < 3; i++) |
| 41 | for(j = 0; j <= i; j++){ | |
| 42 | gsl_matrix_set(matrix,i,j,mat[i][j]); | |
| 43 | gsl_matrix_set(matrix,j,i,mat[i][j]); | |
| 44 | } | |
| 45 | /* | |
| 46 | for(i = 0; i < 3; i++) | |
| 47 | { | |
| 48 | for(j = 0; j < 3; j++) | |
| 49 | printf("%g ", gsl_matrix_get(matrix, i, j)); | |
| 50 | printf("\n"); | |
| 51 | } | |
| 52 | printf("\n"); | |
| 53 | */ | |
| 54 | gsl_eigen_symm(matrix, eig, workspace); | |
| 55 | gsl_sort_vector(eig); | |
| 56 | /* printf("** eigenvalues \n"); | |
| 57 | for(i = 0; i < 3; i++){ | |
| 58 | printf("%3d %25.17e\n", i, gsl_vector_get(eig, i)); | |
| 59 | } */ | |
| 60 | *l1=gsl_vector_get(eig,0); | |
| 61 | *l2=gsl_vector_get(eig,1); | |
| 62 | *l3=gsl_vector_get(eig,2); | |
| 63 | ||
| 64 | gsl_eigen_symm_free(workspace); | |
| 65 | gsl_matrix_free(matrix); | |
| 66 | gsl_vector_free(eig); | |
| 67 | } | |
| 68 | ||
| 69 | ts_ulong get_epoch(){ | |
| 70 | time_t currTime; | |
| 71 | currTime = time(NULL); | |
| 72 | return (ts_ulong)currTime; | |
| 73 | } | |
| 74 | ||
| 75 | ts_bool get_area_volume(ts_vesicle *vesicle, ts_double *area, ts_double *volume){ | |
| 76 | ts_uint i; | |
| 77 | *volume=0.0; | |
| 78 | *area=0.0; | |
| 79 | for(i=0;i<vesicle->tlist->n;i++){ | |
| 80 | *volume+=vesicle->tlist->tria[i]->volume; | |
| 81 | *area+=vesicle->tlist->tria[i]->area; | |
| 82 | } | |
| 83 | return TS_SUCCESS; | |
| 84 | } | |
