Trisurf Monte Carlo simulator
efb615c065452b58799e32735647b0559698163f..3197854e89eb9ffe5ff3137e274d16806c2028a5
2021-07-27 Samo Penic
Found one bug in the code when calculating mprod. Fixed.
319785 diff | tree
2021-04-19 Samo Penic
Work in progress
51f7fd diff | tree
2 files modified
30 ■■■■ changed files
src/energy.c 22 ●●●● patch | view | raw | blame | history
src/tape 8 ●●●● patch | view | raw | blame | history
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==0){
    printf("Eigenvalues: %f, %f, %f\n", eigenval[0], eigenval[1], eigenval[2] );
//    exit(0);
    }
    vtx->energy=(pow(eigenval[0]+eigenval[1],2))*Av;
/*    if(vtx->idx==1){
    printf("Eigenvalues: %f, %f, %f\n", eigenval[0], eigenval[1], eigenval[2] );
    exit(0);
    }
*/
    vtx->energy=0.5*(eigenval[0]+eigenval[1])*(eigenval[0]+eigenval[1])*Av;
    gsl_matrix_free(gsl_Sv);
    gsl_vector_free(Sv_eigen);
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)