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