Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2012-06-07 523bf18206f550a315c6c17e5a0a253381b0f8bf
commit | author | age
88f451 1 #include<math.h>
SP 2 #include<stdlib.h>
3 #include "general.h"
4 #include "sh.h"
5
6 /* Gives you legendre polynomials. Taken from NR, p. 254 */
7 ts_double plgndr(ts_int l, ts_int m, ts_float x){
8     ts_double fact, pll, pmm, pmmp1, somx2;
9     ts_int i,ll;
10
11 #ifdef TS_DOUBLE_DOUBLE
12     if(m<0 || m>l || fabs(x)>1.0)
13         fatal("Bad arguments in routine plgndr",1);
14 #endif
15 #ifdef TS_DOUBLE_FLOAT
16     if(m<0 || m>l || fabsf(x)>1.0)
17         fatal("Bad arguments in routine plgndr",1);
18 #endif
19 #ifdef TS_DOUBLE_LONGDOUBLE
20     if(m<0 || m>l || fabsl(x)>1.0)
21         fatal("Bad arguments in routine plgndr",1);
22 #endif
23     pmm=1.0;
24     if (m>0) {
25 #ifdef TS_DOUBLE_DOUBLE
26         somx2=sqrt((1.0-x)*(1.0+x));
27 #endif
28 #ifdef TS_DOUBLE_FLOAT
29         somx2=sqrtf((1.0-x)*(1.0+x));
30 #endif
31 #ifdef TS_DOUBLE_LONGDOUBLE
32         somx2=sqrtl((1.0-x)*(1.0+x));
33 #endif
34         fact=1.0;
35         for (i=1; i<=m;i++){
36             pmm *= -fact*somx2;
37             fact +=2.0;
38         }
39     }
40
41     if (l == m) return pmm;
42     else {
43         pmmp1=x*(2*m+1)*pmm;
44         if(l==(m+1)) return(pmmp1);
45         else {
46             pll=0; /* so it can not be uninitialized */
47             for(ll=m+2;ll<=l;ll++){
48                 pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m);
49                 pmm=pmmp1;
50                 pmmp1=pll;
51             }
52             return(pll);
53         }
54     }
55 }
56
57
523bf1 58
SP 59 ts_bool precomputeShCoeff(ts_spharm *sph){
60     ts_uint i,j;
61     for(i=0;i<sph->l;i++){
62         sph->co[i][i]=sqrt((2.0*i+1.0)/2.0/M_PI);
63         for(j=0;j<i-1;j++){
64             
65         }
66     }
67     return TS_SUCCESS;
68
69 }
70
71
88f451 72 /*Computes Y(l,m,theta,fi) (Miha's definition that is different from common definition for  factor srqt(1/(2*pi)) */
SP 73 ts_double shY(ts_int l,ts_int m,ts_double theta,ts_double fi){
74     ts_double fac1, fac2, K;
75     int i;
76
77     if(l<0 || m>l || m<-l)
78         fatal("Error using shY function!",1);
79
80     fac1=1.0;
af3bad 81     for(i=1; i<=l-abs(m);i++){
88f451 82         fac1 *= i;
SP 83     }
84     fac2=1.0;
af3bad 85     for(i=1; i<=l+abs(m);i++){
88f451 86         fac2 *= i;
SP 87     }
88
89     if(m==0){
90         K=sqrt(1.0/(2.0*M_PI));
91     }
92     else if (m>0) {
93         K=sqrt(1.0/(M_PI))*cos(m*fi);
94     } 
95     else {
96         //K=pow(-1.0,abs(m))*sqrt(1.0/(2.0*M_PI))*cos(m*fi);
97         if(abs(m)%2==0)
af3bad 98         K=sqrt(1.0/(M_PI))*cos(m*fi);
88f451 99         else
af3bad 100         K=-sqrt(1.0/(M_PI))*cos(m*fi);
88f451 101     }
SP 102     
103     return K*sqrt((2.0*l+1.0)/2.0*fac1/fac2)*plgndr(l,abs(m),cos(theta));    
104 }
523bf1 105
SP 106
107 /* Function transforms coordinates from cartesian to spherical coordinates
108  * (r,phi, theta). */
109 ts_bool *cart2sph(ts_coord *coord, ts_double x, ts_double y, ts_double z){
110     coord->coord_type=TS_COORD_SPHERICAL;
111 #ifdef TS_DOUBLE_DOUBLE
112     coord->e1=sqrt(x*x+y*y+z*z);
113     if(z==0) coord->e3=M_PI/2.0;
114     else coord->e3=atan(sqrt(x*x+y*y)/z);
115     coord->e2=atan2(y,x);
116 #endif
117 #ifdef TS_DOUBLE_FLOAT
118     coord->e1=sqrtf(x*x+y*y+z*z);
119     if(z==0) coord->e3=M_PI/2.0;
120     else coord->e3=atanf(sqrtf(x*x+y*y)/z);
121     coord->e2=atan2f(y,x);
122 #endif
123 #ifdef TS_DOUBLE_LONGDOUBLE
124     coord->e1=sqrtl(x*x+y*y+z*z);
125     if(z==0) coord->e3=M_PI/2.0;
126     else coord->e3=atanl(sqrtl(x*x+y*y)/z);
127     coord->e2=atan2l(y,x);
128 #endif
129
130     return TS_SUCCESS;
131 }
132
133 /* Function returns radius of the sphere with the same volume as vesicle (r0) */
134 ts_double getR0(ts_vesicle *vesicle){
135     ts_double r0;
136  #ifdef TS_DOUBLE_DOUBLE
137    r0=pow(vesicle->volume*3.0/4.0/M_PI,1.0/3.0);
138 #endif
139 #ifdef TS_DOUBLE_FLOAT
140    r0=powf(vesicle->volume*3.0/4.0/M_PI,1.0/3.0);
141 #endif
142 #ifdef TS_DOUBLE_LONGDOUBLE
143    r0=powl(vesicle->volume*3.0/4.0/M_PI,1.0/3.0);
144 #endif
145     return r0;
146 }
147
148
149 ts_bool preparationSh(ts_vesicle *vesicle, ts_double r0){
150 //TODO: before calling or during the call calculate area of each triangle! Can
151 //be also done after vertexmove and bondflip //
152     ts_uint i,j;
153     ts_vertex **vtx=vesicle->vlist->vtx;
154     ts_vertex *cvtx;
155     ts_triangle *ctri;
156     ts_double centroid[3];
157     ts_double r;
158     for (i=0;  i<vesicle->vlist->n; i++){
159         cvtx=vtx[i];
160         //cvtx->projArea=4.0*M_PI/1447.0*(cvtx->x*cvtx->x+cvtx->y*cvtx->y+cvtx->z*cvtx->z)/r0/r0;
161         cvtx->projArea=0.0;
162
163         /* go over all triangles that have a common vertex i */
164         for(j=0; j<cvtx->tristar_no; j++){
165             ctri=cvtx->tristar[j];
166             centroid[0]=(ctri->vertex[0]->x + ctri->vertex[1]->x + ctri->vertex[2]->x)/3.0;
167             centroid[1]=(ctri->vertex[0]->y + ctri->vertex[1]->y + ctri->vertex[2]->y)/3.0;
168             centroid[2]=(ctri->vertex[0]->z + ctri->vertex[1]->z + ctri->vertex[2]->z)/3.0;
169         /* calculating projArea+= area(triangle)*cos(theta) */
170 #ifdef TS_DOUBLE_DOUBLE
171             cvtx->projArea = cvtx->projArea + ctri->area*(-centroid[0]*ctri->xnorm - centroid[1]*ctri->ynorm - centroid[2]*ctri->znorm)/ sqrt(centroid[0]*centroid[0]+centroid[1]*centroid[1]+centroid[2]*centroid[2]);
172 #endif
173 #ifdef TS_DOUBLE_FLOAT
174             cvtx->projArea = cvtx->projArea + ctri->area*(-centroid[0]*ctri->xnorm - centroid[1]*ctri->ynorm - centroid[2]*ctri->znorm)/ sqrtf(centroid[0]*centroid[0]+centroid[1]*centroid[1]+centroid[2]*centroid[2]);
175 #endif
176 #ifdef TS_DOUBLE_LONGDOUBLE
177             cvtx->projArea = cvtx->projArea + ctri->area*(-centroid[0]*ctri->xnorm - centroid[1]*ctri->ynorm - centroid[2]*ctri->znorm)/ sqrtl(centroid[0]*centroid[0]+centroid[1]*centroid[1]+centroid[2]*centroid[2]);
178 #endif
179         }
180
181     cvtx->projArea=cvtx->projArea/3.0;
182         //we dont store spherical coordinates of vertex, so we have to calculate
183         //r(i) at this point.
184 #ifdef TS_DOUBLE_DOUBLE
185     r=sqrt(cvtx->x*cvtx->x+cvtx->y*cvtx->y+cvtx->z*cvtx->z);
186 #endif
187 #ifdef TS_DOUBLE_FLOAT
188     r=sqrtf(cvtx->x*cvtx->x+cvtx->y*cvtx->y+cvtx->z*cvtx->z);
189 #endif
190 #ifdef TS_DOUBLE_LONGDOUBLE
191     r=sqrtl(cvtx->x*cvtx->x+cvtx->y*cvtx->y+cvtx->z*cvtx->z);
192 #endif
193     cvtx->relR=(r-r0)/r0;
194     cvtx->solAngle=cvtx->projArea/cvtx->relR * cvtx->projArea/cvtx->relR;
195     }
196     return TS_SUCCESS;
197 }
198
199
200
201 ts_bool calculateYlmi(ts_vesicle *vesicle){
202     ts_uint i,j,k;
203     ts_spharm *sph=vesicle->sphHarmonics;
204     ts_coord *coord=(ts_coord *)malloc(sizeof(ts_coord));
205     ts_double fi, theta;
206     for(k=0;k<vesicle->vlist->n;k++){
207         sph->Ylmi[0][0][k]=sqrt(1.0/4.0/M_PI);
208         cart2sph(coord,vesicle->vlist->vtx[k]->x, vesicle->vlist->vtx[k]->y, vesicle->vlist->vtx[k]->z);
209         fi=coord->e2;
210         theta=coord->e3; 
211         for(i=0; i<sph->l; i++){
212             for(j=0;j<i;j++){
213                 sph->Ylmi[i][j][k]=sph->co[i][j]*cos((j-i-1)*fi)*pow(-1,j-i-1)*plgndr(i,abs(j-i-1),cos(theta));
214             }
215                 sph->Ylmi[i][j+1][k]=sph->co[i][j+1]*plgndr(i,0,cos(theta));
216             for(j=sph->l;j<2*i;j++){
217                 sph->Ylmi[i][j][k]=sph->co[i][j]*sin((j-i-1)*fi)*plgndr(i,j-i-1,cos(theta));
218             }
219         }
220
221     }
222     free(coord);
223     return TS_SUCCESS;
224 }
225
226
227
228 ts_bool calculateUlm(ts_vesicle *vesicle){
229     ts_uint i,j,k;
230     ts_vertex *cvtx;
231     for(i=0;i<vesicle->sphHarmonics->l;i++){
232         for(j=0;j<2*i;j++) vesicle->sphHarmonics->ulm[i][j]=0.0;
233     }
234
235 //TODO: call calculateYlmi !!!
236
237
238     for(k=0;k<vesicle->vlist->n; k++){
239         cvtx=vesicle->vlist->vtx[k];
240         for(i=0;i<vesicle->sphHarmonics->l;i++){
241             for(j=0;j<2*i;j++){
242                 vesicle->sphHarmonics->ulm[i][j]+= cvtx->solAngle*cvtx->relR*vesicle->sphHarmonics->Ylmi[i][j][k];
243             }
244
245         }
246     }
247
248     return TS_SUCCESS;
249 }