diff --git a/envi/EnviClose.h b/envi/EnviClose.h
new file mode 100644
index 0000000..bb21076
--- /dev/null
+++ b/envi/EnviClose.h
@@ -0,0 +1,58 @@
+#ifndef ENVICLOSE_H
+#define ENVICLOSE_H
+
+
+#include <fstream>
+#include <stdio.h>
+#include "EnviHeader.h"
+#include "EnviFid.h"
+#include "EnviOpen.h"
+
+// this is to close all the opened files and to create a header file for the written files
+
+static void EnviClose(EnviFidClass &fid, int dimX=0)
+{
+
+    float* zeroPtr = (float*) malloc(sizeof(float) * fid.header.nF);
+    memset(zeroPtr, 0, sizeof(float) * fid.header.nF);
+
+	if(fid.mode=="r")
+	{
+		fclose(fid.file);
+	}
+
+	else
+	{
+        int totalPixels = fid.header.sX;
+        int dimY = totalPixels/dimX;
+
+        if(totalPixels % dimX > 0)
+        {
+            dimY++;
+
+            int padding = dimY*dimX-totalPixels;
+
+            for(int p=0; p < padding; p++)
+                fwrite(zeroPtr, sizeof(float), fid.header.nF, fid.file);
+        }
+
+        //write the calculated values to the header object
+        fid.header.sX = dimX;
+        fid.header.sY = dimY;
+
+
+        //write the header file to disk
+        //TODO
+
+
+        //close the binary file
+		fclose(fid.file);
+
+	}
+
+
+
+
+}
+#endif
+
diff --git a/envi/EnviFid.h b/envi/EnviFid.h
new file mode 100644
index 0000000..57d9a4f
--- /dev/null
+++ b/envi/EnviFid.h
@@ -0,0 +1,47 @@
+#ifndef ENVIFID_H
+#define ENVIFID_H
+
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <ctime>
+#include <sstream>
+#include "EnviHeader.h"
+using namespace std;
+
+class EnviFidClass
+{
+public:
+    //binary file name
+    string fileName;
+
+    //header structure
+	EnviHeaderClass header;
+
+	//valid flag (is the current file valid)
+	int isValid;
+
+	//file mask (NULL if the whole file is to be read)
+	char* mask;
+
+	//current position (in pixels) in the file
+	int filePos;
+
+	//read/write mode ("r" or "w")
+	string mode;
+
+	//create a variable to store a pointer to the binary file
+	FILE * file;
+
+    //constructor, set default values
+	EnviFidClass()
+	{
+        isValid = 0;
+        filePos = 0;
+    }
+
+
+};
+
+
+#endif
diff --git a/envi/EnviHeader.h b/envi/EnviHeader.h
new file mode 100644
index 0000000..b0c93e8
--- /dev/null
+++ b/envi/EnviHeader.h
@@ -0,0 +1,126 @@
+#ifndef ENVIHEADER_H
+#define ENVIHEADER_H
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <ctime>
+#include <sstream>
+using namespace std;
+
+enum InterleaveType {bip, bil, bsq};
+class EnviHeaderClass	//"EnviheaderClass" will be used to store data into variables sX, sY, etc.
+{
+public:
+	unsigned int sX;				//samples = "sX"
+	unsigned int sY;				//lines = "sY"
+	unsigned int nF;				//bands = "nF"
+	unsigned int headersize;		//header offset = "headersize"
+	unsigned int datatype;		//data type = "datatype"
+	InterleaveType type;
+	string units;		//wavelength units = "units" (Usully in Wavenumber or Wavelength)
+	int LoadHeader(string filename)
+	{
+        ifstream inFile;
+        //Load the header file in the same directory, for future note, might want to modify code to either "cin" the file name of have a feature in the GUI which take into account different fileneames
+        inFile.open (filename.c_str());
+
+        //If or not there is a file
+        if (!inFile)
+        {
+            cout<<"Error"<<endl;
+            return 1;
+        }
+
+        string line;
+        getline(inFile,line);
+
+        //If the 1st line of text file does not read ENVI, spit out an error.
+        if (line.compare("ENVI"))
+        {
+            cout<<"Error, not an ENVI file"<<endl;
+            return 1;
+        }
+
+        while(inFile) //Loops through entire file
+        {
+                getline(inFile,line);
+                if(line.find("samples")==0)
+                {
+                    sX=getNumber(line);				//Calls getNumber function for "samples", stores in header.sX
+                }
+                if(line.find("lines")==0)
+                {
+                    sY=getNumber(line);				//Calls getNumber function for "lines", stores in header.sY
+                }
+                if(line.find("bands")==0)
+                {
+                    nF=getNumber(line);				//Calls getNumber function for "bands", stores in header.nF
+                }
+                if(line.find("header offset")==0)
+                {
+                    headersize=getNumber(line);		//Calls getNumber function for "header offset", stores in header.headersize
+                }
+                if(line.find("data type")==0)
+                {
+                    datatype=getNumber(line);		//Calls getNumber function for "data type", stores in header.datatype
+                }
+
+                line.find("interleave");
+                if(line.find("interleave")==0)
+                {
+                    string typeStr;
+                    typeStr=getString(line);
+                    if (typeStr.find("bip")==0)
+                        type=bip;
+                    else if(typeStr.find("bsq")==0)
+                        type=bsq;
+                    else if(typeStr.find("bil")==0)
+                        type=bil;
+                    else
+                    {
+                        //throw an error and exit
+                        cout<<"Type not recognized"<<endl;
+                        return 1;
+                    }
+
+                }
+                if(line.find("wavelength units")==0)
+                {
+                    //Calls getString function for "wavelength units", stores in header.units
+                    units=getString(line);
+                }
+            }
+        return 0;
+    }
+
+
+private:
+    //retrieves a numerical parameter from a line of text
+	int getNumber(string line)
+	{
+        int result=0;
+        int position;
+        position=line.find("="); //Finds the position of the equal sign
+
+        string numStr; //Creates a string for storing data called "numStr"
+        numStr=line.substr(position+2); //Moves from the position of "=" 2 spaces to the right...
+        istringstream convert(numStr);  //...and converts the string into numbers using "istringstream"
+
+        convert>>result; //converts result from previous line
+
+        return result;
+    }
+    //retrieves a string parameter from a line of text
+    string getString(string line)
+    {
+        int position;
+        position=line.find("="); //Finds the position of the equal sign
+
+        string word; //Creates a string for storing word data called "word"
+        word=line.substr(position+2); //Moves from the position of "=" 2 spaces to the right
+                                      //Note: We do not need to convert word using "istringstream"
+
+        return word;
+    }
+};
+#endif
diff --git a/envi/EnviOpen.h b/envi/EnviOpen.h
new file mode 100644
index 0000000..7bfb1ee
--- /dev/null
+++ b/envi/EnviOpen.h
@@ -0,0 +1,87 @@
+#ifndef ENVIOPEN_H
+#define ENVIOPEN_H
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <ctime>
+#include <sstream>
+#include <stdio.h>
+#include "EnviHeader.h"
+#include "EnviFid.h"
+using namespace std;
+
+
+
+//return value is incorrect
+static EnviFidClass EnviOpen(string filename, char* mask = NULL, string mode = "r", int nBands = 0, int sX = 0)
+{
+	//create a file ID object
+	EnviFidClass fid;
+
+    //set the file read/write mode
+    fid.mode = mode;
+
+    //set the mask (if specified)
+	fid.mask = mask;
+
+	// load the header file only while reading
+	if(mode=="r")
+	{
+        //load the header file
+        int error = fid.header.LoadHeader(filename);
+
+        //if the LoadHeader function is not zero, the header is invalid
+        if (error != 0)
+        {
+            fid.isValid = 0;
+            return fid;
+        }
+
+	}
+
+	//default header to be used in case of writing a file
+	else
+	{
+
+		fid.header.datatype= 4;
+		fid.header.nF= nBands;
+		fid.header.sX=0;
+		fid.header.sY=1;
+		fid.header.headersize=0;
+	}
+    //load the binary file associated with this header
+        //get the binary file name
+
+	 //find the binary file name
+
+    //locate the position of the extension
+    int pos = filename.find(".hdr");
+    string binaryname = filename.substr(0, pos);
+        //load the binary file using fopen()
+        //store the resulting FILE pointer in the EnviFidClass variable
+	mode += "b";
+	fid.file = fopen (binaryname.c_str(), mode.c_str());
+	if(fid.file == NULL)
+	{
+		cout<<"Error opening binary file."<<endl;
+		fid.isValid = 0;
+		return fid;
+	}
+
+	//set the file name
+	fid.fileName = binaryname;
+
+    fid.isValid = 1;
+	return fid;
+}
+
+static EnviFidClass EnviOpen(EnviFidClass fid, char* mask = NULL, string mode = "r")
+{
+	string filename = fid.fileName;
+	filename += ".hdr";
+
+	return EnviOpen(filename, mask, mode);
+}
+
+#endif
+
diff --git a/envi/EnviRead.h b/envi/EnviRead.h
new file mode 100644
index 0000000..2525f1c
--- /dev/null
+++ b/envi/EnviRead.h
@@ -0,0 +1,65 @@
+#ifndef ENVIREAD_H
+#define ENVIREAD_H
+
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <ctime>
+#include <sstream>
+#include <stdio.h>
+#include "EnviHeader.h"
+#include "EnviFid.h"
+#include "EnviOpen.h"
+using namespace std;
+
+
+static int EnviRead(float* ptr, int nPixels, EnviFidClass &fid)
+{
+    //number of values to read from the file
+    int totalPixels = fid.header.sX * fid.header.sY;
+
+    //ERROR TEST: make sure that fid.file is a valid file pointer
+    if (fid.file == NULL)
+    {
+    cout<<"Binary file isn't open."<<endl;
+    return 0;
+    }
+
+    //---------MASK implementation
+
+
+    int i=0;
+    unsigned int sizeRead;
+    //iterate through every pixel in the image
+    //while we have not encountered the end of the file AND we have not read nPixels pixels
+    while(fid.filePos < totalPixels && i < nPixels)
+    {
+        //if there is no mask OR the current pixel is valid, read it
+        if(fid.mask == NULL || fid.mask[fid.filePos] == 1)
+        {
+            //read the pixel directly from disk
+            sizeRead = fread(&ptr[fid.header.nF * i], sizeof(float), fid.header.nF, fid.file);
+
+            //catch an error if one occurs
+            if (sizeRead != fid.header.nF)
+            {
+                cout<<"Error during read # "<<i<<endl;
+                return 0;
+            }
+            i++;
+        }
+
+
+        //otherwise, skip it
+        else
+            fseek(fid.file, sizeof(float) * fid.header.nF, SEEK_CUR);
+
+        //increment the file position (used to index into the mask)
+        fid.filePos++;
+
+    }
+
+    return i;
+}
+
+#endif
diff --git a/envi/EnviWrite.h b/envi/EnviWrite.h
new file mode 100644
index 0000000..03127ea
--- /dev/null
+++ b/envi/EnviWrite.h
@@ -0,0 +1,71 @@
+#ifndef ENVIWRITE_H
+#define ENVIWRITE_H
+#include <fstream>
+#include <stdio.h>
+#include "EnviHeader.h"
+#include "EnviFid.h"
+#include "EnviOpen.h"
+
+// Take a memory pointer, and copy all of that data to disk.
+static int EnviWrite(float* ptr, int npixels, EnviFidClass &fid)
+{
+    if(fid.mode == "r") return 0;
+
+	//open a binary file for writing
+    if (fid.file == NULL)
+    {
+        cout<<"Binary file isn't open."<<endl;
+        return 1;
+    }
+
+    //allocate an array of zeros equal to the size of a single pixel
+    float* zeroPtr = (float*) malloc(sizeof(float) * fid.header.nF);
+    memset(zeroPtr, 0, sizeof(float) * fid.header.nF);
+
+
+    int j = 0;
+    //for each pixel send to the EnviRead function
+    for(int p = 0; p<npixels; p++)
+    {
+
+        //while the mask is zero
+        while(fid.mask[fid.filePos] == 0)
+        {
+            //increment filePos
+            fid.filePos++;
+            fid.header.sX++;
+            j++;
+
+            //write a sequence of zeros
+            size_t results= fwrite(zeroPtr, sizeof(float), fid.header.nF, fid.file);
+            //check for errors
+            if (results != fid.header.nF)
+            {
+                cout<<"Error during read # "<<j<<endl;
+                return 0;
+            }
+        }
+
+        //write a pixel
+        size_t results= fwrite(&ptr[fid.header.nF * p], sizeof(float), fid.header.nF, fid.file);
+        if (results != fid.header.nF)
+        {
+            cout<<"Error during read # "<<j<<endl;
+            return 0;
+        }
+
+        //increment filePos
+        fid.filePos++;
+        fid.header.sX++;
+        j++;
+    }
+
+    return j;
+
+}
+
+
+
+
+
+#endif
diff --git a/envi/rtsEnvi.h b/envi/rtsEnvi.h
new file mode 100644
index 0000000..16a02f9
--- /dev/null
+++ b/envi/rtsEnvi.h
@@ -0,0 +1,651 @@
+#ifndef RTS_ENVI_H
+#define RTS_ENVI_H
+
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <sstream>
+
+/*	This is a class for reading and writing ENVI binary files. This class will be updated as needed.
+
+What this class CAN currently do:
+	*)  Write band images to a BSQ file.
+	*)  Append band images to a BSQ file.
+
+*/
+
+//information from an ENVI header file
+//A good resource can be found here: http://www.exelisvis.com/docs/enviheaderfiles.html
+struct EnviHeader
+{
+	std::string name;
+
+	std::string description;
+
+	unsigned int samples;	//x-axis
+	unsigned int lines;		//y-axis
+	unsigned int bands;		//spectral axis
+	unsigned int header_offset;		//header offset for binary file (in bytes)
+	std::string file_type;			//should be "ENVI Standard"
+	unsigned int data_type;			//data representation; common value is 4 (32-bit float)
+
+	enum interleaveType {BIP, BIL, BSQ};	//bip = Z,X,Y; bil = X,Z,Y; bsq = X,Y,Z
+	interleaveType interleave;
+
+	std::string sensor_type;		//not really used
+
+	enum endianType {endianLittle, endianBig};
+	endianType byte_order;			//little = least significant bit first (most systems)
+
+	double x_start, y_start;		//coordinates of the upper-left corner of the image
+	std::string wavelength_units;	//stored wavelength units
+	std::string z_plot_titles[2];
+
+	double pixel_size[2];			//pixel size along X and Y
+
+	std::vector<std::string> band_names;	//name for each band in the image
+	std::vector<double> wavelength;		//wavelength for each band
+
+	EnviHeader()
+	{
+        name = "";
+
+		//specify default values for a new or empty ENVI file
+		samples = 0;
+		lines = 0;
+		bands = 0;
+		header_offset = 0;
+		data_type = 4;
+		interleave = BSQ;
+		byte_order = endianLittle;
+		x_start = y_start = 0;
+		pixel_size[0] = pixel_size[1] = 1;
+
+		//strings
+		file_type = "ENVI Standard";
+		sensor_type = "Unknown";
+		wavelength_units = "Unknown";
+		z_plot_titles[0] = z_plot_titles[1] = "Unknown";
+	}
+
+	std::string trim(std::string line)
+	{
+		//trims whitespace from the beginning and end of line
+		int start_i, end_i;
+		for(start_i=0; start_i < line.length(); start_i++)
+			if(line[start_i] != 32)
+			{
+				break;
+			}
+
+		for(end_i = line.length()-1; end_i >= start_i; end_i--)
+			if(line[end_i] != ' ')
+			{
+				break;
+			}
+
+		return line.substr(start_i, end_i - start_i+1);
+	}
+
+
+	std::string get_token(std::string line)
+	{
+		//returns a variable name; in this case we look for the '=' sign
+		size_t i = line.find_first_of('=');
+
+		std::string result;
+		if(i != std::string::npos)
+			result = trim(line.substr(0, i-1));
+
+		return result;
+	}
+
+	std::string get_data_str(std::string line)
+	{
+		size_t i = line.find_first_of('=');
+
+		std::string result;
+		if(i != std::string::npos)
+			result = trim(line.substr(i+1));
+		else
+		{
+			std::cout<<"ENVI Header error - data not found for token: "<<get_token(line)<<std::endl;
+			exit(1);
+		}
+		return result;
+	}
+
+	std::string get_brace_str(std::string token, std::string line, std::ifstream &file)
+	{
+		//this function assembles all of the characters between curly braces
+		//this is how strings are defined within an ENVI file
+
+		std::string result;
+
+		//first, find the starting brace
+		size_t i;
+		do
+		{
+			i = line.find_first_of('{');
+			if(i != std::string::npos)
+				break;
+		}while(file);
+
+		//if i is still npos, we have reached the end of the file without a brace...something is wrong
+		if(i == std::string::npos)
+		{
+			std::cout<<"ENVI Header error - string token declared without being defined: "<<token<<std::endl;
+			exit(1);
+		}
+		line = line.substr(i+1);
+
+		//copy character data into the result string until we find a closing brace
+		while(file)
+		{
+			i = line.find_first_of('}');
+
+
+			if(i != std::string::npos)
+			{
+				result += line.substr(0, i);
+				break;
+			}
+			else
+				result += line;
+
+			getline(file, line);
+		}
+
+		if(i == std::string::npos)
+		{
+			std::cout<<"ENVI Header error - string token declared without a terminating '}': "<<token<<std::endl;
+			exit(1);
+		}
+
+		return trim(result);
+	}
+
+	std::vector<std::string> get_string_seq(std::string token, std::string sequence)
+	{
+		//this function returns a sequence of comma-delimited strings
+		std::vector<std::string> result;
+
+		std::string entry;
+		size_t i;
+		do
+		{
+			i = sequence.find_first_of(',');
+			entry = sequence.substr(0, i);
+			sequence = sequence.substr(i+1);
+			result.push_back(trim(entry));
+		}while(i != std::string::npos);
+
+		return result;
+	}
+
+	std::vector<double> get_double_seq(std::string token, std::string sequence)
+	{
+		//this function returns a sequence of comma-delimited strings
+		std::vector<double> result;
+
+		std::string entry;
+		size_t i;
+		do
+		{
+			i = sequence.find_first_of(',');
+			entry = sequence.substr(0, i);
+			sequence = sequence.substr(i+1);
+			result.push_back(atof(entry.c_str()));
+		}while(i != std::string::npos);
+
+		return result;
+	}
+
+	void load(std::string filename, std::string mode = "r")
+	{
+		name = filename;
+
+
+		//open the header file
+		std::ifstream file(name.c_str());
+
+		if(!file)
+		{
+            //if we are in read mode, throw an exception
+            if(mode == "r")
+            {
+                std::cout<<"ENVI header file not found: "<<name<<std::endl;
+                exit(1);
+            }
+            else
+            {
+                //there is no file, so just return the default structure
+                return;
+            }
+		}
+
+		//the first line should just be "ENVI"
+		std::string line;
+		getline(file, line);
+		if(line != "ENVI")
+		{
+			std::cout<<"The header doesn't appear to be an ENVI file. The first line should be 'ENVI'."<<std::endl;
+			exit(1);
+		}
+
+		//for each line in the file, get the token
+		std::string token;
+		while(file)
+		{
+
+			//get a line
+			getline(file, line);
+
+			//get the token
+			token = get_token(line);
+
+			if(token == "description")
+				description = get_brace_str(token, line, file);
+			else if(token == "band names")
+			{
+				std::string string_sequence = get_brace_str(token, line, file);
+				band_names = get_string_seq(token, string_sequence);
+			}
+			else if(token == "wavelength")
+			{
+				std::string string_sequence = get_brace_str(token, line, file);
+				wavelength = get_double_seq(token, string_sequence);
+			}
+			else if(token == "pixel size")
+			{
+				std::string string_sequence = get_brace_str(token, line, file);
+				std::vector<double> pxsize = get_double_seq(token, string_sequence);
+				pixel_size[0] = pxsize[0];
+				pixel_size[1] = pxsize[1];
+			}
+			else if(token == "z plot titles")
+			{
+				std::string string_sequence = get_brace_str(token, line, file);
+				std::vector<std::string> titles = get_string_seq(token, string_sequence);
+				z_plot_titles[0] = titles[0];
+				z_plot_titles[1] = titles[1];
+			}
+
+			else if(token == "samples")
+				samples = atoi(get_data_str(line).c_str());
+			else if(token == "lines")
+			{
+				lines = atoi(get_data_str(line).c_str());
+				std::cout<<"Number of lines loaded: "<<lines<<std::endl;
+            }
+			else if(token == "bands")
+				bands = atoi(get_data_str(line).c_str());
+			else if(token == "header offset")
+				header_offset = atoi(get_data_str(line).c_str());
+			else if(token == "file type")
+				file_type = get_data_str(line);
+			else if(token == "data type")
+				data_type = atoi(get_data_str(line).c_str());
+			else if(token == "interleave")
+			{
+				std::string interleave_str = get_data_str(line);
+				if(interleave_str == "bip")
+					interleave = BIP;
+				else if(interleave_str == "bil")
+					interleave = BIL;
+				else if(interleave_str == "bsq")
+					interleave = BSQ;
+			}
+			else if(token == "sensor type")
+				sensor_type = get_data_str(line);
+			else if(token == "byte order")
+				byte_order = (endianType)atoi(get_data_str(line).c_str());
+			else if(token == "x start")
+				x_start = atof(get_data_str(line).c_str());
+			else if(token == "y start")
+				y_start = atof(get_data_str(line).c_str());
+			else if(token == "wavelength units")
+				wavelength_units = get_data_str(line);
+		}
+
+		//close the file
+		file.close();
+	}
+
+	void save(std::string filename)
+	{
+		//open a file
+		std::ofstream outfile(filename.c_str());
+
+		//write the ENVI type identifier
+		outfile<<"ENVI"<<std::endl;
+
+		//output all of the data
+		outfile<<"description = {"<<std::endl;
+		outfile<<"  "<<description<<"}"<<std::endl;
+
+		outfile<<"samples = "<<samples<<std::endl;
+		outfile<<"lines = "<<lines<<std::endl;
+		outfile<<"bands = "<<bands<<std::endl;
+		outfile<<"header offset = "<<header_offset<<std::endl;
+		outfile<<"file type = "<<file_type<<std::endl;
+		outfile<<"data type = "<<data_type<<std::endl;
+		outfile<<"interleave = ";
+		if(interleave == BIP)
+			outfile<<"bip";
+		if(interleave == BIL)
+			outfile<<"bil";
+		if(interleave == BSQ)
+			outfile<<"bsq";
+			outfile<<std::endl;
+		outfile<<"sensor type = "<<sensor_type<<std::endl;
+		outfile<<"byte order = "<<byte_order<<std::endl;
+		outfile<<"x start = "<<x_start<<std::endl;
+		outfile<<"y start = "<<y_start<<std::endl;
+		outfile<<"wavelength units = "<<wavelength_units<<std::endl;
+		outfile<<"z plot titles = {";
+			outfile<<z_plot_titles[0]<<", "<<z_plot_titles[1]<<"}"<<std::endl;
+		outfile<<"pixel size = {"<<pixel_size[0]<<", "<<pixel_size[1]<<", units=Meters}"<<std::endl;
+		if(band_names.size() > 0)
+		{
+			outfile<<"band names = {"<<std::endl;
+			for(int i=0; i<band_names.size(); i++)
+			{
+				outfile<<band_names[i];
+				if(i < band_names.size() - 1)
+					outfile<<", ";
+			}
+			outfile<<"}"<<std::endl;
+		}
+		outfile<<"wavelength = {"<<std::endl;
+			for(int i=0; i<wavelength.size()-1; i++)
+				outfile<<wavelength[i]<<", ";
+			outfile<<wavelength.back()<<"}"<<std::endl;
+
+		outfile.close();
+	}
+
+	void save()
+	{
+        //std::cout<<"ENVI Header Name: "<<name<<std::endl;
+		save(name);
+	}
+
+	//returns the size of a single value (in bytes)
+	unsigned int valsize()
+	{
+		switch(data_type)
+		{
+		case 1:			//1 = 8-bit byte
+			return 1;
+		case 2:			//16-bit signed integer
+		case 12:		//16-bit unsigned integer
+			return 2;
+		case 3:			//32-bit signed long integer
+		case 4:			//32-bit floating point
+		case 13:		//32-bit unsigned long integer
+			return 4;
+		case 5:			//64-bit double-precision floating point
+		case 6:			//32-bit complex value
+		case 14:		//64-bit signed long integer
+		case 15:		//64-bit unsigned long integer
+			return 8;
+		case 9:			//64-bit complex value
+			return 16;
+		default:
+			return 0;
+		}
+
+	}
+};		//end EnviHeader
+
+class EnviFile
+{
+	EnviHeader header;
+	FILE* data;
+	std::string mode;
+
+	//a locked file has alread had data written..certain things can't be changed
+	//		accessor functions deal with these issues
+	bool locked;
+
+	void init()
+	{
+		locked = false;
+		data = NULL;
+	}
+
+	void readonly()
+	{
+		if(mode == "r")
+		{
+			std::cout<<"ENVI Error: changes cannot be made to a read-only file."<<std::endl;
+			exit(1);
+		}
+	}
+
+	//this function determines if the current software is capable of the requested change
+	bool caps(EnviHeader new_header)
+	{
+		readonly();		//if the file is read-only, throw an exception
+
+		//we currently only support BSQ
+		if(new_header.interleave == EnviHeader::BIP || new_header.interleave == EnviHeader::BIL)
+		{
+			std::cout<<"This software only supports BSQ.";
+			return false;
+		}
+
+		//if the number of bands has changed
+		if(header.bands != new_header.bands)
+		{
+			//the new header must be a BSQ file
+			if(new_header.interleave != EnviHeader::BSQ)
+			{
+				std::cout<<"ENVI Error: can only add bands to a BSQ file."<<std::endl;
+				return false;
+			}
+		}
+
+		//disallow anything that can't be done when the file is locked
+		if(locked)
+		{
+			if(header.lines != new_header.lines)
+			{
+				std::cout<<"ENVI Error: cannot change the number of lines from "<<header.lines<<" to "<<new_header.lines<<" in a locked BSQ file."<<std::endl;
+				return false;
+			}
+			if(header.samples != new_header.samples)
+			{
+				std::cout<<"ENVI Error: cannot change the number of samples in a locked BSQ file."<<std::endl;
+				return false;
+			}
+			if(header.data_type != new_header.data_type)
+			{
+				std::cout<<"ENVI Error: cannot change the data type of a locked file."<<std::endl;
+				return false;
+			}
+			if(header.byte_order != new_header.byte_order)
+			{
+				std::cout<<"ENVI Error: cannot change the byte order of a locked file."<<std::endl;
+				return false;
+			}
+		}
+
+		return true;
+
+	}
+
+	EnviHeader load_header(std::string header_file, std::string file_mode)
+	{
+		EnviHeader new_header;
+		new_header.load(header_file, file_mode);
+		return new_header;
+	}
+
+	void set_header(EnviHeader new_header)
+	{
+		if(caps(new_header))
+			header = new_header;
+		else
+			exit(1);
+	}
+
+	bool test_exist(std::string filename)
+	{
+        //tests for the existance of the specified binary file
+        FILE* f;
+        long sz = 0;
+        f = fopen(filename.c_str(), "rb");
+        if(f == NULL)
+        {
+            std::cout<<"ENVI Error: error testing file existance."<<std::endl;
+            exit(1);
+        }
+
+        fseek(f, 0, SEEK_END);
+        sz = ftell(f);
+        fclose(f);
+
+        if(sz>0)
+            return true;
+        else
+            return false;
+
+    }
+
+public:
+
+	EnviFile()
+	{
+		init();
+	}
+
+	EnviFile(std::string filename, std::string file_mode = "r")
+	{
+		init();
+		open(filename, filename + ".hdr", file_mode);
+	}
+
+	void open(std::string file_name, std::string header_name, std::string file_mode = "r")
+	{
+		mode = file_mode;
+
+		header.name = header_name;
+
+		//lock the file if we are loading or appending
+		if(mode == "r" || mode == "a" || mode == "r+" || mode == "a+")
+		{
+            //load the ENVI header - failure will cause an exit
+			set_header(load_header(header_name, file_mode));
+
+            if(mode == "r" || test_exist(file_name))
+            {
+                locked = true;
+            }
+
+		}
+		else
+		{
+			header.name = header_name;
+		}
+
+		//open the data file
+		data = fopen(file_name.c_str(), (mode + "b").c_str());
+
+
+	}
+
+	void setDescription(std::string desc)
+	{
+		readonly();		//if the file is read-only, throw an exception
+		header.description = desc;
+	}
+
+	void setZPlotTitles(std::string spectrum, std::string value)
+	{
+		readonly();		//if the file is read-only, throw an exception
+		header.z_plot_titles[0] = spectrum;
+		header.z_plot_titles[1] = value;
+	}
+
+	void setPixelSize(double sx, double sy)
+	{
+		readonly();		//if the file is read-only, throw an exception
+		header.pixel_size[0] = sx;
+		header.pixel_size[1] = sy;
+	}
+
+	void setWavelengthUnits(std::string units)
+	{
+		readonly();		//if the file is read-only, throw an exception
+		header.wavelength_units = units;
+	}
+
+	void setCoordinates(double x, double y)
+	{
+		readonly();		//if the file is read-only, throw an exception
+		header.x_start = x;
+		header.y_start = y;
+	}
+
+	//FUNCTIONS THAT MAY BE DISALLOWED IN A LOCKED FILE
+	void setSamples(unsigned int samples)
+	{
+		EnviHeader newHeader = header;
+		newHeader.samples = samples;
+		set_header(newHeader);
+	}
+
+	void setLines(unsigned int lines)
+	{
+		EnviHeader newHeader = header;
+		newHeader.lines = lines;
+		set_header(newHeader);
+	}
+
+	void setBands(unsigned int bands)
+	{
+		EnviHeader newHeader = header;
+		newHeader.bands = bands;
+		set_header(newHeader);
+	}
+
+	//DATA MANIPULATION
+	void addBand(void* band, unsigned int samples, unsigned int lines, double wavelength, std::string band_name ="")
+	{
+		//add a band to the file
+
+		EnviHeader newHeader = header;
+		newHeader.bands++;
+		newHeader.samples = samples;
+		newHeader.lines = lines;
+		newHeader.wavelength.push_back(wavelength);
+		if(band_name != "")
+			newHeader.band_names.push_back(band_name);
+		set_header(newHeader);
+
+		//copy the data to the file
+		fwrite(band, header.valsize(), samples * lines, data);
+
+		//since data was written, the file must be locked
+		locked = true;
+	}
+
+	//close file
+	void close()
+	{
+        //std::cout<<"ENVI File Mode: "<<mode<<std::endl;
+
+		//close the data stream
+		fclose(data);
+
+		//close the header
+		if(mode != "r")
+			header.save();
+	}
+
+};		//end EnviFile
+
+
+#endif
diff --git a/rts/material.h b/rts/material.h
index d7d763d..c6b08ac 100644
--- a/rts/material.h
+++ b/rts/material.h
@@ -8,7 +8,7 @@
 #include <complex>
 #include <algorithm>
 #include <sstream>
-#include "rts/rtcomplex.h"
+#include "rts/rtsComplex.h"
 
 #define PI  3.14159
 
@@ -49,7 +49,7 @@ namespace rts{
     {
         //wavelength (in microns)
         T lambda;
-        rtcomplex<T> n;
+        rtsComplex<T> n;
     };
 
     template <class T>
@@ -196,6 +196,15 @@ namespace rts{
             dispersion.push_back(ri);
         }
 
+        void add(T lambda, rtsComplex<T> n)
+        {
+            refIndex<T> ri;
+            ri.lambda = lambda;
+            ri.n = n;
+
+            dispersion.push_back(ri);
+        }
+
         //comparison function for sorting
         static bool compare(refIndex<T> a, refIndex<T> b)
         {
@@ -234,9 +243,9 @@ namespace rts{
 
 #ifdef FFTW_AVAILABLE
 			//allocate memory for the FFT
-			rtcomplex<T>* Chi2 = (rtcomplex<T>*)fftw_malloc(sizeof(rtcomplex<T>) * N * pf);
-			rtcomplex<T>* Chi2FFT = (rtcomplex<T>*)fftw_malloc(sizeof(rtcomplex<T>) * N * pf);
-			rtcomplex<T>* Chi1 = (rtcomplex<T>*)fftw_malloc(sizeof(rtcomplex<T>) * N * pf);
+			rtsComplex<T>* Chi2 = (rtsComplex<T>*)fftw_malloc(sizeof(rtsComplex<T>) * N * pf);
+			rtsComplex<T>* Chi2FFT = (rtsComplex<T>*)fftw_malloc(sizeof(rtsComplex<T>) * N * pf);
+			rtsComplex<T>* Chi1 = (rtsComplex<T>*)fftw_malloc(sizeof(rtsComplex<T>) * N * pf);
 
 			//create an FFT plan for the forward and inverse transforms
 			fftw_plan planForward, planInverse;
@@ -267,8 +276,8 @@ namespace rts{
 			}
 
 			//use linear interpolation between the start and end points to pad the spectrum
-			//rtcomplex<T> nMin = dispersion.back();
-			//rtcomplex<T> nMax = dispersion.front();
+			//rtsComplex<T> nMin = dispersion.back();
+			//rtsComplex<T> nMax = dispersion.front();
 			T a;
 			for(int i=N; i<N*pf; i++)
 			{
@@ -282,7 +291,7 @@ namespace rts{
 			fftw_execute(planForward);
 
 			//perform the Hilbert transform in the Fourier domain
-			rtcomplex<T> j(0, 1);
+			rtsComplex<T> j(0, 1);
 			for(int i=0; i<N*pf; i++)
 			{
                 //if w = 0, set the DC component to zero
@@ -375,12 +384,12 @@ namespace rts{
             n0 = n;
         }
 
-        material(std::string filename, std::string format, T scaleA = 1.0)
+        material(std::string filename, std::string format = "", T scaleA = 1.0)
         {
             fromFile(filename, format);
         }
 
-        void fromFile(std::string filename, std::string format, T scaleA = 1.0)
+        void fromFile(std::string filename, std::string format = "", T scaleA = 1.0)
         {
             //clear any previous values
             dispersion.clear();
@@ -402,7 +411,7 @@ namespace rts{
 
         }
 
-        void fromStr(std::string str, std::string format, T scaleA = 1.0)
+        void fromStr(std::string str, std::string format = "", T scaleA = 1.0)
         {
             //create a string stream to process the input data
             std::stringstream ss(str);
@@ -410,9 +419,38 @@ namespace rts{
             //this string will be read line-by-line (where each line is an entry)
             std::string line;
 
-            //create an entry structure (for now just a basic one)
+            //if the format is not provided, see if it is in the file, otherwise use a default
+            if(format == "")
+            {
+                //set a default format of "lambda,n,k"
+                format = "microns,n,k";
+
+                //see if the first line is a comment
+                char c = ss.peek();
+                if(c == '#')
+                {
+                    //get the first line
+                    getline(ss, line);
+                    //look for a bracket, denoting the start of a format string
+                    int istart = line.find('[');
+                    if(istart != std::string::npos)
+                    {
+                        //look for a bracket terminating the format string
+                        int iend = line.find(']');
+                        if(iend != std::string::npos)
+                        {
+                            //read the string between the brackets
+                            format = line.substr(istart+1, iend - istart - 1);
+                        }
+                    }
+                }
+
+            }
+
             entryType<T> entry(format);
 
+            std::cout<<"Loading material with format: "<<format<<std::endl;
+
             T lambda, n, k;
             while(!ss.eof())
             {
@@ -434,24 +472,37 @@ namespace rts{
         }
 
         //convert the material to a string
-        std::string toStr(std::string format = "microns,n,k")
+        std::string toStr(std::string format = "microns,n,k", bool reverse_order = false)
         {
             std::stringstream ss;
 			entryType<T> entry(format);
-            for(unsigned int l=0; l<dispersion.size(); l++)
-            {
-                if(l > 0) ss<<std::endl;
-                ss<<entry.outputEntry(dispersion[l]);
+
+			if(reverse_order)
+			{
+                for(int l=dispersion.size() - 1; l>=0; l--)
+                {
+                    if(l < dispersion.size() - 1) ss<<std::endl;
+                    ss<<entry.outputEntry(dispersion[l]);
+                }
+
+			}
+			else
+			{
+                for(unsigned int l=0; l<dispersion.size(); l++)
+                {
+                    if(l > 0) ss<<std::endl;
+                    ss<<entry.outputEntry(dispersion[l]);
+                }
             }
 
             return ss.str();
         }
 
-        void save(std::string filename, std::string format = "microns,n,k")
+        void save(std::string filename, std::string format = "microns,n,k", bool reverse_order = false)
         {
-            ofstream outfile(filename.c_str());
+            std::ofstream outfile(filename.c_str());
             outfile<<"#material file saved as [" + format + "]"<<std::endl;
-            outfile<<toStr(format)<<std::endl;
+            outfile<<toStr(format, reverse_order)<<std::endl;
 
         }
 
@@ -475,7 +526,7 @@ namespace rts{
 
 		}
 
-		rtcomplex<T> getN(T l)
+		rtsComplex<T> getN(T l)
 		{
             //declare an iterator
             typename std::vector< refIndex<T> >::iterator it;
@@ -499,16 +550,16 @@ namespace rts{
                 //std::cout<<lMin<<"----------"<<lMax<<std::endl;
 
                 T a = (l - lMin) / (lMax - lMin);
-                rtcomplex<T> riMin = (*(it - 1)).n;
-                rtcomplex<T> riMax = (*it).n;
-                rtcomplex<T> interp;
-                interp = rtcomplex<T>(a, 0.0) * riMin + rtcomplex<T>(1.0 - a, 0.0) * riMax;
+                rtsComplex<T> riMin = (*(it - 1)).n;
+                rtsComplex<T> riMax = (*it).n;
+                rtsComplex<T> interp;
+                interp = rtsComplex<T>(a, 0.0) * riMin + rtsComplex<T>(1.0 - a, 0.0) * riMax;
                 return interp;
             }
 
 		}
         //interpolate the given lambda value and return the index of refraction
-        rtcomplex<T> operator()(T l)
+        rtsComplex<T> operator()(T l)
         {
             return getN(l);
         }
diff --git a/rts/point.h b/rts/point.h
deleted file mode 100644
index 2cdf47b..0000000
--- a/rts/point.h
+++ /dev/null
@@ -1,124 +0,0 @@
-#ifndef RTS_POINT_H
-#define RTS_POINT_H
-
-//#include "rts/vector.h"
-#include <string.h>
-
-namespace rts
-{
-
-template <class T, int N>
-struct point
-{
-	T p[N];
-
-	CUDA_CALLABLE point()
-	{
-
-	}
-
-	//efficiency constructor, makes construction easier for 1D-4D vectors
-	CUDA_CALLABLE point(T x, T y = (T)0.0, T z = (T)0.0, T w = (T)0.0)
-	{
-		if(N >= 1)
-			p[0] = x;
-		if(N >= 2)
-			p[1] = y;
-		if(N >= 3)
-			p[2] = z;
-		if(N >= 4)
-			p[3] = w;
-	}
-
-	//arithmetic operators
-	CUDA_CALLABLE rts::point<T, N> operator+(rts::vector<T, N> v)
-	{
-        rts::point<T, N> r;
-
-        //calculate the position of the resulting point
-        for(int i=0; i<N; i++)
-            r.p[i] = p[i] + v.v[i];
-
-        return r;
-	}
-	CUDA_CALLABLE rts::point<T, N> operator-(rts::vector<T, N> v)
-	{
-        rts::point<T, N> r;
-
-        //calculate the position of the resulting point
-        for(int i=0; i<N; i++)
-            r.p[i] = p[i] - v.v[i];
-
-        return r;
-	}
-	CUDA_CALLABLE rts::vector<T, N> operator-(rts::point<T, N> rhs)
-	{
-        rts::vector<T, N> r;
-
-        //calculate the position of the resulting point
-        for(int i=0; i<N; i++)
-            r.v[i] = p[i] - rhs.p[i];
-
-        return r;
-	}
-	CUDA_CALLABLE rts::point<T, N> operator*(T rhs)
-	{
-        rts::point<T, N> r;
-
-        //calculate the position of the resulting point
-        for(int i=0; i<N; i++)
-            r.p[i] = p[i] * rhs;
-
-        return r;
-	}
-
-	CUDA_CALLABLE point(const T(&data)[N])
-	{
-		memcpy(p, data, sizeof(T) * N);
-	}
-
-	std::string toStr()
-	{
-		std::stringstream ss;
-
-		ss<<"(";
-		for(int i=0; i<N; i++)
-		{
-			ss<<p[i];
-			if(i != N-1)
-				ss<<", ";
-		}
-		ss<<")";
-
-		return ss.str();
-	}
-
-	//bracket operator
-	CUDA_CALLABLE T& operator[](int i)
-	{
-        return p[i];
-    }
-
-};
-
-}	//end namespace rts
-
-template <typename T, int N>
-std::ostream& operator<<(std::ostream& os, rts::point<T, N> p)
-{
-    os<<p.toStr();
-    return os;
-}
-
-//arithmetic
-template <typename T, int N>
-CUDA_CALLABLE rts::point<T, N> operator*(T lhs, rts::point<T, N> rhs)
-{
-    rts::point<T, N> r;
-
-    return rhs * lhs;
-}
-
-
-
-#endif
diff --git a/rts/rect.h b/rts/rect.h
index 53e07f3..d2d288d 100644
--- a/rts/rect.h
+++ b/rts/rect.h
@@ -3,8 +3,8 @@
 
 //enable CUDA_CALLABLE macro
 #include "rts/cuda_callable.h"
-#include "rts/vector.h"
-#include "rts/point.h"
+#include "rts/rtsVector.h"
+#include "rts/rtsPoint.h"
 #include <iostream>
 
 namespace rts{
@@ -27,9 +27,9 @@ struct rect
 	T B[N];
 	T C[N];*/
 
-	rts::point<T, N> A;
-	rts::vector<T, N> X;
-	rts::vector<T, N> Y;
+	rts::rtsPoint<T, N> A;
+	rts::rtsVector<T, N> X;
+	rts::rtsVector<T, N> Y;
 
 
 	CUDA_CALLABLE rect()
@@ -37,7 +37,7 @@ struct rect
 
 	}
 
-	CUDA_CALLABLE rect(point<T, N> a, point<T, N> b, point<T, N> c)
+	CUDA_CALLABLE rect(rtsPoint<T, N> a, rtsPoint<T, N> b, rtsPoint<T, N> c)
 	{
 
 		A = a;
@@ -46,27 +46,27 @@ struct rect
 
 	}
 
-	CUDA_CALLABLE rect(rts::point<T, N> pMin, rts::point<T, N> pMax, rts::vector<T, N> normal)
+	CUDA_CALLABLE rect(rts::rtsPoint<T, N> pMin, rts::rtsPoint<T, N> pMax, rts::rtsVector<T, N> normal)
 	{
 
-        //assign the corner point
+        //assign the corner rtsPoint
         A = pMin;
 
         //compute the vector from pMin to pMax
-        rts::vector<T, 3> v0;
+        rts::rtsVector<T, 3> v0;
         v0 = pMax - pMin;
 
         //compute the cross product of A and the plane normal
-        rts::vector<T, 3> v1;
+        rts::rtsVector<T, 3> v1;
         v1 = v0.cross(normal);
 
 
-        //calculate point B
-        rts::point<T, 3> B;
+        //calculate rtsPoint B
+        rts::rtsPoint<T, 3> B;
         B = A + v0 * 0.5 + v1 * 0.5;
 
-        //calculate point C
-        rts::point<T, 3> C;
+        //calculate rtsPoint C
+        rts::rtsPoint<T, 3> C;
         C = A  + v0 * 0.5 - v1 * 0.5;
 
         //calculate X and Y
@@ -78,7 +78,7 @@ struct rect
 
 	}
 
-	/*CUDA_CALLABLE rect(rts::point<T, N> p, rts::vector<T, N> x, rts::vector<T, N> y, T sx, T sy)
+	/*CUDA_CALLABLE rect(rts::rtsPoint<T, N> p, rts::vector<T, N> x, rts::vector<T, N> y, T sx, T sy)
 	{
 		//This constructor creates a rect given a position, orientation, and size
 		//	p	= center position of the rect
@@ -100,16 +100,16 @@ struct rect
 
 	}*/
 
-	CUDA_CALLABLE rts::point<T, N> p(T a, T b)
+	CUDA_CALLABLE rts::rtsPoint<T, N> p(T a, T b)
 	{
-		rts::point<T, N> result;
+		rts::rtsPoint<T, N> result;
 		//given the two parameters a, b = [0 1], returns the position in world space
 		result = A + X * a + Y * b;
 		
 		return result;
 	}
 
-	CUDA_CALLABLE rts::point<T, N> operator()(T a, T b)
+	CUDA_CALLABLE rts::rtsPoint<T, N> operator()(T a, T b)
 	{
 		return p(a, b);
 	}
@@ -126,6 +126,23 @@ struct rect
         return ss.str();
 
 	}
+
+	CUDA_CALLABLE rect<T, N> operator*(T rhs)
+	{
+		//scales the plane by a scalar value
+
+		//compute the center rtsPoint
+		rts::rtsPoint<T, N> c = A + X*0.5 + Y*0.5;
+
+		//create the new rectangle
+		rect<T, N> result;
+		result.X = X * rhs;
+		result.Y = Y * rhs;
+		result.A = c - result.X*0.5 - result.Y*0.5;
+
+		return result;
+
+	}
 };
 
 }	//end namespace rts
diff --git a/rts/rtcomplex.h b/rts/rtcomplex.h
deleted file mode 100644
index daeaf98..0000000
--- a/rts/rtcomplex.h
+++ /dev/null
@@ -1,434 +0,0 @@
-/*RTS Complex number class.  This class is CUDA compatible,
-and can therefore be used in CUDA code and on CUDA devices.
-*/
-
-#ifndef RTS_COMPLEX
-#define RTS_COMPLEX
-
-#include "cuda_callable.h"
-#include <cmath>
-#include <string>
-#include <sstream>
-#include <iostream>
-
-namespace rts
-{
-
-template <class T>
-struct rtcomplex
-{
-    T r, i;
-
-    //default constructor
-    CUDA_CALLABLE rtcomplex()
-    {
-        r = 0.0;
-		i = 0.0;
-    }
-
-	//access methods
-	T real()
-	{
-		return r;
-	}
-
-	T real(T r_val)
-	{
-		r = r_val;
-		return r_val;
-	}
-
-	T imag()
-	{
-		return i;
-	}
-	T imag(T i_val)
-	{
-		i = i_val;
-		return i_val;
-	}
-
-
-
-
-
-    //constructor when given real and imaginary values
-    CUDA_CALLABLE rtcomplex(T r, T i)
-    {
-        this->r = r;
-        this->i = i;
-    }
-
-    //return the current value multiplied by i
-    CUDA_CALLABLE rtcomplex<T> imul()
-    {
-        rtcomplex<T> result;
-        result.r = -i;
-        result.i = r;
-
-        return result;
-    }
-
-	//ARITHMETIC OPERATORS--------------------
-
-    //binary + operator (returns the result of adding two complex values)
-    CUDA_CALLABLE rtcomplex<T> operator+ (const rtcomplex<T> rhs)
-    {
-        rtcomplex<T> result;
-        result.r = r + rhs.r;
-        result.i = i + rhs.i;
-        return result;
-    }
-
-	CUDA_CALLABLE rtcomplex<T> operator+ (const T rhs)
-    {
-        rtcomplex<T> result;
-        result.r = r + rhs;
-        result.i = i;
-        return result;
-    }
-
-    //binary - operator (returns the result of adding two complex values)
-    CUDA_CALLABLE rtcomplex<T> operator- (const rtcomplex<T> rhs)
-    {
-        rtcomplex<T> result;
-        result.r = r - rhs.r;
-        result.i = i - rhs.i;
-        return result;
-    }
-
-    //binary - operator (returns the result of adding two complex values)
-    CUDA_CALLABLE rtcomplex<T> operator- (const T rhs)
-    {
-        rtcomplex<T> result;
-        result.r = r - rhs;
-        result.i = i;
-        return result;
-    }
-
-    //binary MULTIPLICATION operators (returns the result of multiplying complex values)
-    CUDA_CALLABLE rtcomplex<T> operator* (const rtcomplex<T> rhs)
-    {
-        rtcomplex<T> result;
-        result.r = r * rhs.r - i * rhs.i;
-        result.i = r * rhs.i + i * rhs.r;
-        return result;
-    }
-    CUDA_CALLABLE rtcomplex<T> operator* (const T rhs)
-    {
-        return rtcomplex<T>(r * rhs, i * rhs);
-    }
-
-    //binary DIVISION operators (returns the result of dividing complex values)
-    CUDA_CALLABLE rtcomplex<T> operator/ (const rtcomplex<T> rhs)
-    {
-        rtcomplex<T> result;
-        T denom = rhs.r * rhs.r + rhs.i * rhs.i;
-        result.r = (r * rhs.r + i * rhs.i) / denom;
-        result.i = (- r * rhs.i + i * rhs.r) / denom;
-
-        return result;
-    }
-    CUDA_CALLABLE rtcomplex<T> operator/ (const T rhs)
-    {
-        return rtcomplex<T>(r / rhs, i / rhs);
-    }
-
-    //ASSIGNMENT operators-----------------------------------
-    CUDA_CALLABLE rtcomplex<T> & operator=(const rtcomplex<T> &rhs)
-    {
-        //check for self-assignment
-        if(this != &rhs)
-        {
-            this->r = rhs.r;
-            this->i = rhs.i;
-        }
-        return *this;
-    }
-    CUDA_CALLABLE rtcomplex<T> & operator=(const T &rhs)
-    {
-        this->r = rhs;
-        this->i = 0;
-
-		return *this;
-    }
-
-    //arithmetic assignment operators
-    CUDA_CALLABLE rtcomplex<T> operator+=(const rtcomplex<T> &rhs)
-    {
-		*this = *this + rhs;
-        return *this;
-    }
-    CUDA_CALLABLE rtcomplex<T> operator+=(const T &rhs)
-    {
-		*this = *this + rhs;
-        return *this;
-    }
-
-    CUDA_CALLABLE rtcomplex<T> operator*=(const rtcomplex<T> &rhs)
-    {
-		*this = *this * rhs;
-        return *this;
-    }
-	CUDA_CALLABLE rtcomplex<T> operator*=(const T &rhs)
-    {
-		*this = *this * rhs;
-        return *this;
-    }
-	//divide and assign
-	CUDA_CALLABLE rtcomplex<T> operator/=(const rtcomplex<T> &rhs)
-    {
-		*this = *this / rhs;
-        return *this;
-    }
-    CUDA_CALLABLE rtcomplex<T> operator/=(const T &rhs)
-    {
-		*this = *this / rhs;
-        return *this;
-    }
-
-    //absolute value operator (returns the absolute value of the complex number)
-	CUDA_CALLABLE T abs()
-	{
-		return std::sqrt(r * r + i * i);
-	}
-
-	CUDA_CALLABLE rtcomplex<T> log()
-	{
-        rtcomplex<T> result;
-        result.r = std::log(std::sqrt(r * r + i * i));
-        result.i = std::atan2(i, r);
-
-
-        return result;
-	}
-
-	CUDA_CALLABLE rtcomplex<T> exp()
-	{
-        rtcomplex<T> result;
-
-        T e_r = std::exp(r);
-        result.r = e_r * std::cos(i);
-        result.i = e_r * std::sin(i);
-
-        return result;
-	}
-
-	/*CUDA_CALLABLE complex<T> pow(int y)
-	{
-
-        return pow((double)y);
-	}*/
-
-	CUDA_CALLABLE rtcomplex<T> pow(T y)
-	{
-        rtcomplex<T> result;
-
-        result = log() * y;
-
-        return result.exp();
-	}
-
-	CUDA_CALLABLE rtcomplex<T> sqrt()
-	{
-		rtcomplex<T> result;
-
-		//convert to polar coordinates
-		T a = std::sqrt(r*r + i*i);
-		T theta = std::atan2(i, r);
-
-		//find the square root
-		T a_p = std::sqrt(a);
-		T theta_p = theta/2.0;
-
-		//convert back to cartesian coordinates
-		result.r = a_p * std::cos(theta_p);
-		result.i = a_p * std::sin(theta_p);
-
-		return result;
-	}
-
-	std::string toStr()
-	{
-		std::stringstream ss;
-		ss<<"("<<r<<","<<i<<")";
-
-		return ss.str();
-	}
-
-	//COMPARISON operators
-	CUDA_CALLABLE bool operator==(rtcomplex<T> rhs)
-	{
-        if(r == rhs.r && i == rhs.i)
-            return true;
-        return false;
-    }
-
-    CUDA_CALLABLE bool operator==(T rhs)
-	{
-        if(r == rhs && i == (T)0.0)
-            return true;
-        return false;
-    }
-
-	/*//FRIEND functions
-    //unary minus operator (for negating the complex number)
-	template<class A> CUDA_CALLABLE friend complex<A> operator-(const complex<A> &rhs);
-
-	//multiplication by T values when the complex number isn't on the left hand side
-	template<class A> CUDA_CALLABLE friend complex<A> operator*(const A a, const complex<A> b);
-
-	//division by T values when the complex number isn't on the left hand side
-	template<class A> CUDA_CALLABLE friend complex<A> operator/(const A a, const complex<A> b);
-
-	//POW function
-	//template<class A> CUDA_CALLABLE friend complex<A> pow(const complex<A> x, T y);
-	template<class A> CUDA_CALLABLE friend complex<A> pow(const complex<A> x, int y);
-
-	//log function
-	template<class A> CUDA_CALLABLE friend complex<A> log(complex<A> x);
-
-	//exp function
-	template<class A> CUDA_CALLABLE friend complex<A> exp(complex<A> x);
-
-	//sqrt function
-	template<class A> CUDA_CALLABLE friend complex<A> sqrt(complex<A> x);
-
-	//trigonometric functions
-	template<class A> CUDA_CALLABLE friend complex<A> sin(complex<A> x);
-
-	template<class A> CUDA_CALLABLE friend complex<A> cos(complex<A> x);*/
-
-};
-
-//addition
-template<typename T>
-CUDA_CALLABLE static rtcomplex<T> operator+(const double a, const rtcomplex<T> b)
-{
-    return rtcomplex<T>(a + b.r, b.i);
-}
-
-//subtraction with a real value
-template<typename T>
-CUDA_CALLABLE static rtcomplex<T> operator-(const double a, const rtcomplex<T> b)
-{
-    return rtcomplex<T>(a - b.r, -b.i);
-}
-
-//minus sign
-template<typename T>
-CUDA_CALLABLE static rtcomplex<T> operator-(const rtcomplex<T> &rhs)
-{
-    return rtcomplex<T>(-rhs.r, -rhs.i);
-}
-
-//multiply a T value by a complex value
-template<typename T>
-CUDA_CALLABLE static rtcomplex<T> operator*(const double a, const rtcomplex<T> b)
-{
-    return rtcomplex<T>(a * b.r, a * b.i);
-}
-
-//divide a T value by a complex value
-template<typename T>
-CUDA_CALLABLE static rtcomplex<T> operator/(const double a, const rtcomplex<T> b)
-{
-    //return complex<T>(a * b.r, a * b.i);
-    rtcomplex<T> result;
-
-    T denom = b.r * b.r + b.i * b.i;
-
-    result.r = (a * b.r) / denom;
-    result.i = -(a * b.i) / denom;
-
-    return result;
-}
-
-//POW function
-/*template<typename T>
-CUDA_CALLABLE static complex<T> pow(complex<T> x, int y)
-{
-	return x.pow(y);
-}*/
-
-template<typename T>
-CUDA_CALLABLE static rtcomplex<T> pow(rtcomplex<T> x, T y)
-{
-	return x.pow(y);
-}
-
-//log function
-template<typename T>
-CUDA_CALLABLE static rtcomplex<T> log(rtcomplex<T> x)
-{
-	return x.log();
-}
-
-//exp function
-template<typename T>
-CUDA_CALLABLE static rtcomplex<T> exp(rtcomplex<T> x)
-{
-	return x.exp();
-}
-
-//sqrt function
-template<typename T>
-CUDA_CALLABLE static rtcomplex<T> sqrt(rtcomplex<T> x)
-{
-	return x.sqrt();
-}
-
-
-template <typename T>
-CUDA_CALLABLE static T abs(rtcomplex<T> a)
-{
-    return a.abs();
-}
-
-template <typename T>
-CUDA_CALLABLE static T real(rtcomplex<T> a)
-{
-    return a.r;
-}
-
-template <typename T>
-CUDA_CALLABLE static T imag(rtcomplex<T> a)
-{
-    return a.i;
-}
-
-//trigonometric functions
-template<class A>
-CUDA_CALLABLE rtcomplex<A> sin(const rtcomplex<A> x)
-{
-	rtcomplex<A> result;
-	result.r = std::sin(x.r) * std::cosh(x.i);
-	result.i = std::cos(x.r) * std::sinh(x.i);
-
-	return result;
-}
-
-template<class A>
-CUDA_CALLABLE rtcomplex<A> cos(const rtcomplex<A> x)
-{
-	rtcomplex<A> result;
-	result.r = std::cos(x.r) * std::cosh(x.i);
-	result.i = -(std::sin(x.r) * std::sinh(x.i));
-
-	return result;
-}
-
-
-
-}	//end RTS namespace
-
-template<class A>
-std::ostream& operator<<(std::ostream& os, rts::rtcomplex<A> x)
-{
-    os<<x.toStr();
-    return os;
-}
-
-
-
-#endif
diff --git a/rts/rtsComplex.h b/rts/rtsComplex.h
new file mode 100644
index 0000000..1a338f6
--- /dev/null
+++ b/rts/rtsComplex.h
@@ -0,0 +1,435 @@
+/*RTS Complex number class.  This class is CUDA compatible,
+and can therefore be used in CUDA code and on CUDA devices.
+*/
+
+#ifndef RTS_COMPLEX
+#define RTS_COMPLEX
+
+#include "cuda_callable.h"
+#include <cmath>
+#include <string>
+#include <sstream>
+#include <iostream>
+
+namespace rts
+{
+
+template <class T>
+struct rtsComplex
+{
+    T r, i;
+
+    //default constructor
+    CUDA_CALLABLE rtsComplex()
+    {
+        r = 0.0;
+		i = 0.0;
+    }
+
+	//access methods
+	CUDA_CALLABLE T real()
+	{
+		return r;
+	}
+
+	CUDA_CALLABLE T real(T r_val)
+	{
+		r = r_val;
+		return r_val;
+	}
+
+	CUDA_CALLABLE T imag()
+	{
+		return i;
+	}
+	CUDA_CALLABLE T imag(T i_val)
+	{
+		i = i_val;
+		return i_val;
+	}
+
+    //constructor when given real and imaginary values
+    CUDA_CALLABLE rtsComplex(T r, T i)
+    {
+        this->r = r;
+        this->i = i;
+    }
+
+    //return the current value multiplied by i
+    CUDA_CALLABLE rtsComplex<T> imul()
+    {
+        rtsComplex<T> result;
+        result.r = -i;
+        result.i = r;
+
+        return result;
+    }
+
+	//ARITHMETIC OPERATORS--------------------
+
+    //binary + operator (returns the result of adding two complex values)
+    CUDA_CALLABLE rtsComplex<T> operator+ (const rtsComplex<T> rhs)
+    {
+        rtsComplex<T> result;
+        result.r = r + rhs.r;
+        result.i = i + rhs.i;
+        return result;
+    }
+
+	CUDA_CALLABLE rtsComplex<T> operator+ (const T rhs)
+    {
+        rtsComplex<T> result;
+        result.r = r + rhs;
+        result.i = i;
+        return result;
+    }
+
+    //binary - operator (returns the result of adding two complex values)
+    CUDA_CALLABLE rtsComplex<T> operator- (const rtsComplex<T> rhs)
+    {
+        rtsComplex<T> result;
+        result.r = r - rhs.r;
+        result.i = i - rhs.i;
+        return result;
+    }
+
+    //binary - operator (returns the result of adding two complex values)
+    CUDA_CALLABLE rtsComplex<T> operator- (const T rhs)
+    {
+        rtsComplex<T> result;
+        result.r = r - rhs;
+        result.i = i;
+        return result;
+    }
+
+    //binary MULTIPLICATION operators (returns the result of multiplying complex values)
+    CUDA_CALLABLE rtsComplex<T> operator* (const rtsComplex<T> rhs)
+    {
+        rtsComplex<T> result;
+        result.r = r * rhs.r - i * rhs.i;
+        result.i = r * rhs.i + i * rhs.r;
+        return result;
+    }
+    CUDA_CALLABLE rtsComplex<T> operator* (const T rhs)
+    {
+        return rtsComplex<T>(r * rhs, i * rhs);
+    }
+
+    //binary DIVISION operators (returns the result of dividing complex values)
+    CUDA_CALLABLE rtsComplex<T> operator/ (const rtsComplex<T> rhs)
+    {
+        rtsComplex<T> result;
+        T denom = rhs.r * rhs.r + rhs.i * rhs.i;
+        result.r = (r * rhs.r + i * rhs.i) / denom;
+        result.i = (- r * rhs.i + i * rhs.r) / denom;
+
+        return result;
+    }
+    CUDA_CALLABLE rtsComplex<T> operator/ (const T rhs)
+    {
+        return rtsComplex<T>(r / rhs, i / rhs);
+    }
+
+    //ASSIGNMENT operators-----------------------------------
+    CUDA_CALLABLE rtsComplex<T> & operator=(const rtsComplex<T> &rhs)
+    {
+        //check for self-assignment
+        if(this != &rhs)
+        {
+            this->r = rhs.r;
+            this->i = rhs.i;
+        }
+        return *this;
+    }
+    CUDA_CALLABLE rtsComplex<T> & operator=(const T &rhs)
+    {
+        this->r = rhs;
+        this->i = 0;
+
+		return *this;
+    }
+
+    //arithmetic assignment operators
+    CUDA_CALLABLE rtsComplex<T> operator+=(const rtsComplex<T> &rhs)
+    {
+		*this = *this + rhs;
+        return *this;
+    }
+    CUDA_CALLABLE rtsComplex<T> operator+=(const T &rhs)
+    {
+		*this = *this + rhs;
+        return *this;
+    }
+
+    CUDA_CALLABLE rtsComplex<T> operator*=(const rtsComplex<T> &rhs)
+    {
+		*this = *this * rhs;
+        return *this;
+    }
+	CUDA_CALLABLE rtsComplex<T> operator*=(const T &rhs)
+    {
+		*this = *this * rhs;
+        return *this;
+    }
+	//divide and assign
+	CUDA_CALLABLE rtsComplex<T> operator/=(const rtsComplex<T> &rhs)
+    {
+		*this = *this / rhs;
+        return *this;
+    }
+    CUDA_CALLABLE rtsComplex<T> operator/=(const T &rhs)
+    {
+		*this = *this / rhs;
+        return *this;
+    }
+
+    //absolute value operator (returns the absolute value of the complex number)
+	CUDA_CALLABLE T abs()
+	{
+		return std::sqrt(r * r + i * i);
+	}
+
+	CUDA_CALLABLE rtsComplex<T> log()
+	{
+        rtsComplex<T> result;
+        result.r = std::log(std::sqrt(r * r + i * i));
+        result.i = std::atan2(i, r);
+
+
+        return result;
+	}
+
+	CUDA_CALLABLE rtsComplex<T> exp()
+	{
+        rtsComplex<T> result;
+
+        T e_r = std::exp(r);
+        result.r = e_r * std::cos(i);
+        result.i = e_r * std::sin(i);
+
+        return result;
+	}
+
+	/*CUDA_CALLABLE complex<T> pow(int y)
+	{
+
+        return pow((double)y);
+	}*/
+
+	CUDA_CALLABLE rtsComplex<T> pow(T y)
+	{
+        rtsComplex<T> result;
+
+        result = log() * y;
+
+        return result.exp();
+	}
+
+	CUDA_CALLABLE rtsComplex<T> sqrt()
+	{
+		rtsComplex<T> result;
+
+		//convert to polar coordinates
+		T a = std::sqrt(r*r + i*i);
+		T theta = std::atan2(i, r);
+
+		//find the square root
+		T a_p = std::sqrt(a);
+		T theta_p = theta/2.0;
+
+		//convert back to cartesian coordinates
+		result.r = a_p * std::cos(theta_p);
+		result.i = a_p * std::sin(theta_p);
+
+		return result;
+	}
+
+	std::string toStr()
+	{
+		std::stringstream ss;
+		ss<<"("<<r<<","<<i<<")";
+
+		return ss.str();
+	}
+
+	//COMPARISON operators
+	CUDA_CALLABLE bool operator==(rtsComplex<T> rhs)
+	{
+        if(r == rhs.r && i == rhs.i)
+            return true;
+        return false;
+    }
+
+    CUDA_CALLABLE bool operator==(T rhs)
+	{
+        if(r == rhs && i == (T)0.0)
+            return true;
+        return false;
+    }
+
+	/*//FRIEND functions
+    //unary minus operator (for negating the complex number)
+	template<class A> CUDA_CALLABLE friend complex<A> operator-(const complex<A> &rhs);
+
+	//multiplication by T values when the complex number isn't on the left hand side
+	template<class A> CUDA_CALLABLE friend complex<A> operator*(const A a, const complex<A> b);
+
+	//division by T values when the complex number isn't on the left hand side
+	template<class A> CUDA_CALLABLE friend complex<A> operator/(const A a, const complex<A> b);
+
+	//POW function
+	//template<class A> CUDA_CALLABLE friend complex<A> pow(const complex<A> x, T y);
+	template<class A> CUDA_CALLABLE friend complex<A> pow(const complex<A> x, int y);
+
+	//log function
+	template<class A> CUDA_CALLABLE friend complex<A> log(complex<A> x);
+
+	//exp function
+	template<class A> CUDA_CALLABLE friend complex<A> exp(complex<A> x);
+
+	//sqrt function
+	template<class A> CUDA_CALLABLE friend complex<A> sqrt(complex<A> x);
+
+	//trigonometric functions
+	template<class A> CUDA_CALLABLE friend complex<A> sin(complex<A> x);
+
+	template<class A> CUDA_CALLABLE friend complex<A> cos(complex<A> x);*/
+
+};
+
+}	//end RTS namespace
+
+//addition
+template<typename T>
+CUDA_CALLABLE static rts::rtsComplex<T> operator+(const double a, const rts::rtsComplex<T> b)
+{
+    return rts::rtsComplex<T>(a + b.r, b.i);
+}
+
+//subtraction with a real value
+template<typename T>
+CUDA_CALLABLE static rts::rtsComplex<T> operator-(const double a, const rts::rtsComplex<T> b)
+{
+    return rts::rtsComplex<T>(a - b.r, -b.i);
+}
+
+//minus sign
+template<typename T>
+CUDA_CALLABLE static rts::rtsComplex<T> operator-(const rts::rtsComplex<T> &rhs)
+{
+    return rts::rtsComplex<T>(-rhs.r, -rhs.i);
+}
+
+//multiply a T value by a complex value
+template<typename T>
+CUDA_CALLABLE static rts::rtsComplex<T> operator*(const double a, const rts::rtsComplex<T> b)
+{
+    return rts::rtsComplex<T>(a * b.r, a * b.i);
+}
+
+//divide a T value by a complex value
+template<typename T>
+CUDA_CALLABLE static rts::rtsComplex<T> operator/(const double a, const rts::rtsComplex<T> b)
+{
+    //return complex<T>(a * b.r, a * b.i);
+    rts::rtsComplex<T> result;
+
+    T denom = b.r * b.r + b.i * b.i;
+
+    result.r = (a * b.r) / denom;
+    result.i = -(a * b.i) / denom;
+
+    return result;
+}
+
+//POW function
+/*template<typename T>
+CUDA_CALLABLE static complex<T> pow(complex<T> x, int y)
+{
+	return x.pow(y);
+}*/
+
+template<typename T>
+CUDA_CALLABLE static rts::rtsComplex<T> pow(rts::rtsComplex<T> x, T y)
+{
+	return x.pow(y);
+}
+
+//log function
+template<typename T>
+CUDA_CALLABLE static rts::rtsComplex<T> log(rts::rtsComplex<T> x)
+{
+	return x.log();
+}
+
+//exp function
+template<typename T>
+CUDA_CALLABLE static rts::rtsComplex<T> exp(rts::rtsComplex<T> x)
+{
+	return x.exp();
+}
+
+//sqrt function
+template<typename T>
+CUDA_CALLABLE static rts::rtsComplex<T> sqrt(rts::rtsComplex<T> x)
+{
+	return x.sqrt();
+}
+
+
+template <typename T>
+CUDA_CALLABLE static T abs(rts::rtsComplex<T> a)
+{
+    return a.abs();
+}
+
+template <typename T>
+CUDA_CALLABLE static T real(rts::rtsComplex<T> a)
+{
+    return a.r;
+}
+
+//template <typename T>
+CUDA_CALLABLE static float real(float a)
+{
+    return a;
+}
+
+template <typename T>
+CUDA_CALLABLE static T imag(rts::rtsComplex<T> a)
+{
+    return a.i;
+}
+
+//trigonometric functions
+template<class A>
+CUDA_CALLABLE rts::rtsComplex<A> sin(const rts::rtsComplex<A> x)
+{
+	rts::rtsComplex<A> result;
+	result.r = std::sin(x.r) * std::cosh(x.i);
+	result.i = std::cos(x.r) * std::sinh(x.i);
+
+	return result;
+}
+
+template<class A>
+CUDA_CALLABLE rts::rtsComplex<A> cos(const rts::rtsComplex<A> x)
+{
+	rts::rtsComplex<A> result;
+	result.r = std::cos(x.r) * std::cosh(x.i);
+	result.i = -(std::sin(x.r) * std::sinh(x.i));
+
+	return result;
+}
+
+
+template<class A>
+std::ostream& operator<<(std::ostream& os, rts::rtsComplex<A> x)
+{
+    os<<x.toStr();
+    return os;
+}
+
+
+
+#endif
diff --git a/rts/rtsMatrix.h b/rts/rtsMatrix.h
new file mode 100644
index 0000000..06e2a84
--- /dev/null
+++ b/rts/rtsMatrix.h
@@ -0,0 +1,91 @@
+#ifndef RTS_MATRIX_H
+#define RTS_MATRIX_H
+
+//#include "rts/vector.h"
+#include <string.h>
+#include <iostream>
+
+namespace rts
+{
+
+template <class T, int N>
+struct rtsMatrix
+{
+	//the matrix will be stored in column-major order (compatible with OpenGL)
+	T M[N*N];
+
+	rtsMatrix()
+	{
+		for(int r=0; r<N; r++)
+			for(int c=0; c<N; c++)
+				if(r == c)
+					(*this)(r, c) = 1;
+				else
+					(*this)(r, c) = 0;
+	}
+
+	T& operator()(int row, int col)
+	{
+		return M[col * N + row];
+	}
+
+	rtsMatrix<T, N> operator=(T rhs)
+	{
+		int Nsq = N*N;
+		for(int i=0; i<Nsq; i++)
+			M[i] = rhs;
+
+		return *this;
+	}
+
+	/*rtsMatrix<T, N> operator=(rtsMatrix<T, N> rhs)
+	{
+		for(int i=0; i<N; i++)
+			M[i] = rhs.M[i];
+
+		return *this;
+	}*/
+
+	rtsVector<T, N> operator*(rtsVector<T, N> rhs)
+	{
+		rtsVector<T, N> result;
+
+		for(int r=0; r<N; r++)
+			for(int c=0; c<N; c++)
+				result[r] += (*this)(r, c) * rhs[c];
+
+		return result;
+	}
+
+	std::string toStr()
+	{
+		std::stringstream ss;
+
+		for(int r = 0; r < N; r++)
+		{
+			ss<<"| ";
+			for(int c=0; c<N; c++)
+			{
+				ss<<(*this)(r, c)<<" ";
+			}
+			ss<<"|"<<std::endl;
+		}
+
+		return ss.str();
+	}
+
+
+
+
+};
+
+}	//end namespace rts
+
+template <typename T, int N>
+std::ostream& operator<<(std::ostream& os, rts::rtsMatrix<T, N> M)
+{
+    os<<M.toStr();
+    return os;
+}
+
+#endif
diff --git a/rts/rtsPoint.h b/rts/rtsPoint.h
new file mode 100644
index 0000000..68f0701
--- /dev/null
+++ b/rts/rtsPoint.h
@@ -0,0 +1,124 @@
+#ifndef RTS_rtsPoint_H
+#define RTS_rtsPoint_H
+
+//#include "rts/vector.h"
+#include <string.h>
+
+namespace rts
+{
+
+template <class T, int N>
+struct rtsPoint
+{
+	T p[N];
+
+	CUDA_CALLABLE rtsPoint()
+	{
+
+	}
+
+	//efficiency constructor, makes construction easier for 1D-4D vectors
+	CUDA_CALLABLE rtsPoint(T x, T y = (T)0.0, T z = (T)0.0, T w = (T)0.0)
+	{
+		if(N >= 1)
+			p[0] = x;
+		if(N >= 2)
+			p[1] = y;
+		if(N >= 3)
+			p[2] = z;
+		if(N >= 4)
+			p[3] = w;
+	}
+
+	//arithmetic operators
+	CUDA_CALLABLE rts::rtsPoint<T, N> operator+(rts::rtsVector<T, N> v)
+	{
+        rts::rtsPoint<T, N> r;
+
+        //calculate the position of the resulting rtsPoint
+        for(int i=0; i<N; i++)
+            r.p[i] = p[i] + v.v[i];
+
+        return r;
+	}
+	CUDA_CALLABLE rts::rtsPoint<T, N> operator-(rts::rtsVector<T, N> v)
+	{
+        rts::rtsPoint<T, N> r;
+
+        //calculate the position of the resulting rtsPoint
+        for(int i=0; i<N; i++)
+            r.p[i] = p[i] - v.v[i];
+
+        return r;
+	}
+	CUDA_CALLABLE rts::rtsVector<T, N> operator-(rts::rtsPoint<T, N> rhs)
+	{
+        rts::rtsVector<T, N> r;
+
+        //calculate the position of the resulting rtsPoint
+        for(int i=0; i<N; i++)
+            r.v[i] = p[i] - rhs.p[i];
+
+        return r;
+	}
+	CUDA_CALLABLE rts::rtsPoint<T, N> operator*(T rhs)
+	{
+        rts::rtsPoint<T, N> r;
+
+        //calculate the position of the resulting rtsPoint
+        for(int i=0; i<N; i++)
+            r.p[i] = p[i] * rhs;
+
+        return r;
+	}
+
+	CUDA_CALLABLE rtsPoint(const T(&data)[N])
+	{
+		memcpy(p, data, sizeof(T) * N);
+	}
+
+	std::string toStr()
+	{
+		std::stringstream ss;
+
+		ss<<"(";
+		for(int i=0; i<N; i++)
+		{
+			ss<<p[i];
+			if(i != N-1)
+				ss<<", ";
+		}
+		ss<<")";
+
+		return ss.str();
+	}
+
+	//bracket operator
+	CUDA_CALLABLE T& operator[](int i)
+	{
+        return p[i];
+    }
+
+};
+
+}	//end namespace rts
+
+template <typename T, int N>
+std::ostream& operator<<(std::ostream& os, rts::rtsPoint<T, N> p)
+{
+    os<<p.toStr();
+    return os;
+}
+
+//arithmetic
+template <typename T, int N>
+CUDA_CALLABLE rts::rtsPoint<T, N> operator*(T lhs, rts::rtsPoint<T, N> rhs)
+{
+    rts::rtsPoint<T, N> r;
+
+    return rhs * lhs;
+}
+
+
+
+#endif
diff --git a/rts/rtsProgressBar.h b/rts/rtsProgressBar.h
new file mode 100644
index 0000000..f62f7e0
--- /dev/null
+++ b/rts/rtsProgressBar.h
@@ -0,0 +1,33 @@
+#ifndef RTS_PROGRESSBAR_H
+#define RTS_PROGRESSBAR_H
+
+#include <iostream>
+#include <sstream>
+using namespace std;
+
+static void rtsProgressBar(int percent)
+{
+	stringstream bar;
+	static int x = 0;
+	string slash[4];
+	slash[0] = "\\";
+	slash[1] = "-";
+	slash[2] = "/";
+	slash[3] = "|";
+	bar<<"[";
+	for(int i=0; i<40; i++)
+	{
+		if(percent > (float)i/(float)40 *100)
+			bar << "*";
+		else
+			bar << " ";
+	}
+	bar<<"]";
+	cout << "\r"; // carriage return back to beginning of line
+	cout << bar.str() << " " << slash[x] << " " << percent << " %"; // print the bars and percentage
+	x++; // increment to make the slash appear to rotate
+	if(x == 4)
+	x = 0; // reset slash animation
+}
+
+#endif
diff --git a/rts/rtsQuad.h b/rts/rtsQuad.h
new file mode 100644
index 0000000..1e3f24b
--- /dev/null
+++ b/rts/rtsQuad.h
@@ -0,0 +1,158 @@
+#ifndef RTS_RECT_H
+#define RTS_RECT_H
+
+//enable CUDA_CALLABLE macro
+#include "rts/cuda_callable.h"
+#include "rts/rtsVector.h"
+#include "rts/rtsPoint.h"
+#include <iostream>
+
+namespace rts{
+
+//template for a rtsQuadangle class in ND space
+template <class T, int N>
+struct rtsQuad
+{
+	/*
+		C------------------>O
+		^                   ^
+		|                   |
+		Y                   |
+		|                   |
+		|                   |
+		A---------X-------->B
+	*/
+
+	/*T A[N];
+	T B[N];
+	T C[N];*/
+
+	rts::rtsPoint<T, N> A;
+	rts::rtsVector<T, N> X;
+	rts::rtsVector<T, N> Y;
+
+
+	CUDA_CALLABLE rtsQuad()
+	{
+
+	}
+
+	CUDA_CALLABLE rtsQuad(rtsPoint<T, N> a, rtsPoint<T, N> b, rtsPoint<T, N> c)
+	{
+
+		A = a;
+		X = b - a;
+		Y = c - a;
+
+	}
+
+	CUDA_CALLABLE rtsQuad(rts::rtsPoint<T, N> pMin, rts::rtsPoint<T, N> pMax, rts::rtsVector<T, N> normal)
+	{
+
+        //assign the corner rtsPoint
+        A = pMin;
+
+        //compute the vector from pMin to pMax
+        rts::rtsVector<T, 3> v0;
+        v0 = pMax - pMin;
+
+        //compute the cross product of A and the plane normal
+        rts::rtsVector<T, 3> v1;
+        v1 = v0.cross(normal);
+
+
+        //calculate rtsPoint B
+        rts::rtsPoint<T, 3> B;
+        B = A + v0 * 0.5 + v1 * 0.5;
+
+        //calculate rtsPoint C
+        rts::rtsPoint<T, 3> C;
+        C = A  + v0 * 0.5 - v1 * 0.5;
+
+        //calculate X and Y
+        X = B - A;
+        Y = C - A;
+
+
+
+
+	}
+
+	/*CUDA_CALLABLE rtsQuad(rts::rtsPoint<T, N> p, rts::vector<T, N> x, rts::vector<T, N> y, T sx, T sy)
+	{
+		//This constructor creates a rtsQuad given a position, orientation, and size
+		//	p	= center position of the rtsQuad
+		//	x	= x-axis for the rtsQuadangle
+		//	y	= y-axis for the rtsQuadangle
+		//	sx	= size of the rtsQuad along the A-B axis
+		//	sy	= size of the rtsQuad along the A-C axis
+
+		//normalize x and y
+		rts::vector<T, N> nx = x.norm();
+		rts::vector<T, N> ny = y.norm();
+
+		//compute X and Y
+		X = sx * x;
+		Y = sy * y;
+
+		//compute A
+		A = p - 0.5 * X - 0.5 * Y;
+
+	}*/
+
+	CUDA_CALLABLE rts::rtsPoint<T, N> p(T a, T b)
+	{
+		rts::rtsPoint<T, N> result;
+		//given the two parameters a, b = [0 1], returns the position in world space
+		result = A + X * a + Y * b;
+		
+		return result;
+	}
+
+	CUDA_CALLABLE rts::rtsPoint<T, N> operator()(T a, T b)
+	{
+		return p(a, b);
+	}
+
+	std::string toStr()
+	{
+		std::stringstream ss;
+
+		ss<<"A = "<<A<<std::endl;
+		ss<<"B = "<<A + X<<std::endl;
+		ss<<"C = "<<A + X + Y<<std::endl;
+		ss<<"D = "<<A + Y<<std::endl;
+
+        return ss.str();
+
+	}
+
+	CUDA_CALLABLE rtsQuad<T, N> operator*(T rhs)
+	{
+		//scales the plane by a scalar value
+
+		//compute the center rtsPoint
+		rts::rtsPoint<T, N> c = A + X*0.5 + Y*0.5;
+
+		//create the new rtsQuadangle
+		rtsQuad<T, N> result;
+		result.X = X * rhs;
+		result.Y = Y * rhs;
+		result.A = c - result.X*0.5 - result.Y*0.5;
+
+		return result;
+
+	}
+};
+
+}	//end namespace rts
+
+template <typename T, int N>
+std::ostream& operator<<(std::ostream& os, rts::rtsQuad<T, N> R)
+{
+    os<<R.toStr();
+    return os;
+}
+
+
+#endif
\ No newline at end of file
diff --git a/rts/rtsQuaternion.h b/rts/rtsQuaternion.h
new file mode 100644
index 0000000..f6ded4c
--- /dev/null
+++ b/rts/rtsQuaternion.h
@@ -0,0 +1,153 @@
+#include "rtsMatrix.h"
+
+#ifndef RTS_QUATERNION_H
+#define RTS_QUATERNION_H
+
+namespace rts{
+
+template<typename T>
+class rtsQuaternion
+{
+public:
+	T w;
+	T x;
+	T y;
+	T z;
+
+	void normalize();
+	void CreateRotation(T theta, T axis_x, T axis_y, T axis_z);
+	void CreateRotation(T theta, rtsVector<T, 3> axis);
+	rtsQuaternion<T> operator*(rtsQuaternion<T> &rhs);
+	rtsMatrix<T, 3> toMatrix();
+	
+
+	rtsQuaternion();
+	rtsQuaternion(T w, T x, T y, T z);
+
+};
+
+template<typename T>
+void rtsQuaternion<T>::normalize()
+{
+	double length=sqrt(w*w + x*x + y*y + z*z);
+	w=w/length;
+	x=x/length;
+	y=y/length;
+	z=z/length;
+}
+
+template<typename T>
+void rtsQuaternion<T>::CreateRotation(T theta, T axis_x, T axis_y, T axis_z)
+{
+	//assign the given Euler rotation to this quaternion
+	w = cos(theta/2.0);
+	x = axis_x*sin(theta/2.0);
+	y = axis_y*sin(theta/2.0);
+	z = axis_z*sin(theta/2.0);
+}
+
+template<typename T>
+void rtsQuaternion<T>::CreateRotation(T theta, rtsVector<T, 3> axis)
+{
+	CreateRotation(theta, axis[0], axis[1], axis[2]);
+}
+
+template<typename T>
+rtsQuaternion<T> rtsQuaternion<T>::operator *(rtsQuaternion<T> &param)
+{
+	float A, B, C, D, E, F, G, H;
+
+
+	A = (w + x)*(param.w + param.x);
+	B = (z - y)*(param.y - param.z);
+	C = (w - x)*(param.y + param.z); 
+	D = (y + z)*(param.w - param.x);
+	E = (x + z)*(param.x + param.y);
+	F = (x - z)*(param.x - param.y);
+	G = (w + y)*(param.w - param.z);
+	H = (w - y)*(param.w + param.z);
+
+	rtsQuaternion<T> result;
+	result.w = B + (-E - F + G + H) /2;
+	result.x = A - (E + F + G + H)/2; 
+	result.y = C + (E - F + G - H)/2; 
+	result.z = D + (E - F - G + H)/2;
+
+	return result;
+}
+
+template<typename T>
+rtsMatrix<T, 3> rtsQuaternion<T>::toMatrix()
+{
+	rtsMatrix<T, 3> result;
+
+
+    double wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2;
+
+
+    // calculate coefficients
+    x2 = x + x; y2 = y + y;
+    z2 = z + z;
+    xx = x * x2; xy = x * y2; xz = x * z2;
+    yy = y * y2; yz = y * z2; zz = z * z2;
+    wx = w * x2; wy = w * y2; wz = w * z2;
+
+	result(0, 0) = 1.0 - (yy + zz);
+	result(0, 1) = xy - wz;
+    //m[0][0] = 1.0 - (yy + zz); m[1][0] = xy - wz;
+
+	result(0, 2) = xz + wy;
+	//result(0, 3) = 0.0;
+    //m[2][0] = xz + wy; m[3][0] = 0.0;
+
+	result(1, 0) = xy + wz;
+	result(1, 1) = 1.0 - (xx + zz);
+    //m[0][1] = xy + wz; m[1][1] = 1.0 - (xx + zz);
+
+	result(1, 2) = yz - wx;
+	//result(1, 3) = 0.0;
+    //m[2][1] = yz - wx; m[3][1] = 0.0;
+
+	result(2, 0) = xz - wy;
+	result(2, 1) = yz + wx;
+    //m[0][2] = xz - wy; m[1][2] = yz + wx;
+
+	result(2, 2) = 1.0 - (xx + yy);
+	//result(3, 2) = 0.0;
+    //m[2][2] = 1.0 - (xx + yy); m[3][2] = 0.0;
+
+	/*result(3, 0) = 0.0;
+	result(3, 1) = 0.0;
+	result(3, 2) = 0.0;
+	result(3, 3) = 1.0;*/
+    //m[0][3] = 0; m[1][3] = 0;
+    //m[2][3] = 0; m[3][3] = 1;
+	/*
+	double* orientationmatrix=(double*)m;
+	char c;
+
+
+	double* result=new double[16];
+	double* array=(double*)m;
+	for(int i=0; i<16; i++)
+		result[i]=array[i];
+	*/
+
+	return result;
+}
+
+template<typename T>
+rtsQuaternion<T>::rtsQuaternion()
+{
+	w=0.0; x=0.0; y=0.0; z=0.0;
+}
+
+template<typename T>
+rtsQuaternion<T>::rtsQuaternion(T c, T i, T j, T k)
+{
+	w=c;  x=i;  y=j;  z=k;
+}
+
+}	//end rts namespace
+
+#endif
diff --git a/rts/rtsVector.h b/rts/rtsVector.h
new file mode 100644
index 0000000..f9a3957
--- /dev/null
+++ b/rts/rtsVector.h
@@ -0,0 +1,213 @@
+#ifndef RTS_VECTOR_H
+#define RTS_VECTOR_H
+
+#include <iostream>
+//#include "rts/point.h"
+
+namespace rts
+{
+
+
+
+template <class T, int N>
+struct rtsVector
+{
+	T v[N];
+
+	CUDA_CALLABLE rtsVector()
+	{
+		//memset(v, 0, sizeof(T) * N);
+		for(int i=0; i<N; i++)
+			v[i] = 0;
+	}
+
+	//efficiency constructor, makes construction easier for 1D-4D vectors
+	CUDA_CALLABLE rtsVector(T x, T y = (T)0.0, T z = (T)0.0, T w = (T)0.0)
+	{
+		if(N >= 1)
+			v[0] = x;
+		if(N >= 2)
+			v[1] = y;
+		if(N >= 3)
+			v[2] = z;
+		if(N >= 4)
+			v[3] = w;
+	}
+
+	CUDA_CALLABLE rtsVector(const T(&data)[N])
+	{
+		memcpy(v, data, sizeof(T) * N);
+	}
+
+	CUDA_CALLABLE T len()
+	{
+        //compute and return the vector length
+        T sum_sq = (T)0;
+        for(int i=0; i<N; i++)
+        {
+            sum_sq += v[i] * v[i];
+        }
+        return std::sqrt(sum_sq);
+
+	}
+
+	CUDA_CALLABLE rtsVector<T, N> cart2sph()
+	{
+		//convert the vector from cartesian to spherical coordinates
+		//x, y, z -> r, theta, phi (where theta = 0 to 2*pi)
+
+		rtsVector<T, N> sph;
+		sph[0] = std::sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
+		sph[1] = std::atan2(v[1], v[0]);
+		sph[2] = std::acos(v[2] / sph[0]);
+
+		return sph;
+	}
+
+	CUDA_CALLABLE rtsVector<T, N> sph2cart()
+	{
+		//convert the vector from cartesian to spherical coordinates
+		//r, theta, phi -> x, y, z (where theta = 0 to 2*pi)
+
+		rtsVector<T, N> cart;
+		cart[0] = v[0] * std::cos(v[1]) * std::sin(v[2]);
+		cart[1] = v[0] * std::sin(v[1]) * std::sin(v[2]);
+		cart[2] = v[0] * std::cos(v[2]);
+
+		return cart;
+	}
+
+	CUDA_CALLABLE rtsVector<T, N> norm()
+	{
+        //compute and return the vector norm
+        rtsVector<T, N> result;
+
+        //compute the vector length
+        T l = len();
+
+        //normalize
+        for(int i=0; i<N; i++)
+        {
+            result.v[i] = v[i] / l;
+        }
+
+        return result;
+	}
+
+	CUDA_CALLABLE rtsVector<T, 3> cross(rtsVector<T, 3> rhs)
+	{
+		rtsVector<T, 3> result;
+
+		//compute the cross product (only valid for 3D vectors)
+		result[0] = v[1] * rhs[2] - v[2] * rhs[1];
+		result[1] = v[2] * rhs[0] - v[0] * rhs[2];
+		result[2] = v[0] * rhs[1] - v[1] * rhs[0];
+
+		return result;
+	}
+
+    CUDA_CALLABLE T dot(rtsVector<T, N> rhs)
+    {
+        T result = (T)0;
+
+        for(int i=0; i<N; i++)
+            result += v[i] * rhs.v[i];
+
+        return result;
+
+    }
+
+	//arithmetic
+	CUDA_CALLABLE rtsVector<T, N> operator+(rtsVector<T, N> rhs)
+	{
+        rtsVector<T, N> result;
+
+        for(int i=0; i<N; i++)
+            result.v[i] = v[i] + rhs.v[i];
+
+        return result;
+	}
+	CUDA_CALLABLE rtsVector<T, N> operator-(rtsVector<T, N> rhs)
+	{
+        rtsVector<T, N> result;
+
+        for(int i=0; i<N; i++)
+            result.v[i] = v[i] - rhs.v[i];
+
+        return result;
+	}
+	CUDA_CALLABLE rtsVector<T, N> operator*(T rhs)
+	{
+        rtsVector<T, N> result;
+
+        for(int i=0; i<N; i++)
+            result.v[i] = v[i] * rhs;
+
+        return result;
+	}
+	CUDA_CALLABLE rtsVector<T, N> operator/(T rhs)
+	{
+        rtsVector<T, N> result;
+
+        for(int i=0; i<N; i++)
+            result.v[i] = v[i] / rhs;
+
+        return result;
+	}
+
+	std::string toStr()
+	{
+		std::stringstream ss;
+
+		ss<<"[";
+		for(int i=0; i<N; i++)
+		{
+			ss<<v[i];
+			if(i != N-1)
+				ss<<", ";
+		}
+		ss<<"]";
+
+		return ss.str();
+	}
+
+	//bracket operator
+	CUDA_CALLABLE T& operator[](int i)
+	{
+        return v[i];
+    }
+
+};
+
+
+}	//end namespace rts
+
+template <typename T, int N>
+std::ostream& operator<<(std::ostream& os, rts::rtsVector<T, N> v)
+{
+    os<<v.toStr();
+    return os;
+}
+
+//arithmetic operators
+template <typename T, int N>
+CUDA_CALLABLE rts::rtsVector<T, N> operator-(rts::rtsVector<T, N> v)
+{
+    rts::rtsVector<T, N> r;
+
+    //negate the vector
+    for(int i=0; i<N; i++)
+        r.v[i] = -v.v[i];
+
+    return r;
+}
+
+template <typename T, int N>
+CUDA_CALLABLE rts::rtsVector<T, N> operator*(T lhs, rts::rtsVector<T, N> rhs)
+{
+    rts::rtsVector<T, N> r;
+
+    return rhs * lhs;
+}
+
+#endif
diff --git a/rts/sbessel.h b/rts/sbessel.h
index 46e0d6f..b432408 100644
--- a/rts/sbessel.h
+++ b/rts/sbessel.h
@@ -5,10 +5,11 @@
 
 namespace rts{
 
-#define RTS_BESSEL_CONVERGENCE_FLOAT	0.00045
+#define RTS_BESSEL_CONVERGENCE_FLOAT	0.0145
+#define RTS_BESSEL_MAXIMUM_FLOAT		-1e33
 
 template <typename T>
-CUDA_CALLABLE void sbesselj(int n, rtcomplex<T> x, rtcomplex<T>* j)
+CUDA_CALLABLE void sbesselj(int n, rtsComplex<T> x, rtsComplex<T>* j)
 {
     //compute the first bessel function
     if(n >= 0)
@@ -34,7 +35,7 @@ CUDA_CALLABLE void sbesselj(int n, rtcomplex<T> x, rtcomplex<T>* j)
 }
 
 template <typename T>
-CUDA_CALLABLE void sbessely(int n, rtcomplex<T> x, rtcomplex<T>* y)
+CUDA_CALLABLE void sbessely(int n, rtsComplex<T> x, rtsComplex<T>* y)
 {
     //compute the first bessel function
     if(n >= 0)
@@ -54,14 +55,14 @@ CUDA_CALLABLE void sbessely(int n, rtcomplex<T> x, rtcomplex<T>* y)
 
 //spherical Hankel functions of the first kind
 template <typename T>
-CUDA_CALLABLE void sbesselh1(int n, rtcomplex<T> x, rtcomplex<T>* h)
+CUDA_CALLABLE void sbesselh1(int n, rtsComplex<T> x, rtsComplex<T>* h)
 {
     //compute j_0 and j_1
-    rtcomplex<T> j[2];
+    rtsComplex<T> j[2];
     sbesselj(1, x, j);
 
     //compute y_0 and y_1
-    rtcomplex<T> y[2];
+    rtsComplex<T> y[2];
     sbessely(1, x, y);
 
     //compute the first-order Hhankel function
@@ -82,55 +83,60 @@ CUDA_CALLABLE void sbesselh1(int n, rtcomplex<T> x, rtcomplex<T>* h)
 template <typename T>
 CUDA_CALLABLE void init_sbesselj(T x, T* j)
 {
-	/*if(x < EPSILON_FLOAT)
-	{
-		j[0] = (T)1;
-		j[1] = (T)0;
-		return;
-	}*/
-	
 	//compute the first 2 bessel functions
-	j[0] = std::sin(x) / x;
-
-	j[1] = j[0] / x - std::cos(x) / x;
+	j[0] = sin(x) / x;
 
-	//deal with the degenerate case of x = 0
-	/*if(isnan(j[0]))
-	{
-		j[0] = (T)1;
-		j[1] = (T)0;
-	}*/
-	
+	j[1] = j[0] / x - cos(x) / x;
 }
 
 template <typename T>
-CUDA_CALLABLE void shift_sbesselj(int n, T x, T* j)//, T stability = 1.4)
+CUDA_CALLABLE void init_sbessely(T x, T* y)
 {
+	//compute the first 2 bessel functions
+	y[0] = -cos(x) / x;
 
-	/*if(x < EPSILON_FLOAT)
-	{
-		j[0] = j[1];
-		j[1] = (T)0;
-		return;
-	}*/
+	y[1] = y[0] / x - sin(x) / x;
+}
 
-	T jnew;
+template <typename T>
+CUDA_CALLABLE void shift_sbesselj(int n, T x, T* b)//, T stability = 1.4)
+{
+
+	T bnew;
 
 	//compute the next (order n) Bessel function
-	jnew = ((2 * n - 1) * j[1])/x - j[0];
+	bnew = ((2 * n - 1) * b[1])/x - b[0];
 
 	//if(n > stability*x)
-	if(n > x)
-		if(jnew < RTS_BESSEL_CONVERGENCE_FLOAT)
-			jnew = (T)0;
+	if(n > real(x))
+		if(real(bnew) < RTS_BESSEL_CONVERGENCE_FLOAT)
+			bnew = 0.0;
+
+	//shift and add the new value to the array
+	b[0] = b[1];
+	b[1] = bnew;
+}
+
+template <typename T>
+CUDA_CALLABLE void shift_sbessely(int n, T x, T* b)//, T stability = 1.4)
+{
+
+	T bnew;
+
+	//compute the next (order n) Bessel function
+	bnew = ((2 * n - 1) * b[1])/x - b[0];
+
+	if(bnew < RTS_BESSEL_MAXIMUM_FLOAT ||
+	   (n > x && bnew > 0))
+	{
+		bnew = (T)0;
+		b[1] = (T)0;
+	}
 
-	//deal with the degenerate case of x = 0
-	//if(isnan(jnew))
-	//	jnew = (T)0;
 
 	//shift and add the new value to the array
-	j[0] = j[1];
-	j[1] = jnew;
+	b[0] = b[1];
+	b[1] = bnew;
 }
 
 
diff --git a/rts/vector.h b/rts/vector.h
deleted file mode 100644
index 776d9d6..0000000
--- a/rts/vector.h
+++ /dev/null
@@ -1,211 +0,0 @@
-#ifndef RTS_VECTOR_H
-#define RTS_VECTOR_H
-
-#include <iostream>
-//#include "rts/point.h"
-
-namespace rts
-{
-
-
-
-template <class T, int N>
-struct vector
-{
-	T v[N];
-
-	CUDA_CALLABLE vector()
-	{
-		//memset(v, 0, sizeof(T) * N);
-	}
-
-	//efficiency constructor, makes construction easier for 1D-4D vectors
-	CUDA_CALLABLE vector(T x, T y = (T)0.0, T z = (T)0.0, T w = (T)0.0)
-	{
-		if(N >= 1)
-			v[0] = x;
-		if(N >= 2)
-			v[1] = y;
-		if(N >= 3)
-			v[2] = z;
-		if(N >= 4)
-			v[3] = w;
-	}
-
-	CUDA_CALLABLE vector(const T(&data)[N])
-	{
-		memcpy(v, data, sizeof(T) * N);
-	}
-
-	CUDA_CALLABLE T len()
-	{
-        //compute and return the vector length
-        T sum_sq = (T)0;
-        for(int i=0; i<N; i++)
-        {
-            sum_sq += v[i] * v[i];
-        }
-        return std::sqrt(sum_sq);
-
-	}
-
-	CUDA_CALLABLE vector<T, N> cart2sph()
-	{
-		//convert the vector from cartesian to spherical coordinates
-		//x, y, z -> r, theta, phi
-
-		vector<T, N> sph;
-		sph[0] = std::sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
-		sph[1] = std::atan2(v[1], v[0]);
-		sph[2] = std::acos(v[2] / sph[0]);
-
-		return sph;
-	}
-
-	CUDA_CALLABLE vector<T, N> sph2cart()
-	{
-		//convert the vector from cartesian to spherical coordinates
-		//r, theta, phi -> x, y, z
-
-		vector<T, N> cart;
-		cart[0] = v[0] * std::cos(v[1]) * std::sin(v[2]);
-		cart[1] = v[0] * std::sin(v[1]) * std::sin(v[2]);
-		cart[2] = v[0] * std::cos(v[2]);
-
-		return cart;
-	}
-
-	CUDA_CALLABLE vector<T, N> norm()
-	{
-        //compute and return the vector norm
-        vector<T, N> result;
-
-        //compute the vector length
-        T l = len();
-
-        //normalize
-        for(int i=0; i<N; i++)
-        {
-            result.v[i] = v[i] / l;
-        }
-
-        return result;
-	}
-
-	CUDA_CALLABLE vector<T, 3> cross(vector<T, 3> rhs)
-	{
-		vector<T, 3> result;
-
-		//compute the cross product (only valid for 3D vectors)
-		result[0] = v[1] * rhs[2] - v[2] * rhs[1];
-		result[1] = v[2] * rhs[0] - v[0] * rhs[2];
-		result[2] = v[0] * rhs[1] - v[1] * rhs[0];
-
-		return result;
-	}
-    
-    CUDA_CALLABLE T dot(vector<T, N> rhs)
-    {
-        T result = (T)0;
-
-        for(int i=0; i<N; i++)
-            result += v[i] * rhs.v[i];
-
-        return result;
-    
-    }
-
-	//arithmetic
-	CUDA_CALLABLE vector<T, N> operator+(vector<T, N> rhs)
-	{
-        vector<T, N> result;
-
-        for(int i=0; i<N; i++)
-            result.v[i] = v[i] + rhs.v[i];
-
-        return result;
-	}
-	CUDA_CALLABLE vector<T, N> operator-(vector<T, N> rhs)
-	{
-        vector<T, N> result;
-
-        for(int i=0; i<N; i++)
-            result.v[i] = v[i] - rhs.v[i];
-
-        return result;
-	}
-	CUDA_CALLABLE vector<T, N> operator*(T rhs)
-	{
-        vector<T, N> result;
-
-        for(int i=0; i<N; i++)
-            result.v[i] = v[i] * rhs;
-
-        return result;
-	}
-	CUDA_CALLABLE vector<T, N> operator/(T rhs)
-	{
-        vector<T, N> result;
-
-        for(int i=0; i<N; i++)
-            result.v[i] = v[i] / rhs;
-
-        return result;
-	}
-
-	std::string toStr()
-	{
-		std::stringstream ss;
-
-		ss<<"[";
-		for(int i=0; i<N; i++)
-		{
-			ss<<v[i];
-			if(i != N-1)
-				ss<<", ";
-		}
-		ss<<"]";
-
-		return ss.str();
-	}
-
-	//bracket operator
-	CUDA_CALLABLE T& operator[](int i)
-	{
-        return v[i];
-    }
-
-};
-
-
-}	//end namespace rts
-
-template <typename T, int N>
-std::ostream& operator<<(std::ostream& os, rts::vector<T, N> v)
-{
-    os<<v.toStr();
-    return os;
-}
-
-//arithmetic operators
-template <typename T, int N>
-CUDA_CALLABLE rts::vector<T, N> operator-(rts::vector<T, N> v)
-{
-    rts::vector<T, N> r;
-
-    //negate the vector
-    for(int i=0; i<N; i++)
-        r.v[i] = -v.v[i];
-
-    return r;
-}
-
-template <typename T, int N>
-CUDA_CALLABLE rts::vector<T, N> operator*(T lhs, rts::vector<T, N> rhs)
-{
-    rts::vector<T, N> r;
-
-    return rhs * lhs;
-}
-
-#endif
--
libgit2 0.21.4