From f06f5f9ad1e42a10d6a9367513e66c0c06295de2 Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@fe.uni-lj.si>
Date: Wed, 23 Apr 2014 12:18:03 +0000
Subject: [PATCH] Calculating Kc. Needs debugging.
---
src/shcomplex.c | 33 +++++++++++++++++++++++++++++++++
1 files changed, 33 insertions(+), 0 deletions(-)
diff --git a/src/shcomplex.c b/src/shcomplex.c
index d9488ca..9350fe5 100644
--- a/src/shcomplex.c
+++ b/src/shcomplex.c
@@ -3,6 +3,10 @@
#include<gsl/gsl_complex.h>
#include<gsl/gsl_complex_math.h>
#include<gsl/gsl_sf_legendre.h>
+
+#include<gsl/gsl_matrix.h>
+#include<gsl/gsl_vector.h>
+#include<gsl/gsl_linalg.h>
#include "general.h"
#include "sh.h"
#include "shcomplex.h"
@@ -125,3 +129,32 @@
return TS_SUCCESS;
}
+
+ts_bool calculateKc(ts_vesicle *vesicle){
+
+ ts_int i;
+ gsl_matrix *A=gsl_matrix_alloc(vesicle->sphHarmonics->l,2);
+ gsl_vector *tau=gsl_vector_alloc(2);
+ gsl_vector *b=gsl_vector_alloc(vesicle->sphHarmonics->l);
+ gsl_vector *x=gsl_vector_alloc(2);
+ gsl_vector *res=gsl_vector_alloc(vesicle->sphHarmonics->l);
+
+ //solving (A^T*A)*x=A^T*UlmSqAvg
+ //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]);
+ }
+ 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));
+
+ gsl_matrix_free(A);
+ gsl_vector_free(tau);
+ gsl_vector_free(b);
+ gsl_vector_free(x);
+ gsl_vector_free(res);
+ return TS_SUCCESS;
+}
--
Gitblit v1.9.3