From 88f45143117f3ae973be6e3ea5dac8063a7b8ece Mon Sep 17 00:00:00 2001
From: Samo Penic <samo@amalthea>
Date: Sat, 19 Mar 2011 18:35:08 +0000
Subject: [PATCH] Added spherical harmonics routine.I've got a feeling that some harmonics are not calculated correctly!
---
src/Makefile.am | 7 +
src/shdiscover.c | 80 ++++++++++++++++++++
src/sh.h | 9 ++
src/sh.c | 90 ++++++++++++++++++++++
4 files changed, 184 insertions(+), 2 deletions(-)
diff --git a/src/Makefile.am b/src/Makefile.am
index 25fc76f..66eec22 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,5 +1,8 @@
trisurfdir=../
trisurf_PROGRAMS = trisurf
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 main.c
-trisurf_LDFLAGS = -lm -lconfuse
-trisurf_CFLAGS = -Wall -Werror
+#trisurf_LDFLAGS = -lm -lconfuse
+shdiscoverdir=../
+shdiscover_PROGRAMS= shdiscover
+shdiscover_SOURCES= general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c energy.c sh.c shdiscover.c
+AM_CFLAGS = -Wall -Werror
diff --git a/src/sh.c b/src/sh.c
new file mode 100644
index 0000000..b494b25
--- /dev/null
+++ b/src/sh.c
@@ -0,0 +1,90 @@
+#include<math.h>
+#include<stdlib.h>
+#include "general.h"
+#include "sh.h"
+
+/* Gives you legendre polynomials. Taken from NR, p. 254 */
+ts_double plgndr(ts_int l, ts_int m, ts_float x){
+ ts_double fact, pll, pmm, pmmp1, somx2;
+ ts_int i,ll;
+
+#ifdef TS_DOUBLE_DOUBLE
+ if(m<0 || m>l || fabs(x)>1.0)
+ fatal("Bad arguments in routine plgndr",1);
+#endif
+#ifdef TS_DOUBLE_FLOAT
+ if(m<0 || m>l || fabsf(x)>1.0)
+ fatal("Bad arguments in routine plgndr",1);
+#endif
+#ifdef TS_DOUBLE_LONGDOUBLE
+ if(m<0 || m>l || fabsl(x)>1.0)
+ fatal("Bad arguments in routine plgndr",1);
+#endif
+ pmm=1.0;
+ if (m>0) {
+#ifdef TS_DOUBLE_DOUBLE
+ somx2=sqrt((1.0-x)*(1.0+x));
+#endif
+#ifdef TS_DOUBLE_FLOAT
+ somx2=sqrtf((1.0-x)*(1.0+x));
+#endif
+#ifdef TS_DOUBLE_LONGDOUBLE
+ somx2=sqrtl((1.0-x)*(1.0+x));
+#endif
+ fact=1.0;
+ for (i=1; i<=m;i++){
+ pmm *= -fact*somx2;
+ fact +=2.0;
+ }
+ }
+
+ if (l == m) return pmm;
+ else {
+ pmmp1=x*(2*m+1)*pmm;
+ if(l==(m+1)) return(pmmp1);
+ else {
+ pll=0; /* so it can not be uninitialized */
+ for(ll=m+2;ll<=l;ll++){
+ pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m);
+ pmm=pmmp1;
+ pmmp1=pll;
+ }
+ return(pll);
+ }
+ }
+}
+
+
+/*Computes Y(l,m,theta,fi) (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;
+
+ if(l<0 || m>l || m<-l)
+ fatal("Error using shY function!",1);
+
+ fac1=1.0;
+ for(i=1; i<=l-m;i++){
+ fac1 *= i;
+ }
+ fac2=1.0;
+ for(i=1; i<=l+m;i++){
+ fac2 *= i;
+ }
+
+ if(m==0){
+ K=sqrt(1.0/(2.0*M_PI));
+ }
+ else if (m>0) {
+ K=sqrt(1.0/(M_PI))*cos(m*fi);
+ }
+ else {
+ //K=pow(-1.0,abs(m))*sqrt(1.0/(2.0*M_PI))*cos(m*fi);
+ if(abs(m)%2==0)
+ K=sqrt(1.0/(M_PI))*sin(m*fi);
+ else
+ K=-sqrt(1.0/(M_PI))*sin(m*fi);
+ }
+
+ return K*sqrt((2.0*l+1.0)/2.0*fac1/fac2)*plgndr(l,abs(m),cos(theta));
+}
diff --git a/src/sh.h b/src/sh.h
new file mode 100644
index 0000000..1255f52
--- /dev/null
+++ b/src/sh.h
@@ -0,0 +1,9 @@
+#ifndef _H_SH
+#define _H_SH
+#include "general.h"
+
+ts_double plgndr(ts_int l, ts_int m, ts_float x);
+ts_double shY(ts_int l,ts_int m,ts_double theta,ts_double fi);
+
+
+#endif
diff --git a/src/shdiscover.c b/src/shdiscover.c
new file mode 100644
index 0000000..9783ef6
--- /dev/null
+++ b/src/shdiscover.c
@@ -0,0 +1,80 @@
+#include "general.h"
+#include "vertex.h"
+#include "initial_distribution.h"
+#include "io.h"
+#include "vesicle.h"
+#include "sh.h"
+#include "frame.c"
+#include <math.h>
+#include <stdlib.h>
+int main(int argc, char *argv[]){
+
+ts_fprintf(stdout,"SHdiscover was called with %d coefficients!\n",argc-1);
+ts_uint n,i,j,l;
+ts_int m;
+ts_double fi,theta,r,Y;
+ts_vesicle *vesicle=initial_distribution_dipyramid(17,60,60,60,0.15);
+ts_vertex_list *vlist=vesicle->vlist;
+centermass(vesicle);
+ts_fprintf(stdout,"Vesicle has a CenterMass in %f,%f,%f\n",vesicle->cm[0],vesicle->cm[1], vesicle->cm[2]);
+
+n=vlist->n;
+
+ts_fprintf(stdout,"Tests\n");
+ts_fprintf(stdout,"P(0,0,0.5)=%f (%f)\n",plgndr(0,0,0.5),1.0);
+ts_fprintf(stdout,"P(1,0,0.5)=%f (%f)\n",plgndr(1,0,0.5),0.5);
+ts_fprintf(stdout,"P(2,0,0.5)=%f (%f)\n",plgndr(2,0,0.5),0.5*(3*0.5*0.5-1));
+
+ts_fprintf(stdout,"Y(0,0,pi/6,pi/4)=%f (%f)\n",shY(0,0,M_PI/6,M_PI/4),sqrt(1/(4*M_PI)));
+ts_fprintf(stdout,"Y(1,0,pi/6,pi/4)=%f (%f)\n",shY(1,0,M_PI/6,M_PI/4),sqrt(3/(4*M_PI))*cos(M_PI/6));
+ts_fprintf(stdout,"Y(1,0,4*pi/6,6*pi/4)=%f (%f)\n",shY(1,0,4*M_PI/6,6*M_PI/4),sqrt(3/(4*M_PI))*cos(4*M_PI/6));
+ts_fprintf(stdout,"Y(1,1,pi/6,pi/4)=%f (%f)\n",shY(1,1,M_PI/6,M_PI/4),-sqrt(3/(8*M_PI))*sin(M_PI/6)*sin(M_PI/4));
+ts_fprintf(stdout,"Y(2,0,pi/6,pi/4)=%f (%f)\n",shY(2,0,M_PI/6,M_PI/4),sqrt(5/(4*M_PI))*(3.0/2.0*cos(M_PI/6)*cos(M_PI/6)-1.0/2.0));
+ts_fprintf(stdout,"Y(2,-2,pi/6,pi/4)=%f (0)\n",shY(2,-2,M_PI/6,M_PI/4));
+ts_fprintf(stdout,"Y(2,2,pi/6,pi/4)=%f (0)\n",shY(2,2,M_PI/6,M_PI/4));
+
+ for(j=1;j<argc;j++){
+ l=(int)sqrt(j-1); /* determine l from dataline */
+ m=j-1-l*(l+1); /* determine m from dataline */
+ ts_fprintf(stdout,"l=%d, m=%d, u=%s\n",l,m,argv[j]);
+ }
+
+/*we calculate new position of each vertex of vesicle */
+for(i=0;i<n;i++){
+ fi=atan2(vlist->vtx[i]->data->y, vlist->vtx[i]->data->x);
+/* theta=atan2(
+ sqrt(vlist->vtx[i]->data->x*vlist->vtx[i]->data->x +
+ vlist->vtx[i]->data->y*vlist->vtx[i]->data->y),
+ vlist->vtx[i]->data->z
+ ); */
+ theta=acos(
+ vlist->vtx[i]->data->z /
+ sqrt(vlist->vtx[i]->data->x*vlist->vtx[i]->data->x +
+ vlist->vtx[i]->data->y*vlist->vtx[i]->data->y+
+ vlist->vtx[i]->data->z*vlist->vtx[i]->data->z)
+
+ );
+
+
+
+ r=0.0;
+ for(j=1;j<argc;j++){
+ l=(int)sqrt(j-1); /* determine l from dataline */
+ m=j-1-l*(l+1); /* determine m from dataline */
+ Y=shY(l,m,theta,fi);
+ r+=fabs(atof(argv[j])*Y);
+ /*ts_fprintf(stdout,"l=%d, m=%d, u=%s\n",l,m,argv[j]);*/
+ }
+
+ vlist->vtx[i]->data->z=fabs(r)*cos(theta);
+ vlist->vtx[i]->data->x=fabs(r)*sin(theta)*cos(fi);
+ vlist->vtx[i]->data->y=fabs(r)*sin(theta)*sin(fi);
+}
+
+write_vertex_xml_file(vesicle,0);
+write_master_xml_file("test.pvd");
+
+
+vesicle_free(vesicle);
+return 0;
+}
--
Gitblit v1.9.3