spharmonics.h 3.32 KB
``````#ifndef STIM_SPH_HARMONICS
#define STIM_SPH_HARMONICS

#include <stim/math/vector.h>
#include <boost/math/special_functions/spherical_harmonic.hpp>
#include <vector>

#define PI 3.14159
#define WIRE_SCALE 1.001
namespace stim{

template<class T>
class spharmonics{

protected:

std::vector<T> 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<stim::vec<double> > 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<<C[i]<<'\t';

m++;			//increment m

//if we're in a new tier, increment l and set m = -l
if(m > l){
l++;
m = -l;

ss<<std::endl;

}
}

return ss.str();

}

/// Returns the value of the function at the coordinate (theta, phi)

/// @param theta = [0, 2pi]
/// @param phi = [0, pi]
double operator()(double theta, double phi){

double fx = 0;

int l = 0;
int m = 0;
for(unsigned int i = 0; i < C.size(); i++){
fx += C[i] * SH(l, m, theta, phi);
m++;
if(m > l){
l++;
m = -l;
}

}

return fx;
}

};		//end class sph_harmonics

}

#endif
``````