Trisurf Monte Carlo simulator
Miha
2012-07-13 60bc6aa6b61edad39b664939667a1616711d580c
commit | author | age
262607 1 #include<stdio.h>
SP 2 #include<math.h>
3 #include "general.h"
4 #include "vertex.h"
5 #include "bond.h"
6 #include "triangle.h"
7 #include "cell.h"
8 #include "vesicle.h"
9 #include "io.h"
10 #include "initial_distribution.h"
11 #include "frame.h"
12 #include "timestep.h"
13 #include "sh.h"
14
15 /** Entrance function to the program
16   * @param argv is a number of parameters used in program call (including the program name
17   * @param argc is a pointer to strings (character arrays) which holds the arguments
18   * @returns returns 0 on success, any other number on fail.
19 */
20 ts_bool saveAvgUlm2(ts_vesicle *vesicle);
21 int main(int argv, char *argc[]){
a63f17 22 ts_uint i,j,k;
262607 23 ts_vesicle *vesicle;
SP 24 ts_double r0;
25 vesicle=initial_distribution_dipyramid(17,60,60,60,0.15);
26 //parsetape(vesicle,&i);
27
59c7f5 28 //similar to nmax in fortran code
M 29 ts_uint nmax;
30
262607 31 //these four must come from parsetype!
SP 32 vesicle->dmax=1.67*1.67;
33 vesicle->stepsize=0.15;
34 vesicle->clist->max_occupancy=8;
35 vesicle->bending_rigidity=25.0;
36 //fprintf(stderr,"xk=%f",vesicle->bending_rigidity);
37
919555 38     centermass(vesicle);
a63f17 39 cell_occupation(vesicle);
SP 40
41 //test if the structure is internally organized into cells correctly 
42 ts_uint cind;
43 for(i=0;i<vesicle->vlist->n;i++){
44     cind=vertex_self_avoidance(vesicle, vesicle->vlist->vtx[i]);
45
46     if(vesicle->clist->cell[cind]==vesicle->vlist->vtx[i]->cell){
47         //fprintf(stdout,"(T) Idx match!\n");
48     } else {
49         fprintf(stderr,"(T) ***** Idx don't match!\n");
50
51     }
52 }
53 //end test
262607 54 vesicle->sphHarmonics=sph_init(vesicle->vlist, 21);
SP 55
56 vesicle_volume(vesicle);
57 r0=getR0(vesicle);
58
59 preparationSh(vesicle,r0);
60 calculateYlmi(vesicle);
61 calculateUlm(vesicle);
62
60bc6a 63 //preloop:
M 64 for(i=0;i<1000;i++){
65     cell_occupation(vesicle);
66     for(j=0;j<1000;j++){
67         single_timestep(vesicle);
68     }    
69     centermass(vesicle);
70     fprintf(stderr, "Preloop %d completed.\n",i+1);
71 }
72
59c7f5 73 nmax=1000;
M 74 for(i=0;i<nmax;i++){
672ae4 75     for(j=0;j<200;j++){
a63f17 76         cell_occupation(vesicle);
SP 77         for(k=0;k<5;k++){
262607 78         single_timestep(vesicle);
a63f17 79         }
SP 80         centermass(vesicle);
262607 81     }    
919555 82     vesicle_volume(vesicle);
SP 83     r0=getR0(vesicle);
262607 84
919555 85     preparationSh(vesicle,r0);
SP 86     calculateYlmi(vesicle);
87     calculateUlm(vesicle);
262607 88
919555 89     storeUlm2(vesicle);
SP 90     saveAvgUlm2(vesicle);
262607 91
SP 92     write_vertex_xml_file(vesicle,i);
59c7f5 93     fprintf(stderr, "Loop %d out of %d completed.\n",i+1,nmax);
a63f17 94
262607 95 }
a63f17 96
262607 97 write_master_xml_file("test.pvd");
SP 98 write_dout_fcompat_file(vesicle,"dout");
99 vesicle_free(vesicle);
100
101 return 0; //program finished perfectly ok. We return 0.
102 }
103
104
105
106 ts_bool saveAvgUlm2(ts_vesicle *vesicle){
107
108     FILE *fh;
109     
110     fh=fopen("sph2out.dat", "w");
111     if(fh==NULL){
112         err("Cannot open file %s for writing");
113         return TS_FAIL;
114     }
115
116     ts_spharm *sph=vesicle->sphHarmonics;
117     ts_int i,j;
118     fprintf(fh,"l,\tm,\tulm^2avg\n");
119     for(i=0;i<sph->l;i++){
120             for(j=0;j<2*i+1;j++){
121         fprintf(fh,"%d,\t%d,\t%e\n", i, j-i, sph->sumUlm2[i][j]/(ts_double)sph->N);
122
123             }
443ba1 124     fprintf(fh,"\n");
262607 125     }
SP 126     fclose(fh);
127     return TS_SUCCESS;
128 }