Commit 9bfae4d0ec70ed9fb1a34fd3a2b0d14ecc396bea

Authored by Pavel Govyadinov
1 parent 761ebaa9

extended gl_spharmonics

stim/math/spharmonics.h
1 1 #ifndef STIM_SPH_HARMONICS
2 2 #define STIM_SPH_HARMONICS
3 3  
  4 +#include <complex>
4 5 #include <stim/math/vector.h>
5 6 #include <boost/math/special_functions/spherical_harmonic.hpp>
  7 +#include <stim/math/constants.h>
6 8 #include <vector>
7 9  
8 10 #define WIRE_SCALE 1.001
... ... @@ -11,16 +13,20 @@ namespace stim{
11 13 template<class T>
12 14 class spharmonics{
13 15  
  16 +public:
  17 + std::vector<T> C; //list of SH coefficients
  18 +
14 19 protected:
15 20  
16   - std::vector<T> C; //list of SH coefficients
17 21  
18 22 unsigned int mcN; //number of Monte-Carlo samples
19 23  
20 24 //calculate the value of the SH basis function (l, m) at (theta, phi)
21 25 //here, theta = [0, PI], phi = [0, 2*PI]
22 26 double SH(int l, int m, double theta, double phi){
23   - return boost::math::spherical_harmonic_r(l, m, phi, theta);
  27 + std::complex<T> result = boost::math::spherical_harmonic(l, m, phi, theta);
  28 + return result.imag() + result.real();
  29 +// return boost::math::spherical_harmonic_i(l, m, phi, theta);
24 30 }
25 31  
26 32 unsigned int coeff_1d(unsigned int l, int m){
... ... @@ -55,6 +61,16 @@ public:
55 61 C[c] = value;
56 62 }
57 63  
  64 + unsigned int getSize() const{
  65 + return C.size();
  66 + }
  67 +
  68 + std::vector<T> getC() const{
  69 + return C;
  70 + }
  71 +
  72 +
  73 +
58 74 /// Initialize Monte-Carlo sampling of a function using N spherical harmonics coefficients
59 75  
60 76 /// @param N is the number of spherical harmonics coefficients used to represent the user function
... ... @@ -167,7 +183,7 @@ public:
167 183  
168 184 return fx;
169 185 }
170   -
  186 +/*
171 187 //overload arithmetic operations
172 188 spharmonics<T> operator*(T rhs) const {
173 189 spharmonics<T> result(C.size()); //create a new spherical harmonics object
... ... @@ -199,7 +215,7 @@ public:
199 215 spharmonics<T> operator-(spharmonics<T> rhs) {
200 216 return (*this) + (rhs * (T)(-1));
201 217 }
202   -
  218 +*/
203 219 }; //end class sph_harmonics
204 220  
205 221  
... ...
stim/visualization/gl_spharmonics.h
... ... @@ -190,6 +190,40 @@ public:
190 190 gen_function();
191 191 }
192 192  
  193 + //overload arithmetic operations
  194 + gl_spharmonics<T> operator*(T rhs) const {
  195 + gl_spharmonics<T> result(C.size()); //create a new sph erical harmonics object
  196 + for (size_t c = 0; c < C.size(); c++) //for each coefficient
  197 + result.C[c] = C[c] * rhs; //calculat e the factor and store the result in the new spharmonics object
  198 + return result;
  199 + }
  200 +
  201 + gl_spharmonics<T> operator+(gl_spharmonics<T> rhs) {
  202 + size_t low = std::min(C.size(), rhs.C.size()); //store th e number of coefficients in the lowest object
  203 + size_t high = std::max(C.size(), rhs.C.size()); //store th e number of coefficients in the result
  204 + bool rhs_lowest = false; //true if rhs has the lowest number of coefficients
  205 + if (rhs.C.size() < C.size()) rhs_lowest = true; //if rhs has a low er number of coefficients, set the flag
  206 +
  207 + gl_spharmonics<T> result(high); //create a new object
  208 + size_t c;
  209 + for (c = 0; c < low; c++) //perform the first batch of additions
  210 + result.C[c] = C[c] + rhs.C[c]; // perform the addition
  211 +
  212 + for (c = low; c < high; c++) {
  213 + if (rhs_lowest)
  214 + result.C[c] = C[c];
  215 + else
  216 + result.C[c] = rhs.C[c];
  217 + }
  218 + return result;
  219 + }
  220 +
  221 + gl_spharmonics<T> operator-(gl_spharmonics<T> rhs) {
  222 + return (*this) + (rhs * (T)(-1));
  223 + }
  224 +
  225 +
  226 +
193 227 }; //end gl_spharmonics
194 228  
195 229  
... ...