#include<stdio.h>
|
#include<stdlib.h>
|
#include<math.h>
|
#include "general.h"
|
#include "sh.h"
|
#include "shreal.h"
|
#include "shcomplex.h"
|
#include "string.h"
|
#include "io.h"
|
#include<gsl/gsl_sf_legendre.h>
|
|
/* here we just initialize missing memory space and pointers that are not in complex sh initialization and are used to calculate real spherical harmonics */
|
ts_bool shreal_init(ts_spharm *sph, ts_uint l){
|
ts_uint j;
|
sph->ulmReal=(ts_double **)calloc(l,sizeof(ts_double *));
|
for(j=0;j<l;j++){
|
sph->ulmReal[j]=(ts_double *)calloc(2*j+1,sizeof(ts_double));
|
}
|
|
sph->sumUlm2Real=(ts_double **)calloc(l,sizeof(ts_double *));
|
for(j=0;j<l;j++){
|
sph->sumUlm2Real[j]=(ts_double *)calloc(2*j+1,sizeof(ts_double));
|
}
|
|
|
sph->sumUlm2Old=(ts_double **)calloc(l,sizeof(ts_double *));
|
for(j=0;j<l;j++){
|
sph->sumUlm2Old[j]=(ts_double *)calloc(2*j+1,sizeof(ts_double));
|
}
|
return TS_SUCCESS;
|
}
|
|
/* frees memory initialized in function above */
|
ts_bool shreal_free(ts_spharm *sph){
|
ts_uint i;
|
for(i=0;i<sph->l;i++){
|
if(sph->ulmReal[i]!=NULL) free(sph->ulmReal[i]);
|
if(sph->sumUlm2Real[i]!=NULL) free(sph->sumUlm2Real[i]);
|
if(sph->sumUlm2Old[i]!=NULL) free(sph->sumUlm2Old[i]);
|
}
|
return TS_SUCCESS;
|
}
|
|
ts_double factorial(ts_int n){
|
ts_int i;
|
ts_double fac=0.0;
|
if(n<0) fatal("Error calculation factorial",10);
|
for(i=1;i<=n;i++){
|
fac=fac+(ts_double)i;
|
}
|
if(n==0) fac=1.0;
|
return fac;
|
}
|
|
ts_double ZlmiCoeff(ts_int l, ts_int m){
|
ts_double coeff=(ts_double)(2*l+1)/2.0*M_PI*(factorial(l-m)/factorial(l+m));
|
return sqrt(coeff);
|
}
|
|
ts_bool calculateUlmReal(ts_vesicle *vesicle){
|
ts_int i,j,k;
|
ts_vertex *cvtx;
|
ts_double Zlmi;
|
ts_int m;
|
ts_coord *coord=(ts_coord *)malloc(sizeof(ts_coord));
|
ts_double phi, theta;
|
for(i=0;i<vesicle->sphHarmonics->l;i++){
|
for(j=0;j<2*i+1;j++) vesicle->sphHarmonics->ulmReal[i][j]=0.0;
|
}
|
|
|
for(k=0;k<vesicle->vlist->n; k++){
|
cvtx=vesicle->vlist->vtx[k];
|
cart2sph(coord,cvtx->x, cvtx->y, cvtx->z);
|
phi=coord->e2;
|
theta=coord->e3;
|
for(i=0;i<vesicle->sphHarmonics->l;i++){
|
for(j=0;j<2*i+1;j++){
|
m=(j-i);
|
if(m<0){ //negative m
|
// gsl_sf_legendre_sphPlm(l,-m,cos(coord.e3)))
|
Zlmi=gsl_sf_legendre_sphPlm(i,-m,cos(theta))*cos(m*phi)*sqrt(2);
|
} else if(m==0){ //m=0
|
Zlmi=gsl_sf_legendre_sphPlm(i,0,cos(theta));
|
} else { //positive m
|
Zlmi=gsl_sf_legendre_sphPlm(i,m,cos(theta))*sin(m*phi)*sqrt(2);
|
// Zlmi=ZlmiCoeff(i,m)*plgndr(i,m,cos(theta))*cos(m*phi);
|
}
|
vesicle->sphHarmonics->ulmReal[i][j]+= cvtx->solAngle*cvtx->relR*Zlmi;
|
}
|
|
}
|
}
|
return TS_SUCCESS;
|
}
|
|
|
ts_bool storeUlm2Real(ts_vesicle *vesicle){
|
|
ts_spharm *sph=vesicle->sphHarmonics;
|
ts_int i,j;
|
for(i=0;i<sph->l;i++){
|
for(j=0;j<2*i+1;j++){
|
/* DEBUG fprintf(stderr,"sph->sumUlm2[%d][%d]=%e\n",i,j,sph->ulm[i][j]* sph->ulm[i][j]); */
|
sph->sumUlm2Real[i][j]+=sph->ulmReal[i][j]* sph->ulmReal[i][j];
|
}
|
}
|
sph->N++;
|
return TS_SUCCESS;
|
}
|
|
|
ts_bool saveAvgUlm2Real(ts_vesicle *vesicle){
|
|
FILE *fh;
|
char filename[10000];
|
strcpy(filename, command_line_args.path);
|
strcat(filename, "sph2outReal.dat");
|
fh=fopen(filename, "w");
|
if(fh==NULL){
|
err("Cannot open file %s for writing");
|
return TS_FAIL;
|
}
|
|
ts_spharm *sph=vesicle->sphHarmonics;
|
ts_int i,j;
|
fprintf(fh,"l,\tm,\tulm^2avg\n");
|
for(i=0;i<sph->l;i++){
|
for(j=0;j<2*i+1;j++){
|
fprintf(fh,"%d,\t%d,\t%e\n", i, j-i, sph->sumUlm2Real[i][j]/(ts_double)sph->N);
|
|
}
|
fprintf(fh,"\n");
|
}
|
fclose(fh);
|
return TS_SUCCESS;
|
}
|
