diff --git a/stim/math/random.h b/stim/math/random.h new file mode 100644 index 0000000..b79c31d --- /dev/null +++ b/stim/math/random.h @@ -0,0 +1,114 @@ +#ifndef STIM_RANDOM +#define STIM_RANDOM + +#include +#include +#include +#include +#include + +namespace stim{ + +template +class Random{ +protected: + void init(int seed = 0) + { + if(seed <= 0) + srand(time(NULL)); + else if(seed > 0) + srand(time(seed)); + else + std::cout << "Error: Unknown value: in STIM::RANDOM" << std::endl; + } + +public: + /// Default Constructor + Random(){ + init(-1); + } + + /// Constructor from a seed. + /// A positive seed sets, 0 or negative yeilds the + Random(int seed){ + init(seed); + } + + ///Returns a random number uniformly sampled between 0 and 1 + static + T uniformRandom() + { + return ( (T)(rand()))/( (T)(RAND_MAX)); ///generates a random number between 0 and 1 using the uniform distribution. + } + + ///Returns a random number from a normal distribution between 0 to 1. + static + T normalRandom() + { + T u1 = uniformRandom(); + T u2 = uniformRandom(); + return cos(2.0*atan(1.0)*u2)*sqrt(-1.0*log(u1)); ///generate a random number using the normal distribution between 0 and 1. + } + ///Return a random vec3 each value between 0 and 1 from a uniform distribution. + static + stim::vec3 uniformRandVector() + { + stim::vec3 r(uniformRandom(), uniformRandom(), 1.0); ///generate a random vector using the uniform distribution between 0 and 1. + return r; + } + ///Return a random vec3, each value between 0 and 1 from a normal distribution. + static + stim::vec3 normalRandVector() + { + stim::vec3 r(normalRandom(), normalRandom(), 1.0); ///generate a random vector using the normal distribution between 0 and 1. + return r; + } + + ///place num_samples of samples on the surface of a sphere of radius r. + ///returns an std::vector of vec3's in cartisian coordinates. + static std::vector > + sample_sphere(unsigned int num_samples, T radius = 1, T solidAngle = stim::TAU) + { + std::cout << "did this" << std::endl; + T PHI[2], Z[2], range; ///Range of angles in cylinderical coordinates + PHI[0] = solidAngle/2; ///project the solid angle into spherical coords + PHI[1] = asin(0); /// + Z[0] = cos(PHI[0]); ///project the z into spherical coordinates + Z[1] = cos(PHI[1]); /// + range = Z[0] - Z[1]; ///the range of all possible z values. + + T z, theta, phi; /// temporary individual + + std::vector > samples; + + //srand(100); ///set random seed + + for(int i = 0; i < num_samples; i++) ///for each sample + { + z = uniformRandom()*range + Z[1]; ///find a random z based on the solid angle + theta = uniformRandom() * stim::TAU; ///find theta + phi = acos(z); ///project into spherical coord phi + stim::vec3 sph(radius, theta, phi); ///assume spherical + stim::vec3 cart = sph.sph2cart(); ///conver to cartesisn + samples.push_back(cart); ///push into list + } + +///THIS IS DEBUGGING CODE, UNCOMMENT TO CHECK WHETHER THE SURFACE IS WELL SAMPLED! +/* + std::stringstream name; + for(int i = 0; i < num_samples; i++) + name << samples[i].str() << std::endl; + + + std::ofstream outFile; + outFile.open("Sampled Surface.txt"); + outFile << name.str().c_str(); +*/ + + return samples; ///return full list. + } +}; + +} + +#endif diff --git a/stim/math/spharmonics.h b/stim/math/spharmonics.h index 607ac43..266af7d 100644 --- a/stim/math/spharmonics.h +++ b/stim/math/spharmonics.h @@ -5,6 +5,7 @@ #include #include #include +#include #include #define WIRE_SCALE 1.001 @@ -26,7 +27,6 @@ protected: double SH(int l, int m, double theta, double phi){ std::complex result = boost::math::spherical_harmonic(l, m, phi, theta); return result.imag() + result.real(); -// return boost::math::spherical_harmonic_i(l, m, phi, theta); } unsigned int coeff_1d(unsigned int l, int m){ @@ -117,10 +117,10 @@ public: } /// 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 @@ -134,6 +134,42 @@ public: mcEnd(); } + /// Generates a PDF describing the probability distribution of points on a spherical surface + /// @param sph_pts is a list of points in cartesian coordinates + /// @param l is the maximum degree of the spherical harmonic function + /// @param m is the maximum order + /// @param c is the centroid of the points in sph_pts. DEFAULT 0,0,0 + /// @param n is the number of points of the surface of the sphere used to create the PDF. DEFAULT 1000 + void pdf(std::vector > sph_pts, unsigned int l, int m, stim::vec3 c = stim::vec3(0,0,0), unsigned int n = 1000) + { + std::vector weights; ///the weight at each point on the surface of the sphere. +// weights.resize(n); + unsigned int nP = sph_pts.size(); + std::vector > sphere = stim::Random::sample_sphere(n, 1.0, stim::TAU); + for(int i = 0; i < n; i++) + { + double val = 0; + for(int j = 0; j < nP; j++) + { + stim::vec3 temp = sph_pts[j] - c; + if(temp.dot(sphere[i]) > 0) + val += pow(temp.dot(sphere[i]),4); + } + weights.push_back(val); + } + + + mcBegin(l, m); //begin spherical harmonic sampling + + for(unsigned int i = 0; i < n; i++) + { + stim::vec3 sph = sphere[i].cart2sph(); + mcSample(sph[1], sph[2], weights[i]); + } + + mcEnd(); + } + std::string str(){ std::stringstream ss; -- libgit2 0.21.4