From 3096998157b37ad4d6c1af9ce32e657407795d5d Mon Sep 17 00:00:00 2001 From: Samo Penic <samo.penic@fe.uni-lj.si> Date: Thu, 07 Jun 2012 13:43:50 +0000 Subject: [PATCH] Merge branch 'Ylmi2sh' --- src/sh.c | 79 +++++++++++++++++++++++++++++++++++++-- 1 files changed, 74 insertions(+), 5 deletions(-) diff --git a/src/sh.c b/src/sh.c index 055a3c0..7cbb04f 100644 --- a/src/sh.c +++ b/src/sh.c @@ -3,6 +3,66 @@ #include "general.h" #include "sh.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)); + + /* 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++){ + sph->Ylmi[i][j]=(ts_double *)calloc(vlist->n,sizeof(ts_double)); + } + } + + /* lets initialize ulm */ + sph->ulm=(ts_double **)calloc(l,sizeof(ts_double *)); + for(j=0;j<l;j++){ + sph->ulm[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)); + } + + /* Calculate coefficients that will remain constant during all the simulation */ + precomputeShCoeff(sph); + + return sph; +} + + +ts_bool sph_free(ts_spharm *sph){ + int i,j; + for(i=0;i<sph->l;i++){ + if(sph->ulm[i]!=NULL) free(sph->ulm[i]); + if(sph->co[i]!=NULL) free(sph->co[i]); + } + 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++){ + if(sph->Ylmi[i][j]!=NULL) free (sph->Ylmi[i][j]); + } + free(sph->Ylmi[i]); + } + } + free(sph->Ylmi); + } + + free(sph); + return TS_SUCCESS; +} + /* 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; @@ -57,12 +117,19 @@ ts_bool precomputeShCoeff(ts_spharm *sph){ - ts_uint i,j; + ts_int i,j,al,am; + ts_double **co=sph->co; for(i=0;i<sph->l;i++){ - sph->co[i][i]=sqrt((2.0*i+1.0)/2.0/M_PI); - for(j=0;j<i-1;j++){ - + al=i+1; + 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)); + 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][i+1]=sqrt((2.0*al+1.0)/4.0/M_PI); } return TS_SUCCESS; @@ -203,9 +270,11 @@ ts_spharm *sph=vesicle->sphHarmonics; ts_coord *coord=(ts_coord *)malloc(sizeof(ts_coord)); ts_double fi, theta; + ts_vertex *cvtx; for(k=0;k<vesicle->vlist->n;k++){ + cvtx=vesicle->vlist->vtx[k]; sph->Ylmi[0][0][k]=sqrt(1.0/4.0/M_PI); - cart2sph(coord,vesicle->vlist->vtx[k]->x, vesicle->vlist->vtx[k]->y, vesicle->vlist->vtx[k]->z); + cart2sph(coord,cvtx->x, cvtx->y, cvtx->z); fi=coord->e2; theta=coord->e3; for(i=0; i<sph->l; i++){ -- Gitblit v1.9.3