Trisurf Monte Carlo simulator
Samo Penic
2012-02-23 e016c401c98b6a33392b551b2c40fccded21efe3
commit | author | age
d7639a 1 #include<stdlib.h>
SP 2 #include "general.h"
3 #include "energy.h"
4 #include "vertex.h"
5 #include<math.h>
6 #include<stdio.h>
7 ts_bool mean_curvature_and_energy(ts_vesicle *vesicle){
8
f74313 9     ts_uint i;
d7639a 10     
f74313 11     ts_vertex_list *vlist=vesicle->vlist;
SP 12     ts_vertex **vtx=vlist->vtx;
d7639a 13
SP 14     for(i=0;i<vlist->n;i++){
f74313 15         energy_vertex(vtx[i]);
b01cc1 16         
d7639a 17     }
SP 18
19     return TS_SUCCESS;
20 }
21
22
23 inline ts_bool energy_vertex(ts_vertex *vtx){
24 //    ts_vertex *vtx=&vlist->vertex[n]-1; // Caution! 0 Indexed value!
25 //    ts_triangle *tristar=vtx->tristar-1;
f74313 26     ts_vertex_data *data=vtx->data;
d7639a 27     ts_uint jj;
SP 28     ts_uint jjp,jjm;
29     ts_vertex *j,*jp, *jm;
30     ts_triangle *jt;
31     ts_double s=0,xh=0,yh=0,zh=0,txn=0,tyn=0,tzn=0;
32     ts_double x1,x2,x3,ctp,ctm,tot,xlen;
33     ts_double h,ht;
f74313 34     for(jj=1; jj<=data->neigh_no;jj++){
d7639a 35         jjp=jj+1;
f74313 36         if(jjp>data->neigh_no) jjp=1;
d7639a 37         jjm=jj-1;
f74313 38         if(jjm<1) jjm=data->neigh_no;
SP 39         j=data->neigh[jj-1];
40         jp=data->neigh[jjp-1];
41         jm=data->neigh[jjm-1];
b01cc1 42 //        printf("tristar_no=%u, neigh_no=%u, jj=%u\n",data->tristar_no,data->neigh_no,jj);
f74313 43         jt=data->tristar[jj-1];
SP 44         x1=vtx_distance_sq(vtx,jp); //shouldn't be zero!
45         x2=vtx_distance_sq(j,jp); // shouldn't be zero!
46         x3=(j->data->x-jp->data->x)*(data->x-jp->data->x)+
47            (j->data->y-jp->data->y)*(data->y-jp->data->y)+
48            (j->data->z-jp->data->z)*(data->z-jp->data->z);
d7639a 49         
SP 50 #ifdef TS_DOUBLE_DOUBLE
51         ctp=x3/sqrt(x1*x2-x3*x3);
52 #endif
53 #ifdef TS_DOUBLE_FLOAT
54         ctp=x3/sqrtf(x1*x2-x3*x3);
55 #endif
56 #ifdef TS_DOUBLE_LONGDOUBLE
57         ctp=x3/sqrtl(x1*x2-x3*x3);
58 #endif
f74313 59         x1=vtx_distance_sq(vtx,jm);
SP 60         x2=vtx_distance_sq(j,jm);
61         x3=(j->data->x-jm->data->x)*(data->x-jm->data->x)+
62            (j->data->y-jm->data->y)*(data->y-jm->data->y)+
63            (j->data->z-jm->data->z)*(data->z-jm->data->z);
d7639a 64 #ifdef TS_DOUBLE_DOUBLE
SP 65         ctm=x3/sqrt(x1*x2-x3*x3);
66 #endif
67 #ifdef TS_DOUBLE_FLOAT
68         ctm=x3/sqrtf(x1*x2-x3*x3);
69 #endif
70 #ifdef TS_DOUBLE_LONGDOUBLE
71         ctm=x3/sqrtl(x1*x2-x3*x3);
72 #endif
73         tot=ctp+ctm;
74         tot=0.5*tot;
f74313 75         xlen=vtx_distance_sq(j,vtx);
d7639a 76 #ifdef  TS_DOUBLE_DOUBLE 
e016c4 77         data->bond[jj-1]->bond_length=sqrt(xlen); 
d7639a 78 #endif
SP 79 #ifdef  TS_DOUBLE_FLOAT
e016c4 80         data->bond[jj-1]->bond_length=sqrtf(xlen); 
d7639a 81 #endif
SP 82 #ifdef  TS_DOUBLE_LONGDOUBLE 
e016c4 83         data->bond[jj-1]->bond_length=sqrtl(xlen); 
d7639a 84 #endif
SP 85
e016c4 86         data->bond[jj-1]->bond_length_dual=tot*data->bond[jj-1]->bond_length;
d7639a 87
SP 88         s+=tot*xlen;
f74313 89         xh+=tot*(j->data->x - data->x);
SP 90         yh+=tot*(j->data->y - data->y);
91         zh+=tot*(j->data->z - data->z);
41a035 92         txn+=jt->xnorm;
SP 93         tyn+=jt->ynorm;
94         tzn+=jt->znorm;
d7639a 95     }
SP 96     
97     h=xh*xh+yh*yh+zh*zh;
98     ht=txn*xh+tyn*yh + tzn*zh;
99     s=s/4.0; 
100 #ifdef TS_DOUBLE_DOUBLE
101     if(ht>=0.0) {
f74313 102         data->curvature=sqrt(h);
d7639a 103     } else {
f74313 104         data->curvature=-sqrt(h);
d7639a 105     }
SP 106 #endif
107 #ifdef TS_DOUBLE_FLOAT
108     if(ht>=0.0) {
f74313 109         data->curvature=sqrtf(h);
d7639a 110     } else {
f74313 111         data->curvature=-sqrtf(h);
d7639a 112     }
SP 113 #endif
114 #ifdef TS_DOUBLE_LONGDOUBLE
115     if(ht>=0.0) {
f74313 116         data->curvature=sqrtl(h);
d7639a 117     } else {
f74313 118         data->curvature=-sqrtl(h);
d7639a 119     }
SP 120 #endif
dac2e5 121 // What is vtx->data->c?????????????? Here it is 0!
SP 122 // c is forced curvature energy for each vertex. Should be set to zero for
123 // norman circumstances.
f74313 124     data->energy=0.5*s*(data->curvature/s-data->c)*(data->curvature/s-data->c);
d7639a 125
SP 126     return TS_SUCCESS;
127 }