Trisurf Monte Carlo simulator
Samo Penic
2012-07-13 a63f1719d2c7fd2c69accc0eb3eb038af50e555e
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
28 //these four must come from parsetype!
29 vesicle->dmax=1.67*1.67;
30 vesicle->stepsize=0.15;
31 vesicle->clist->max_occupancy=8;
32 vesicle->bending_rigidity=25.0;
33 //fprintf(stderr,"xk=%f",vesicle->bending_rigidity);
34
919555 35     centermass(vesicle);
a63f17 36 cell_occupation(vesicle);
SP 37
38 //test if the structure is internally organized into cells correctly 
39 ts_uint cind;
40 for(i=0;i<vesicle->vlist->n;i++){
41     cind=vertex_self_avoidance(vesicle, vesicle->vlist->vtx[i]);
42
43     if(vesicle->clist->cell[cind]==vesicle->vlist->vtx[i]->cell){
44         //fprintf(stdout,"(T) Idx match!\n");
45     } else {
46         fprintf(stderr,"(T) ***** Idx don't match!\n");
47
48     }
49 }
50 //end test
262607 51 vesicle->sphHarmonics=sph_init(vesicle->vlist, 21);
SP 52
53 vesicle_volume(vesicle);
54 r0=getR0(vesicle);
55
56 preparationSh(vesicle,r0);
57 calculateYlmi(vesicle);
58 calculateUlm(vesicle);
59
60
1ad6d1 61 for(i=0;i<1;i++){
a63f17 62     for(j=0;j<20;j++){
SP 63         cell_occupation(vesicle);
64         for(k=0;k<5;k++){
262607 65         single_timestep(vesicle);
a63f17 66         }
SP 67         centermass(vesicle);
262607 68     }    
919555 69     vesicle_volume(vesicle);
SP 70     r0=getR0(vesicle);
262607 71
919555 72     preparationSh(vesicle,r0);
SP 73     calculateYlmi(vesicle);
74     calculateUlm(vesicle);
262607 75
919555 76     storeUlm2(vesicle);
SP 77     saveAvgUlm2(vesicle);
262607 78
SP 79     write_vertex_xml_file(vesicle,i);
80     fprintf(stderr, "Loop %d completed.\n",i+1);
a63f17 81
262607 82 }
a63f17 83
262607 84 write_master_xml_file("test.pvd");
SP 85 write_dout_fcompat_file(vesicle,"dout");
86 vesicle_free(vesicle);
87
88 return 0; //program finished perfectly ok. We return 0.
89 }
90
91
92
93 ts_bool saveAvgUlm2(ts_vesicle *vesicle){
94
95     FILE *fh;
96     
97     fh=fopen("sph2out.dat", "w");
98     if(fh==NULL){
99         err("Cannot open file %s for writing");
100         return TS_FAIL;
101     }
102
103     ts_spharm *sph=vesicle->sphHarmonics;
104     ts_int i,j;
105     fprintf(fh,"l,\tm,\tulm^2avg\n");
106     for(i=0;i<sph->l;i++){
107             for(j=0;j<2*i+1;j++){
108         fprintf(fh,"%d,\t%d,\t%e\n", i, j-i, sph->sumUlm2[i][j]/(ts_double)sph->N);
109
110             }
443ba1 111     fprintf(fh,"\n");
262607 112     }
SP 113     fclose(fh);
114     return TS_SUCCESS;
115 }