From 262607715dced9c01b742ac2127af65750a3be3b Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@gmail.com>
Date: Wed, 27 Jun 2012 14:53:51 +0000
Subject: [PATCH] Added program that calculates square spherical harmonics averages and stores them into file
---
src/Makefile.am | 3 +
src/sh.h | 2
src/spherical_trisurf.c | 95 +++++++++++++++++++++++++++++++
src/general.h | 2
src/sh.c | 25 ++++++++
5 files changed, 125 insertions(+), 2 deletions(-)
diff --git a/src/Makefile.am b/src/Makefile.am
index 0ba3c2a..d1ae2ab 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -9,3 +9,6 @@
co_testdir=../
co_test_PROGRAMS=co_test
co_test_SOURCES= general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c energy.c sh.c co_test.c frame.c
+spherical_trisurfdir=../
+spherical_trisurf_PROGRAMS = spherical_trisurf
+spherical_trisurf_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c spherical_trisurf.c sh.c
diff --git a/src/general.h b/src/general.h
index 87c6174..6263862 100644
--- a/src/general.h
+++ b/src/general.h
@@ -217,6 +217,8 @@
typedef struct {
ts_uint l;
ts_double **ulm;
+ ts_double **sumUlm2;
+ ts_uint N;
ts_double **co;
ts_double ***Ylmi;
} ts_spharm;
diff --git a/src/sh.c b/src/sh.c
index ea050b9..31df74f 100644
--- a/src/sh.c
+++ b/src/sh.c
@@ -9,7 +9,7 @@
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<l;i++){
@@ -25,6 +25,11 @@
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 */
//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!
@@ -46,6 +51,7 @@
int i,j;
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]);
@@ -354,3 +360,20 @@
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++){
+ sph->sumUlm2[i][j]+=sph->ulm[i][j]* sph->ulm[i][j];
+ }
+}
+ sph->N++;
+return TS_SUCCESS;
+}
diff --git a/src/sh.h b/src/sh.h
index 93bed39..7853139 100644
--- a/src/sh.h
+++ b/src/sh.h
@@ -1,7 +1,7 @@
#ifndef _H_SH
#define _H_SH
#include "general.h"
-
+ts_bool storeUlm2(ts_vesicle *vesicle);
ts_double plgndr(ts_int l, ts_int m, ts_double x);
ts_double shY(ts_int l,ts_int m,ts_double theta,ts_double fi);
diff --git a/src/spherical_trisurf.c b/src/spherical_trisurf.c
new file mode 100644
index 0000000..f81e45d
--- /dev/null
+++ b/src/spherical_trisurf.c
@@ -0,0 +1,95 @@
+#include<stdio.h>
+#include<math.h>
+#include "general.h"
+#include "vertex.h"
+#include "bond.h"
+#include "triangle.h"
+#include "cell.h"
+#include "vesicle.h"
+#include "io.h"
+#include "initial_distribution.h"
+#include "frame.h"
+#include "timestep.h"
+#include "sh.h"
+
+/** Entrance function to the program
+ * @param argv is a number of parameters used in program call (including the program name
+ * @param argc is a pointer to strings (character arrays) which holds the arguments
+ * @returns returns 0 on success, any other number on fail.
+*/
+ts_bool saveAvgUlm2(ts_vesicle *vesicle);
+int main(int argv, char *argc[]){
+ts_uint i,j;
+ts_vesicle *vesicle;
+ts_double r0;
+vesicle=initial_distribution_dipyramid(17,60,60,60,0.15);
+//parsetape(vesicle,&i);
+
+//these four must come from parsetype!
+vesicle->dmax=1.67*1.67;
+vesicle->stepsize=0.15;
+vesicle->clist->max_occupancy=8;
+vesicle->bending_rigidity=25.0;
+//fprintf(stderr,"xk=%f",vesicle->bending_rigidity);
+
+vesicle->sphHarmonics=sph_init(vesicle->vlist, 21);
+
+vesicle_volume(vesicle);
+r0=getR0(vesicle);
+
+preparationSh(vesicle,r0);
+calculateYlmi(vesicle);
+calculateUlm(vesicle);
+
+
+
+for(i=0;i<10;i++){
+ centermass(vesicle);
+ cell_occupation(vesicle);
+ for(j=0;j<1000;j++){
+ single_timestep(vesicle);
+ }
+vesicle_volume(vesicle);
+r0=getR0(vesicle);
+
+preparationSh(vesicle,r0);
+calculateYlmi(vesicle);
+calculateUlm(vesicle);
+
+storeUlm2(vesicle);
+saveAvgUlm2(vesicle);
+
+ write_vertex_xml_file(vesicle,i);
+ fprintf(stderr, "Loop %d completed.\n",i+1);
+}
+write_master_xml_file("test.pvd");
+write_dout_fcompat_file(vesicle,"dout");
+vesicle_free(vesicle);
+
+return 0; //program finished perfectly ok. We return 0.
+}
+
+
+
+ts_bool saveAvgUlm2(ts_vesicle *vesicle){
+
+ FILE *fh;
+
+ fh=fopen("sph2out.dat", "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);
+
+ }
+ }
+ fclose(fh);
+ return TS_SUCCESS;
+}
--
Gitblit v1.9.3