diff --git a/stim/math/spharmonics.h b/stim/math/spharmonics.h
index defb374..607ac43 100644
--- a/stim/math/spharmonics.h
+++ b/stim/math/spharmonics.h
@@ -1,8 +1,10 @@
 #ifndef STIM_SPH_HARMONICS
 #define STIM_SPH_HARMONICS
 
+#include <complex>
 #include <stim/math/vector.h>
 #include <boost/math/special_functions/spherical_harmonic.hpp>
+#include <stim/math/constants.h>
 #include <vector>
 
 #define WIRE_SCALE 1.001
@@ -11,16 +13,20 @@ namespace stim{
 template<class T>
 class spharmonics{
 
+public:
+	std::vector<T> C;	//list of SH coefficients
+
 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);
+		std::complex<T> 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){
@@ -55,6 +61,16 @@ public:
 		C[c] = value;
 	}
 
+	unsigned int getSize() const{
+		return C.size();
+	}
+
+	std::vector<T> getC() const{
+		return C;
+	}
+
+	
+
 	/// 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
@@ -167,7 +183,7 @@ public:
 
 		return fx;
 	}
-
+/*
 	//overload arithmetic operations
 	spharmonics<T> operator*(T rhs) const {
 		spharmonics<T> result(C.size());				//create a new spherical harmonics object
@@ -199,7 +215,7 @@ public:
 	spharmonics<T> operator-(spharmonics<T> rhs) {
 		return (*this) + (rhs * (T)(-1));
 	}
-
+*/
 };		//end class sph_harmonics
 
 
diff --git a/stim/visualization/gl_spharmonics.h b/stim/visualization/gl_spharmonics.h
index 37b9a2e..921cc50 100644
--- a/stim/visualization/gl_spharmonics.h
+++ b/stim/visualization/gl_spharmonics.h
@@ -190,6 +190,40 @@ public:
 		gen_function();
 	}
 
+	//overload arithmetic operations
+	gl_spharmonics<T> operator*(T rhs) const {
+		gl_spharmonics<T> result(C.size());                                //create a new sph    erical harmonics object
+		for (size_t c = 0; c < C.size(); c++)                   //for each coefficient
+			result.C[c] = C[c] * rhs;                                       //calculat    e the factor and store the result in the new spharmonics object
+		return result;
+	}
+
+	gl_spharmonics<T> operator+(gl_spharmonics<T> rhs) {
+		size_t low = std::min(C.size(), rhs.C.size());                          //store th    e number of coefficients in the lowest object
+		size_t high = std::max(C.size(), rhs.C.size());                         //store th    e number of coefficients in the result
+		bool rhs_lowest = false;                                                                    //true if rhs has the lowest number of coefficients
+		if (rhs.C.size() < C.size()) rhs_lowest = true;                 //if rhs has a low    er number of coefficients, set the flag
+
+		gl_spharmonics<T> result(high);                                                                //create a new object
+		size_t c;
+		for (c = 0; c < low; c++)                                               //perform     the first batch of additions
+			result.C[c] = C[c] + rhs.C[c];                                          //    perform the addition
+
+		for (c = low; c < high; c++) {
+			if (rhs_lowest)
+				result.C[c] = C[c];
+			else
+				result.C[c] = rhs.C[c];
+		}
+		return result;
+	}
+
+	gl_spharmonics<T> operator-(gl_spharmonics<T> rhs) {
+		return (*this) + (rhs * (T)(-1));
+	}
+
+
+
 };		//end gl_spharmonics
 
 
--
libgit2 0.21.4