From df111bb0cf32bbe48718a5b77f8b9d9a7400d496 Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@gmail.com>
Date: Thu, 15 Aug 2019 06:52:41 +0000
Subject: [PATCH] Fix in printing the ulms
---
src/shcomplex.c | 67 +++++++++++++++++++++++++--------
1 files changed, 51 insertions(+), 16 deletions(-)
diff --git a/src/shcomplex.c b/src/shcomplex.c
index 9350fe5..4090346 100644
--- a/src/shcomplex.c
+++ b/src/shcomplex.c
@@ -1,5 +1,7 @@
+/* vim: set ts=4 sts=4 sw=4 noet : */
#include<math.h>
#include<stdlib.h>
+#include<string.h>
#include<gsl/gsl_complex.h>
#include<gsl/gsl_complex_math.h>
#include<gsl/gsl_sf_legendre.h>
@@ -68,7 +70,7 @@
if(sph->co != NULL) free(sph->co);
if(sph->ulm !=NULL) free(sph->ulm);
if(sph->ulmComplex !=NULL) free(sph->ulmComplex);
-
+ if(sph->sumUlm2 !=NULL) free(sph->sumUlm2);
if(sph->Ylmi!=NULL) {
for(i=0;i<sph->l;i++){
if(sph->Ylmi[i]!=NULL){
@@ -116,6 +118,27 @@
return TS_SUCCESS;
}
+char *Ulm2Complex2String(ts_vesicle *vesicle){
+ ts_int i,j;
+ char *strng=(char *)calloc(5000, sizeof(char));
+ char tmpstrng[255];
+ for(i=0;i<vesicle->sphHarmonics->l;i++){
+ for(j=i;j<2*i+1;j++){
+ sprintf(tmpstrng,"%e ", gsl_complex_abs2(vesicle->sphHarmonics->ulmComplex[i][j]));
+ strcat(strng,tmpstrng);
+ }
+ }
+ //strcat(strng,"\n");
+
+ return strng;
+}
+
+ts_bool freeUlm2String(char *strng){
+ free(strng);
+ return TS_SUCCESS;
+}
+
+
ts_bool storeUlmComplex2(ts_vesicle *vesicle){
ts_spharm *sph=vesicle->sphHarmonics;
@@ -130,31 +153,43 @@
}
-ts_bool calculateKc(ts_vesicle *vesicle){
-
- ts_int i;
- gsl_matrix *A=gsl_matrix_alloc(vesicle->sphHarmonics->l,2);
+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);
gsl_vector *tau=gsl_vector_alloc(2);
- gsl_vector *b=gsl_vector_alloc(vesicle->sphHarmonics->l);
+ gsl_vector *b=gsl_vector_alloc(max-min);
gsl_vector *x=gsl_vector_alloc(2);
- gsl_vector *res=gsl_vector_alloc(vesicle->sphHarmonics->l);
+ gsl_vector *res=gsl_vector_alloc(max-min);
- //solving (A^T*A)*x=A^T*UlmSqAvg
+ //solving (A^T*A)*x=A^T*b
//fill the data for matrix A and vector b
- for(i=1;i<=vesicle->sphHarmonics->l;i++){
- gsl_matrix_set(A, i-1,0,(ts_double)((i-1)*(i+2)));
- gsl_matrix_set(A, i-1,1,(ts_double)((i-1)*(i+2)*(i+1)*i));
- gsl_vector_set(b,i-1,(ts_double)vesicle->sphHarmonics->N/vesicle->sphHarmonics->sumUlm2[i-1][(i-1)*2]);
+ 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));
+ bval=0.0;
+ //average for m from 0..l (only positive m's)
+ for(j=0;j<=i;j++){
+ bval+=vesicle->sphHarmonics->sumUlm2[i][(j+i)];
+ }
+ 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,"b[2]=%e\n",gsl_vector_get(b,1));
+// fprintf(stderr,"b[2]=%e\n",gsl_vector_get(b,1));
gsl_linalg_QR_decomp(A,tau);
gsl_linalg_QR_lssolve(A,tau,b,x,res);
- fprintf(stderr,"kc=%e\n",gsl_vector_get(x,1));
-
+// fprintf(stderr,"kc=%e\n",gsl_vector_get(x,1));
+ retval=gsl_vector_get(x,1);
gsl_matrix_free(A);
gsl_vector_free(tau);
gsl_vector_free(b);
gsl_vector_free(x);
gsl_vector_free(res);
- return TS_SUCCESS;
+
+ return retval;
}
--
Gitblit v1.9.3