From 3197854e89eb9ffe5ff3137e274d16806c2028a5 Mon Sep 17 00:00:00 2001 From: Samo Penic <samo.penic@gmail.com> Date: Tue, 27 Jul 2021 12:51:09 +0000 Subject: [PATCH] Found one bug in the code when calculating mprod. Fixed. --- src/tape | 8 ++++---- src/energy.c | 18 +++++++++--------- 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/src/energy.c b/src/energy.c index 170c971..85a4fb4 100644 --- a/src/energy.c +++ b/src/energy.c @@ -166,9 +166,9 @@ } */ for(jj=0;jj<vtx->neigh_no;jj++){ - edge_vector_x[jj]=vtx->neigh[jj]->x-vtx->x; - edge_vector_y[jj]=vtx->neigh[jj]->y-vtx->y; - edge_vector_z[jj]=vtx->neigh[jj]->z-vtx->z; + edge_vector_x[jj]=(vtx->neigh[jj]->x-vtx->x); + edge_vector_y[jj]=(vtx->neigh[jj]->y-vtx->y); + edge_vector_z[jj]=(vtx->neigh[jj]->z-vtx->z); //Here we calculate normalized edge vector @@ -239,10 +239,9 @@ edge_binormal_z[jj]=(edge_normal_x[jj]*edge_vector_y[jj])-(edge_normal_y[jj]*edge_vector_x[jj]); - //mprod[jj]=it->x*(k->y*edge_vector_z[jj]-edge_vector_y[jj]*k->z)-it->y*(k->x*edge_vector_z[jj]-k->z*edge_vector_x[jj])+it->z*(k->x*edge_vector_y[jj]-k->y*edge_vector_x[jj]); - mprod[jj]=lm->xnorm*(lp->ynorm*edge_vector_z[jj]-lp->znorm*edge_vector_y[jj]) - lm->ynorm*(lp->xnorm*edge_vector_z[jj]-lp->znorm*edge_vector_z[jj])+ lm->znorm*(lp->xnorm*edge_vector_y[jj]-lp->ynorm*edge_vector_x[jj]); + mprod[jj]=lm->xnorm*(lp->ynorm*edge_vector_z[jj]-lp->znorm*edge_vector_y[jj]) - lm->ynorm*(lp->xnorm*edge_vector_z[jj]-lp->znorm*edge_vector_x[jj])+ lm->znorm*(lp->xnorm*edge_vector_y[jj]-lp->ynorm*edge_vector_x[jj]); - phi[jj]=copysign(acos(lm->xnorm*lp->xnorm+lm->ynorm*lp->ynorm+lm->znorm*lp->znorm-1e-10),-mprod[jj])+M_PI; + phi[jj]=-copysign(acos(lm->xnorm*lp->xnorm+lm->ynorm*lp->ynorm+lm->znorm*lp->znorm),mprod[jj])+M_PI; /* if(vtx->idx==0){ printf("Angle PHI vertex %d (angle %d): %f\n",vtx->idx,jj,phi[jj]); } @@ -307,12 +306,13 @@ eigenval[2]= gsl_vector_get(Sv_eigen, 2); qsort(eigenval, 3, sizeof(ts_double), cmpfunc); - if(vtx->idx==1){ + +/* if(vtx->idx==1){ printf("Eigenvalues: %f, %f, %f\n", eigenval[0], eigenval[1], eigenval[2] ); exit(0); } - - vtx->energy=(pow(eigenval[0]+eigenval[1],2))*Av; +*/ + vtx->energy=0.5*(eigenval[0]+eigenval[1])*(eigenval[0]+eigenval[1])*Av; gsl_matrix_free(gsl_Sv); gsl_vector_free(Sv_eigen); diff --git a/src/tape b/src/tape index 41bd4da..8de15ff 100644 --- a/src/tape +++ b/src/tape @@ -53,9 +53,9 @@ R_nucleusY=0 R_nucleusZ=0 ####### Cell definitions ############ -nxmax=100 -nymax=100 -nzmax=100 +nxmax=150 +nymax=150 +nzmax=150 ####### Program Control ############ @@ -92,7 +92,7 @@ #number of vertices with spontaneous curvature (integer) -number_of_vertices_with_c0=100 +number_of_vertices_with_c0=0 #c0/2 is spontaneous curvature. c0 is used as (c1+c1-c0)^2 in energy term (float) c0=0.5 #energy of attraction of vertices with spontaneous curvature (float, positive value for attraction) -- Gitblit v1.9.3