From fda1ab6babed79842534b3a21a6ee96bc26f9d93 Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@fe.uni-lj.si>
Date: Tue, 16 Dec 2014 14:47:56 +0000
Subject: [PATCH] Real spherical harmonics are programmed. Inefficiently, but for testing purposes they works.
---
src/shcomplex.c | 44 +++++++++++++++++++++++++++++++++++++-------
1 files changed, 37 insertions(+), 7 deletions(-)
diff --git a/src/shcomplex.c b/src/shcomplex.c
index 7afe27e..8df133c 100644
--- a/src/shcomplex.c
+++ b/src/shcomplex.c
@@ -9,8 +9,10 @@
#include<gsl/gsl_linalg.h>
#include "general.h"
#include "sh.h"
+#include "shreal.h"
#include "shcomplex.h"
-
+#include "string.h"
+#include "io.h"
ts_spharm *complex_sph_init(ts_vertex_list *vlist, ts_uint l){
ts_uint j,i;
@@ -51,13 +53,14 @@
/* Calculate coefficients that will remain constant during all the simulation */
precomputeShCoeff(sph);
-
+ shreal_init(sph,l); //this is for TESTING only!
return sph;
}
ts_bool complex_sph_free(ts_spharm *sph){
int i,j;
if(sph==NULL) return TS_FAIL;
+ shreal_free(sph); //this is for TESTING only!
for(i=0;i<sph->l;i++){
if(sph->ulm[i]!=NULL) free(sph->ulm[i]);
if(sph->ulmComplex[i]!=NULL) free(sph->ulmComplex[i]);
@@ -130,9 +133,9 @@
}
-ts_double calculateKc(ts_vesicle *vesicle){
- ts_int min=4;
- ts_int max=vesicle->sphHarmonics->l-6;
+ts_double calculateKc(ts_vesicle *vesicle, ts_int lmin, ts_int lmax){
+ ts_int min=lmin;
+ ts_int max=lmax; //vesicle->sphHarmonics->l-3;
ts_long i,j;
ts_double retval, bval;
gsl_matrix *A=gsl_matrix_alloc(max-min,2);
@@ -146,7 +149,7 @@
for(i=min;i<max;i++){
gsl_matrix_set(A, i-min,0,(ts_double)((i-1)*(i+2)));
gsl_matrix_set(A, i-min,1,(ts_double)((i-1)*(i+2)*(i+1)*i));
- fprintf(stderr,"%e %e\n", gsl_matrix_get(A,i-min,0), gsl_matrix_get(A,i-min,1));
+// fprintf(stderr,"%e %e\n", gsl_matrix_get(A,i-min,0), gsl_matrix_get(A,i-min,1));
bval=0.0;
//average for m from 0..l (only positive m's)
for(j=0;j<=i;j++){
@@ -155,7 +158,7 @@
bval=bval/(ts_double)vesicle->sphHarmonics->N/(ts_double)(i+1);
gsl_vector_set(b,i-min,1.0/bval);
- fprintf(stderr,"%e\n", 1.0/gsl_vector_get(b,i-min));
+// fprintf(stderr,"%e\n", 1.0/gsl_vector_get(b,i-min));
}
// fprintf(stderr,"b[2]=%e\n",gsl_vector_get(b,1));
gsl_linalg_QR_decomp(A,tau);
@@ -170,3 +173,30 @@
return retval;
}
+
+
+ts_bool saveAvgUlm2Complex(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