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