From 0ee39cda7ef53b70ed466d0277d22913eddee689 Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@fe.uni-lj.si>
Date: Wed, 23 Apr 2014 17:20:36 +0000
Subject: [PATCH] in debuggingg ptrocedure
---
src/timestep.c | 8 ++++----
src/shcomplex.c | 34 ++++++++++++++++++++++------------
src/shcomplex.h | 2 +-
3 files changed, 27 insertions(+), 17 deletions(-)
diff --git a/src/shcomplex.c b/src/shcomplex.c
index 9350fe5..1c97600 100644
--- a/src/shcomplex.c
+++ b/src/shcomplex.c
@@ -130,31 +130,41 @@
}
-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)];
+ }
+ 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));
+// 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;
}
diff --git a/src/shcomplex.h b/src/shcomplex.h
index 8968b31..72f8030 100644
--- a/src/shcomplex.h
+++ b/src/shcomplex.h
@@ -5,5 +5,5 @@
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
diff --git a/src/timestep.c b/src/timestep.c
index 415be74..d5fd583 100644
--- a/src/timestep.c
+++ b/src/timestep.c
@@ -15,7 +15,7 @@
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;
@@ -24,7 +24,7 @@
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);
@@ -56,7 +56,7 @@
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++){
@@ -71,7 +71,7 @@
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);
--
Gitblit v1.9.3