Trisurf Monte Carlo simulator
Samo Penic
2014-03-05 719c9febac2eaff9483fda487b57684afbb59bb2
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
fedf2b 23 inline ts_bool bond_energy(ts_bond *bond,ts_poly *poly){
304510 24 //TODO: This value to be changed and implemented in data structure:
M 25     ts_double d_relaxed=1.0;
26     bond->energy=poly->k*pow(bond->bond_length-d_relaxed,2);
fedf2b 27     return TS_SUCCESS;
M 28 };
29
d7639a 30
SP 31 inline ts_bool energy_vertex(ts_vertex *vtx){
32 //    ts_vertex *vtx=&vlist->vertex[n]-1; // Caution! 0 Indexed value!
33 //    ts_triangle *tristar=vtx->tristar-1;
8f6a69 34     //ts_vertex_data *data=vtx->data;
d7639a 35     ts_uint jj;
SP 36     ts_uint jjp,jjm;
37     ts_vertex *j,*jp, *jm;
38     ts_triangle *jt;
a63f17 39     ts_double s=0.0,xh=0.0,yh=0.0,zh=0.0,txn=0.0,tyn=0.0,tzn=0.0;
d7639a 40     ts_double x1,x2,x3,ctp,ctm,tot,xlen;
SP 41     ts_double h,ht;
20ae83 42     for(jj=1; jj<=vtx->neigh->n;jj++){
d7639a 43         jjp=jj+1;
20ae83 44         if(jjp>vtx->neigh->n) jjp=1;
d7639a 45         jjm=jj-1;
20ae83 46         if(jjm<1) jjm=vtx->neigh->n;
SP 47         j=vtx->neigh->vtx[jj-1];
48         jp=vtx->neigh->vtx[jjp-1];
49         jm=vtx->neigh->vtx[jjm-1];
b01cc1 50 //        printf("tristar_no=%u, neigh_no=%u, jj=%u\n",data->tristar_no,data->neigh_no,jj);
8f6a69 51         jt=vtx->tristar[jj-1];
f74313 52         x1=vtx_distance_sq(vtx,jp); //shouldn't be zero!
SP 53         x2=vtx_distance_sq(j,jp); // shouldn't be zero!
8f6a69 54         x3=(j->x-jp->x)*(vtx->x-jp->x)+
SP 55            (j->y-jp->y)*(vtx->y-jp->y)+
56            (j->z-jp->z)*(vtx->z-jp->z);
d7639a 57         
SP 58 #ifdef TS_DOUBLE_DOUBLE
59         ctp=x3/sqrt(x1*x2-x3*x3);
60 #endif
61 #ifdef TS_DOUBLE_FLOAT
62         ctp=x3/sqrtf(x1*x2-x3*x3);
63 #endif
64 #ifdef TS_DOUBLE_LONGDOUBLE
65         ctp=x3/sqrtl(x1*x2-x3*x3);
66 #endif
f74313 67         x1=vtx_distance_sq(vtx,jm);
SP 68         x2=vtx_distance_sq(j,jm);
8f6a69 69         x3=(j->x-jm->x)*(vtx->x-jm->x)+
SP 70            (j->y-jm->y)*(vtx->y-jm->y)+
71            (j->z-jm->z)*(vtx->z-jm->z);
d7639a 72 #ifdef TS_DOUBLE_DOUBLE
SP 73         ctm=x3/sqrt(x1*x2-x3*x3);
74 #endif
75 #ifdef TS_DOUBLE_FLOAT
76         ctm=x3/sqrtf(x1*x2-x3*x3);
77 #endif
78 #ifdef TS_DOUBLE_LONGDOUBLE
79         ctm=x3/sqrtl(x1*x2-x3*x3);
80 #endif
81         tot=ctp+ctm;
82         tot=0.5*tot;
a63f17 83
f74313 84         xlen=vtx_distance_sq(j,vtx);
a63f17 85 /*
d7639a 86 #ifdef  TS_DOUBLE_DOUBLE 
8f6a69 87         vtx->bond[jj-1]->bond_length=sqrt(xlen); 
d7639a 88 #endif
SP 89 #ifdef  TS_DOUBLE_FLOAT
8f6a69 90         vtx->bond[jj-1]->bond_length=sqrtf(xlen); 
d7639a 91 #endif
SP 92 #ifdef  TS_DOUBLE_LONGDOUBLE 
8f6a69 93         vtx->bond[jj-1]->bond_length=sqrtl(xlen); 
d7639a 94 #endif
SP 95
8f6a69 96         vtx->bond[jj-1]->bond_length_dual=tot*vtx->bond[jj-1]->bond_length;
a63f17 97 */
d7639a 98         s+=tot*xlen;
8f6a69 99         xh+=tot*(j->x - vtx->x);
SP 100         yh+=tot*(j->y - vtx->y);
101         zh+=tot*(j->z - vtx->z);
41a035 102         txn+=jt->xnorm;
SP 103         tyn+=jt->ynorm;
104         tzn+=jt->znorm;
d7639a 105     }
SP 106     
107     h=xh*xh+yh*yh+zh*zh;
108     ht=txn*xh+tyn*yh + tzn*zh;
109     s=s/4.0; 
110 #ifdef TS_DOUBLE_DOUBLE
111     if(ht>=0.0) {
8f6a69 112         vtx->curvature=sqrt(h);
d7639a 113     } else {
8f6a69 114         vtx->curvature=-sqrt(h);
d7639a 115     }
SP 116 #endif
117 #ifdef TS_DOUBLE_FLOAT
118     if(ht>=0.0) {
8f6a69 119         vtx->curvature=sqrtf(h);
d7639a 120     } else {
8f6a69 121         vtx->curvature=-sqrtf(h);
d7639a 122     }
SP 123 #endif
124 #ifdef TS_DOUBLE_LONGDOUBLE
125     if(ht>=0.0) {
8f6a69 126         vtx->curvature=sqrtl(h);
d7639a 127     } else {
8f6a69 128         vtx->curvature=-sqrtl(h);
d7639a 129     }
SP 130 #endif
8f6a69 131 // What is vtx->c?????????????? Here it is 0!
dac2e5 132 // c is forced curvature energy for each vertex. Should be set to zero for
8f6a69 133 // normal circumstances.
SP 134     vtx->energy=0.5*s*(vtx->curvature/s-vtx->c)*(vtx->curvature/s-vtx->c);
d7639a 135
SP 136     return TS_SUCCESS;
137 }