| 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(2*i+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(2*j+1,sizeof(ts_double)); | |
| 35 | sph->ulmComplex[j]=(gsl_complex *)calloc(2*j+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(2*j+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 2*j+2 instead of 2*j+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(2*j+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)m*coord.e2),cvtx->solAngle*cvtx->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)m*coord.e2),cvtx->solAngle*cvtx->relR*pow(-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^T*A)*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 | } |
