#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;
}