#ifndef STIM_SPH_HARMONICS #define STIM_SPH_HARMONICS #include #include #include #define PI 3.14159 #define WIRE_SCALE 1.001 namespace stim{ template class spharmonics{ protected: std::vector C; //list of SH coefficients unsigned int mcN; //number of Monte-Carlo samples //calculate the value of the SH basis function (l, m) at (theta, phi) //here, theta = [0, PI], phi = [0, 2*PI] double SH(int l, int m, double theta, double phi){ return boost::math::spherical_harmonic_r(l, m, phi, theta); } unsigned int coeff_1d(unsigned int l, int m){ return pow(l + 1, 2) - (l - m) - 1; } public: void push(double c){ C.push_back(c); } void resize(unsigned int n){ C.resize(n); } void setc(unsigned int l, int m, T value){ unsigned int c = coeff_1d(l, m); C[c] = value; } void setc(unsigned int c, T value){ C[c] = value; } /// Initialize Monte-Carlo sampling of a function using N spherical harmonics coefficients /// @param N is the number of spherical harmonics coefficients used to represent the user function void mcBegin(unsigned int coefficients){ C.resize(coefficients, 0); mcN = 0; } void mcBegin(unsigned int l, int m){ unsigned int c = pow(l + 1, 2) - (l - m); mcBegin(c); } void mcSample(double theta, double phi, double val){ int l, m; double sh; l = m = 0; for(unsigned int i = 0; i < C.size(); i++){ sh = SH(l, m, theta, phi); C[i] += sh * val; m++; //increment m //if we're in a new tier, increment l and set m = -l if(m > l){ l++; m = -l; } } //end for all coefficients //increment the number of samples mcN++; } //end mcSample() void mcEnd(){ //divide all coefficients by the number of samples for(unsigned int i = 0; i < C.size(); i++) C[i] /= mcN; } /// Generates a PDF describing the probability distribution of points on a spherical surface /// @param sph_pts is a list of points in spherical coordinates (theta, phi) where theta = [0, 2pi] and phi = [0, pi] /// @param l is the maximum degree of the spherical harmonic function /// @param m is the maximum order void pdf(std::vector > sph_pts, unsigned int l, int m){ mcBegin( l, m ); //begin spherical harmonic sampling unsigned int nP = sph_pts.size(); for(unsigned int p = 0; p < nP; p++){ mcSample(sph_pts[p][1], sph_pts[p][2], 1.0); } mcEnd(); } std::string str(){ std::stringstream ss; int l, m; l = m = 0; for(unsigned int i = 0; i < C.size(); i++){ ss< l){ l++; m = -l; ss< l){ l++; m = -l; } } return fx; } }; //end class sph_harmonics } #endif