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