| | |
| | | } |
| | | |
| | | |
| | | ts_bool calculateKc(ts_vesicle *vesicle){ |
| | | ts_double calculateKc(ts_vesicle *vesicle){ |
| | | |
| | | ts_int i; |
| | | gsl_matrix *A=gsl_matrix_alloc(vesicle->sphHarmonics->l,2); |
| | | ts_int i,j; |
| | | ts_double retval, bval; |
| | | gsl_matrix *A=gsl_matrix_alloc(vesicle->sphHarmonics->l-1,2); |
| | | gsl_vector *tau=gsl_vector_alloc(2); |
| | | gsl_vector *b=gsl_vector_alloc(vesicle->sphHarmonics->l); |
| | | gsl_vector *b=gsl_vector_alloc(vesicle->sphHarmonics->l-1); |
| | | gsl_vector *x=gsl_vector_alloc(2); |
| | | gsl_vector *res=gsl_vector_alloc(vesicle->sphHarmonics->l); |
| | | gsl_vector *res=gsl_vector_alloc(vesicle->sphHarmonics->l-1); |
| | | |
| | | //solving (A^T*A)*x=A^T*UlmSqAvg |
| | | //solving (A^T*A)*x=A^T*b |
| | | //fill the data for matrix A and vector b |
| | | for(i=1;i<=vesicle->sphHarmonics->l;i++){ |
| | | for(i=1;i<vesicle->sphHarmonics->l;i++){ |
| | | gsl_matrix_set(A, i-1,0,(ts_double)((i-1)*(i+2))); |
| | | gsl_matrix_set(A, i-1,1,(ts_double)((i-1)*(i+2)*(i+1)*i)); |
| | | gsl_vector_set(b,i-1,(ts_double)vesicle->sphHarmonics->N/vesicle->sphHarmonics->sumUlm2[i-1][(i-1)*2]); |
| | | |
| | | bval=0; |
| | | //average for m from 0..l (only positive m's) |
| | | for(j=0;j<i;j++){ |
| | | bval+=vesicle->sphHarmonics->sumUlm2[i][(j+i)]; |
| | | } |
| | | fprintf(stderr,"b[2]=%e\n",gsl_vector_get(b,1)); |
| | | bval=bval/(ts_double)vesicle->sphHarmonics->N/(ts_double)i; |
| | | |
| | | gsl_vector_set(b,i-1,1.0/bval); |
| | | } |
| | | // fprintf(stderr,"b[2]=%e\n",gsl_vector_get(b,1)); |
| | | gsl_linalg_QR_decomp(A,tau); |
| | | gsl_linalg_QR_lssolve(A,tau,b,x,res); |
| | | fprintf(stderr,"kc=%e\n",gsl_vector_get(x,1)); |
| | | |
| | | // fprintf(stderr,"kc=%e\n",gsl_vector_get(x,1)); |
| | | retval=gsl_vector_get(x,1); |
| | | gsl_matrix_free(A); |
| | | gsl_vector_free(tau); |
| | | gsl_vector_free(b); |
| | | gsl_vector_free(x); |
| | | gsl_vector_free(res); |
| | | return TS_SUCCESS; |
| | | |
| | | return retval; |
| | | } |
| | |
| | | ts_spharm *complex_sph_init(ts_vertex_list *vlist, ts_uint l); |
| | | ts_bool complex_sph_free(ts_spharm *sph); |
| | | ts_bool calculateUlmComplex(ts_vesicle *vesicle); |
| | | ts_bool calculateKc(ts_vesicle *vesicle); |
| | | ts_double calculateKc(ts_vesicle *vesicle); |
| | | #endif |
| | |
| | | |
| | | ts_bool run_simulation(ts_vesicle *vesicle, ts_uint mcsweeps, ts_uint inititer, ts_uint iterations, ts_uint start_iteration){ |
| | | ts_uint i, j,k; |
| | | ts_double r0; |
| | | ts_double r0,kc; |
| | | ts_double l1,l2,l3,volume=0.0,area=0.0,vmsr,bfsr, vmsrt, bfsrt; |
| | | ts_ulong epochtime; |
| | | FILE *fd1; |
| | |
| | | if(fd==NULL){ |
| | | fatal("Cannot open statistics.csv file for writing",1); |
| | | } |
| | | fprintf(fd, "Epoch OuterLoop VertexMoveSucessRate BondFlipSuccessRate Volume Area lamdba1 lambda2 lambda3\n"); |
| | | fprintf(fd, "Epoch OuterLoop VertexMoveSucessRate BondFlipSuccessRate Volume Area lamdba1 lambda2 lambda3 Kc\n"); |
| | | centermass(vesicle); |
| | | cell_occupation(vesicle); |
| | | if(start_iteration<inititer) ts_fprintf(stdout, "Starting simulation (first %d x %d MC sweeps will not be recorded on disk)\n", inititer, mcsweeps); |
| | |
| | | calculateUlmComplex(vesicle); |
| | | storeUlmComplex2(vesicle); |
| | | saveAvgUlm2(vesicle); |
| | | calculateKc(vesicle); |
| | | kc=calculateKc(vesicle); |
| | | fd1=fopen("state.dat","w"); |
| | | fprintf(fd1,"%e %e\n",vesicle->volume, getR0(vesicle)); |
| | | for(k=0;k<vesicle->vlist->n;k++){ |
| | |
| | | fclose(fd1); |
| | | } |
| | | |
| | | fprintf(fd, "%lu %u %e %e %1.16e %1.16e %1.16e %1.16e %1.16e\n",epochtime,i,vmsr,bfsr,volume, area,l1,l2,l3); |
| | | fprintf(fd, "%lu %u %e %e %1.16e %1.16e %1.16e %1.16e %1.16e %1.16e\n",epochtime,i,vmsr,bfsr,volume, area,l1,l2,l3,kc); |
| | | fflush(fd); |
| | | // sprintf(filename,"timestep-%05d.pov",i-inititer); |
| | | // write_pov_file(vesicle,filename); |