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/energy.c | 22 +++++++++++-----------
1 files changed, 11 insertions(+), 11 deletions(-)
diff --git a/src/energy.c b/src/energy.c
index efe5911..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==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);
--
Gitblit v1.9.3