trisurf-ng.git - Gitblit

Samo Penic

2016-07-09 9e5a9c0f110d3ed02751c9d44a3716caa6010a70

commit \| author \| age
7f6076	1	/* vim: set ts=4 sts=4 sw=4 noet : */
459ff9	2	#include<math.h>
SP	3	#include<stdlib.h>
	4	#include<gsl/gsl_complex.h>
	5	#include<gsl/gsl_complex_math.h>
	6	#include<gsl/gsl_sf_legendre.h>
f06f5f	7
SP	8	#include<gsl/gsl_matrix.h>
	9	#include<gsl/gsl_vector.h>
	10	#include<gsl/gsl_linalg.h>
459ff9	11	#include "general.h"
SP	12	#include "sh.h"
	13	#include "shcomplex.h"
	14
	15
	16	ts_spharm complex_sph_init(ts_vertex_list vlist, ts_uint l){
	17	ts_uint j,i;
	18	ts_spharm sph=(ts_spharm )malloc(sizeof(ts_spharm));
	19
	20	sph->N=0;
	21	/* lets initialize Ylm for each vertex. */
	22	sph->Ylmi=(ts_double *)calloc(l,sizeof(ts_double ));
	23	for(i=0;i<l;i++){
	24	sph->Ylmi[i]=(ts_double *)calloc(2i+1,sizeof(ts_double *));
	25	for(j=0;j<(2*i+1);j++){
	26	sph->Ylmi[i][j]=(ts_double *)calloc(vlist->n,sizeof(ts_double));
	27	}
	28	}
	29
	30	/* lets initialize ulm */
	31	sph->ulm=(ts_double *)calloc(l,sizeof(ts_double ));
	32	sph->ulmComplex=(gsl_complex *)calloc(l,sizeof(gsl_complex ));
	33	for(j=0;j<l;j++){
	34	sph->ulm[j]=(ts_double )calloc(2j+1,sizeof(ts_double));
	35	sph->ulmComplex[j]=(gsl_complex )calloc(2j+1,sizeof(gsl_complex));
	36	}
	37
	38	/* lets initialize sum of Ulm2 */
	39	sph->sumUlm2=(ts_double *)calloc(l,sizeof(ts_double ));
	40	for(j=0;j<l;j++){
	41	sph->sumUlm2[j]=(ts_double )calloc(2j+1,sizeof(ts_double));
	42	}
	43
	44	/* lets initialize co */
	45	//NOTE: C is has zero based indexing. Code is imported from fortran and to comply with original indexes we actually generate one index more. Also second dimension is 2j+2 instead of 2j+2. elements starting with 0 are useles and should be ignored!
	46	sph->co=(ts_double *)calloc(l+1,sizeof(ts_double ));
	47	for(j=0;j<=l;j++){
	48	sph->co[j]=(ts_double )calloc(2j+2,sizeof(ts_double));
	49	}
	50
	51	sph->l=l;
	52
	53	/* Calculate coefficients that will remain constant during all the simulation */
	54	precomputeShCoeff(sph);
	55
	56	return sph;
	57	}
	58
	59	ts_bool complex_sph_free(ts_spharm *sph){
	60	int i,j;
	61	if(sph==NULL) return TS_FAIL;
	62	for(i=0;i<sph->l;i++){
	63	if(sph->ulm[i]!=NULL) free(sph->ulm[i]);
	64	if(sph->ulmComplex[i]!=NULL) free(sph->ulmComplex[i]);
	65	if(sph->sumUlm2[i]!=NULL) free(sph->sumUlm2[i]);
	66	if(sph->co[i]!=NULL) free(sph->co[i]);
	67	}
	68	if(sph->co[sph->l]!=NULL) free(sph->co[sph->l]);
	69	if(sph->co != NULL) free(sph->co);
	70	if(sph->ulm !=NULL) free(sph->ulm);
	71	if(sph->ulmComplex !=NULL) free(sph->ulmComplex);
701026	72	if(sph->sumUlm2 !=NULL) free(sph->sumUlm2);
459ff9	73	if(sph->Ylmi!=NULL) {
SP	74	for(i=0;i<sph->l;i++){
	75	if(sph->Ylmi[i]!=NULL){
	76	for(j=0;j<i*2+1;j++){
	77	if(sph->Ylmi[i][j]!=NULL) free (sph->Ylmi[i][j]);
	78	}
	79	free(sph->Ylmi[i]);
	80	}
	81	}
	82	free(sph->Ylmi);
	83	}
	84
	85	free(sph);
	86	return TS_SUCCESS;
	87	}
	88
	89
	90	ts_bool calculateUlmComplex(ts_vesicle *vesicle){
	91	ts_int i,j,k,m,l;
	92	ts_vertex *cvtx;
	93	ts_coord coord;
	94	/* set all values to zero */
	95	for(i=0;i<vesicle->sphHarmonics->l;i++){
	96	for(j=0;j<2*i+1;j++) GSL_SET_COMPLEX(&(vesicle->sphHarmonics->ulmComplex[i][j]),0.0,0.0);
	97	}
	98
	99	for(k=0;k<vesicle->vlist->n; k++){
	100	cvtx=vesicle->vlist->vtx[k];
	101	cart2sph(&coord,cvtx->x,cvtx->y,cvtx->z);
	102	for(i=0;i<vesicle->sphHarmonics->l;i++){
	103	for(j=0;j<2*i+1;j++){
	104	m=j-i;
	105	l=i;
	106	if(m>=0){
	107	// fprintf(stderr, "Racunam za l=%d, m=%d\n", l,m);
	108	vesicle->sphHarmonics->ulmComplex[i][j]=gsl_complex_add(vesicle->sphHarmonics->ulmComplex[i][j], gsl_complex_conjugate(gsl_complex_mul_real(gsl_complex_polar(1.0,(ts_double)mcoord.e2),cvtx->solAnglecvtx->relR*gsl_sf_legendre_sphPlm(l,m,cos(coord.e3)))) );
	109	} else {
	110	// fprintf(stderr, "Racunam za l=%d, abs(m=%d)\n", l,m);
	111	vesicle->sphHarmonics->ulmComplex[i][j]=gsl_complex_add(vesicle->sphHarmonics->ulmComplex[i][j], gsl_complex_conjugate(gsl_complex_mul_real(gsl_complex_polar(1.0,(ts_double)mcoord.e2),cvtx->solAnglecvtx->relRpow(-1,m)gsl_sf_legendre_sphPlm(l,-m,cos(coord.e3)))) );
	112
	113	}
	114	}
	115	}
	116	}
	117	return TS_SUCCESS;
	118	}
	119
	120	ts_bool storeUlmComplex2(ts_vesicle *vesicle){
	121
	122	ts_spharm *sph=vesicle->sphHarmonics;
	123	ts_int i,j;
	124	for(i=0;i<sph->l;i++){
	125	for(j=0;j<2*i+1;j++){
	126	sph->sumUlm2[i][j]+=gsl_complex_abs2(sph->ulmComplex[i][j]);
	127	}
	128	}
	129	sph->N++;
	130	return TS_SUCCESS;
	131	}
	132
f06f5f	133
22cdfd	134	ts_double calculateKc(ts_vesicle *vesicle, ts_int lmin, ts_int lmax){
SP	135	ts_int min=lmin;
	136	ts_int max=lmax; //vesicle->sphHarmonics->l-3;
9f2ad6	137	ts_long i,j;
0ee39c	138	ts_double retval, bval;
9f2ad6	139	gsl_matrix *A=gsl_matrix_alloc(max-min,2);
f06f5f	140	gsl_vector *tau=gsl_vector_alloc(2);
9f2ad6	141	gsl_vector *b=gsl_vector_alloc(max-min);
f06f5f	142	gsl_vector *x=gsl_vector_alloc(2);
9f2ad6	143	gsl_vector *res=gsl_vector_alloc(max-min);
f06f5f	144
0ee39c	145	//solving (A^TA)x=A^T*b
f06f5f	146	//fill the data for matrix A and vector b
9f2ad6	147	for(i=min;i<max;i++){
SP	148	gsl_matrix_set(A, i-min,0,(ts_double)((i-1)*(i+2)));
	149	gsl_matrix_set(A, i-min,1,(ts_double)((i-1)(i+2)(i+1)*i));
22cdfd	150	// fprintf(stderr,"%e %e\n", gsl_matrix_get(A,i-min,0), gsl_matrix_get(A,i-min,1));
9f2ad6	151	bval=0.0;
0ee39c	152	//average for m from 0..l (only positive m's)
9f2ad6	153	for(j=0;j<=i;j++){
0ee39c	154	bval+=vesicle->sphHarmonics->sumUlm2[i][(j+i)];
SP	155	}
9f2ad6	156	bval=bval/(ts_double)vesicle->sphHarmonics->N/(ts_double)(i+1);
0ee39c	157
9f2ad6	158	gsl_vector_set(b,i-min,1.0/bval);
22cdfd	159	// fprintf(stderr,"%e\n", 1.0/gsl_vector_get(b,i-min));
f06f5f	160	}
0ee39c	161	// fprintf(stderr,"b[2]=%e\n",gsl_vector_get(b,1));
f06f5f	162	gsl_linalg_QR_decomp(A,tau);
SP	163	gsl_linalg_QR_lssolve(A,tau,b,x,res);
0ee39c	164	// fprintf(stderr,"kc=%e\n",gsl_vector_get(x,1));
SP	165	retval=gsl_vector_get(x,1);
f06f5f	166	gsl_matrix_free(A);
SP	167	gsl_vector_free(tau);
	168	gsl_vector_free(b);
	169	gsl_vector_free(x);
	170	gsl_vector_free(res);
0ee39c	171
SP	172	return retval;
f06f5f	173	}