From bd826de2f539f2e48c8c01d2d7f9f34c7e97104a Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@gmail.com>
Date: Fri, 13 May 2016 07:43:27 +0000
Subject: [PATCH] Fix in trisurf output, inhibiting print of successful reconstruction. Multiple fixes and improvements in python module. Added symlinking of tapes into the running directories and dumping tapes from snapshots into tape files.
---
src/sh.c | 172 +++++++++++++++++++++++++++++++++++++++++++++++----------
1 files changed, 141 insertions(+), 31 deletions(-)
diff --git a/src/sh.c b/src/sh.c
index 7cbb04f..3f64939 100644
--- a/src/sh.c
+++ b/src/sh.c
@@ -1,19 +1,22 @@
+/* vim: set ts=4 sts=4 sw=4 noet : */
#include<math.h>
#include<stdlib.h>
#include "general.h"
#include "sh.h"
-
+#include "io.h"
+#include <string.h>
ts_spharm *sph_init(ts_vertex_list *vlist, ts_uint l){
ts_uint j,i;
ts_spharm *sph=(ts_spharm *)malloc(sizeof(ts_spharm));
+ sph->N=0;
/* lets initialize Ylm for each vertex. */
sph->Ylmi=(ts_double ***)calloc(l,sizeof(ts_double **));
- for(i=0;i<vlist->n;i++){
- sph->Ylmi[i]=(ts_double **)calloc(2*l+1,sizeof(ts_double *));
- for(j=0;j<l;j++){
+ for(i=0;i<l;i++){
+ sph->Ylmi[i]=(ts_double **)calloc(2*i+1,sizeof(ts_double *));
+ for(j=0;j<(2*i+1);j++){
sph->Ylmi[i][j]=(ts_double *)calloc(vlist->n,sizeof(ts_double));
}
}
@@ -24,33 +27,44 @@
sph->ulm[j]=(ts_double *)calloc(2*j+1,sizeof(ts_double));
}
+ /* lets initialize sum of Ulm2 */
+ sph->sumUlm2=(ts_double **)calloc(l,sizeof(ts_double *));
+ for(j=0;j<l;j++){
+ sph->sumUlm2[j]=(ts_double *)calloc(2*j+1,sizeof(ts_double));
+ }
/* lets initialize co */
- sph->co=(ts_double **)calloc(l,sizeof(ts_double *));
- for(j=0;j<l;j++){
- sph->co[j]=(ts_double *)calloc(2*j+1,sizeof(ts_double));
+//NOTE: C is has zero based indexing. Code is imported from fortran and to comply with original indexes we actually generate one index more. Also second dimension is 2*j+2 instead of 2*j+2. elements starting with 0 are useles and should be ignored!
+ sph->co=(ts_double **)calloc(l+1,sizeof(ts_double *));
+ for(j=0;j<=l;j++){
+ sph->co[j]=(ts_double *)calloc(2*j+2,sizeof(ts_double));
}
-
- /* Calculate coefficients that will remain constant during all the simulation */
- precomputeShCoeff(sph);
+ sph->l=l;
+
+ /* Calculate coefficients that will remain constant during all the simulation */
+ precomputeShCoeff(sph);
+
return sph;
}
ts_bool sph_free(ts_spharm *sph){
int i,j;
+ if(sph==NULL) return TS_FAIL;
for(i=0;i<sph->l;i++){
if(sph->ulm[i]!=NULL) free(sph->ulm[i]);
+ if(sph->sumUlm2[i]!=NULL) free(sph->sumUlm2[i]);
if(sph->co[i]!=NULL) free(sph->co[i]);
}
+ if(sph->co[sph->l]!=NULL) free(sph->co[sph->l]);
if(sph->co != NULL) free(sph->co);
if(sph->ulm !=NULL) free(sph->ulm);
if(sph->Ylmi!=NULL) {
for(i=0;i<sph->l;i++){
if(sph->Ylmi[i]!=NULL){
- for(j=0;j<sph->l*2+1;j++){
+ for(j=0;j<i*2+1;j++){
if(sph->Ylmi[i][j]!=NULL) free (sph->Ylmi[i][j]);
}
free(sph->Ylmi[i]);
@@ -64,7 +78,7 @@
}
/* Gives you legendre polynomials. Taken from NR, p. 254 */
-ts_double plgndr(ts_int l, ts_int m, ts_float x){
+ts_double plgndr(ts_int l, ts_int m, ts_double x){
ts_double fact, pll, pmm, pmmp1, somx2;
ts_int i,ll;
@@ -115,20 +129,22 @@
}
+/** @brief: Precomputes coefficients that are required for spherical harmonics computations.
+*/
ts_bool precomputeShCoeff(ts_spharm *sph){
ts_int i,j,al,am;
ts_double **co=sph->co;
- for(i=0;i<sph->l;i++){
- al=i+1;
+ for(i=1;i<=sph->l;i++){
+ al=i;
sph->co[i][i+1]=sqrt((2.0*al+1.0)/2.0/M_PI);
- for(j=0;j<al;j++){
- am=j+1;
- sph->co[i][i+1+j]=co[i][i+j]*sqrt(1.0/(al-am+1)/(al+am));
+ for(j=1;j<=i-1;j++){
+ am=j;
+ sph->co[i][i+1+j]=co[i][i+j]*sqrt(1.0/(al-am+1.0)/(al+am));
sph->co[i][i+1-j]=co[i][i+1+j];
}
- co[i][2*i]=co[i][2*i]*sqrt(1.0/(2.0*al));
- co[i][0]=co[i][2*i+1];
+ co[i][2*i+1]=co[i][2*i]*sqrt(1.0/(2.0*al));
+ co[i][1]=co[i][2*i+1];
co[i][i+1]=sqrt((2.0*al+1.0)/4.0/M_PI);
}
return TS_SUCCESS;
@@ -136,7 +152,16 @@
}
-/*Computes Y(l,m,theta,fi) (Miha's definition that is different from common definition for factor srqt(1/(2*pi)) */
+/** @brief: Computes Y(l,m,theta,fi)
+ *
+ * Function calculates Y^l_m for vertex with given (\theta, \fi) coordinates in
+ * spherical coordinate system.
+ * @param l is an ts_int argument.
+ * @param m is an ts_int argument.
+ * @param theta is ts_double argument.
+ * @param fi is a ts_double argument.
+ *
+ * (Miha's definition that is different from common definition for factor srqt(1/(2*pi)) */
ts_double shY(ts_int l,ts_int m,ts_double theta,ts_double fi){
ts_double fac1, fac2, K;
int i;
@@ -167,7 +192,7 @@
K=-sqrt(1.0/(M_PI))*cos(m*fi);
}
- return K*sqrt((2.0*l+1.0)/2.0*fac1/fac2)*plgndr(l,abs(m),cos(theta));
+ return K*sqrt((2.0*l+1.0)/2.0*(ts_double)(fac1/fac2))*plgndr(l,abs(m),cos(theta));
}
@@ -178,7 +203,7 @@
#ifdef TS_DOUBLE_DOUBLE
coord->e1=sqrt(x*x+y*y+z*z);
if(z==0) coord->e3=M_PI/2.0;
- else coord->e3=atan(sqrt(x*x+y*y)/z);
+ else coord->e3=atan2(sqrt(x*x+y*y),z);
coord->e2=atan2(y,x);
#endif
#ifdef TS_DOUBLE_FLOAT
@@ -196,6 +221,23 @@
return TS_SUCCESS;
}
+
+
+ts_bool sph2cart(ts_coord *coord){
+ coord->coord_type=TS_COORD_CARTESIAN;
+ ts_double x,y,z;
+
+ x=coord->e1*cos(coord->e2)*sin(coord->e3);
+ y=coord->e1*sin(coord->e2)*sin(coord->e3);
+ z=coord->e1*cos(coord->e3);
+
+ coord->e1=x;
+ coord->e2=y;
+ coord->e3=z;
+
+ return TS_SUCCESS;
+}
+
/* Function returns radius of the sphere with the same volume as vesicle (r0) */
ts_double getR0(ts_vesicle *vesicle){
@@ -216,6 +258,7 @@
ts_bool preparationSh(ts_vesicle *vesicle, ts_double r0){
//TODO: before calling or during the call calculate area of each triangle! Can
//be also done after vertexmove and bondflip //
+//DONE: in energy calculation! //
ts_uint i,j;
ts_vertex **vtx=vesicle->vlist->vtx;
ts_vertex *cvtx;
@@ -258,7 +301,7 @@
r=sqrtl(cvtx->x*cvtx->x+cvtx->y*cvtx->y+cvtx->z*cvtx->z);
#endif
cvtx->relR=(r-r0)/r0;
- cvtx->solAngle=cvtx->projArea/cvtx->relR * cvtx->projArea/cvtx->relR;
+ cvtx->solAngle=cvtx->projArea/r/r;
}
return TS_SUCCESS;
}
@@ -266,10 +309,11 @@
ts_bool calculateYlmi(ts_vesicle *vesicle){
- ts_uint i,j,k;
+ ts_int i,j,k;
ts_spharm *sph=vesicle->sphHarmonics;
ts_coord *coord=(ts_coord *)malloc(sizeof(ts_coord));
ts_double fi, theta;
+ ts_int m;
ts_vertex *cvtx;
for(k=0;k<vesicle->vlist->n;k++){
cvtx=vesicle->vlist->vtx[k];
@@ -277,13 +321,34 @@
cart2sph(coord,cvtx->x, cvtx->y, cvtx->z);
fi=coord->e2;
theta=coord->e3;
- for(i=0; i<sph->l; i++){
+ for(i=1; i<sph->l; i++){
for(j=0;j<i;j++){
- sph->Ylmi[i][j][k]=sph->co[i][j]*cos((j-i-1)*fi)*pow(-1,j-i-1)*plgndr(i,abs(j-i-1),cos(theta));
+ m=j+1;
+//Nastudiraj!!!!!
+ sph->Ylmi[i][j][k]=sph->co[i][m]*cos((m-i-1)*fi)*pow(-1,m-i-1)*plgndr(i,abs(m-i-1),cos(theta));
+ if(i==2 && j==0){
+ /* fprintf(stderr," **** vtx %d ****\n", k+1);
+ fprintf(stderr,"m-i-1 =%d\n",m-i-1);
+ fprintf(stderr,"fi =%e\n",fi);
+ fprintf(stderr,"(m-i-1)*fi =%e\n",((ts_double)(m-i-1))*fi);
+ fprintf(stderr,"-2*fi =%e\n",-2*fi);
+ fprintf(stderr,"m =%d\n",m);
+
+ fprintf(stderr," cos(m-i-1)=%e\n",cos((m-i-1)*fi));
+ fprintf(stderr," cos(-2*fi)=%e\n",cos(-2*fi));
+ fprintf(stderr," sph->co[i][m]=%e\n",sph->co[i][m]);
+ fprintf(stderr," plgndr(i,abs(m-i-1),cos(theta))=%e\n",plgndr(i,abs(m-i-1),cos(theta)));
+*/
+ }
}
- sph->Ylmi[i][j+1][k]=sph->co[i][j+1]*plgndr(i,0,cos(theta));
- for(j=sph->l;j<2*i;j++){
- sph->Ylmi[i][j][k]=sph->co[i][j]*sin((j-i-1)*fi)*plgndr(i,j-i-1,cos(theta));
+//Nastudiraj!!!!!
+ j=i;
+ m=j+1;
+ sph->Ylmi[i][j][k]=sph->co[i][m]*plgndr(i,0,cos(theta));
+ for(j=i+1;j<2*i+1;j++){
+ m=j+1;
+//Nastudiraj!!!!!
+ sph->Ylmi[i][j][k]=sph->co[i][m]*sin((m-i-1)*fi)*plgndr(i,m-i-1,cos(theta));
}
}
@@ -298,7 +363,7 @@
ts_uint i,j,k;
ts_vertex *cvtx;
for(i=0;i<vesicle->sphHarmonics->l;i++){
- for(j=0;j<2*i;j++) vesicle->sphHarmonics->ulm[i][j]=0.0;
+ for(j=0;j<2*i+1;j++) vesicle->sphHarmonics->ulm[i][j]=0.0;
}
//TODO: call calculateYlmi !!!
@@ -307,7 +372,7 @@
for(k=0;k<vesicle->vlist->n; k++){
cvtx=vesicle->vlist->vtx[k];
for(i=0;i<vesicle->sphHarmonics->l;i++){
- for(j=0;j<2*i;j++){
+ for(j=0;j<2*i+1;j++){
vesicle->sphHarmonics->ulm[i][j]+= cvtx->solAngle*cvtx->relR*vesicle->sphHarmonics->Ylmi[i][j][k];
}
@@ -316,3 +381,48 @@
return TS_SUCCESS;
}
+
+
+
+
+
+ts_bool storeUlm2(ts_vesicle *vesicle){
+
+ts_spharm *sph=vesicle->sphHarmonics;
+ts_int i,j;
+for(i=0;i<sph->l;i++){
+ for(j=0;j<2*i+1;j++){
+ /* DEBUG fprintf(stderr,"sph->sumUlm2[%d][%d]=%e\n",i,j,sph->ulm[i][j]* sph->ulm[i][j]); */
+ sph->sumUlm2[i][j]+=sph->ulm[i][j]* sph->ulm[i][j];
+ }
+}
+ sph->N++;
+return TS_SUCCESS;
+}
+
+
+ts_bool saveAvgUlm2(ts_vesicle *vesicle){
+
+ FILE *fh;
+ char filename[10000];
+ strcpy(filename, command_line_args.path);
+ strcat(filename, "sph2out.dat");
+ fh=fopen(filename, "w");
+ if(fh==NULL){
+ err("Cannot open file %s for writing");
+ return TS_FAIL;
+ }
+
+ ts_spharm *sph=vesicle->sphHarmonics;
+ ts_int i,j;
+ fprintf(fh,"l,\tm,\tulm^2avg\n");
+ for(i=0;i<sph->l;i++){
+ for(j=0;j<2*i+1;j++){
+ fprintf(fh,"%d,\t%d,\t%e\n", i, j-i, sph->sumUlm2[i][j]/(ts_double)sph->N);
+
+ }
+ fprintf(fh,"\n");
+ }
+ fclose(fh);
+ return TS_SUCCESS;
+}
--
Gitblit v1.9.3