| commit | author | age | ||
| fda1ab | 1 | #include<stdio.h> |
| SP | 2 | #include<stdlib.h> |
| 3 | #include<math.h> | |
| 4 | #include "general.h" | |
| 5 | #include "sh.h" | |
| 6 | #include "shreal.h" | |
| 7 | #include "shcomplex.h" | |
| 8 | #include "string.h" | |
| 9 | #include "io.h" | |
| 10 | #include<gsl/gsl_sf_legendre.h> | |
| 11 | ||
| 12 | /* here we just initialize missing memory space and pointers that are not in complex sh initialization and are used to calculate real spherical harmonics */ | |
| 13 | ts_bool shreal_init(ts_spharm *sph, ts_uint l){ | |
| 14 | ts_uint j; | |
| 15 | sph->ulmReal=(ts_double **)calloc(l,sizeof(ts_double *)); | |
| 16 | for(j=0;j<l;j++){ | |
| 17 | sph->ulmReal[j]=(ts_double *)calloc(2*j+1,sizeof(ts_double)); | |
| 18 | } | |
| 19 | ||
| 20 | sph->sumUlm2Real=(ts_double **)calloc(l,sizeof(ts_double *)); | |
| 21 | for(j=0;j<l;j++){ | |
| 22 | sph->sumUlm2Real[j]=(ts_double *)calloc(2*j+1,sizeof(ts_double)); | |
| 23 | } | |
| 24 | ||
| 25 | ||
| 26 | sph->sumUlm2Old=(ts_double **)calloc(l,sizeof(ts_double *)); | |
| 27 | for(j=0;j<l;j++){ | |
| 28 | sph->sumUlm2Old[j]=(ts_double *)calloc(2*j+1,sizeof(ts_double)); | |
| 29 | } | |
| 30 | return TS_SUCCESS; | |
| 31 | } | |
| 32 | ||
| 33 | /* frees memory initialized in function above */ | |
| 34 | ts_bool shreal_free(ts_spharm *sph){ | |
| 35 | ts_uint i; | |
| 36 | for(i=0;i<sph->l;i++){ | |
| 37 | if(sph->ulmReal[i]!=NULL) free(sph->ulmReal[i]); | |
| 38 | if(sph->sumUlm2Real[i]!=NULL) free(sph->sumUlm2Real[i]); | |
| 39 | if(sph->sumUlm2Old[i]!=NULL) free(sph->sumUlm2Old[i]); | |
| 40 | } | |
| 41 | return TS_SUCCESS; | |
| 42 | } | |
| 43 | ||
| 44 | ts_double factorial(ts_int n){ | |
| 45 | ts_int i; | |
| 46 | ts_double fac=0.0; | |
| 47 | if(n<0) fatal("Error calculation factorial",10); | |
| 48 | for(i=1;i<=n;i++){ | |
| 49 | fac=fac+(ts_double)i; | |
| 50 | } | |
| 51 | if(n==0) fac=1.0; | |
| 52 | return fac; | |
| 53 | } | |
| 54 | ||
| 55 | ts_double ZlmiCoeff(ts_int l, ts_int m){ | |
| 56 | ts_double coeff=(ts_double)(2*l+1)/2.0*M_PI*(factorial(l-m)/factorial(l+m)); | |
| 57 | return sqrt(coeff); | |
| 58 | } | |
| 59 | ||
| 60 | ts_bool calculateUlmReal(ts_vesicle *vesicle){ | |
| 61 | ts_int i,j,k; | |
| 62 | ts_vertex *cvtx; | |
| 63 | ts_double Zlmi; | |
| 64 | ts_int m; | |
| 65 | ts_coord *coord=(ts_coord *)malloc(sizeof(ts_coord)); | |
| 66 | ts_double phi, theta; | |
| 67 | for(i=0;i<vesicle->sphHarmonics->l;i++){ | |
| 68 | for(j=0;j<2*i+1;j++) vesicle->sphHarmonics->ulmReal[i][j]=0.0; | |
| 69 | } | |
| 70 | ||
| 71 | ||
| 72 | for(k=0;k<vesicle->vlist->n; k++){ | |
| 73 | cvtx=vesicle->vlist->vtx[k]; | |
| 74 | cart2sph(coord,cvtx->x, cvtx->y, cvtx->z); | |
| 75 | phi=coord->e2; | |
| 76 | theta=coord->e3; | |
| 77 | for(i=0;i<vesicle->sphHarmonics->l;i++){ | |
| 78 | for(j=0;j<2*i+1;j++){ | |
| 79 | m=(j-i); | |
| 80 | if(m<0){ //negative m | |
| 81 | // gsl_sf_legendre_sphPlm(l,-m,cos(coord.e3))) | |
| 82 | Zlmi=gsl_sf_legendre_sphPlm(i,-m,cos(theta))*cos(m*phi)*sqrt(2); | |
| 83 | } else if(m==0){ //m=0 | |
| 84 | Zlmi=gsl_sf_legendre_sphPlm(i,0,cos(theta)); | |
| 85 | } else { //positive m | |
| 86 | Zlmi=gsl_sf_legendre_sphPlm(i,m,cos(theta))*sin(m*phi)*sqrt(2); | |
| 87 | // Zlmi=ZlmiCoeff(i,m)*plgndr(i,m,cos(theta))*cos(m*phi); | |
| 88 | } | |
| 89 | vesicle->sphHarmonics->ulmReal[i][j]+= cvtx->solAngle*cvtx->relR*Zlmi; | |
| 90 | } | |
| 91 | ||
| 92 | } | |
| 93 | } | |
| 94 | return TS_SUCCESS; | |
| 95 | } | |
| 96 | ||
| 97 | ||
| 98 | ts_bool storeUlm2Real(ts_vesicle *vesicle){ | |
| 99 | ||
| 100 | ts_spharm *sph=vesicle->sphHarmonics; | |
| 101 | ts_int i,j; | |
| 102 | for(i=0;i<sph->l;i++){ | |
| 103 | for(j=0;j<2*i+1;j++){ | |
| 104 | /* DEBUG fprintf(stderr,"sph->sumUlm2[%d][%d]=%e\n",i,j,sph->ulm[i][j]* sph->ulm[i][j]); */ | |
| 105 | sph->sumUlm2Real[i][j]+=sph->ulmReal[i][j]* sph->ulmReal[i][j]; | |
| 106 | } | |
| 107 | } | |
| 108 | sph->N++; | |
| 109 | return TS_SUCCESS; | |
| 110 | } | |
| 111 | ||
| 112 | ||
| 113 | ts_bool saveAvgUlm2Real(ts_vesicle *vesicle){ | |
| 114 | ||
| 115 | FILE *fh; | |
| 116 | char filename[10000]; | |
| 117 | strcpy(filename, command_line_args.path); | |
| 118 | strcat(filename, "sph2outReal.dat"); | |
| 119 | fh=fopen(filename, "w"); | |
| 120 | if(fh==NULL){ | |
| 121 | err("Cannot open file %s for writing"); | |
| 122 | return TS_FAIL; | |
| 123 | } | |
| 124 | ||
| 125 | ts_spharm *sph=vesicle->sphHarmonics; | |
| 126 | ts_int i,j; | |
| 127 | fprintf(fh,"l,\tm,\tulm^2avg\n"); | |
| 128 | for(i=0;i<sph->l;i++){ | |
| 129 | for(j=0;j<2*i+1;j++){ | |
| 130 | fprintf(fh,"%d,\t%d,\t%e\n", i, j-i, sph->sumUlm2Real[i][j]/(ts_double)sph->N); | |
| 131 | ||
| 132 | } | |
| 133 | fprintf(fh,"\n"); | |
| 134 | } | |
| 135 | fclose(fh); | |
| 136 | return TS_SUCCESS; | |
| 137 | } | |
