Trisurf Monte Carlo simulator
Samo Penic
2014-04-25 1665aaef46e1bf85eb99f709bb596ab4fcc07f59
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