| | |
| | | /* 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> |
| | |
| | | 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){ |
| | |
| | | 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; |
| | |
| | | } |
| | | |
| | | |
| | | ts_double calculateKc(ts_vesicle *vesicle){ |
| | | |
| | | ts_int i,j; |
| | | 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(vesicle->sphHarmonics->l-1,2); |
| | | 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-1); |
| | | gsl_vector *b=gsl_vector_alloc(max-min); |
| | | gsl_vector *x=gsl_vector_alloc(2); |
| | | gsl_vector *res=gsl_vector_alloc(vesicle->sphHarmonics->l-1); |
| | | gsl_vector *res=gsl_vector_alloc(max-min); |
| | | |
| | | //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)); |
| | | |
| | | bval=0; |
| | | 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++){ |
| | | for(j=0;j<=i;j++){ |
| | | bval+=vesicle->sphHarmonics->sumUlm2[i][(j+i)]; |
| | | } |
| | | bval=bval/(ts_double)vesicle->sphHarmonics->N/(ts_double)i; |
| | | bval=bval/(ts_double)vesicle->sphHarmonics->N/(ts_double)(i+1); |
| | | |
| | | gsl_vector_set(b,i-1,1.0/bval); |
| | | 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)); |
| | | gsl_linalg_QR_decomp(A,tau); |