diff --git a/math/function.h b/math/function.h
index 431cf4e..11655be 100644
--- a/math/function.h
+++ b/math/function.h
@@ -5,145 +5,166 @@
 
 namespace rts{
 
-template <class X, class Y>
+//template class for a one-dimensional function
+template <class Tx, class Ty>
 class function
 {
-	//datapoint class for storing function points
-	struct dataPoint
-	{
-		X x;
-		Y y;
-	};
-
-	//function data
-	std::vector<dataPoint> f;
+	std::vector<Tx> X;
+	std::vector<Ty> Y;
 
 	//comparison function for searching lambda
-    static bool findCeiling(dataPoint a, dataPoint b)
-    {
-        return (a.x > b.x);
-    }
+	static bool findCeiling(Tx ax, Tx bx){
+	    return (ax > bx);
+	}
 
+	//process a string to extract a function (generally delimited by tabs, commas, spaces, etc.)
+	void process_string(std::string s){
+		std::stringstream ss(s);
 
-public:
-	function()
-	{
-		//insert(0, 0);
-	}
+		//std::string test;
+		//std::getline(ss, test);
+		//std::cout<<test;
+		//exit(1);
 
-	Y linear(X x) const
-	{
-		if(f.size() == 0)	return (Y)0;	//return zero if the function is empty
-		//declare an iterator
-		typedef typename std::vector< dataPoint >::iterator f_iter;
-        f_iter it;
-
-		//dataPoint s;
-		//s.x = x;
-
-        //it = search(f.begin(), f.end(), &s, &s + 1, &function<X, Y>::findCeiling);
-        unsigned int i;
-        for(i = 0; i<f.size(); i++)
-        {
-        	if(f[i].x > x)
-        		break;
-        }
+		Tx x;
+		Ty y;
+		std::string line;
 
-        //if the wavelength is past the end of the list, return the back
-        if(i == f.size())
-            return f.back().y;
-        //if the wavelength is before the beginning of the list, return the front
-        else if(i == 0)
-            return f.front().y;
-        //otherwise interpolate
-        else
-        {
-            X xMax = f[i].x;
-            X xMin = f[i - 1].x;
-            //std::cout<<lMin<<"----------"<<lMax<<std::endl;
-
-            X a = (x - xMin) / (xMax - xMin);
-            Y riMin = f[i - 1].y;
-            Y riMax = f[i].y;
-            Y interp;
-            interp = riMax * a + riMin * (1.0 - a);
-            return interp;
-        }
+		while(!ss.eof()){
+			
+			std::getline(ss, line);
+			if(line[0] == '#') continue;
+
+			std::stringstream lstream(line);
+
+			lstream>>x;			//read the x value
+			lstream>>y;			//read the y value
+
+			//std::cout<<x<<", "<<y<<std::endl;
+
+			if(ss.eof()) break;
+			insert(x, y);			//insert the read value into the function
+			
+		}	
 	}
 
-	///add a data point to a function
-	void insert(X x, Y y)
+
+public:
+
+	//linear interpolation
+	Ty linear(Tx x) const
 	{
-		dataPoint s;
-		s.x = x;
-		s.y = y;
+		if(X.size() == 0)	return (Ty)0;	//return zero if the function is empty
 
-		if(f.size() == 0 || f.back().x < x)
-			return f.push_back(s);
+		unsigned int N = X.size();			//number of sample points
 
 		//declare an iterator
-        typename std::vector< dataPoint >::iterator it;
+		typedef typename std::vector< Tx >::iterator f_iter;	//declare an iterator
+		f_iter it;
+
+		//find the first X-coordinate that is greater than x
+		unsigned int i;
+		for(i = 0; i<N; i++){
+			if(X[i] > x)
+				break;
+		}
+		//i currently holds the ceiling
+
+		//if the wavelength is past the end of the list, return the last sample point
+		if(i == N) return Y[N];
+		//if the wavelength is before the beginning of the list, return the front
+		else if(i == 0) return Y[0];
+		//otherwise interpolate
+		else{
+		    Tx xMax = X[i];
+		    Tx xMin = X[i-1];
+
+		    Tx a = (x - xMin) / (xMax - xMin);
+		    Ty riMin = Y[i - 1];
+		    Ty riMax = Y[i];
+		    Ty interp = riMax * a + riMin * (1 - a);
+		    return interp;
+		}
+	}
+
+	///add a data point to a function
+	void insert(Tx x, Ty y)
+	{
+		unsigned int N = X.size();		//number of sample points
 
-		
+		if(N == 0 || X[N-1] < x){
+			X.push_back(x);
+			Y.push_back(y);
+			return;
+		}
 
-        it = search(f.begin(), f.end(), &s, &s + 1, &function<X, Y>::findCeiling);
+		//declare an iterator and search for the x value
+        typename std::vector< Tx >::iterator it;		
+        it = search(X.begin(), X.end(), &x, &x + 1, &function<Tx, Ty>::findCeiling);
 
         //if the function value is past the end of the vector, add it to the back
-        if(it == f.end())
-            return f.push_back(s);
+        if(*it == N){
+        	X.push_back(x);
+        	Y.push_back(y);
+        }
         //otherwise add the value at the iterator position
-        else
-		{
-			f.insert(it, s);
+        else{
+			X.insert(it, x);
+			Y.insert(Y.begin() + *it, y);
 		}
 
 	}
 
-	X getX(unsigned int i) const
-	{
-		return f[i].x;
+	Tx getX(unsigned int i) const{
+		return X[i];
 	}
 
-	Y getY(unsigned int i) const
-	{
-		return f[i].y;
+	Ty getY(unsigned int i) const{
+		return Y[i];
 	}
 
 	///get the number of data points in the function
-	unsigned int getN() const
-	{
-		return f.size();
+	unsigned int getN() const{
+		return X.size();
 	}
 
 	//look up an indexed component
-	dataPoint operator[](int i) const
-	{
-		return f[i];
+	Ty operator[](int i) const{
+		if(i <= X.size()){
+			std::cout<<"ERROR: accessing non-existing sample point in 'function'"<<std::endl;
+			exit(1);
+		}
+		return Y[i];
 	}
 
 	///linear interpolation
-	Y operator()(X x) const
-	{
+	Tx operator()(Tx x) const{
 		return linear(x);
 	}
 
-	function<X, Y> operator+(Y r) const
-	{
-		function<X, Y> result;
+	//add a constant to the function
+	function<Tx, Ty> operator+(Ty r) const{
+
+		function<Tx, Ty> result;
 
-		//add r to every point in f
-		for(int i=0; i<f.size(); i++)
-		{
-			result.f.push_back(f[i]);
-			result.f[i].y += r;
+		//if there are points in the function
+		if(X.size() > 0){
+			//add r to every point in f
+			for(unsigned int i=0; i<X.size(); i++){
+				result.X.push_back(X[i]);
+				result.Y.push_back(Y[i] + r);
+			}
+		}
+		else{
+			result = r;
 		}
 
 		return result;
 	}
 
-	function<X, Y> & operator= (const Y & rhs)
-	{		
-		f.clear();
+	function<Tx, Ty> & operator= (const Ty & rhs){		
+		X.clear();
+		Y.clear();
 		if(rhs != 0)			//if the RHS is zero, just clear, otherwise add one value of RHS
 			insert(0, rhs);
 
@@ -151,6 +172,29 @@ public:
 	}
 
 
+
+	std::string str(){
+		stringstream ss;
+
+		unsigned int N = X.size();		//number of sample points
+		//output each function value
+		for(unsigned int i = 0; i<N; i++){
+			ss<<X[i]<<", "<<Y[i]<<std::endl;
+		}
+
+		return ss.str();
+
+	}
+
+	void load(std::string filename){
+		std::ifstream t(filename.c_str());
+		std::string str((std::istreambuf_iterator<char>(t)),
+		std::istreambuf_iterator<char>());
+
+		process_string(str);
+	}
+
+
 };
 
 }	//end namespace rts
--
libgit2 0.21.4