Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2014-12-16 fda1ab6babed79842534b3a21a6ee96bc26f9d93
commit | author | age
459ff9 1 #include<math.h>
SP 2 #include<stdlib.h>
3 #include<gsl/gsl_complex.h>
4 #include<gsl/gsl_complex_math.h>
5 #include<gsl/gsl_sf_legendre.h>
f06f5f 6
SP 7 #include<gsl/gsl_matrix.h>
8 #include<gsl/gsl_vector.h>
9 #include<gsl/gsl_linalg.h>
459ff9 10 #include "general.h"
SP 11 #include "sh.h"
fda1ab 12 #include "shreal.h"
459ff9 13 #include "shcomplex.h"
fda1ab 14 #include "string.h"
SP 15 #include "io.h"
459ff9 16
SP 17 ts_spharm *complex_sph_init(ts_vertex_list *vlist, ts_uint l){
18     ts_uint j,i;
19     ts_spharm *sph=(ts_spharm *)malloc(sizeof(ts_spharm));
20
21     sph->N=0;
22     /* lets initialize Ylm for each vertex. */
23     sph->Ylmi=(ts_double ***)calloc(l,sizeof(ts_double **));
24     for(i=0;i<l;i++){
25             sph->Ylmi[i]=(ts_double **)calloc(2*i+1,sizeof(ts_double *));
26             for(j=0;j<(2*i+1);j++){
27                 sph->Ylmi[i][j]=(ts_double *)calloc(vlist->n,sizeof(ts_double));
28             }
29     }
30         
31     /* lets initialize ulm */
32     sph->ulm=(ts_double **)calloc(l,sizeof(ts_double *));
33     sph->ulmComplex=(gsl_complex **)calloc(l,sizeof(gsl_complex *));
34     for(j=0;j<l;j++){
35         sph->ulm[j]=(ts_double *)calloc(2*j+1,sizeof(ts_double));
36         sph->ulmComplex[j]=(gsl_complex *)calloc(2*j+1,sizeof(gsl_complex));
37     }
38
39     /* lets initialize sum of Ulm2 */
40     sph->sumUlm2=(ts_double **)calloc(l,sizeof(ts_double *));
41     for(j=0;j<l;j++){
42         sph->sumUlm2[j]=(ts_double *)calloc(2*j+1,sizeof(ts_double));
43     }
44
45     /* lets initialize co */
46 //NOTE: C is has zero based indexing. Code is imported from fortran and to comply with original indexes we actually generate one index more. Also second dimension is 2*j+2 instead of 2*j+2. elements starting with 0 are useles and should be ignored!
47     sph->co=(ts_double **)calloc(l+1,sizeof(ts_double *));
48     for(j=0;j<=l;j++){
49         sph->co[j]=(ts_double *)calloc(2*j+2,sizeof(ts_double));
50     }
51
52     sph->l=l;   
53
54     /* Calculate coefficients that will remain constant during all the simulation */ 
55    precomputeShCoeff(sph);
fda1ab 56    shreal_init(sph,l); //this is for TESTING only! 
459ff9 57     return sph;
SP 58 }
59
60 ts_bool complex_sph_free(ts_spharm *sph){
61     int i,j;
62     if(sph==NULL) return TS_FAIL;
fda1ab 63     shreal_free(sph); //this is for TESTING only!
459ff9 64     for(i=0;i<sph->l;i++){
SP 65         if(sph->ulm[i]!=NULL) free(sph->ulm[i]);
66         if(sph->ulmComplex[i]!=NULL) free(sph->ulmComplex[i]);
67         if(sph->sumUlm2[i]!=NULL) free(sph->sumUlm2[i]);
68         if(sph->co[i]!=NULL) free(sph->co[i]);
69     }
70         if(sph->co[sph->l]!=NULL) free(sph->co[sph->l]);
71     if(sph->co != NULL) free(sph->co);
72     if(sph->ulm !=NULL) free(sph->ulm);
73     if(sph->ulmComplex !=NULL) free(sph->ulmComplex);
74
75         if(sph->Ylmi!=NULL) {
76             for(i=0;i<sph->l;i++){
77                 if(sph->Ylmi[i]!=NULL){
78                     for(j=0;j<i*2+1;j++){
79                         if(sph->Ylmi[i][j]!=NULL) free (sph->Ylmi[i][j]);
80                     }
81                     free(sph->Ylmi[i]);
82                 }
83             }
84             free(sph->Ylmi);
85         }
86
87     free(sph);
88     return TS_SUCCESS;
89 }
90
91
92 ts_bool calculateUlmComplex(ts_vesicle *vesicle){
93     ts_int i,j,k,m,l;
94     ts_vertex *cvtx;
95     ts_coord coord;
96 /* set all values to zero */
97     for(i=0;i<vesicle->sphHarmonics->l;i++){
98         for(j=0;j<2*i+1;j++) GSL_SET_COMPLEX(&(vesicle->sphHarmonics->ulmComplex[i][j]),0.0,0.0);
99     }
100
101     for(k=0;k<vesicle->vlist->n; k++){
102         cvtx=vesicle->vlist->vtx[k];
103     cart2sph(&coord,cvtx->x,cvtx->y,cvtx->z);
104         for(i=0;i<vesicle->sphHarmonics->l;i++){
105             for(j=0;j<2*i+1;j++){
106         m=j-i;
107         l=i;
108         if(m>=0){    
109     //    fprintf(stderr, "Racunam za l=%d, m=%d\n", l,m);
110                 vesicle->sphHarmonics->ulmComplex[i][j]=gsl_complex_add(vesicle->sphHarmonics->ulmComplex[i][j], gsl_complex_conjugate(gsl_complex_mul_real(gsl_complex_polar(1.0,(ts_double)m*coord.e2),cvtx->solAngle*cvtx->relR*gsl_sf_legendre_sphPlm(l,m,cos(coord.e3)))) );
111         } else {
112     //    fprintf(stderr, "Racunam za l=%d, abs(m=%d)\n", l,m);
113                 vesicle->sphHarmonics->ulmComplex[i][j]=gsl_complex_add(vesicle->sphHarmonics->ulmComplex[i][j], gsl_complex_conjugate(gsl_complex_mul_real(gsl_complex_polar(1.0,(ts_double)m*coord.e2),cvtx->solAngle*cvtx->relR*pow(-1,m)*gsl_sf_legendre_sphPlm(l,-m,cos(coord.e3)))) );
114
115         }
116             }
117         }
118     }
119     return TS_SUCCESS;
120 }
121
122 ts_bool storeUlmComplex2(ts_vesicle *vesicle){
123
124     ts_spharm *sph=vesicle->sphHarmonics;
125     ts_int i,j;
126     for(i=0;i<sph->l;i++){
127             for(j=0;j<2*i+1;j++){
128                 sph->sumUlm2[i][j]+=gsl_complex_abs2(sph->ulmComplex[i][j]);
129             }
130     }
131     sph->N++;
132     return TS_SUCCESS;
133 }
134
f06f5f 135
22cdfd 136 ts_double calculateKc(ts_vesicle *vesicle, ts_int lmin, ts_int lmax){
SP 137     ts_int min=lmin;
138     ts_int max=lmax; //vesicle->sphHarmonics->l-3;
9f2ad6 139     ts_long i,j;
0ee39c 140     ts_double retval, bval;
9f2ad6 141     gsl_matrix *A=gsl_matrix_alloc(max-min,2);
f06f5f 142     gsl_vector *tau=gsl_vector_alloc(2);
9f2ad6 143     gsl_vector *b=gsl_vector_alloc(max-min);
f06f5f 144     gsl_vector *x=gsl_vector_alloc(2);
9f2ad6 145     gsl_vector *res=gsl_vector_alloc(max-min);
f06f5f 146
0ee39c 147     //solving (A^T*A)*x=A^T*b
f06f5f 148     //fill the data for matrix A and vector b
9f2ad6 149     for(i=min;i<max;i++){
SP 150             gsl_matrix_set(A, i-min,0,(ts_double)((i-1)*(i+2)));
151             gsl_matrix_set(A, i-min,1,(ts_double)((i-1)*(i+2)*(i+1)*i));
22cdfd 152 //            fprintf(stderr,"%e %e\n", gsl_matrix_get(A,i-min,0), gsl_matrix_get(A,i-min,1));
9f2ad6 153             bval=0.0;
0ee39c 154             //average for m from 0..l (only positive m's)
9f2ad6 155             for(j=0;j<=i;j++){
0ee39c 156                 bval+=vesicle->sphHarmonics->sumUlm2[i][(j+i)];
SP 157             }
9f2ad6 158                 bval=bval/(ts_double)vesicle->sphHarmonics->N/(ts_double)(i+1);
0ee39c 159
9f2ad6 160             gsl_vector_set(b,i-min,1.0/bval);
22cdfd 161 //            fprintf(stderr,"%e\n", 1.0/gsl_vector_get(b,i-min));
f06f5f 162     }
0ee39c 163 //    fprintf(stderr,"b[2]=%e\n",gsl_vector_get(b,1));
f06f5f 164     gsl_linalg_QR_decomp(A,tau);
SP 165     gsl_linalg_QR_lssolve(A,tau,b,x,res);
0ee39c 166 //    fprintf(stderr,"kc=%e\n",gsl_vector_get(x,1));
SP 167     retval=gsl_vector_get(x,1);
f06f5f 168     gsl_matrix_free(A);
SP 169     gsl_vector_free(tau);
170     gsl_vector_free(b);
171     gsl_vector_free(x);
172     gsl_vector_free(res);
0ee39c 173     
SP 174     return retval;
f06f5f 175 }
fda1ab 176
SP 177
178 ts_bool saveAvgUlm2Complex(ts_vesicle *vesicle){
179
180     FILE *fh;
181     char filename[10000];
182     strcpy(filename, command_line_args.path);
183     strcat(filename, "sph2out.dat");
184     fh=fopen(filename, "w");
185     if(fh==NULL){
186         err("Cannot open file %s for writing");
187         return TS_FAIL;
188     }
189
190     ts_spharm *sph=vesicle->sphHarmonics;
191     ts_int i,j;
192     fprintf(fh,"l,\tm,\tulm^2avg\n");
193     for(i=0;i<sph->l;i++){
194             for(j=0;j<2*i+1;j++){
195         fprintf(fh,"%d,\t%d,\t%e\n", i, j-i, sph->sumUlm2[i][j]/(ts_double)sph->N);
196
197             }
198     fprintf(fh,"\n");
199     }
200     fclose(fh);
201     return TS_SUCCESS;
202 }