diff --git a/lib/poems/POEMSChain.h b/lib/poems/POEMSChain.h
index 9f0db59333..8baa293dbd 100644
--- a/lib/poems/POEMSChain.h
+++ b/lib/poems/POEMSChain.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: PoemsChain.h                                            *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,7 +11,7 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
@@ -19,56 +19,59 @@
 #define POEMSCHAIN_H_
 
 #include "poemslist.h"
+#include <iostream>
 
+namespace POEMS {
 struct ChildRingData {
-	List<int> * childRing;
-	int entranceNodeId;
+   List<int> * childRing;
+   int entranceNodeId;
 };
 
 struct POEMSChain{
-	~POEMSChain(){
-		for(int i = 0; i < childChains.GetNumElements(); i++)
-		{
-			delete childChains(i);
-		}
-        listOfNodes.DeleteValues();
-	}
-	//void printTreeStructure(int tabs);
-	//void getTreeAsList(List<int> * temp);
-	List<int> listOfNodes;
-	List<POEMSChain> childChains;
-	POEMSChain * parentChain;
-	List<ChildRingData> childRings;
-	
-	
-	void printTreeStructure(int tabs){
-		for(int i = 0; i < tabs; i++)
-		{
-			cout << "\t";
-		}
-		cout << "Chain: ";
-		for(int i = 0; i < listOfNodes.GetNumElements(); i++)
-		{
-			cout << *(listOfNodes(i)) << " ";
-		}
-		cout << endl;
-		for(int i = 0; i < childChains.GetNumElements(); i++)
-		{
-			childChains(i)->printTreeStructure(tabs + 1);
-		}
-	}
-	void getTreeAsList(List<int> * temp)
-	{
-		for(int i = 0; i < listOfNodes.GetNumElements(); i++)
-		{
-			int * integer = new int;
-			*integer = *(listOfNodes(i));
-			temp->Append(integer);
-		}
-		for(int i = 0; i < childChains.GetNumElements(); i++)
-		{
-			childChains(i)->getTreeAsList(temp);
-		}
-	}
+   ~POEMSChain(){
+ for(int i = 0; i < childChains.GetNumElements(); i++)
+ {
+    delete childChains(i);
+ }
+   listOfNodes.DeleteValues();
+   }
+   //void printTreeStructure(int tabs);
+   //void getTreeAsList(List<int> * temp);
+   List<int> listOfNodes;
+   List<POEMSChain> childChains;
+   POEMSChain * parentChain;
+   List<ChildRingData> childRings;
+
+
+   void printTreeStructure(int tabs){
+ for(int i = 0; i < tabs; i++)
+ {
+    std::cout << "\t";
+ }
+ std::cout << "Chain: ";
+ for(int i = 0; i < listOfNodes.GetNumElements(); i++)
+ {
+    std::cout << *(listOfNodes(i)) << " ";
+ }
+ std::cout << std::endl;
+ for(int i = 0; i < childChains.GetNumElements(); i++)
+ {
+    childChains(i)->printTreeStructure(tabs + 1);
+ }
+   }
+   void getTreeAsList(List<int> * temp)
+   {
+ for(int i = 0; i < listOfNodes.GetNumElements(); i++)
+ {
+    int * integer = new int;
+    *integer = *(listOfNodes(i));
+    temp->Append(integer);
+ }
+ for(int i = 0; i < childChains.GetNumElements(); i++)
+ {
+    childChains(i)->getTreeAsList(temp);
+ }
+   }
 };
+}
 #endif
diff --git a/lib/poems/SystemProcessor.h b/lib/poems/SystemProcessor.h
index 3be168c34d..a6001d96d0 100644
--- a/lib/poems/SystemProcessor.h
+++ b/lib/poems/SystemProcessor.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: SystemProcessor.h                                       *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,59 +11,61 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
 #ifndef _SYS_PROCESSOR_H_
 #define _SYS_PROCESSOR_H_
+#include <iostream>
 #include "poemslist.h"
 #include "poemstree.h"
 #include "POEMSChain.h"
 
+namespace POEMS {
 
 struct POEMSNode {
-	List<POEMSNode> links;
-	List<bool> taken;
-	int idNumber;
-	bool visited;
-	
-	~POEMSNode(){
-		for(int i = 0; i < taken.GetNumElements(); i++)
-		{
-			delete taken(i);
-		}
-	};	
+        List<POEMSNode> links;
+        List<bool> taken;
+        int idNumber;
+        bool visited;
+
+        ~POEMSNode(){
+                for(int i = 0; i < taken.GetNumElements(); i++)
+                {
+                        delete taken(i);
+                }
+        };
 };
 
 
 class SystemProcessor{
 private:
-	Tree nodes;
-	static void POEMSNodeDelete_cb(void *node) {
-		delete (POEMSNode *) node;
-	}
-	List<POEMSChain> headsOfSystems;
-	List<List<int> > ringsInSystem;
-	POEMSNode * findSingleLink(TreeNode * aNode);
-	POEMSChain * AddNewChain(POEMSNode * currentNode);
-	bool setLinkVisited(POEMSNode * firstNode, POEMSNode * secondNode);
+        Tree nodes;
+        static void POEMSNodeDelete_cb(void *node) {
+                delete (POEMSNode *) node;
+        }
+        List<POEMSChain> headsOfSystems;
+        List<List<int> > ringsInSystem;
+        POEMSNode * findSingleLink(TreeNode * aNode);
+        POEMSChain * AddNewChain(POEMSNode * currentNode);
+        bool setLinkVisited(POEMSNode * firstNode, POEMSNode * secondNode);
 public:
-	SystemProcessor(void);
-	
-	~SystemProcessor(void)	{
-		headsOfSystems.DeleteValues();
-		for(int i = 0; i < ringsInSystem.GetNumElements(); i++)
-		{
-			for(int k = 0; k < ringsInSystem(i)->GetNumElements(); i++)
-			{
-				delete (*ringsInSystem(i))(k);
-			}
-		}
-	};
-	void processArray(int** links, int numLinks);
-	List<POEMSChain> * getSystemData();
-	int getNumberOfHeadChains();
+        SystemProcessor(void);
+
+        ~SystemProcessor(void)  {
+                headsOfSystems.DeleteValues();
+                for(int i = 0; i < ringsInSystem.GetNumElements(); i++)
+                {
+                        for(int k = 0; k < ringsInSystem(i)->GetNumElements(); i++)
+                        {
+                                delete (*ringsInSystem(i))(k);
+                        }
+                }
+        };
+        void processArray(int** links, int numLinks);
+        List<POEMSChain> * getSystemData();
+        int getNumberOfHeadChains();
 };
 
 SystemProcessor::SystemProcessor(void){
@@ -73,145 +75,145 @@ SystemProcessor::SystemProcessor(void){
 
 void SystemProcessor::processArray(int** links, int numLinks)
 {
-	bool * false_var;					//holds the value false; needed because a constant cannot be put into a list; the list requires a
-										//reference.
-	for(int i = 0; i < numLinks; i++)	//go through all the links in the input array
-	{
-		if(!nodes.Find(links[i][0]))	//if the first node in the pair is not found in the storage tree
-		{
-			POEMSNode * newNode = new POEMSNode;	//make a new node
-//			forDeletion.Append(newNode);
-			newNode->idNumber = links[i][0];		//set its ID to the value
-			newNode->visited = false;				//set it to be unvisited
-			nodes.Insert(links[i][0], links[i][0], (void *) newNode);	//and add it to the tree storage structure
-		}
-		if(!nodes.Find(links[i][1]))				//repeat process for the other half of each link
-		{
-			POEMSNode * newNode = new POEMSNode;
-//			forDeletion.Append(newNode);
-			newNode->idNumber = links[i][1];
-			newNode->visited = false;
-			nodes.Insert(links[i][1], links[i][1], (void *) newNode);
-		}
-		POEMSNode * firstNode = (POEMSNode *)nodes.Find(links[i][0]);	//now that we are sure both nodes exist,
-		POEMSNode * secondNode = (POEMSNode *)nodes.Find(links[i][1]);	//we can get both of them out of the tree
-		firstNode->links.Append(secondNode);							//and add the link from the first to the second...
-		false_var = new bool;
-		*false_var = false;												//make a new false boolean to note that the link between these two
-		firstNode->taken.Append(false_var);								//has not already been taken, and append it to the taken list
-		secondNode->links.Append(firstNode);							//repeat process for link from second node to first
-		false_var = new bool;
-		*false_var = false;
-		secondNode->taken.Append(false_var);
-	}	
-	
-	TreeNode * temp = nodes.GetRoot();							//get the root node of the node storage tree
-	POEMSNode * currentNode;
-	do
-	{
-		currentNode = findSingleLink(temp);						//find the start of the next available chain
-		if(currentNode != NULL)
-		{
-			headsOfSystems.Append(AddNewChain(currentNode));							//and add it to the headsOfSystems list of chains
-		}
-	} 
-	while(currentNode != NULL);									//repeat this until all chains have been added		
+        bool * false_var;                                       //holds the value false; needed because a constant cannot be put into a list; the list requires a
+                                                                                //reference.
+        for(int i = 0; i < numLinks; i++)       //go through all the links in the input array
+        {
+                if(!nodes.Find(links[i][0]))    //if the first node in the pair is not found in the storage tree
+                {
+                        POEMSNode * newNode = new POEMSNode;    //make a new node
+//                      forDeletion.Append(newNode);
+                        newNode->idNumber = links[i][0];                //set its ID to the value
+                        newNode->visited = false;                               //set it to be unvisited
+                        nodes.Insert(links[i][0], links[i][0], (void *) newNode);       //and add it to the tree storage structure
+                }
+                if(!nodes.Find(links[i][1]))                            //repeat process for the other half of each link
+                {
+                        POEMSNode * newNode = new POEMSNode;
+//                      forDeletion.Append(newNode);
+                        newNode->idNumber = links[i][1];
+                        newNode->visited = false;
+                        nodes.Insert(links[i][1], links[i][1], (void *) newNode);
+                }
+                POEMSNode * firstNode = (POEMSNode *)nodes.Find(links[i][0]);   //now that we are sure both nodes exist,
+                POEMSNode * secondNode = (POEMSNode *)nodes.Find(links[i][1]);  //we can get both of them out of the tree
+                firstNode->links.Append(secondNode);                                                    //and add the link from the first to the second...
+                false_var = new bool;
+                *false_var = false;                                                                                             //make a new false boolean to note that the link between these two
+                firstNode->taken.Append(false_var);                                                             //has not already been taken, and append it to the taken list
+                secondNode->links.Append(firstNode);                                                    //repeat process for link from second node to first
+                false_var = new bool;
+                *false_var = false;
+                secondNode->taken.Append(false_var);
+        }
+
+        TreeNode * temp = nodes.GetRoot();                                                      //get the root node of the node storage tree
+        POEMSNode * currentNode;
+        do
+        {
+                currentNode = findSingleLink(temp);                                             //find the start of the next available chain
+                if(currentNode != NULL)
+                {
+                        headsOfSystems.Append(AddNewChain(currentNode));                                                        //and add it to the headsOfSystems list of chains
+                }
+        }
+        while(currentNode != NULL);                                                                     //repeat this until all chains have been added
 }
 
 POEMSChain * SystemProcessor::AddNewChain(POEMSNode * currentNode){
-	if(currentNode == NULL)	//Termination condition; if the currentNode is null, then return null
-	{
-		return NULL;
-	}
-	int * tmp;
-	POEMSNode * nextNode = NULL;	//nextNode stores the proposed next node to add to the chain.  this will be checked to make sure no backtracking is occuring before being assigned as the current node.
-	POEMSChain * newChain = new POEMSChain;	//make a new POEMSChain object.  This will be the object returned
+        if(currentNode == NULL) //Termination condition; if the currentNode is null, then return null
+        {
+                return NULL;
+        }
+        int * tmp;
+        POEMSNode * nextNode = NULL;    //nextNode stores the proposed next node to add to the chain.  this will be checked to make sure no backtracking is occuring before being assigned as the current node.
+        POEMSChain * newChain = new POEMSChain; //make a new POEMSChain object.  This will be the object returned
 
-	if(currentNode->links.GetNumElements() == 0)	//if we have no links from this node, then the whole chain is only one node.  Add this node to the chain and return it; mark node as visited for future reference
-	{
-		currentNode->visited = true;
-		tmp = new int;
-		*tmp = currentNode->idNumber;
-		newChain->listOfNodes.Append(tmp);
-		return newChain;
-	}
-	while(currentNode->links.GetNumElements() <= 2)	//we go until we get to a node that branches, or both branches have already been taken both branches can already be taken if a loop with no spurs is found in the input data
-	{
-		currentNode->visited = true;
-		tmp = new int;
-		*tmp = currentNode->idNumber;
-		newChain->listOfNodes.Append(tmp);	//append the current node to the chain & mark as visited
-		//cout << "Appending node " << currentNode->idNumber << " to chain" << endl;
-		nextNode = currentNode->links.GetHeadElement()->value;	//the next node is the first or second value stored in the links array
-																//of the current node.  We get the first value...
-		if(!setLinkVisited(currentNode, nextNode))					//...and see if it points back to where we came from. If it does...
-		{														//either way, we set this link as visited
-			if(currentNode->links.GetNumElements() == 1)		//if it does, then if that is the only link to this node, we're done with the chain, so append the chain to the list and return the newly created chain
-			{
-//				headsOfSystems.Append(newChain);
-				return newChain;
-			}
-			nextNode = currentNode->links.GetHeadElement()->next->value;//follow the other link if there is one, so we go down the chain
-			if(!setLinkVisited(currentNode, nextNode))				//mark link as followed, so we know not to backtrack			
-			{
-	//			headsOfSystems.Append(newChain);
-				return newChain;								//This condition, where no branches have occurred but both links have already
-																//been taken can only occur in a loop with no spurs; add this loop to the
-																//system (currently added as a chain for consistency), and return.
-			}
-		}
-		currentNode = nextNode;									//set the current node to be the next node in the chain
-	}
-	currentNode->visited = true;
-	tmp = new int;
-	*tmp = currentNode->idNumber;
-	newChain->listOfNodes.Append(tmp);		//append the last node before branch (node shared jointly with branch chains)
-																//re-mark as visited, just to make sure
-	ListElement<POEMSNode> * tempNode = currentNode->links.GetHeadElement();	//go through all of the links, one at a time that branch
-	POEMSChain * tempChain = NULL;								//temporary variable to hold data 
-	while(tempNode != NULL)										//when we have followed all links, stop
-	{
-		if(setLinkVisited(tempNode->value, currentNode))		//dont backtrack, or create closed loops
-		{
-			tempChain = AddNewChain(tempNode->value);			//Add a new chain created out of the next node down that link
-			tempChain->parentChain = newChain;					//set the parent to be this chain
-			newChain->childChains.Append(tempChain);			//append the chain to this chain's list of child chains
-		}
-		tempNode = tempNode->next;								//go to process the next chain
-	}
-	//headsOfSystems.Append(newChain);							//append this chain to the system list		
-	return newChain;	
+        if(currentNode->links.GetNumElements() == 0)    //if we have no links from this node, then the whole chain is only one node.  Add this node to the chain and return it; mark node as visited for future reference
+        {
+                currentNode->visited = true;
+                tmp = new int;
+                *tmp = currentNode->idNumber;
+                newChain->listOfNodes.Append(tmp);
+                return newChain;
+        }
+        while(currentNode->links.GetNumElements() <= 2) //we go until we get to a node that branches, or both branches have already been taken both branches can already be taken if a loop with no spurs is found in the input data
+        {
+                currentNode->visited = true;
+                tmp = new int;
+                *tmp = currentNode->idNumber;
+                newChain->listOfNodes.Append(tmp);      //append the current node to the chain & mark as visited
+                //std::cout << "Appending node " << currentNode->idNumber << " to chain" << std::endl;
+                nextNode = currentNode->links.GetHeadElement()->value;  //the next node is the first or second value stored in the links array
+                                                                                                                                //of the current node.  We get the first value...
+                if(!setLinkVisited(currentNode, nextNode))                                      //...and see if it points back to where we came from. If it does...
+                {                                                                                                               //either way, we set this link as visited
+                        if(currentNode->links.GetNumElements() == 1)            //if it does, then if that is the only link to this node, we're done with the chain, so append the chain to the list and return the newly created chain
+                        {
+//                              headsOfSystems.Append(newChain);
+                                return newChain;
+                        }
+                        nextNode = currentNode->links.GetHeadElement()->next->value;//follow the other link if there is one, so we go down the chain
+                        if(!setLinkVisited(currentNode, nextNode))                              //mark link as followed, so we know not to backtrack
+                        {
+        //                      headsOfSystems.Append(newChain);
+                                return newChain;                                                                //This condition, where no branches have occurred but both links have already
+                                                                                                                                //been taken can only occur in a loop with no spurs; add this loop to the
+                                                                                                                                //system (currently added as a chain for consistency), and return.
+                        }
+                }
+                currentNode = nextNode;                                                                 //set the current node to be the next node in the chain
+        }
+        currentNode->visited = true;
+        tmp = new int;
+        *tmp = currentNode->idNumber;
+        newChain->listOfNodes.Append(tmp);              //append the last node before branch (node shared jointly with branch chains)
+                                                                                                                                //re-mark as visited, just to make sure
+        ListElement<POEMSNode> * tempNode = currentNode->links.GetHeadElement();        //go through all of the links, one at a time that branch
+        POEMSChain * tempChain = NULL;                                                          //temporary variable to hold data
+        while(tempNode != NULL)                                                                         //when we have followed all links, stop
+        {
+                if(setLinkVisited(tempNode->value, currentNode))                //dont backtrack, or create closed loops
+                {
+                        tempChain = AddNewChain(tempNode->value);                       //Add a new chain created out of the next node down that link
+                        tempChain->parentChain = newChain;                                      //set the parent to be this chain
+                        newChain->childChains.Append(tempChain);                        //append the chain to this chain's list of child chains
+                }
+                tempNode = tempNode->next;                                                              //go to process the next chain
+        }
+        //headsOfSystems.Append(newChain);                                                      //append this chain to the system list
+        return newChain;
 }
 
 POEMSNode * SystemProcessor::findSingleLink(TreeNode * aNode)
 //This function takes the root of a search tree containing POEMSNodes and returns a POEMSNode corresponding to the start of a chain in the
 //system.  It finds a node that has not been visited before, and only has one link; this node will be used as the head of the chain.
 {
-	if(aNode == NULL)	
-	{
-		return NULL;
-	}
-	POEMSNode * returnVal =  (POEMSNode *)aNode->GetAuxData();	//get the poemsnode data out of the treenode
-	POEMSNode * detectLoneLoops = NULL;							//is used to handle a loop that has no protruding chains
-	if(returnVal->visited == false)
-	{
-		detectLoneLoops = returnVal;							//if we find any node that has not been visited yet, save it
-	}
-	if(returnVal->links.GetNumElements() == 1 && returnVal->visited == false)	//see if it has one element and hasnt been visited already
-	{
-		return returnVal;														//return the node is it meets this criteria
-	}
-	returnVal = findSingleLink(aNode->Left());									//otherwise, check the left subtree
-	if(returnVal == NULL)														//and if we find nothing...
-	{
-		returnVal = findSingleLink(aNode->Right());								//check the right subtree
-	}
-	if(returnVal == NULL)														//if we could not find any chains
-	{
-		returnVal = detectLoneLoops;											//see if we found any nodes at all that havent been processed
-	}
-	return returnVal;															//return what we find (will be NULL if no new chains are 
-																				//found)
+        if(aNode == NULL)
+        {
+                return NULL;
+        }
+        POEMSNode * returnVal =  (POEMSNode *)aNode->GetAuxData();      //get the poemsnode data out of the treenode
+        POEMSNode * detectLoneLoops = NULL;                                                     //is used to handle a loop that has no protruding chains
+        if(returnVal->visited == false)
+        {
+                detectLoneLoops = returnVal;                                                    //if we find any node that has not been visited yet, save it
+        }
+        if(returnVal->links.GetNumElements() == 1 && returnVal->visited == false)       //see if it has one element and hasnt been visited already
+        {
+                return returnVal;                                                                                                               //return the node is it meets this criteria
+        }
+        returnVal = findSingleLink(aNode->Left());                                                                      //otherwise, check the left subtree
+        if(returnVal == NULL)                                                                                                           //and if we find nothing...
+        {
+                returnVal = findSingleLink(aNode->Right());                                                             //check the right subtree
+        }
+        if(returnVal == NULL)                                                                                                           //if we could not find any chains
+        {
+                returnVal = detectLoneLoops;                                                                                    //see if we found any nodes at all that havent been processed
+        }
+        return returnVal;                                                                                                                       //return what we find (will be NULL if no new chains are
+                                                                                                                                                                //found)
 }
 
 bool SystemProcessor::setLinkVisited(POEMSNode * firstNode, POEMSNode * secondNode)
@@ -223,65 +225,66 @@ bool SystemProcessor::setLinkVisited(POEMSNode * firstNode, POEMSNode * secondNo
 //value for that particular link.  Because each link is represented twice, (once at each node in the link), both of the boolean values need
 //to be set in the event that the link has to be set as visited.
 {
-	//cout << "Checking link between nodes " << firstNode->idNumber << " and " << secondNode->idNumber << "... ";
-	ListElement<POEMSNode> * tmp = firstNode->links.GetHeadElement();	//get the head element of the list of pointers for node 1
-	ListElement<bool> * tmp2 = firstNode->taken.GetHeadElement();		//get the head element of the list of bool isVisited flags for node 1
-	while(tmp->value != NULL || tmp2->value != NULL)					//go through untill we reach the end of the lists
-	{
-		if(tmp->value == secondNode)							//if we find the link to the other node
-		{
-			if(*(tmp2->value) == true)							//if the link has already been visited
-			{
-				 //cout << "visited already" << endl;
-				return false;									//return false to indicate that the link has been visited before this attempt
-			}
-			else												//otherwise, visit it
-			{
-				*tmp2->value = true;
-			}
-			break;
-		}
-		tmp = tmp->next;										//go check next link
-		tmp2 = tmp2->next;
-	}
+        //std::cout << "Checking link between nodes " << firstNode->idNumber << " and " << secondNode->idNumber << "... ";
+        ListElement<POEMSNode> * tmp = firstNode->links.GetHeadElement();       //get the head element of the list of pointers for node 1
+        ListElement<bool> * tmp2 = firstNode->taken.GetHeadElement();           //get the head element of the list of bool isVisited flags for node 1
+        while(tmp->value != NULL || tmp2->value != NULL)                                        //go through untill we reach the end of the lists
+        {
+                if(tmp->value == secondNode)                                                    //if we find the link to the other node
+                {
+                        if(*(tmp2->value) == true)                                                      //if the link has already been visited
+                        {
+                                 //std::cout << "visited already" << std::endl;
+                                return false;                                                                   //return false to indicate that the link has been visited before this attempt
+                        }
+                        else                                                                                            //otherwise, visit it
+                        {
+                                *tmp2->value = true;
+                        }
+                        break;
+                }
+                tmp = tmp->next;                                                                                //go check next link
+                tmp2 = tmp2->next;
+        }
 
-	tmp = secondNode->links.GetHeadElement();			//now, if the link was unvisited, we need to go set the other node's list such that
-														//it also knows this link is being visited
-	tmp2 = secondNode->taken.GetHeadElement();
-	while(tmp->value != NULL || tmp2->value != NULL)	//go through the list
-	{
-		if(tmp->value == firstNode)						//if we find the link
-		{
-			if(*(tmp2->value) == true)					//and it has already been visited, then signal an error; this shouldnt ever happen
-			{
-				cout << "Error in parsing structure! Should never reach this condition! \n" << 
-						"Record of visited links out of synch between two adjacent nodes.\n";
-				return false;
-			}
-			else
-			{
-				*tmp2->value = true;					//set the appropriate value to true to indicate this link has been visited
-			}
-			break;
-		}
-		tmp = tmp->next;
-		tmp2 = tmp2->next;
-	}
-	//cout << "not visited" << endl;
-	return true;										//return true to indicate that this is the first time the link has been visited
+        tmp = secondNode->links.GetHeadElement();                       //now, if the link was unvisited, we need to go set the other node's list such that
+                                                                                                                //it also knows this link is being visited
+        tmp2 = secondNode->taken.GetHeadElement();
+        while(tmp->value != NULL || tmp2->value != NULL)        //go through the list
+        {
+                if(tmp->value == firstNode)                                             //if we find the link
+                {
+                        if(*(tmp2->value) == true)                                      //and it has already been visited, then signal an error; this shouldnt ever happen
+                        {
+                                std::cout << "Error in parsing structure! Should never reach this condition! \n" <<
+                                                "Record of visited links out of synch between two adjacent nodes.\n";
+                                return false;
+                        }
+                        else
+                        {
+                                *tmp2->value = true;                                    //set the appropriate value to true to indicate this link has been visited
+                        }
+                        break;
+                }
+                tmp = tmp->next;
+                tmp2 = tmp2->next;
+        }
+        //std::cout << "not visited" << std::endl;
+        return true;                                                                            //return true to indicate that this is the first time the link has been visited
 }
 
-List<POEMSChain> * SystemProcessor::getSystemData(void)	//Gets the list of POEMSChains that comprise the system.  Might eventually only
-														//return chains linked to the reference plane, but currently returns every chain
-														//in the system.
+List<POEMSChain> * SystemProcessor::getSystemData(void) //Gets the list of POEMSChains that comprise the system.  Might eventually only
+                                                                                                                //return chains linked to the reference plane, but currently returns every chain
+                                                                                                                //in the system.
 {
-	return &headsOfSystems;
+        return &headsOfSystems;
 }
 
 int SystemProcessor::getNumberOfHeadChains(void) //This function isnt implemented yet, and might be taken out entirely; this was a holdover
-												//from when I intended to return an array of chain pointers, rather than a list of chains
-												//It will probably be deleted once I finish figuring out exactly what needs to be returned
+                                                                                                //from when I intended to return an array of chain pointers, rather than a list of chains
+                                                                                                //It will probably be deleted once I finish figuring out exactly what needs to be returned
 {
-	return 0;
+        return 0;
+}
 }
 #endif
diff --git a/lib/poems/body.cpp b/lib/poems/body.cpp
index 487d5e1ab7..6825941b93 100644
--- a/lib/poems/body.cpp
+++ b/lib/poems/body.cpp
@@ -26,9 +26,8 @@
 #include "rigidbody.h"
 #include "vect3.h"
 
-class Joint;
-
 using namespace std;
+using namespace POEMS;
 
 Body::Body()
 {
@@ -132,7 +131,7 @@ void Body::AddPoint(Point* point){
 // global body functions
 //
 
-Body* NewBody(int type){
+Body* POEMS::NewBody(int type){
   switch( BodyType(type) )
     {
       case INERTIALFRAME :  // The inertial reference frame
diff --git a/lib/poems/body.h b/lib/poems/body.h
index 06a9c72417..e2394b40aa 100644
--- a/lib/poems/body.h
+++ b/lib/poems/body.h
@@ -25,6 +25,7 @@
 #include "mat3x3.h"
 #include "vect3.h"
 
+namespace POEMS {
 // emumerated type
 enum BodyType {
   INERTIALFRAME = 0,
@@ -75,5 +76,5 @@ public:
 
 // global body functions
 Body* NewBody(int type);
-
+}
 #endif
diff --git a/lib/poems/body23joint.cpp b/lib/poems/body23joint.cpp
index 7aa761fe2a..e3ef32b2a9 100644
--- a/lib/poems/body23joint.cpp
+++ b/lib/poems/body23joint.cpp
@@ -27,6 +27,9 @@
 #include "vect3.h"
 #include "virtualmatrix.h"
 
+using namespace std;
+using namespace POEMS;
+
 Body23Joint::Body23Joint(){
   DimQandU(4,2);
 }
diff --git a/lib/poems/body23joint.h b/lib/poems/body23joint.h
index c47cdfb9de..43a727c815 100644
--- a/lib/poems/body23joint.h
+++ b/lib/poems/body23joint.h
@@ -22,6 +22,7 @@
 #include "joint.h"
 #include "matrix.h"
 
+namespace POEMS {
 class Body23Joint : public Joint  {
   Matrix const_sP;
 public:
@@ -38,5 +39,5 @@ public:
   void ForwardKinematics();
   void BackwardKinematics();
 };
-
+}
 #endif
diff --git a/lib/poems/colmatmap.cpp b/lib/poems/colmatmap.cpp
index 1954e7ff15..34a95feb03 100644
--- a/lib/poems/colmatmap.cpp
+++ b/lib/poems/colmatmap.cpp
@@ -21,6 +21,8 @@
 #include <cstdlib>
 
 using namespace std;
+using namespace POEMS;
+
 
 ColMatMap::ColMatMap(){
 	numrows = 0;
diff --git a/lib/poems/colmatmap.h b/lib/poems/colmatmap.h
index 6d84870046..5a12cf487f 100644
--- a/lib/poems/colmatmap.h
+++ b/lib/poems/colmatmap.h
@@ -24,6 +24,7 @@
 #include "virtualcolmatrix.h"
 #include "virtualmatrix.h"
 
+namespace POEMS {
 class ColMatrix;
 
 class ColMatMap : public VirtualColMatrix  {
@@ -64,5 +65,5 @@ public:
         friend void FastCKRK5(ColMatMap& X, ColMatrix& Xi, ColMatrix* f, double* c, double dt);
         friend void FastFRK5(ColMatMap& X, ColMatrix& Xi, ColMatrix* f, double* c, double dt);
 };
-
+}
 #endif
diff --git a/lib/poems/colmatrix.cpp b/lib/poems/colmatrix.cpp
index e7ae98a368..0cc666645a 100644
--- a/lib/poems/colmatrix.cpp
+++ b/lib/poems/colmatrix.cpp
@@ -21,6 +21,8 @@
 #include <cstdlib>
 
 using namespace std;
+using namespace POEMS;
+
 
 ColMatrix::ColMatrix(){
 	numrows = 0;
diff --git a/lib/poems/colmatrix.h b/lib/poems/colmatrix.h
index 70f13356e2..bcfa86843a 100644
--- a/lib/poems/colmatrix.h
+++ b/lib/poems/colmatrix.h
@@ -23,12 +23,14 @@
 #include "virtualcolmatrix.h"
 #include "virtualmatrix.h"
 
+namespace POEMS {
 class Matrix;
 class Vect6;
 class Mat3x3;
 class Vect3;
 
 class ColMatrix : public VirtualColMatrix  {
+protected:
   double* elements;
 public:
   ColMatrix();
@@ -58,28 +60,28 @@ public:
   ColMatrix& operator=(const VirtualMatrix& A); // overloaded =
   ColMatrix& operator*=(double b);
 
-        void Abs();
-        void BasicMax(double& value, int& index);
-        void BasicMin(double& value, int& index);
+  void Abs();
+  void BasicMax(double& value, int& index);
+  void BasicMin(double& value, int& index);
 
   // fast matrix operations
-		  friend void FastQuaternions(ColMatrix& q, Mat3x3& C);
-		  friend void FastInvQuaternions(Mat3x3& C, ColMatrix& q);
-		  friend void FastQuaternionDerivatives(ColMatrix& q, ColMatrix& omega, ColMatrix& qdot);
-		  friend void FastTMult(Matrix& A, Vect6& B, ColMatrix& C);
-		  friend void FastMult(Matrix& A, ColMatrix& B, Vect6& C);
-		  friend void FastAssign(ColMatrix& A, ColMatrix& C);
+  friend void FastQuaternions(ColMatrix& q, Mat3x3& C);
+  friend void FastInvQuaternions(Mat3x3& C, ColMatrix& q);
+  friend void FastQuaternionDerivatives(ColMatrix& q, ColMatrix& omega, ColMatrix& qdot);
+  friend void FastTMult(Matrix& A, Vect6& B, ColMatrix& C);
+  friend void FastMult(Matrix& A, ColMatrix& B, Vect6& C);
+  friend void FastAssign(ColMatrix& A, ColMatrix& C);
 
-		  friend void FastMult(Mat3x3& A, ColMatrix& B, Vect3& C);
-		  friend void FastMult(Mat3x3& A, Vect3& B, ColMatrix& C);
-		  friend void FastAssign(ColMatrix&A, Vect3& C);
-  
-		  friend void EP_Derivatives(ColMatrix& q, ColMatrix& u, ColMatrix& qdot);
-		  friend void EP_Transformation(ColMatrix& q, Mat3x3& C);
-		  friend void EP_FromTransformation(ColMatrix& q, Mat3x3& C);
-		  friend void EP_Normalize(ColMatrix& q);
-		  friend void EPdotdot_udot(ColMatrix& Audot, ColMatrix& Aqdot, ColMatrix& Aq,ColMatrix& Aqddot);
-		  friend void qdot_to_u(ColMatrix& q, ColMatrix& u, ColMatrix& qdot);
+  friend void FastMult(Mat3x3& A, ColMatrix& B, Vect3& C);
+  friend void FastMult(Mat3x3& A, Vect3& B, ColMatrix& C);
+  friend void FastAssign(ColMatrix&A, Vect3& C);
+
+  friend void EP_Derivatives(ColMatrix& q, ColMatrix& u, ColMatrix& qdot);
+  friend void EP_Transformation(ColMatrix& q, Mat3x3& C);
+  friend void EP_FromTransformation(ColMatrix& q, Mat3x3& C);
+  friend void EP_Normalize(ColMatrix& q);
+  friend void EPdotdot_udot(ColMatrix& Audot, ColMatrix& Aqdot, ColMatrix& Aq,ColMatrix& Aqddot);
+  friend void qdot_to_u(ColMatrix& q, ColMatrix& u, ColMatrix& qdot);
 };
-
+}
 #endif
diff --git a/lib/poems/defines.h b/lib/poems/defines.h
index 455bc298d5..877bc0fa0b 100644
--- a/lib/poems/defines.h
+++ b/lib/poems/defines.h
@@ -15,13 +15,14 @@
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
-#ifndef _DEFINES_H_
-#define _DEFINES_H_
+#ifndef POEMS_DEFINES_H
+#define POEMS_DEFINES_H
 
+namespace POEMS {
 enum SolverType {
   ONSOLVER = 0,
   PARTICLESOLVER = 1
 };
-
+}
 #endif
 
diff --git a/lib/poems/eulerparameters.cpp b/lib/poems/eulerparameters.cpp
index f698c40b3a..bb8d676796 100644
--- a/lib/poems/eulerparameters.cpp
+++ b/lib/poems/eulerparameters.cpp
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: eulerparameters.cpp                                     *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,7 +11,7 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
@@ -23,8 +23,10 @@
 #include "mat3x3.h"
 
 using namespace std;
+using namespace POEMS;
 
-void EP_Derivatives(ColMatrix& q, ColMatrix& u, ColMatrix& qdot){  
+
+void POEMS::EP_Derivatives(ColMatrix& q, ColMatrix& u, ColMatrix& qdot){
   EP_Normalize(q);
   int num=u.GetNumRows();
   if (3<num){
@@ -32,17 +34,17 @@ void EP_Derivatives(ColMatrix& q, ColMatrix& u, ColMatrix& qdot){
 		  qdot.elements[i]=u.elements[i-1];
 	  }
   }
-  
-  qdot.elements[0] = 0.5 *(q.elements[3]*u.elements[0] - q.elements[2]*u.elements[1] + q.elements[1]*u.elements[2]);  
+
+  qdot.elements[0] = 0.5 *(q.elements[3]*u.elements[0] - q.elements[2]*u.elements[1] + q.elements[1]*u.elements[2]);
   qdot.elements[1] = 0.5 *(q.elements[2]*u.elements[0] + q.elements[3]*u.elements[1] - q.elements[0]*u.elements[2]);
   qdot.elements[2] = 0.5 *(-q.elements[1]*u.elements[0] + q.elements[0]*u.elements[1] + q.elements[3]*u.elements[2]);
   qdot.elements[3] = -0.5 *(q.elements[0]*u.elements[0] + q.elements[1]*u.elements[1] + q.elements[2]*u.elements[2]);
-    
+
   }
 
-void EP_Transformation(ColMatrix& q, Mat3x3& C){
+void POEMS::EP_Transformation(ColMatrix& q, Mat3x3& C){
   EP_Normalize(q);
-  
+
   double q11 = q.elements[0]*q.elements[0];
   double q22 = q.elements[1]*q.elements[1];
   double q33 = q.elements[2]*q.elements[2];
@@ -69,7 +71,7 @@ void EP_Transformation(ColMatrix& q, Mat3x3& C){
   C.elements[2][1] = 2*(q23 + q14);
 }
 
-void EP_FromTransformation(ColMatrix& q, Mat3x3& C){
+void POEMS::EP_FromTransformation(ColMatrix& q, Mat3x3& C){
   double b[4];
 
   // condition indicators
@@ -117,8 +119,8 @@ void EP_FromTransformation(ColMatrix& q, Mat3x3& C){
   EP_Normalize(q);
 }
 
-void EP_Normalize(ColMatrix& q){
-  double one = 1.0/sqrt(q.elements[0]*q.elements[0] + q.elements[1]*q.elements[1] + q.elements[2]*q.elements[2] + q.elements[3]*q.elements[3]);  
+void POEMS:: EP_Normalize(ColMatrix& q){
+  double one = 1.0/sqrt(q.elements[0]*q.elements[0] + q.elements[1]*q.elements[1] + q.elements[2]*q.elements[2] + q.elements[3]*q.elements[3]);
   q.elements[0] = one*q.elements[0];
   q.elements[1] = one*q.elements[1];
   q.elements[2] = one*q.elements[2];
@@ -128,29 +130,29 @@ void EP_Normalize(ColMatrix& q){
 
 
 
-void EPdotdot_udot(ColMatrix& Audot, ColMatrix& Aqdot, ColMatrix& Aq, ColMatrix& Aqddot){  
+void POEMS::EPdotdot_udot(ColMatrix& Audot, ColMatrix& Aqdot, ColMatrix& Aq, ColMatrix& Aqddot){
   int num=Audot.GetNumRows();
   if (3<num){
 	  for (int i=4; i<=num; i++){
 		  Aqddot.elements[i]=Audot.elements[i-1];
 	  }
   }
-    
+
   double AA;
   AA=Aqdot.elements[0]*Aqdot.elements[0]+Aqdot.elements[1]*Aqdot.elements[1]+Aqdot.elements[2]*Aqdot.elements[2]+Aqdot.elements[3]*Aqdot.elements[3];
-  
+
   Aqddot.elements[0] = 0.5 *(Aq.elements[3]*Audot.elements[0] - Aq.elements[2]*Audot.elements[1] + Aq.elements[1]*Audot.elements[2]-2*Aq.elements[0]*AA);
-  
+
   Aqddot.elements[1] = 0.5 *(Aq.elements[2]*Audot.elements[0] + Aq.elements[3]*Audot.elements[1] - Aq.elements[0]*Audot.elements[2]-2*Aq.elements[1]*AA);
-  
+
   Aqddot.elements[2] = 0.5 *(-Aq.elements[1]*Audot.elements[0] + Aq.elements[0]*Audot.elements[1] + Aq.elements[3]*Audot.elements[2]-2*Aq.elements[2]*AA);
-  
+
   Aqddot.elements[3] = -0.5 *(Aq.elements[0]*Audot.elements[0] + Aq.elements[1]*Audot.elements[1] + Aq.elements[2]*Audot.elements[2]+2*Aq.elements[3]*AA);
-  
+
 }
 
 
-void qdot_to_u(ColMatrix& q, ColMatrix& u, ColMatrix& qdot){
+void POEMS::qdot_to_u(ColMatrix& q, ColMatrix& u, ColMatrix& qdot){
 	EP_Normalize(q);
 	int num=qdot.GetNumRows();
       	if (4<num){
@@ -161,7 +163,7 @@ void qdot_to_u(ColMatrix& q, ColMatrix& u, ColMatrix& qdot){
       	u.elements[0]=2*(q.elements[3]*qdot.elements[0]+q.elements[2]*qdot.elements[1]-q.elements[1]*qdot.elements[2]-q.elements[0]*qdot.elements[3]);
       	u.elements[1]=2*(-q.elements[2]*qdot.elements[0]+q.elements[3]*qdot.elements[1]+q.elements[0]*qdot.elements[2]-q.elements[1]*qdot.elements[3]);
       	u.elements[2]=2*(q.elements[1]*qdot.elements[0]-q.elements[0]*qdot.elements[1]+q.elements[3]*qdot.elements[2]-q.elements[2]*qdot.elements[3]);
-  
+
 }
 
 
diff --git a/lib/poems/eulerparameters.h b/lib/poems/eulerparameters.h
index 2e69a85fb9..d331b47872 100644
--- a/lib/poems/eulerparameters.h
+++ b/lib/poems/eulerparameters.h
@@ -18,6 +18,7 @@
 #ifndef EULERPARAMETERS_H
 #define EULERPARAMETERS_H
 
+namespace POEMS {
 class ColMatrix;
 class Mat3x3;
 
@@ -32,6 +33,6 @@ void EP_Normalize(ColMatrix& q);
 void EPdotdot_udot(ColMatrix& Audot, ColMatrix& Aqdot, ColMatrix& Aq,ColMatrix& Aqddot);
 
 void qdot_to_u(ColMatrix& q, ColMatrix& u, ColMatrix& qdot);
-
+}
 #endif
 
diff --git a/lib/poems/fastmatrixops.cpp b/lib/poems/fastmatrixops.cpp
index e6c81f0659..db0dc3236d 100644
--- a/lib/poems/fastmatrixops.cpp
+++ b/lib/poems/fastmatrixops.cpp
@@ -29,12 +29,14 @@
 #include "vect6.h"
 
 using namespace std;
+using namespace POEMS;
+
 
 //
 // Cross Product (friend of Vect3)
 //
 
-void FastCross(Vect3& a, Vect3& b, Vect3& c){
+void POEMS::FastCross(Vect3& a, Vect3& b, Vect3& c){
   c.elements[0] = a.elements[1]*b.elements[2] - a.elements[2]*b.elements[1];
   c.elements[1] = a.elements[2]*b.elements[0] - a.elements[0]*b.elements[2];
   c.elements[2] = a.elements[0]*b.elements[1] - a.elements[1]*b.elements[0];
@@ -44,7 +46,7 @@ void FastCross(Vect3& a, Vect3& b, Vect3& c){
 //  Simple Rotation (friend of Vect3 and Mat3x3)
 //
 
-void FastSimpleRotation(Vect3& v, double q, Mat3x3& C){
+void POEMS::FastSimpleRotation(Vect3& v, double q, Mat3x3& C){
   // intermediate quantities
   double cq = cos(q);
   double sq = sin(q);
@@ -69,7 +71,7 @@ void FastSimpleRotation(Vect3& v, double q, Mat3x3& C){
 // Quaternion Functions
 //
 
-void FastQuaternions(ColMatrix& q, Mat3x3& C){
+void POEMS::FastQuaternions(ColMatrix& q, Mat3x3& C){
   double* e = q.elements;
 
   // normalize the quaternions
@@ -94,7 +96,7 @@ void FastQuaternions(ColMatrix& q, Mat3x3& C){
   C.elements[2][1] = 2 * (e[2]*e[3] + e[0]*e[1]);
 }
 
-void FastQuaternionDerivatives(ColMatrix& q, ColMatrix& omega, ColMatrix& qdot){
+void POEMS::FastQuaternionDerivatives(ColMatrix& q, ColMatrix& omega, ColMatrix& qdot){
   double* w = omega.elements;
   double* e = q.elements;
 
@@ -104,7 +106,7 @@ void FastQuaternionDerivatives(ColMatrix& q, ColMatrix& omega, ColMatrix& qdot){
   qdot.elements[3] = 0.5 * ( w[2]*e[0] + w[1]*e[1] - w[0]*e[2]);
 }
 
-void FastInvQuaternions(Mat3x3& C, ColMatrix& q){
+void POEMS::FastInvQuaternions(Mat3x3& C, ColMatrix& q){
 }
 
 //
@@ -121,7 +123,7 @@ void FastInvQuaternions(Mat3x3& C, ColMatrix& q){
 //
 
 // friend of Matrix
-void FastLDLT(Matrix& A, Matrix& C){ // C is the LD of the LDL^T decomposition of A (SPD)
+void POEMS::FastLDLT(Matrix& A, Matrix& C){ // C is the LD of the LDL^T decomposition of A (SPD)
   double Lv;
   int n = A.numrows;
 
@@ -143,7 +145,7 @@ void FastLDLT(Matrix& A, Matrix& C){ // C is the LD of the LDL^T decomposition o
 
 
 // friend of Mat6x6
-void FastLDLT(Mat6x6& A, Mat6x6& C){ // C is the LD of the LDL^T decomposition of A (SPD)
+void POEMS::FastLDLT(Mat6x6& A, Mat6x6& C){ // C is the LD of the LDL^T decomposition of A (SPD)
   double v[6];
   double Lv;
 
@@ -165,7 +167,7 @@ void FastLDLT(Mat6x6& A, Mat6x6& C){ // C is the LD of the LDL^T decomposition o
 }
 
 // friend of Matrix
-void FastLDLTSubs(Matrix& LD, Matrix& B, Matrix& C){
+void POEMS::FastLDLTSubs(Matrix& LD, Matrix& B, Matrix& C){
   int n = B.numrows;
   int c = B.numcols;
   double temp;
@@ -190,7 +192,7 @@ void FastLDLTSubs(Matrix& LD, Matrix& B, Matrix& C){
 }
 
 // friend of Matrix
-void FastLDLTSubsLH(Matrix& B, Matrix& LD, Matrix& C){
+void POEMS::FastLDLTSubsLH(Matrix& B, Matrix& LD, Matrix& C){
   int n = B.numcols;
   int c = B.numrows;
   double temp;
@@ -215,7 +217,7 @@ void FastLDLTSubsLH(Matrix& B, Matrix& LD, Matrix& C){
 }
 
 // friend of Mat6x6
-void FastLDLTSubs(Mat6x6& LD, Mat6x6& B, Mat6x6& C){
+void POEMS::FastLDLTSubs(Mat6x6& LD, Mat6x6& B, Mat6x6& C){
   double temp;
 
   for(int k=0;k<6;k++){
@@ -238,7 +240,7 @@ void FastLDLTSubs(Mat6x6& LD, Mat6x6& B, Mat6x6& C){
 }
 
 // friend of Mat6x6 & Vect6
-void FastLDLTSubs(Mat6x6& LD, Vect6& B, Vect6& C){
+void POEMS::FastLDLTSubs(Mat6x6& LD, Vect6& B, Vect6& C){
   double temp;
 
   for(int i=0;i<6;i++){
@@ -259,7 +261,7 @@ void FastLDLTSubs(Mat6x6& LD, Vect6& B, Vect6& C){
 }
 
 // friend of Matrix
-void FastLU(Matrix& A, Matrix& LU, int *indx){ // LU is the LU decomposition of A
+void POEMS::FastLU(Matrix& A, Matrix& LU, int *indx){ // LU is the LU decomposition of A
   int i,imax=0,j,k;
   int n = A.numrows;
   double big, dum, sum, temp;
@@ -309,7 +311,7 @@ void FastLU(Matrix& A, Matrix& LU, int *indx){ // LU is the LU decomposition of
 }
 
 // friend of Mat3x3
-void FastLU(Mat3x3& A, Mat3x3& LU, int *indx){ // LU is the LU decomposition of A
+void POEMS::FastLU(Mat3x3& A, Mat3x3& LU, int *indx){ // LU is the LU decomposition of A
   int i,imax=0,j,k;
   double big, dum, sum, temp;
   double vv[10000];
@@ -357,7 +359,7 @@ void FastLU(Mat3x3& A, Mat3x3& LU, int *indx){ // LU is the LU decomposition of
 }
 
 // friend of Mat4x4
-void FastLU(Mat4x4& A, Mat4x4& LU, int *indx){ // LU is the LU decomposition of A
+void POEMS::FastLU(Mat4x4& A, Mat4x4& LU, int *indx){ // LU is the LU decomposition of A
   int i,imax=0,j,k;
   double big, dum, sum, temp;
   double vv[10000];
@@ -405,7 +407,7 @@ void FastLU(Mat4x4& A, Mat4x4& LU, int *indx){ // LU is the LU decomposition of
 }
 
 // friend of Mat6x6
-void FastLU(Mat6x6& A, Mat6x6& LU, int *indx){ // LU is the LU decomposition of A
+void POEMS::FastLU(Mat6x6& A, Mat6x6& LU, int *indx){ // LU is the LU decomposition of A
   int i,imax=0,j,k;
   double big, dum, sum, temp;
   double vv[10000];
@@ -453,7 +455,7 @@ void FastLU(Mat6x6& A, Mat6x6& LU, int *indx){ // LU is the LU decomposition of
 }
 
 // friend of Matrix
-void FastLUSubs(Matrix& LU, Matrix& B, Matrix& C, int *indx){ // Appropriate Forward and Back Substitution	
+void POEMS::FastLUSubs(Matrix& LU, Matrix& B, Matrix& C, int *indx){ // Appropriate Forward and Back Substitution	
   int i,ip,j,k;
   int n = B.numrows;
   int c = B.numcols;
@@ -477,7 +479,7 @@ void FastLUSubs(Matrix& LU, Matrix& B, Matrix& C, int *indx){ // Appropriate For
 }
 
 // friend of Matrix and Mat3x3
-void FastLUSubs(Mat3x3& LU, Matrix& B, Matrix& C, int *indx){ // Appropriate Forward and Back Substitution	
+void POEMS::FastLUSubs(Mat3x3& LU, Matrix& B, Matrix& C, int *indx){ // Appropriate Forward and Back Substitution	
   int i,ip,j,k;
   int n = B.numrows;
   int c = B.numcols;
@@ -501,7 +503,7 @@ void FastLUSubs(Mat3x3& LU, Matrix& B, Matrix& C, int *indx){ // Appropriate For
 }
 
 // friend of Matrix and Mat4x4
-void FastLUSubs(Mat4x4& LU, Matrix& B, Matrix& C, int *indx){ // Appropriate Forward and Back Substitution
+void POEMS::FastLUSubs(Mat4x4& LU, Matrix& B, Matrix& C, int *indx){ // Appropriate Forward and Back Substitution
   int i,ip,j,k;
   int n = B.numrows;
   int c = B.numcols;
@@ -525,7 +527,7 @@ void FastLUSubs(Mat4x4& LU, Matrix& B, Matrix& C, int *indx){ // Appropriate For
 }
 
 // friend of Matrix and Mat6x6
-void FastLUSubs(Mat6x6& LU, Matrix& B, Matrix& C, int *indx){ // Appropriate Forward and Back Substitution	
+void POEMS::FastLUSubs(Mat6x6& LU, Matrix& B, Matrix& C, int *indx){ // Appropriate Forward and Back Substitution	
   int i,ip,j,k;
   int n = B.numrows;
   int c = B.numcols;
@@ -551,7 +553,7 @@ void FastLUSubs(Mat6x6& LU, Matrix& B, Matrix& C, int *indx){ // Appropriate For
 
 // The following LUSubsLH routine is incomplete at the moment.
 // friend of Matrix
-void FastLUSubsLH(Matrix& LU, Matrix& B, Matrix& C, int *indx){ // Appropriate Forward and Back Substitution	
+void POEMS::FastLUSubsLH(Matrix& LU, Matrix& B, Matrix& C, int *indx){ // Appropriate Forward and Back Substitution	
   int i,ip,j,k;
   int n = B.numcols;
   int c = B.numrows;
@@ -582,13 +584,13 @@ void FastLUSubsLH(Matrix& LU, Matrix& B, Matrix& C, int *indx){ // Appropriate F
 // Triple sum
 //
 
-void FastTripleSum(Vect3& a, Vect3& b, Vect3& c, Vect3& d){ // d = a+b+c
+void POEMS::FastTripleSum(Vect3& a, Vect3& b, Vect3& c, Vect3& d){ // d = a+b+c
   d.elements[0] = a.elements[0]+b.elements[0]+c.elements[0];
   d.elements[1] = a.elements[1]+b.elements[1]+c.elements[1];
   d.elements[2] = a.elements[2]+b.elements[2]+c.elements[2];
 }
 
-void FastTripleSumPPM(Vect3& a, Vect3& b, Vect3& c, Vect3& d){ // d = a+b-c
+void POEMS::FastTripleSumPPM(Vect3& a, Vect3& b, Vect3& c, Vect3& d){ // d = a+b-c
   d.elements[0] = a.elements[0]+b.elements[0]-c.elements[0];
   d.elements[1] = a.elements[1]+b.elements[1]-c.elements[1];
   d.elements[2] = a.elements[2]+b.elements[2]-c.elements[2];
@@ -599,7 +601,7 @@ void FastTripleSumPPM(Vect3& a, Vect3& b, Vect3& c, Vect3& d){ // d = a+b-c
 //
 
 // friend of matrix
-void FastMult(Matrix& A, Matrix& B, Matrix& C){  // C = A*B
+void POEMS::FastMult(Matrix& A, Matrix& B, Matrix& C){  // C = A*B
   // assumes dimensions are already correct!
   int r = A.numrows;
   int ca = A.numcols;
@@ -615,7 +617,7 @@ void FastMult(Matrix& A, Matrix& B, Matrix& C){  // C = A*B
 }
 
 // friend of matrix
-void FastTMult(Matrix& A, Matrix& B, Matrix& C){  // C = A*B
+void POEMS::FastTMult(Matrix& A, Matrix& B, Matrix& C){  // C = A*B
   // assumes dimensions are already correct!
   int r = A.numcols;
   int ca = A.numrows;
@@ -631,14 +633,14 @@ void FastTMult(Matrix& A, Matrix& B, Matrix& C){  // C = A*B
 }
 
 // friend of Mat3x3 & Vect3
-void FastMult(Mat3x3& A, Vect3& B, Vect3& C){  // C = A*B
+void POEMS::FastMult(Mat3x3& A, Vect3& B, Vect3& C){  // C = A*B
   C.elements[0] = A.elements[0][0]*B.elements[0] + A.elements[0][1]*B.elements[1] + A.elements[0][2]*B.elements[2];
   C.elements[1] = A.elements[1][0]*B.elements[0] + A.elements[1][1]*B.elements[1] + A.elements[1][2]*B.elements[2];
   C.elements[2] = A.elements[2][0]*B.elements[0] + A.elements[2][1]*B.elements[1] + A.elements[2][2]*B.elements[2];
 }
 
 // friend of Mat3x3, ColMatrix, & Vect3
-void FastMult(Mat3x3& A, ColMatrix& B, Vect3& C){  // C = A*B
+void POEMS::FastMult(Mat3x3& A, ColMatrix& B, Vect3& C){  // C = A*B
   C.elements[0] = A.elements[0][0]*B.elements[0] + A.elements[0][1]*B.elements[1] + A.elements[0][2]*B.elements[2];
   C.elements[1] = A.elements[1][0]*B.elements[0] + A.elements[1][1]*B.elements[1] + A.elements[1][2]*B.elements[2];
   C.elements[2] = A.elements[2][0]*B.elements[0] + A.elements[2][1]*B.elements[1] + A.elements[2][2]*B.elements[2];
@@ -646,42 +648,42 @@ void FastMult(Mat3x3& A, ColMatrix& B, Vect3& C){  // C = A*B
 
 
 // friend of Mat3x3, ColMatrix, & Vect3
-void FastMult(Mat3x3& A, Vect3& B, ColMatrix& C){  // C = A*B
+void POEMS::FastMult(Mat3x3& A, Vect3& B, ColMatrix& C){  // C = A*B
   C.elements[0] = A.elements[0][0]*B.elements[0] + A.elements[0][1]*B.elements[1] + A.elements[0][2]*B.elements[2];
   C.elements[1] = A.elements[1][0]*B.elements[0] + A.elements[1][1]*B.elements[1] + A.elements[1][2]*B.elements[2];
   C.elements[2] = A.elements[2][0]*B.elements[0] + A.elements[2][1]*B.elements[1] + A.elements[2][2]*B.elements[2];
 }
 
 // friend of Mat3x3 & Vect3
-void FastTMult(Mat3x3& A, Vect3& B, Vect3& C){  // C = A^T*B
+void POEMS::FastTMult(Mat3x3& A, Vect3& B, Vect3& C){  // C = A^T*B
   C.elements[0] = A.elements[0][0]*B.elements[0] + A.elements[1][0]*B.elements[1] + A.elements[2][0]*B.elements[2];
   C.elements[1] = A.elements[0][1]*B.elements[0] + A.elements[1][1]*B.elements[1] + A.elements[2][1]*B.elements[2];
   C.elements[2] = A.elements[0][2]*B.elements[0] + A.elements[1][2]*B.elements[1] + A.elements[2][2]*B.elements[2];
 }
 
 // friend of Mat3x3 & Vect3
-void FastNegMult(Mat3x3& A, Vect3& B, Vect3& C){  // C = -A*B
+void POEMS::FastNegMult(Mat3x3& A, Vect3& B, Vect3& C){  // C = -A*B
   C.elements[0] = -A.elements[0][0]*B.elements[0] - A.elements[0][1]*B.elements[1] - A.elements[0][2]*B.elements[2];
   C.elements[1] = -A.elements[1][0]*B.elements[0] - A.elements[1][1]*B.elements[1] - A.elements[1][2]*B.elements[2];
   C.elements[2] = -A.elements[2][0]*B.elements[0] - A.elements[2][1]*B.elements[1] - A.elements[2][2]*B.elements[2];
 }
 
 // friend of Mat3x3 & Vect3
-void FastNegTMult(Mat3x3& A, Vect3& B, Vect3& C){  // C = -A^T*B
+void POEMS::FastNegTMult(Mat3x3& A, Vect3& B, Vect3& C){  // C = -A^T*B
   C.elements[0] = -A.elements[0][0]*B.elements[0] - A.elements[1][0]*B.elements[1] - A.elements[2][0]*B.elements[2];
   C.elements[1] = -A.elements[0][1]*B.elements[0] - A.elements[1][1]*B.elements[1] - A.elements[2][1]*B.elements[2];
   C.elements[2] = -A.elements[0][2]*B.elements[0] - A.elements[1][2]*B.elements[1] - A.elements[2][2]*B.elements[2];
 }
 
 // friend of Vect3
-void FastMult(double a, Vect3& B, Vect3& C){  // C = a*B
+void POEMS::FastMult(double a, Vect3& B, Vect3& C){  // C = a*B
   C.elements[0] = a*B.elements[0];
   C.elements[1] = a*B.elements[1];
   C.elements[2] = a*B.elements[2];
 }
 
 // friend of Mat4x4 & Vect4
-void FastMult(Mat4x4& A, Vect4& B, Vect4& C){  // C = A*B
+void POEMS::FastMult(Mat4x4& A, Vect4& B, Vect4& C){  // C = A*B
   C.elements[0] = A.elements[0][0]*B.elements[0] + A.elements[0][1]*B.elements[1] + A.elements[0][2]*B.elements[2] + A.elements[0][3]*B.elements[3];
   C.elements[1] = A.elements[1][0]*B.elements[0] + A.elements[1][1]*B.elements[1] + A.elements[1][2]*B.elements[2] + A.elements[1][3]*B.elements[3];
   C.elements[2] = A.elements[2][0]*B.elements[0] + A.elements[2][1]*B.elements[1] + A.elements[2][2]*B.elements[2] + A.elements[2][3]*B.elements[3];
@@ -689,7 +691,7 @@ void FastMult(Mat4x4& A, Vect4& B, Vect4& C){  // C = A*B
 }
 
 // friend of Mat4x4 & Vect4
-void FastTMult(Mat4x4& A, Vect4& B, Vect4& C){  // C = A^T*B
+void POEMS::FastTMult(Mat4x4& A, Vect4& B, Vect4& C){  // C = A^T*B
   C.elements[0] = A.elements[0][0]*B.elements[0] + A.elements[1][0]*B.elements[1] + A.elements[2][0]*B.elements[2] + A.elements[3][0]*B.elements[3];
   C.elements[1] = A.elements[0][1]*B.elements[0] + A.elements[1][1]*B.elements[1] + A.elements[2][1]*B.elements[2] + A.elements[3][1]*B.elements[3];
   C.elements[2] = A.elements[0][2]*B.elements[0] + A.elements[1][2]*B.elements[1] + A.elements[2][2]*B.elements[2] + A.elements[3][2]*B.elements[3];
@@ -697,7 +699,7 @@ void FastTMult(Mat4x4& A, Vect4& B, Vect4& C){  // C = A^T*B
 }
 
 // friend of Mat4x4 & Vect4
-void FastNegMult(Mat4x4& A, Vect4& B, Vect4& C){  // C = -A*B
+void POEMS::FastNegMult(Mat4x4& A, Vect4& B, Vect4& C){  // C = -A*B
   C.elements[0] = -A.elements[0][0]*B.elements[0] - A.elements[0][1]*B.elements[1] - A.elements[0][2]*B.elements[2] - A.elements[0][3]*B.elements[3];
   C.elements[1] = -A.elements[1][0]*B.elements[0] - A.elements[1][1]*B.elements[1] - A.elements[1][2]*B.elements[2] - A.elements[1][3]*B.elements[3];
   C.elements[2] = -A.elements[2][0]*B.elements[0] - A.elements[2][1]*B.elements[1] - A.elements[2][2]*B.elements[2] - A.elements[2][3]*B.elements[3];
@@ -705,7 +707,7 @@ void FastNegMult(Mat4x4& A, Vect4& B, Vect4& C){  // C = -A*B
 }
 
 // friend of Mat4x4 & Vect4
-void FastNegTMult(Mat4x4& A, Vect4& B, Vect4& C){  // C = -A^T*B
+void POEMS::FastNegTMult(Mat4x4& A, Vect4& B, Vect4& C){  // C = -A^T*B
   C.elements[0] = -A.elements[0][0]*B.elements[0] - A.elements[1][0]*B.elements[1] - A.elements[2][0]*B.elements[2] - A.elements[3][0]*B.elements[3];
   C.elements[1] = -A.elements[0][1]*B.elements[0] - A.elements[1][1]*B.elements[1] - A.elements[2][1]*B.elements[2] - A.elements[3][1]*B.elements[3];
   C.elements[2] = -A.elements[0][2]*B.elements[0] - A.elements[1][2]*B.elements[1] - A.elements[2][2]*B.elements[2] - A.elements[3][2]*B.elements[3];
@@ -713,7 +715,7 @@ void FastNegTMult(Mat4x4& A, Vect4& B, Vect4& C){  // C = -A^T*B
 }
 
 // friend of Vect4
-void FastMult(double a, Vect4& B, Vect4& C){  // C = a*B
+void POEMS::FastMult(double a, Vect4& B, Vect4& C){  // C = a*B
   C.elements[0] = a*B.elements[0];
   C.elements[1] = a*B.elements[1];
   C.elements[2] = a*B.elements[2];
@@ -721,7 +723,7 @@ void FastMult(double a, Vect4& B, Vect4& C){  // C = a*B
 }
 
 // friend of Matrix & Mat6x6
-void FastMultT(Matrix& A, Matrix& B, Mat6x6& C){  // C = A*B^T
+void POEMS::FastMultT(Matrix& A, Matrix& B, Mat6x6& C){  // C = A*B^T
   int i,j,k,n;
   n = A.numcols;
 
@@ -734,7 +736,7 @@ void FastMultT(Matrix& A, Matrix& B, Mat6x6& C){  // C = A*B^T
 }
 
 // friend Matrix, Vect6, ColMatrix
-void FastMult(Matrix& A, ColMatrix& B, Vect6& C){
+void POEMS::FastMult(Matrix& A, ColMatrix& B, Vect6& C){
   int ca = A.numcols;
 
   int i,k;
@@ -746,7 +748,7 @@ void FastMult(Matrix& A, ColMatrix& B, Vect6& C){
 }
 
 // friend of Matrix & Mat6x6
-void FastMult(Mat6x6& A, Matrix& B, Matrix& C){  // C = A*B
+void POEMS::FastMult(Mat6x6& A, Matrix& B, Matrix& C){  // C = A*B
   // assumes dimensions are already correct!
   int cb = B.numcols;
 
@@ -760,7 +762,7 @@ void FastMult(Mat6x6& A, Matrix& B, Matrix& C){  // C = A*B
 }
 
 // friend Matrix, Vect6, ColMatrix
-void FastTMult(Matrix& A, Vect6& B, ColMatrix& C){  // C = A^T*B
+void POEMS::FastTMult(Matrix& A, Vect6& B, ColMatrix& C){  // C = A^T*B
   int n = C.numrows;
   int i,k;
   for(i=0;i<n;i++){
@@ -771,7 +773,7 @@ void FastTMult(Matrix& A, Vect6& B, ColMatrix& C){  // C = A^T*B
 }
 
 // friend of Mat3x3
-void FastMult(Mat3x3& A, Mat3x3& B, Mat3x3& C){  // C = A*B
+void POEMS::FastMult(Mat3x3& A, Mat3x3& B, Mat3x3& C){  // C = A*B
   C.elements[0][0] = A.elements[0][0]*B.elements[0][0] + A.elements[0][1]*B.elements[1][0] + A.elements[0][2]*B.elements[2][0];
   C.elements[0][1] = A.elements[0][0]*B.elements[0][1] + A.elements[0][1]*B.elements[1][1] + A.elements[0][2]*B.elements[2][1];
   C.elements[0][2] = A.elements[0][0]*B.elements[0][2] + A.elements[0][1]*B.elements[1][2] + A.elements[0][2]*B.elements[2][2];
@@ -786,7 +788,7 @@ void FastMult(Mat3x3& A, Mat3x3& B, Mat3x3& C){  // C = A*B
 }
 
 // friend of Mat3x3
-void FastMultT(Mat3x3& A, Mat3x3& B, Mat3x3& C){  // C = A*B^T
+void POEMS::FastMultT(Mat3x3& A, Mat3x3& B, Mat3x3& C){  // C = A*B^T
   C.elements[0][0] = A.elements[0][0]*B.elements[0][0] + A.elements[0][1]*B.elements[0][1] + A.elements[0][2]*B.elements[0][2];
   C.elements[0][1] = A.elements[0][0]*B.elements[1][0] + A.elements[0][1]*B.elements[1][1] + A.elements[0][2]*B.elements[1][2];
   C.elements[0][2] = A.elements[0][0]*B.elements[2][0] + A.elements[0][1]*B.elements[2][1] + A.elements[0][2]*B.elements[2][2];
@@ -801,7 +803,7 @@ void FastMultT(Mat3x3& A, Mat3x3& B, Mat3x3& C){  // C = A*B^T
 }
 
 // friend of Mat4x4
-void FastMult(Mat4x4& A, Mat4x4& B, Mat4x4& C){  // C = A*B
+void POEMS::FastMult(Mat4x4& A, Mat4x4& B, Mat4x4& C){  // C = A*B
   C.elements[0][0] = A.elements[0][0]*B.elements[0][0] + A.elements[0][1]*B.elements[1][0] + A.elements[0][2]*B.elements[2][0] + A.elements[0][3]*B.elements[3][0];
   C.elements[0][1] = A.elements[0][0]*B.elements[0][1] + A.elements[0][1]*B.elements[1][1] + A.elements[0][2]*B.elements[2][1] + A.elements[0][3]*B.elements[3][1];
   C.elements[0][2] = A.elements[0][0]*B.elements[0][2] + A.elements[0][1]*B.elements[1][2] + A.elements[0][2]*B.elements[2][2] + A.elements[0][3]*B.elements[3][2];
@@ -824,7 +826,7 @@ void FastMult(Mat4x4& A, Mat4x4& B, Mat4x4& C){  // C = A*B
 }
 
 // friend of Mat4x4
-void FastMultT(Mat4x4& A, Mat4x4& B, Mat4x4& C){  // C = A*B^T
+void POEMS::FastMultT(Mat4x4& A, Mat4x4& B, Mat4x4& C){  // C = A*B^T
   C.elements[0][0] = A.elements[0][0]*B.elements[0][0] + A.elements[0][1]*B.elements[0][1] + A.elements[0][2]*B.elements[0][2] + A.elements[0][3]*B.elements[0][3];
   C.elements[0][1] = A.elements[0][0]*B.elements[1][0] + A.elements[0][1]*B.elements[1][1] + A.elements[0][2]*B.elements[1][2] + A.elements[0][3]*B.elements[1][3];
   C.elements[0][2] = A.elements[0][0]*B.elements[2][0] + A.elements[0][1]*B.elements[2][1] + A.elements[0][2]*B.elements[2][2] + A.elements[0][3]*B.elements[2][3];
@@ -848,7 +850,7 @@ void FastMultT(Mat4x4& A, Mat4x4& B, Mat4x4& C){  // C = A*B^T
 }
 
 // friend of Mat6x6
-void FastMult(Mat6x6& A, Mat6x6& B, Mat6x6& C){  // C = A*B
+void POEMS::FastMult(Mat6x6& A, Mat6x6& B, Mat6x6& C){  // C = A*B
   int i,j,k;
   for(i=0;i<6;i++)
     for(j=0;j<6;j++){
@@ -859,7 +861,7 @@ void FastMult(Mat6x6& A, Mat6x6& B, Mat6x6& C){  // C = A*B
 }
 
 // friend of Mat6x6
-void FastMultT(Mat6x6& A, Mat6x6& B, Mat6x6& C){  // C = A*B
+void POEMS::FastMultT(Mat6x6& A, Mat6x6& B, Mat6x6& C){  // C = A*B
   int i,j,k;
   for(i=0;i<6;i++)
     for(j=0;j<6;j++){
@@ -870,7 +872,7 @@ void FastMultT(Mat6x6& A, Mat6x6& B, Mat6x6& C){  // C = A*B
 }
 
 // friend of Mat6x6
-void FastTMult(Mat6x6& A, Mat6x6& B, Mat6x6& C){  // C = A^T*B
+void POEMS::FastTMult(Mat6x6& A, Mat6x6& B, Mat6x6& C){  // C = A^T*B
   int i,j,k;
   for(i=0;i<6;i++)
     for(j=0;j<6;j++){
@@ -881,13 +883,13 @@ void FastTMult(Mat6x6& A, Mat6x6& B, Mat6x6& C){  // C = A^T*B
 }
 
 // friend of Mat6x6 & Vect6
-void FastMult(Mat6x6& A, Vect6& B, Vect6& C){  // C = A*B
+void POEMS::FastMult(Mat6x6& A, Vect6& B, Vect6& C){  // C = A*B
   for(int i=0;i<6;i++)
     C.elements[i] = A.elements[i][0]*B.elements[0] + A.elements[i][1]*B.elements[1] + A.elements[i][2]*B.elements[2] + A.elements[i][3]*B.elements[3] + A.elements[i][4]*B.elements[4] + A.elements[i][5]*B.elements[5];
 }
 
 // friend of Mat6x6 & Vect6
-void FastTMult(Mat6x6& A, Vect6& B, Vect6& C){  // C = A^T*B
+void POEMS::FastTMult(Mat6x6& A, Vect6& B, Vect6& C){  // C = A^T*B
   for(int i=0;i<6;i++)
     C.elements[i] = A.elements[0][i]*B.elements[0] + A.elements[1][i]*B.elements[1] + A.elements[2][i]*B.elements[2] + A.elements[3][i]*B.elements[3] + A.elements[4][i]*B.elements[4] + A.elements[5][i]*B.elements[5];
 }
@@ -897,21 +899,21 @@ void FastTMult(Mat6x6& A, Vect6& B, Vect6& C){  // C = A^T*B
 //
 
 // friend of Vect3
-void FastAdd(Vect3& A, Vect3& B, Vect3& C){ // C = A+B
+void POEMS::FastAdd(Vect3& A, Vect3& B, Vect3& C){ // C = A+B
   C.elements[0] = A.elements[0] + B.elements[0];
   C.elements[1] = A.elements[1] + B.elements[1];
   C.elements[2] = A.elements[2] + B.elements[2];
 }
 
 // friend of Vect4
-void FastAdd(Vect4& A, Vect4& B, Vect4& C){ // C = A+B
+void POEMS::FastAdd(Vect4& A, Vect4& B, Vect4& C){ // C = A+B
   C.elements[0] = A.elements[0] + B.elements[0];
   C.elements[1] = A.elements[1] + B.elements[1];
   C.elements[2] = A.elements[2] + B.elements[2];
   C.elements[3] = A.elements[3] + B.elements[3];
 }
 
-void FastAdd(Mat6x6& A, Mat6x6& B, Mat6x6& C){  // C = A+B
+void POEMS::FastAdd(Mat6x6& A, Mat6x6& B, Mat6x6& C){  // C = A+B
   int i,j;
   for(i=0;i<6;i++)
     for(j=0;j<6;j++)
@@ -919,7 +921,7 @@ void FastAdd(Mat6x6& A, Mat6x6& B, Mat6x6& C){  // C = A+B
 }
 
 // friend of Vect6
-void FastAdd(Vect6& A, Vect6& B, Vect6& C){ // C = A-B
+void POEMS::FastAdd(Vect6& A, Vect6& B, Vect6& C){ // C = A-B
   C.elements[0] = A.elements[0] + B.elements[0];
   C.elements[1] = A.elements[1] + B.elements[1];
   C.elements[2] = A.elements[2] + B.elements[2];
@@ -933,21 +935,21 @@ void FastAdd(Vect6& A, Vect6& B, Vect6& C){ // C = A-B
 //
 
 // friend of Vect3
-void FastSubt(Vect3& A, Vect3& B, Vect3& C){ // C = A-B
+void POEMS::FastSubt(Vect3& A, Vect3& B, Vect3& C){ // C = A-B
   C.elements[0] = A.elements[0] - B.elements[0];
   C.elements[1] = A.elements[1] - B.elements[1];
   C.elements[2] = A.elements[2] - B.elements[2];
 }
 
 // friend of Vect4
-void FastSubt(Vect4& A, Vect4& B, Vect4& C){ // C = A-B
+void POEMS::FastSubt(Vect4& A, Vect4& B, Vect4& C){ // C = A-B
   C.elements[0] = A.elements[0] - B.elements[0];
   C.elements[1] = A.elements[1] - B.elements[1];
   C.elements[2] = A.elements[2] - B.elements[2];
   C.elements[3] = A.elements[3] - B.elements[3];
 }
 
-void FastSubt(Mat6x6& A, Mat6x6& B, Mat6x6& C){  // C = A-B
+void POEMS::FastSubt(Mat6x6& A, Mat6x6& B, Mat6x6& C){  // C = A-B
   int i,j;
   for(i=0;i<6;i++)
     for(j=0;j<6;j++)
@@ -955,7 +957,7 @@ void FastSubt(Mat6x6& A, Mat6x6& B, Mat6x6& C){  // C = A-B
 }
 
 // friend of Vect6
-void FastSubt(Vect6& A, Vect6& B, Vect6& C){ // C = A-B
+void POEMS::FastSubt(Vect6& A, Vect6& B, Vect6& C){ // C = A-B
   C.elements[0] = A.elements[0] - B.elements[0];
   C.elements[1] = A.elements[1] - B.elements[1];
   C.elements[2] = A.elements[2] - B.elements[2];
@@ -965,33 +967,33 @@ void FastSubt(Vect6& A, Vect6& B, Vect6& C){ // C = A-B
 }
 
 // friend of ColMatMap
-void FastAssign(ColMatMap& A, ColMatMap& C){
+void POEMS::FastAssign(ColMatMap& A, ColMatMap& C){
   for(int i=0;i<C.numrows;i++)
     *(C.elements[i]) = *(A.elements[i]);
 }
 
 // friend of ColMatrix
-void FastAssign(ColMatrix& A, ColMatrix& C){ //C = A
+void POEMS::FastAssign(ColMatrix& A, ColMatrix& C){ //C = A
   for(int i=0;i<C.numrows;i++)
     C.elements[i] = A.elements[i];
 }
 
 // friend of Vect3
-void FastAssign(Vect3& A, Vect3& C){ //C = A
+void POEMS::FastAssign(Vect3& A, Vect3& C){ //C = A
     C.elements[0] = A.elements[0];
     C.elements[1] = A.elements[1];
     C.elements[2] = A.elements[2];
 }
 
 // friend of Vect3 & ColMatrix
-void FastAssign(ColMatrix& A, Vect3& C){ //C = A
+void POEMS::FastAssign(ColMatrix& A, Vect3& C){ //C = A
     C.elements[0] = A.elements[0];
     C.elements[1] = A.elements[1];
     C.elements[2] = A.elements[2];
 }
 
 // friend of Vect4
-void FastAssign(Vect4& A, Vect4& C){ //C = A
+void POEMS::FastAssign(Vect4& A, Vect4& C){ //C = A
     C.elements[0] = A.elements[0];
     C.elements[1] = A.elements[1];
     C.elements[2] = A.elements[2];
@@ -999,7 +1001,7 @@ void FastAssign(Vect4& A, Vect4& C){ //C = A
 }
 
 // friend of Mat3x3
-void FastAssignT(Mat3x3& A, Mat3x3& C){
+void POEMS::FastAssignT(Mat3x3& A, Mat3x3& C){
   C.elements[0][0] = A.elements[0][0];
   C.elements[1][1] = A.elements[1][1];
   C.elements[2][2] = A.elements[2][2];
@@ -1015,7 +1017,7 @@ void FastAssignT(Mat3x3& A, Mat3x3& C){
 }
 
 // friend of Mat4x4
-void FastAssignT(Mat4x4& A, Mat4x4& C){
+void POEMS::FastAssignT(Mat4x4& A, Mat4x4& C){
   C.elements[0][0] = A.elements[0][0];
   C.elements[1][1] = A.elements[1][1];
   C.elements[2][2] = A.elements[2][2];
diff --git a/lib/poems/fastmatrixops.h b/lib/poems/fastmatrixops.h
index b2784a6e53..4c3c1353bd 100644
--- a/lib/poems/fastmatrixops.h
+++ b/lib/poems/fastmatrixops.h
@@ -18,6 +18,7 @@
 #ifndef FASTMATRIXOPS_H
 #define FASTMATRIXOPS_H
 
+namespace POEMS {
 class ColMatMap;
 class ColMatrix;
 class Mat3x3;
@@ -102,5 +103,5 @@ void FastAssign(ColMatrix&A, Vect3& C);
 void FastAssign(Vect4& A, Vect4& C);          // C = A
 void FastAssignT(Mat3x3& A, Mat3x3& C);       // C = A^T
 void FastAssignT(Mat4x4& A, Mat4x4& C);       // C = A^T
-
+}
 #endif
diff --git a/lib/poems/fixedpoint.cpp b/lib/poems/fixedpoint.cpp
index 238ddfed70..c28b160c18 100644
--- a/lib/poems/fixedpoint.cpp
+++ b/lib/poems/fixedpoint.cpp
@@ -24,6 +24,8 @@
 #include "virtualmatrix.h"
 
 using namespace std;
+using namespace POEMS;
+
 
 FixedPoint::FixedPoint(){
 }
diff --git a/lib/poems/fixedpoint.h b/lib/poems/fixedpoint.h
index fb8f23b0ae..5da1d6d69c 100644
--- a/lib/poems/fixedpoint.h
+++ b/lib/poems/fixedpoint.h
@@ -24,6 +24,7 @@
 #include "point.h"
 #include "vect3.h"
 
+namespace POEMS {
 class FixedPoint : public Point  {
 public:
   FixedPoint();
@@ -35,5 +36,5 @@ public:
   bool ReadInPointData(std::istream& in);
   void WriteOutPointData(std::ostream& out);
 };
-
+}
 #endif
diff --git a/lib/poems/freebodyjoint.cpp b/lib/poems/freebodyjoint.cpp
index e6c59340f3..0022ebd9ff 100644
--- a/lib/poems/freebodyjoint.cpp
+++ b/lib/poems/freebodyjoint.cpp
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: freebodyjoint.cpp                                       *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,10 +11,10 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
- 
+
 #include "freebodyjoint.h"
 
 #include <iostream>
@@ -28,10 +28,12 @@
 #include "mat6x6.h"
 #include "vect3.h"
 #include "virtualmatrix.h"
-   
+
+using namespace POEMS;
+
 
 FreeBodyJoint::FreeBodyJoint(){
-  DimQandU(7,6);  
+  DimQandU(7,6);
 }
 
 FreeBodyJoint::~FreeBodyJoint(){
@@ -49,112 +51,112 @@ void FreeBodyJoint::WriteOutJointData(std::ostream& out){
 }
 
 void FreeBodyJoint::ComputeLocalTransform(){
-  Mat3x3 ko_C_k;  
-  EP_Transformation(q, ko_C_k);   
-  FastMult(pk_C_ko,ko_C_k,pk_C_k);  
+  Mat3x3 ko_C_k;
+  EP_Transformation(q, ko_C_k);
+  FastMult(pk_C_ko,ko_C_k,pk_C_k);
 }
 
 Matrix FreeBodyJoint::GetForward_sP(){
   Mat6x6 sP;
-  //sP.Identity();      
-  
-  sP.Zeros();  
+  //sP.Identity();
+
+  sP.Zeros();
   Mat3x3 temp0=T(pk_C_k);
   for(int i=1;i<4;i++){
 	  sP(i,i)=1.0;
 	  for(int j=1;j<4;j++){
 		  sP(3+i,3+j)=temp0(i,j);
 		  }
-	  } 	  	  
-  return sP;    
+	  } 	  	
+  return sP;
 }
 
 Matrix FreeBodyJoint::GetBackward_sP(){
   Mat6x6 sP;
-  sP.Identity();        
+  sP.Identity();
   sP =-1.0*sP;
-  cout<<"Did I come here in "<<endl;
-  return sP;    
+  std::cout<<"Did I come here in "<<std::endl;
+  return sP;
 }
 
 
 void FreeBodyJoint::UpdateForward_sP( Matrix& sP){
   // do nothing
 }
-  
+
 void FreeBodyJoint::UpdateBackward_sP( Matrix& sP){
   // do nothing
-}  
-  
+}
+
 void FreeBodyJoint::ForwardKinematics(){
- //cout<<"Check in freebody "<<q<<" "<<qdot<<endl;
+ //std::cout<<"Check in freebody "<<q<<" "<<qdot<<std::endl;
   EP_Normalize(q);
 
   // COMMENT STEP1: CALCULATE ORIENTATIONS
-  ComputeForwardTransforms();  
-  
-    
+  ComputeForwardTransforms();
+
+
   //COMMENT STEP2: CALCULATE POSITION VECTORS
-  Vect3 result1, result2, result3, result4; 
-  
+  Vect3 result1, result2, result3, result4;
+
   result1.BasicSet(0,q.BasicGet(4));
   result1.BasicSet(1,q.BasicGet(5));
   result1.BasicSet(2,q.BasicGet(6));
-  
-  FastAssign(result1,r12);    
+
+  FastAssign(result1,r12);
   FastNegMult(k_C_pk,r12,r21);
-  
-  FastAssign(r12,body2->r);  
-    
-  //COMMENT STEP3: CALCULATE QDOT  
+
+  FastAssign(r12,body2->r);
+
+  //COMMENT STEP3: CALCULATE QDOT
   qdot_to_u(q, u, qdot);
-  
-  
+
+
   Vect3 WN;
   WN.BasicSet(0,u.BasicGet(0));
   WN.BasicSet(1,u.BasicGet(1));
   WN.BasicSet(2,u.BasicGet(2));
-    
-  Vect3 VN; 
+
+  Vect3 VN;
   VN.BasicSet(0,u.BasicGet(3));
   VN.BasicSet(1,u.BasicGet(4));
-  VN.BasicSet(2,u.BasicGet(5));  
-    
-  FastAssign(WN,body2->omega_k);  
- 
-  Vect3 pk_w_k;
-  FastMult(body2->n_C_k,WN,pk_w_k);      
-  FastAssign(pk_w_k,body2->omega);
-  
+  VN.BasicSet(2,u.BasicGet(5));
+
+  FastAssign(WN,body2->omega_k);
+
+  Vect3 pk_w_k;
+  FastMult(body2->n_C_k,WN,pk_w_k);
+  FastAssign(pk_w_k,body2->omega);
+
+
+
+  //COMMENT STEP5: CALCULATE VELOCITES
+  FastAssign(VN,body2->v);
+  FastTMult(body2->n_C_k,body2->v,body2->v_k);
+
+
+  //CALCULATE KE
 
-  
-  //COMMENT STEP5: CALCULATE VELOCITES        
-  FastAssign(VN,body2->v);    
-  FastTMult(body2->n_C_k,body2->v,body2->v_k);        
-  
-  
-  //CALCULATE KE          
- 
   Matrix tempke;
-  tempke = T(body2->v)*(body2->v);    
+  tempke = T(body2->v)*(body2->v);
   double ke = 0.0;
-  ke = body2->mass*tempke(1,1);  
+  ke = body2->mass*tempke(1,1);
   FastMult(body2->inertia,body2->omega_k,result1);
-  tempke= T(body2->omega_k)*result1;  
+  tempke= T(body2->omega_k)*result1;
   ke = 0.5*ke + 0.5*tempke(1,1);
   body2->KE=ke;
- 
-  
+
+
   //COMMENT STEP6: CALCULATE STATE EXPLICIT ANGULAR ACCELERATIONS
   body2->alpha_t.Zeros();
-    
-    
+
+
   //COMMENT STEP7: CALCULATE STATE EXPLICIT ACCELERATIONS
-  body2->a_t.Zeros();   
-  
+  body2->a_t.Zeros();
+
   }
 
 void FreeBodyJoint::BackwardKinematics(){
-cout<<"Did I come here "<<endl;
-         
+std::cout<<"Did I come here "<<std::endl;
+
 }
diff --git a/lib/poems/freebodyjoint.h b/lib/poems/freebodyjoint.h
index 10574847e1..8b4bec2cbb 100644
--- a/lib/poems/freebodyjoint.h
+++ b/lib/poems/freebodyjoint.h
@@ -23,12 +23,12 @@
 #include "joint.h"
 #include "matrix.h"
 
-
+namespace POEMS {
 class FreeBodyJoint : public Joint{
-public: 
+public:
   FreeBodyJoint();
   ~FreeBodyJoint();
-  
+
   JointType GetType();
   bool ReadInJointData(std::istream& in);
   void WriteOutJointData(std::ostream& out);
@@ -40,5 +40,5 @@ public:
   void ForwardKinematics();
   void BackwardKinematics();
 };
-
+}
 #endif
diff --git a/lib/poems/inertialframe.cpp b/lib/poems/inertialframe.cpp
index 98ec820e7f..19ca7cdb0b 100644
--- a/lib/poems/inertialframe.cpp
+++ b/lib/poems/inertialframe.cpp
@@ -21,6 +21,8 @@
 #include "virtualmatrix.h"
 
 using namespace std;
+using namespace POEMS;
+
 
 InertialFrame::InertialFrame(){
   gravity.Zeros();
diff --git a/lib/poems/inertialframe.h b/lib/poems/inertialframe.h
index 3620641537..ae987cdc56 100644
--- a/lib/poems/inertialframe.h
+++ b/lib/poems/inertialframe.h
@@ -24,7 +24,7 @@
 #include "body.h"
 #include "vect3.h"
 
-
+namespace POEMS {
 class InertialFrame : public Body  {
   Vect3 gravity;
 public:
@@ -36,5 +36,5 @@ public:
   bool ReadInBodyData(std::istream& in);
   void WriteOutBodyData(std::ostream& out);
 };
-
+}
 #endif
diff --git a/lib/poems/joint.cpp b/lib/poems/joint.cpp
index ff24d5ebf5..3e356337cb 100644
--- a/lib/poems/joint.cpp
+++ b/lib/poems/joint.cpp
@@ -36,6 +36,8 @@
 #include "virtualmatrix.h"
 
 using namespace std;
+using namespace POEMS;
+
 
 Joint::Joint(){
   body1 = body2 = 0;
@@ -245,7 +247,7 @@ void Joint::ComputeBackwardGlobalTransform(){
 // global joint functions
 //
 
-Joint* NewJoint(int type){
+Joint* POEMS::NewJoint(int type){
   switch( JointType(type) )
   {
     case FREEBODYJOINT : return new FreeBodyJoint;
diff --git a/lib/poems/joint.h b/lib/poems/joint.h
index 958805c4b5..fa9b6826a7 100644
--- a/lib/poems/joint.h
+++ b/lib/poems/joint.h
@@ -25,6 +25,7 @@
 #include "matrix.h"
 #include "vect3.h"
 
+namespace POEMS {
 class VirtualMatrix;
 
 enum JointType {
@@ -124,5 +125,5 @@ public:
 
 // global joint functions
 Joint* NewJoint(int type);
-
+}
 #endif
diff --git a/lib/poems/mat3x3.cpp b/lib/poems/mat3x3.cpp
index 3d9f1adcca..8a5098ff59 100644
--- a/lib/poems/mat3x3.cpp
+++ b/lib/poems/mat3x3.cpp
@@ -21,6 +21,8 @@
 #include "mat3x3.h"
 
 using namespace std;
+using namespace POEMS;
+
 
 Mat3x3::Mat3x3(){
   numrows = numcols = 3;
diff --git a/lib/poems/mat3x3.h b/lib/poems/mat3x3.h
index 18abeadddc..bc38e6d393 100644
--- a/lib/poems/mat3x3.h
+++ b/lib/poems/mat3x3.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: mat3x3.h                                                *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,7 +11,7 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
@@ -21,6 +21,7 @@
 #include <iostream>
 #include "virtualmatrix.h"
 
+namespace POEMS {
 class Vect3;
 class Mat6x6;
 class ColMatrix;
@@ -67,16 +68,15 @@ public:
   friend void FastMult(Mat3x3& A, Mat3x3& B, Mat3x3& C);
   friend void FastMultT(Mat3x3& A, Mat3x3& B, Mat3x3& C);
   friend void FastAssignT(Mat3x3& A, Mat3x3& C);
-  friend void FastMult(Mat3x3& A, Vect3& B, ColMatrix& C);  
-  
+  friend void FastMult(Mat3x3& A, Vect3& B, ColMatrix& C);
+
   friend void OnPopulateSC(Vect3& gamma, Mat3x3& C, Mat6x6& SC);
   friend void OnPopulateSI(Mat3x3& inertia, double mass, Mat6x6& sI);
 
   friend void FastMult(Mat3x3& A, ColMatrix& B, Vect3& C);
-  
+
   friend void EP_Transformation(ColMatrix& q, Mat3x3& C);
   friend void EP_FromTransformation(ColMatrix& q, Mat3x3& C);
-  
 };
-
+}
 #endif
diff --git a/lib/poems/mat4x4.cpp b/lib/poems/mat4x4.cpp
index a6a587bc03..9a0e7b8d56 100644
--- a/lib/poems/mat4x4.cpp
+++ b/lib/poems/mat4x4.cpp
@@ -21,6 +21,8 @@
 #include <iostream>
 
 using namespace std;
+using namespace POEMS;
+
 
 Mat4x4::Mat4x4(){
   numrows = numcols = 4;
diff --git a/lib/poems/mat4x4.h b/lib/poems/mat4x4.h
index a6ae232473..3f6d97ba15 100644
--- a/lib/poems/mat4x4.h
+++ b/lib/poems/mat4x4.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: mat4x4.h                                                *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,7 +11,7 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
@@ -21,17 +21,18 @@
 #include <iostream>
 #include "virtualmatrix.h"
 
+namespace POEMS {
 class Vect4;
 class Matrix;
 
 class Mat4x4 : public VirtualMatrix  {
   double elements[4][4];
-public: 
+public:
   Mat4x4();
   ~Mat4x4();
   Mat4x4(const Mat4x4& A);  // copy constructor
   Mat4x4(const VirtualMatrix& A);  // copy constructor
-  
+
   double& operator_2int (int i, int j); // array access
   double Get_2int(int i, int j) const;
   void Set_2int(int i, int j, double value);
@@ -64,5 +65,5 @@ public:
   friend void FastMultT(Mat4x4& A, Mat4x4& B, Mat4x4& C);
   friend void FastAssignT(Mat4x4& A, Mat4x4& C);
 };
-
+}
 #endif
diff --git a/lib/poems/mat6x6.cpp b/lib/poems/mat6x6.cpp
index fdc917e03d..e9dde2fe87 100644
--- a/lib/poems/mat6x6.cpp
+++ b/lib/poems/mat6x6.cpp
@@ -21,6 +21,8 @@
 #include <iostream>
 
 using namespace std;
+using namespace POEMS;
+
 
 Mat6x6::Mat6x6(){
   numrows = numcols = 6;
diff --git a/lib/poems/mat6x6.h b/lib/poems/mat6x6.h
index 5eaa4f803b..a611bc8283 100644
--- a/lib/poems/mat6x6.h
+++ b/lib/poems/mat6x6.h
@@ -21,6 +21,7 @@
 
 #include "virtualmatrix.h"
 
+namespace POEMS {
 class Matrix;
 class Mat3x3;
 class Vect6;
@@ -73,5 +74,5 @@ public:
   friend void OnPopulateSC(Vect3& gamma, Mat3x3& C, Mat6x6& SC);
   friend void OnPopulateSI(Mat3x3& inertia, double mass, Mat6x6& sI);
 };
-
+}
 #endif
diff --git a/lib/poems/matrices.h b/lib/poems/matrices.h
index 7a7e2c4b58..15f0d3ad23 100644
--- a/lib/poems/matrices.h
+++ b/lib/poems/matrices.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: matrices.h                                              *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,7 +11,7 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
@@ -20,23 +20,14 @@
 #define MATRICES_H
 
 #include "matrix.h"
-
 #include "colmatrix.h"
-
 #include "rowmatrix.h"
-
 #include "mat3x3.h"
-
 #include "vect3.h"
-
 #include "mat4x4.h"
-
 #include "vect4.h"
-
 #include "mat6x6.h"
-
 #include "vect6.h"
-
 #include "colmatmap.h"
 
 #endif
diff --git a/lib/poems/matrix.cpp b/lib/poems/matrix.cpp
index f05d6f815c..f7b0c631f2 100644
--- a/lib/poems/matrix.cpp
+++ b/lib/poems/matrix.cpp
@@ -21,6 +21,8 @@
 #include <cstdlib>
 
 using namespace std;
+using namespace POEMS;
+
 
 Matrix::Matrix(){
   numrows = numcols = 0;
diff --git a/lib/poems/matrix.h b/lib/poems/matrix.h
index 8bd8110760..c2ed92873c 100644
--- a/lib/poems/matrix.h
+++ b/lib/poems/matrix.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: matrix.h                                                *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,7 +11,7 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
@@ -22,6 +22,7 @@
 
 #include "virtualmatrix.h"
 
+namespace POEMS {
 class Mat3x3;
 class Mat4x4;
 class Mat6x6;
@@ -73,7 +74,6 @@ public:
   friend void FastMult(Mat6x6& A, Matrix& B, Matrix& C);
   friend void FastMult(Matrix& A, ColMatrix& B, Vect6& C);
   friend void FastMultT(Matrix& A, Matrix& B, Mat6x6& C);
-  
 };
-
+}
 #endif
diff --git a/lib/poems/matrixfun.cpp b/lib/poems/matrixfun.cpp
index 99fbae4ca1..037e986a99 100644
--- a/lib/poems/matrixfun.cpp
+++ b/lib/poems/matrixfun.cpp
@@ -27,12 +27,14 @@
 #include "virtualrowmatrix.h"
 
 using namespace std;
+using namespace POEMS;
+
 
 //
 //  Create a new matrix
 //
 
-VirtualMatrix* NewMatrix(int type){
+VirtualMatrix* POEMS::NewMatrix(int type){
   switch( MatrixType(type) )
     {
       case MATRIX : return new Matrix;
@@ -50,7 +52,7 @@ VirtualMatrix* NewMatrix(int type){
 // Transpose
 //
 
-Matrix T(const VirtualMatrix& A){
+Matrix POEMS::T(const VirtualMatrix& A){
 	int numrows = A.GetNumRows();
 	int numcols = A.GetNumCols();
 	Matrix C(numcols,numrows);
@@ -60,8 +62,8 @@ Matrix T(const VirtualMatrix& A){
 	return C;
 }
 
-Mat3x3 T(const Mat3x3& A){
-	Mat3x3 C;
+Mat3x3 POEMS::T(const Mat3x3& A){
+        Mat3x3 C;
 	C.elements[0][0] = A.elements[0][0];
 	C.elements[1][1] = A.elements[1][1];
 	C.elements[2][2] = A.elements[2][2];
@@ -76,7 +78,7 @@ Mat3x3 T(const Mat3x3& A){
 	return C;
 }
 
-Mat6x6 T(const Mat6x6& A){
+Mat6x6 POEMS::T(const Mat6x6& A){
 	Mat6x6 C;
 	int i,j;
 	for(i=0;i<6;i++)
@@ -86,7 +88,7 @@ Mat6x6 T(const Mat6x6& A){
 	return C;
 }
 
-Matrix T(const Vect3& A){
+Matrix POEMS::T(const Vect3& A){
 	Matrix C(1,3);
 	C.BasicSet(0,0,A.elements[0]);
 	C.BasicSet(0,1,A.elements[1]);
@@ -95,7 +97,7 @@ Matrix T(const Vect3& A){
 	return C;
 }
 
-Matrix T(const Vect6& A){
+Matrix POEMS::T(const Vect6& A){
 	Matrix C(1,6);
 	C.BasicSet(0,0,A.elements[0]);
 	C.BasicSet(0,1,A.elements[1]);
@@ -107,7 +109,7 @@ Matrix T(const Vect6& A){
 	return C;
 }
 
-RowMatrix T(const VirtualColMatrix &A){
+RowMatrix POEMS::T(const VirtualColMatrix &A){
 	int numele = A.GetNumRows();
 	RowMatrix C(numele);
 	for(int i=0;i<numele;i++)
@@ -115,7 +117,7 @@ RowMatrix T(const VirtualColMatrix &A){
 	return C;
 }
 
-ColMatrix T(const VirtualRowMatrix &A){
+ColMatrix POEMS::T(const VirtualRowMatrix &A){
 	int numele = A.GetNumCols();
 	ColMatrix C(numele);
 	for(int i=0;i<numele;i++)
@@ -127,7 +129,7 @@ ColMatrix T(const VirtualRowMatrix &A){
 // Symmetric Inverse
 //
 
-Matrix SymInverse(Matrix &A){
+Matrix POEMS::SymInverse(Matrix &A){
 	int r = A.GetNumRows();
 	Matrix C(r,r);
 	Matrix LD(r,r);
@@ -139,7 +141,7 @@ Matrix SymInverse(Matrix &A){
 	return C;
 }
 
-Mat6x6 SymInverse(Mat6x6 &A){
+Mat6x6 POEMS::SymInverse(Mat6x6 &A){
 	Mat6x6 C;
 	Mat6x6 LD;
 	Mat6x6 I;
@@ -207,7 +209,7 @@ Mat6x6 Inverse(Mat6x6& A){
 // overloaded addition
 //
 
-Matrix operator+ (const VirtualMatrix &A, const VirtualMatrix &B){      // addition
+Matrix POEMS::operator+ (const VirtualMatrix &A, const VirtualMatrix &B){      // addition
 	int Arows,Acols,Brows,Bcols;
 	Arows = A.GetNumRows();
 	Acols = A.GetNumCols();
@@ -232,7 +234,7 @@ Matrix operator+ (const VirtualMatrix &A, const VirtualMatrix &B){      // addit
 // overloaded subtraction
 //
 
-Matrix operator- (const VirtualMatrix &A, const VirtualMatrix &B){      // subtraction
+Matrix POEMS::operator- (const VirtualMatrix &A, const VirtualMatrix &B){      // subtraction
 	int Arows,Acols,Brows,Bcols;
 	Arows = A.GetNumRows();
 	Acols = A.GetNumCols();
@@ -257,7 +259,7 @@ Matrix operator- (const VirtualMatrix &A, const VirtualMatrix &B){      // subtr
 // overloaded matrix multiplication
 //
 
-Matrix operator* (const VirtualMatrix &A, const VirtualMatrix &B){      // multiplication
+Matrix POEMS::operator* (const VirtualMatrix &A, const VirtualMatrix &B){      // multiplication
 	int Arows,Acols,Brows,Bcols;
 	Arows = A.GetNumRows();
 	Acols = A.GetNumCols();
@@ -284,13 +286,13 @@ Matrix operator* (const VirtualMatrix &A, const VirtualMatrix &B){      // multi
 // overloaded scalar multiplication
 //
 
-Matrix operator* (const VirtualMatrix &A, double b){    // multiplication
+Matrix POEMS::operator* (const VirtualMatrix &A, double b){    // multiplication
 	Matrix C = A;
 	C *= b;
 	return C;
 }
 
-Matrix operator* (double b, const VirtualMatrix &A){    // multiplication
+Matrix POEMS::operator* (double b, const VirtualMatrix &A){    // multiplication
 	Matrix C = A;  
 	C *= b;
 	return C;
@@ -300,7 +302,7 @@ Matrix operator* (double b, const VirtualMatrix &A){    // multiplication
 // overloaded negative
 //
 
-Matrix operator- (const VirtualMatrix& A){  // negative
+Matrix POEMS::operator- (const VirtualMatrix& A){  // negative
 	int r = A.GetNumRows();
 	int c = A.GetNumCols();
 	Matrix C(r,c);
@@ -316,7 +318,7 @@ Matrix operator- (const VirtualMatrix& A){  // negative
 // Cross product (friend of Vect3)
 //
 
-Vect3 Cross(Vect3& a, Vect3& b){
+Vect3 POEMS::Cross(Vect3& a, Vect3& b){
 	return CrossMat(a)*b;
 }
 
@@ -324,7 +326,7 @@ Vect3 Cross(Vect3& a, Vect3& b){
 // Cross Matrix (friend of Vect3 & Mat3x3)
 //
 
-Mat3x3 CrossMat(Vect3& a){
+Mat3x3 POEMS::CrossMat(Vect3& a){
 	Mat3x3 C;
 	C.Zeros();
 	C.elements[0][1] = -a.elements[2];
@@ -341,7 +343,7 @@ Mat3x3 CrossMat(Vect3& a){
 // Stack
 //
 
-Matrix Stack(VirtualMatrix& A, VirtualMatrix& B){
+Matrix POEMS::Stack(VirtualMatrix& A, VirtualMatrix& B){
 	int m,na,nb;
 	m = A.GetNumCols();
 	if( m != B.GetNumCols()){
@@ -370,7 +372,7 @@ Matrix Stack(VirtualMatrix& A, VirtualMatrix& B){
 
 //Hstack
 
-Matrix HStack(VirtualMatrix& A, VirtualMatrix& B){
+Matrix POEMS::HStack(VirtualMatrix& A, VirtualMatrix& B){
 	int m,na,nb;
 	m = A.GetNumRows();
 	if( m != B.GetNumRows()){
@@ -400,13 +402,13 @@ Matrix HStack(VirtualMatrix& A, VirtualMatrix& B){
 //
 //
 
-void Set6DAngularVector(Vect6& v6, Vect3& v3){
+void POEMS::Set6DAngularVector(Vect6& v6, Vect3& v3){
 	v6.elements[0] = v3.elements[0];
 	v6.elements[1] = v3.elements[1];
 	v6.elements[2] = v3.elements[2];
 }
 
-void Set6DLinearVector(Vect6& v6, Vect3& v3){
+void POEMS::Set6DLinearVector(Vect6& v6, Vect3& v3){
 	v6.elements[3] = v3.elements[0];
 	v6.elements[4] = v3.elements[1];
 	v6.elements[5] = v3.elements[2];
diff --git a/lib/poems/matrixfun.h b/lib/poems/matrixfun.h
index e9721eca7b..a34882140c 100644
--- a/lib/poems/matrixfun.h
+++ b/lib/poems/matrixfun.h
@@ -27,6 +27,7 @@
 #include "vect3.h"
 #include "vect6.h"
 
+namespace POEMS {
 class VirtualColMatrix;
 class VirtualMatrix;
 class VirtualRowMatrix;
@@ -55,22 +56,22 @@ Mat4x4 Inverse(Mat4x4& A);
 Mat6x6 Inverse(Mat6x6& A);
 
 // overloaded addition
-Matrix operator+ (const VirtualMatrix &A, const VirtualMatrix &B);	// addition
-//Mat3x3 operator+ (const Mat3x3 &A, const Mat3x3 &B);	// addition
-//Matrix operator+ (const VirtualMatrix &A, const VirtualMatrix &B);	// addition
+Matrix operator+ (const VirtualMatrix &A, const VirtualMatrix &B);      // addition
+//Mat3x3 operator+ (const Mat3x3 &A, const Mat3x3 &B);  // addition
+//Matrix operator+ (const VirtualMatrix &A, const VirtualMatrix &B);    // addition
 
 // overloaded subtraction
-Matrix operator- (const VirtualMatrix &A, const VirtualMatrix &B);	// subtraction
+Matrix operator- (const VirtualMatrix &A, const VirtualMatrix &B);      // subtraction
 
 // overloaded matrix multiplication
 Matrix operator* (const VirtualMatrix &A, const VirtualMatrix &B);  // multiplication
 
 // overloaded scalar-matrix multiplication
-Matrix operator* (const VirtualMatrix &A, double b);			// overloaded *
-Matrix operator* (double b, const VirtualMatrix &A);			// overloaded *
+Matrix operator* (const VirtualMatrix &A, double b);                    // overloaded *
+Matrix operator* (double b, const VirtualMatrix &A);                    // overloaded *
 
 // overloaded negative
-Matrix operator- (const VirtualMatrix &A);	// negative
+Matrix operator- (const VirtualMatrix &A);      // negative
 
 Vect3 Cross(Vect3& a, Vect3& b);
 Mat3x3 CrossMat(Vect3& a);
@@ -80,5 +81,5 @@ Matrix HStack(VirtualMatrix& A, VirtualMatrix& B);
 
 void Set6DAngularVector(Vect6& v6, Vect3& v3);
 void Set6DLinearVector(Vect6& v6, Vect3& v3);
-
+}
 #endif
diff --git a/lib/poems/mixedjoint.cpp b/lib/poems/mixedjoint.cpp
index ef4ccd25df..1107131545 100644
--- a/lib/poems/mixedjoint.cpp
+++ b/lib/poems/mixedjoint.cpp
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: mixedjoint.cpp 	                                       *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,10 +11,10 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
- 
+
 #include "mixedjoint.h"
 
 #include <iostream>
@@ -29,9 +29,11 @@
 #include "vect3.h"
 #include "virtualmatrix.h"
 
-   
+using namespace POEMS;
+using namespace std;
 
-MixedJoint::MixedJoint(){		 
+
+MixedJoint::MixedJoint(){		
 }
 
 MixedJoint::~MixedJoint(){
@@ -54,23 +56,23 @@ void MixedJoint::SetsP(Matrix& sPr, Vect6& temp_dofs, int i, int j){
 	numrots = i;
 	numtrans = j;
 	if (numrots < 2)
-		DimQandU(numrots+numtrans,numrots+numtrans);  
-	else 
-		DimQandU((4+numtrans),(numrots+numtrans)); 
+		DimQandU(numrots+numtrans,numrots+numtrans);
+	else
+		DimQandU((4+numtrans),(numrots+numtrans));
 	cout<<"Check "<<4+numtrans<<" "<<numrots+numtrans<<" "<<i<<" "<<j<<endl;
 }
 
 void MixedJoint::ComputeLocalTransform(){
-  Mat3x3 ko_C_k;  
-  EP_Transformation(q, ko_C_k);   
-  FastMult(pk_C_ko,ko_C_k,pk_C_k);  
+  Mat3x3 ko_C_k;
+  EP_Transformation(q, ko_C_k);
+  FastMult(pk_C_ko,ko_C_k,pk_C_k);
 }
 
 Matrix MixedJoint::GetForward_sP(){
   Mat6x6 temp_sP;
   Matrix sP;
-    
-  temp_sP.Zeros();  
+
+  temp_sP.Zeros();
   Mat3x3 temp0=T(pk_C_k);
   for(int i=1;i<4;i++){
 	  temp_sP(i,i)=1.0;
@@ -78,37 +80,37 @@ Matrix MixedJoint::GetForward_sP(){
 		  temp_sP(3+i,3+j)=temp0(i,j);
 	  }
   }
-  sP = temp_sP*const_sP;  
-  return sP;    
+  sP = temp_sP*const_sP;
+  return sP;
 }
 
 Matrix MixedJoint::GetBackward_sP(){
   Mat6x6 sP;
-  sP.Identity();        
+  sP.Identity();
   sP =-1.0*sP;
   cout<<"Did I come here in "<<endl;
-  return sP;    
+  return sP;
 }
 
 
 void MixedJoint::UpdateForward_sP( Matrix& sP){
   // do nothing
 }
-  
+
 void MixedJoint::UpdateBackward_sP( Matrix& sP){
   // do nothing
-}  
-  
+}
+
 void MixedJoint::ForwardKinematics(){	
 	if(numrots > 1)
 		EP_Normalize(q);
 
   // COMMENT STEP1: CALCULATE ORIENTATIONS
-  ComputeForwardTransforms();  
-  
-    
+  ComputeForwardTransforms();
+
+
   //COMMENT STEP2: CALCULATE POSITION VECTORS
-  Vect3 result1, result2, result3, result4;   
+  Vect3 result1, result2, result3, result4;
   result1.Zeros();
   for (int k=0; k<3; k++){
 	  if( dofs(3+k) != 0.0 ){
@@ -116,32 +118,32 @@ void MixedJoint::ForwardKinematics(){
 			  result1.BasicSet(k,q.BasicGet(4 + k));
 		  else
 			  result1.BasicSet(k,q.BasicGet(numrots + k));
-	  }  
+	  }
   }
-    
-  
-  
-  FastAssign(result1,r12);     // r12 in parents basis  i.e. Newtonian  
-  FastNegMult(k_C_pk,r12,r21);	    // r21 in body basis           
-  
-  FastAssign(r12,body2->r);  // This is right     
-    
-  //COMMENT STEP3: CALCULATE QDOT  
+
+
+
+  FastAssign(result1,r12);     // r12 in parents basis  i.e. Newtonian
+  FastNegMult(k_C_pk,r12,r21);	    // r21 in body basis
+
+  FastAssign(r12,body2->r);  // This is right
+
+  //COMMENT STEP3: CALCULATE QDOT
   int pp = 0;
   if (numrots > 1){
 	  ColMatrix temp_u(3+numtrans);
 	  qdot_to_u(q,temp_u,qdot);
 	  for (int k=1;k<=6;k++){
 		  if(dofs(k) != 0.0){
-			  u.BasicSet(pp,temp_u.BasicGet(k-1));			  
+			  u.BasicSet(pp,temp_u.BasicGet(k-1));			
 			  pp = pp+1;
-		  }		  
+		  }		
 	  }
   }
   else u = qdot;
 
-  
-  
+
+
   Vect3 WN; WN.Zeros();
   int p = 0;
   for (int k=0;k<3;k++){
@@ -151,51 +153,51 @@ void MixedJoint::ForwardKinematics(){
 	  }
   }// WN is in body basis
 
-  
-  Vect3 VN; VN.Zeros();  
+
+  Vect3 VN; VN.Zeros();
   for (int k=0;k<3;k++){
 	  if( dofs(3+k+1) != 0.0 ) {
 		  VN.BasicSet(k,u.BasicGet(p));
 		  p=p+1;
 	  }
   }// VN is the vector of translational velocities in Newtonian basis
-    
-  FastAssign(WN,body2->omega_k);  
 
-  // cout<<"Angular Velocity "<<WN<<endl;  
+  FastAssign(WN,body2->omega_k);
+
+  // cout<<"Angular Velocity "<<WN<<endl;
   Vect3 pk_w_k;
-  FastMult(body2->n_C_k,WN,pk_w_k);      
+  FastMult(body2->n_C_k,WN,pk_w_k);
   FastAssign(pk_w_k,body2->omega);
-  
-  
-  
-  //COMMENT STEP5: CALCULATE VELOCITES        
-  FastAssign(VN,body2->v);    
-  FastTMult(body2->n_C_k,body2->v,body2->v_k);        
-  
-  
-  //CALCULATE KE          
- 
+
+
+
+  //COMMENT STEP5: CALCULATE VELOCITES
+  FastAssign(VN,body2->v);
+  FastTMult(body2->n_C_k,body2->v,body2->v_k);
+
+
+  //CALCULATE KE
+
   Matrix tempke;
-  tempke = T(body2->v)*(body2->v);    
+  tempke = T(body2->v)*(body2->v);
   double ke = 0.0;
-  ke = body2->mass*tempke(1,1);    
+  ke = body2->mass*tempke(1,1);
   FastMult(body2->inertia,body2->omega_k,result1);
-  tempke= T(body2->omega_k)*result1;  
+  tempke= T(body2->omega_k)*result1;
   ke = 0.5*ke + 0.5*tempke(1,1);
   body2->KE=ke;
-  
-  
+
+
   //COMMENT STEP6: CALCULATE STATE EXPLICIT ANGULAR ACCELERATIONS
   body2->alpha_t.Zeros();
-      
-    
+
+
   //COMMENT STEP7: CALCULATE STATE EXPLICIT ACCELERATIONS
-  body2->a_t.Zeros();   
-  
+  body2->a_t.Zeros();
+
   }
 
 void MixedJoint::BackwardKinematics(){
 cout<<"Did I come here "<<endl;
-         
+
 }
diff --git a/lib/poems/mixedjoint.h b/lib/poems/mixedjoint.h
index 417f689523..b83bef4a9c 100644
--- a/lib/poems/mixedjoint.h
+++ b/lib/poems/mixedjoint.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: mixedjoint.h		                                       *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,7 +11,7 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
@@ -24,16 +24,16 @@
 #include "matrix.h"
 #include "vect6.h"
 
-
+namespace POEMS {
 class MixedJoint : public Joint{
   Matrix const_sP;
   int numrots;
   int numtrans;
   Vect6 dofs;
-public: 
+public:
   MixedJoint();
   ~MixedJoint();
-  
+
   JointType GetType();
   bool ReadInJointData(std::istream& in);
   void WriteOutJointData(std::ostream& out);
@@ -46,5 +46,5 @@ public:
   void ForwardKinematics();
   void BackwardKinematics();
 };
-
+}
 #endif
diff --git a/lib/poems/norm.cpp b/lib/poems/norm.cpp
index 6b8a56fa29..6c983dd06f 100644
--- a/lib/poems/norm.cpp
+++ b/lib/poems/norm.cpp
@@ -24,6 +24,8 @@
 #include "vect4.h"
 #include "vect6.h"
 
+using namespace POEMS;
+
 double Magnitude(ColMatrix& A){
 	double G;
 	G = 0;
diff --git a/lib/poems/norm.h b/lib/poems/norm.h
index b3ff285dea..808be609e8 100644
--- a/lib/poems/norm.h
+++ b/lib/poems/norm.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: norm.h                                                  *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,13 +11,14 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
 #ifndef NORM_H
 #define NORM_H
 
+namespace POEMS {
 class ColMatrix;
 class RowMatrix;
 class Vect3;
@@ -29,5 +30,5 @@ double Magnitude(RowMatrix& A);
 double Magnitude(Vect3& A);
 double Magnitude(Vect4& A);
 double Magnitude(Vect6& A);
-
+}
 #endif
diff --git a/lib/poems/onbody.cpp b/lib/poems/onbody.cpp
index bf843be161..7cea2d2f25 100644
--- a/lib/poems/onbody.cpp
+++ b/lib/poems/onbody.cpp
@@ -30,6 +30,8 @@
 #include "colmatrix.h"
 
 using namespace std;
+using namespace POEMS;
+
 
 
 OnBody::OnBody(){
diff --git a/lib/poems/onbody.h b/lib/poems/onbody.h
index 20190dfee8..df82709ad5 100644
--- a/lib/poems/onbody.h
+++ b/lib/poems/onbody.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: onbody.h                                              *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,7 +11,7 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
@@ -26,6 +26,7 @@
 #include "mat3x3.h"
 #include "vect3.h"
 
+namespace POEMS {
 class Body;
 class ColMatrix;
 class InertialFrame;
@@ -37,65 +38,64 @@ enum Direction {
 };
 
 class OnBody {
-       Body* system_body;
-       Joint* system_joint;
-       OnBody* parent;
-       List<OnBody> children;
+  Body* system_body;
+  Joint* system_joint;
+  OnBody* parent;
+  List<OnBody> children;
 
-       Direction joint_dir;
-       void (Joint::*kinfun)();          // kinematics function
-       void (Joint::*updatesP)(Matrix&); // sP update function
-       Vect3* gamma;                     // pointer to gamma vector
-       Mat3x3* pk_C_k;                   // pointer to transformation
+  Direction joint_dir;
+  void (Joint::*kinfun)();          // kinematics function
+  void (Joint::*updatesP)(Matrix&); // sP update function
+  Vect3* gamma;                     // pointer to gamma vector
+  Mat3x3* pk_C_k;                   // pointer to transformation
 
-  
-       Mat6x6 sI;      // spatial inertias
-       Mat6x6 sIhat;   // recursive spatial inertias
-       Mat6x6 sSC;     // spatial shift
-       Mat6x6 sT;      // spatial triangularization
 
-       Vect6 sF;       // spatial forces
-       Vect6 sFhat;    // recursive spatial forces
-       Vect6 sAhat;    // recursive spatial acceleration
+  Mat6x6 sI;      // spatial inertias
+  Mat6x6 sIhat;   // recursive spatial inertias
+  Mat6x6 sSC;     // spatial shift
+  Mat6x6 sT;      // spatial triangularization
 
-       Matrix sP;      // spatial partial velocities
-       Matrix sM;      // triangularized mass matrix diagonal elements
-       Matrix sMinv;   // inverse of sM
-       Matrix sPsMinv;
-       Matrix sIhatsP;
+  Vect6 sF;       // spatial forces
+  Vect6 sFhat;    // recursive spatial forces
+  Vect6 sAhat;    // recursive spatial acceleration
 
-       // states and state derivatives
-       ColMatrix* q;
-       ColMatrix* u;
-       ColMatrix* qdot;
-       ColMatrix* udot;
-       ColMatrix* qdotdot;
-	  
-       ColMatrix* r;       
-       ColMatrix* acc;
-       ColMatrix* ang;
+  Matrix sP;      // spatial partial velocities
+  Matrix sM;      // triangularized mass matrix diagonal elements
+  Matrix sMinv;   // inverse of sM
+  Matrix sPsMinv;
+  Matrix sIhatsP;
+
+  // states and state derivatives
+  ColMatrix* q;
+  ColMatrix* u;
+  ColMatrix* qdot;
+  ColMatrix* udot;
+  ColMatrix* qdotdot;
+
+  ColMatrix* r;
+  ColMatrix* acc;
+  ColMatrix* ang;
 
   // friend classes
   friend class OnSolver;
-  
 
-public: 
-       OnBody();
-       ~OnBody();
-       int RecursiveSetup(InertialFrame* basebody);
-       int RecursiveSetup(int ID, OnBody* parentbody, Joint* sys_joint);
-       void RecursiveKinematics();
-       void RecursiveTriangularization();
-       void RecursiveForwardSubstitution();
-       Mat3x3 GetN_C_K();
-       Vect3 LocalCart();
-       int GetBodyID();
-       void CalculateAcceleration();       
-       void Setup();
-       void SetupInertialFrame();
-       void LocalKinematics();       
-		 void LocalTriangularization(Vect3& Torque, Vect3& Force);
-       void LocalForwardSubstitution();
+public:
+  OnBody();
+  ~OnBody();
+  int RecursiveSetup(InertialFrame* basebody);
+  int RecursiveSetup(int ID, OnBody* parentbody, Joint* sys_joint);
+  void RecursiveKinematics();
+  void RecursiveTriangularization();
+  void RecursiveForwardSubstitution();
+  Mat3x3 GetN_C_K();
+  Vect3 LocalCart();
+  int GetBodyID();
+  void CalculateAcceleration();
+  void Setup();
+  void SetupInertialFrame();
+  void LocalKinematics();
+  void LocalTriangularization(Vect3& Torque, Vect3& Force);
+  void LocalForwardSubstitution();
 };
-
+}
 #endif
diff --git a/lib/poems/onfunctions.cpp b/lib/poems/onfunctions.cpp
index 1d91dbf378..cfc0c60936 100644
--- a/lib/poems/onfunctions.cpp
+++ b/lib/poems/onfunctions.cpp
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: onfunction.cpp                                          *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,10 +11,10 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
- 
+
 
 #include "onfunctions.h"
 
@@ -25,9 +25,11 @@
 #include "virtualmatrix.h"
 
 using namespace std;
+using namespace POEMS;
+
 
 // friend of Vect3 & Vect6
-void OnPopulateSVect(Vect3& angular, Vect3& linear, Vect6& sV){
+void POEMS::OnPopulateSVect(Vect3& angular, Vect3& linear, Vect6& sV){
   sV.elements[0] = angular.elements[0];
   sV.elements[1] = angular.elements[1];
   sV.elements[2] = angular.elements[2];
@@ -37,38 +39,38 @@ void OnPopulateSVect(Vect3& angular, Vect3& linear, Vect6& sV){
 }
 
 // friend of Vect3, Mat3x3, & Mat6x6
-void OnPopulateSC(Vect3& gamma, Mat3x3& C, Mat6x6& SC){
+void POEMS::OnPopulateSC(Vect3& gamma, Mat3x3& C, Mat6x6& SC){
   // the block diagonals
-  
+
   // the gamma cross with transform
   Mat3x3 temp; Mat3x3 temp2;
   SC.Zeros();
-  temp.Zeros();  
+  temp.Zeros();
   temp2.Zeros();
-  //FastTMult(C,gamma,temp);     
+  //FastTMult(C,gamma,temp);
   temp(1,2)= -gamma(3); temp(1,3)= gamma(2); temp(2,1)= gamma(3);
-  temp(2,3)= -gamma(1); temp(3,1)= -gamma(2); temp(3,2)= gamma(1);  
+  temp(2,3)= -gamma(1); temp(3,1)= -gamma(2); temp(3,2)= gamma(1);
   FastMult(temp,C,temp2);
-  
+
   SC(1,4)=temp2(1,1);   SC(2,4)=temp2(2,1);   SC(3,4)=temp2(3,1);
   SC(1,5)=temp2(1,2);   SC(2,5)=temp2(2,2);   SC(3,5)=temp2(3,2);
   SC(1,6)=temp2(1,3);   SC(2,6)=temp2(2,3);   SC(3,6)=temp2(3,3);
-  
+
   SC(1,1)=C(1,1);   SC(2,1)=C(2,1);   SC(3,1)=C(3,1);
   SC(1,2)=C(1,2);   SC(2,2)=C(2,2);   SC(3,2)=C(3,2);
   SC(1,3)=C(1,3);   SC(2,3)=C(2,3);   SC(3,3)=C(3,3);
 
   SC(4,4)=C(1,1);   SC(5,4)=C(2,1);   SC(6,4)=C(3,1);
   SC(4,5)=C(1,2);   SC(5,5)=C(2,2);   SC(6,5)=C(3,2);
-  SC(4,6)=C(1,3);   SC(5,6)=C(2,3);   SC(6,6)=C(3,3);  
-  
+  SC(4,6)=C(1,3);   SC(5,6)=C(2,3);   SC(6,6)=C(3,3);
+
   }
 
 // friend of Mat3x3 & Mat6x6
-void OnPopulateSI(Mat3x3& inertia, double mass, Mat6x6& sI){
-  
-	sI(4,4)=mass; sI(5,5)=mass; sI(6,6)=mass;
-	sI(1,1)=inertia(1,1); sI(1,2)=inertia(1,2); sI(1,3)=inertia(1,3);
-	sI(2,1)=inertia(2,1); sI(2,2)=inertia(2,2); sI(2,3)=inertia(2,3);
-	sI(3,1)=inertia(3,1); sI(3,2)=inertia(3,2); sI(3,3)=inertia(3,3);	
+void POEMS::OnPopulateSI(Mat3x3& inertia, double mass, Mat6x6& sI){
+
+        sI(4,4)=mass; sI(5,5)=mass; sI(6,6)=mass;
+        sI(1,1)=inertia(1,1); sI(1,2)=inertia(1,2); sI(1,3)=inertia(1,3);
+        sI(2,1)=inertia(2,1); sI(2,2)=inertia(2,2); sI(2,3)=inertia(2,3);
+        sI(3,1)=inertia(3,1); sI(3,2)=inertia(3,2); sI(3,3)=inertia(3,3);
 }
diff --git a/lib/poems/onfunctions.h b/lib/poems/onfunctions.h
index 42d0f4ac1c..bb07425b85 100644
--- a/lib/poems/onfunctions.h
+++ b/lib/poems/onfunctions.h
@@ -18,6 +18,7 @@
 #ifndef ONFUNCTIONS_H
 #define ONFUNCTIONS_H
 
+namespace POEMS {
 class Mat3x3;
 class Mat6x6;
 class Vect3;
@@ -30,5 +31,5 @@ void OnPopulateSI(Mat3x3& inertia, double mass, Mat6x6& sI);
 
 void Create_Map(int MM);
 int ICELL(int IX,int IY,int IZ, int MM);
-
+}
 #endif
diff --git a/lib/poems/onsolver.cpp b/lib/poems/onsolver.cpp
index dd94c9c7ad..609ef32613 100644
--- a/lib/poems/onsolver.cpp
+++ b/lib/poems/onsolver.cpp
@@ -33,6 +33,8 @@
 #include "virtualmatrix.h"
 
 using namespace std;
+using namespace POEMS;
+
 
 OnSolver::OnSolver(){
   numbodies = 0;
diff --git a/lib/poems/onsolver.h b/lib/poems/onsolver.h
index b4bbd5e25d..1e850b8fc1 100644
--- a/lib/poems/onsolver.h
+++ b/lib/poems/onsolver.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: onsolver.h                                              *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,7 +11,7 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
@@ -22,32 +22,31 @@
 #include "solver.h"
 #include "onbody.h"
 
+namespace POEMS {
 class ColMatrix;
 class Matrix;
 
-
-  
 class OnSolver : public Solver {
-       OnBody inertialframe;   
-       int numbodies;
-       OnBody** bodyarray;  
-       ColMatrix** q;		 
-       ColMatrix** qdot;
-		 ColMatrix** qdotdot;
-       ColMatrix** u;
-       ColMatrix** udot;
-       
-       
-       
-       void DeleteModel();
-       int CreateTopologyArray(int i, OnBody* body);
-       void CreateStateMatrixMaps();
-	    void GetType();
-public: 
-       OnSolver();
-       ~OnSolver();
-       void CreateModel();
-       void Solve(double time, Matrix& FF);       
-};
+  OnBody inertialframe;
+  int numbodies;
+  OnBody** bodyarray;
+  ColMatrix** q;
+  ColMatrix** qdot;
+                 ColMatrix** qdotdot;
+  ColMatrix** u;
+  ColMatrix** udot;
 
+
+
+  void DeleteModel();
+  int CreateTopologyArray(int i, OnBody* body);
+  void CreateStateMatrixMaps();
+            void GetType();
+public:
+  OnSolver();
+  ~OnSolver();
+  void CreateModel();
+  void Solve(double time, Matrix& FF);
+};
+}
 #endif
diff --git a/lib/poems/particle.cpp b/lib/poems/particle.cpp
index 0d953e00f6..37ec790666 100644
--- a/lib/poems/particle.cpp
+++ b/lib/poems/particle.cpp
@@ -18,6 +18,8 @@
 
 #include "particle.h"
 
+using namespace POEMS;
+
 Particle::Particle(){
 }
 
diff --git a/lib/poems/particle.h b/lib/poems/particle.h
index 0f13252f99..8a7a521451 100644
--- a/lib/poems/particle.h
+++ b/lib/poems/particle.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: particle.h                                              *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,7 +11,7 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
@@ -23,14 +23,14 @@
 
 #include "body.h"
 
-
+namespace POEMS {
 class Particle : public Body  {
-public: 
+public:
   Particle();
   ~Particle();
   BodyType GetType();
   bool ReadInBodyData(std::istream& in);
   void WriteOutBodyData(std::ostream& out);
 };
-
+}
 #endif
diff --git a/lib/poems/poemslist.h b/lib/poems/poemslist.h
index 9882af013a..bd0c1456b4 100644
--- a/lib/poems/poemslist.h
+++ b/lib/poems/poemslist.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: poemslist.h                                             *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,7 +11,7 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
@@ -21,8 +21,7 @@
 #include <iostream>
 #include <cstdlib>
 
-using namespace std;
-
+namespace POEMS {
 template<class T> class ListElement{
 public:
   ListElement<T>* prev;
@@ -51,7 +50,7 @@ public:
   S* operator()(int id);
   void Append(List<S> * listToAppend);
   void DeleteValues();
-  void RemoveElementAndDeleteValue(ListElement<S>* ele);  
+  void RemoveElementAndDeleteValue(ListElement<S>* ele);
   void PrintList();
 };
 
@@ -90,9 +89,9 @@ template<class S> List<S>::~List(){
 
 template<class S> void List<S>::Append(List<S> * listToAppend)
 {
-	tail->next = listToAppend->head;
-	listToAppend->head->prev = tail;
-	tail = listToAppend->tail;
+   tail->next = listToAppend->head;
+   listToAppend->head->prev = tail;
+   tail = listToAppend->tail;
 }
 
 template<class S> int List<S>::GetNumElements(){
@@ -104,17 +103,17 @@ template<class S> ListElement<S>* List<S>::GetHeadElement(){
 }
 
 template<class S> ListElement<S>* List<S>::GetTailElement(){
-	return tail;
+   return tail;
 }
 
 
 template<class S> void List<S>::Remove(ListElement<S>* ele){
   if(!ele){
-    cerr << "ERROR: ListElement to be removed not defined" << endl;
+    std::cerr << "ERROR: ListElement to be removed not defined" << std::endl;
     exit(0);
   }
   if(!numelements){
-    cerr << "ERROR: List is empty" << endl;
+    std::cerr << "ERROR: List is empty" << std::endl;
     exit(0);
   }
 
@@ -135,7 +134,7 @@ template<class S> void List<S>::Remove(ListElement<S>* ele){
 
 template<class S> ListElement<S>* List<S>::Append(S* v){
   if(!v){
-    cerr << "ERROR: cannot add null Body to list" << endl;
+    std::cerr << "ERROR: cannot add null Body to list" << std::endl;
     exit(0);
   }
 
@@ -145,22 +144,22 @@ template<class S> ListElement<S>* List<S>::Append(S* v){
   if(numelements==1)
     head = tail = ele;
   else{
-	  /*
+     /*
     tail->next = ele;
     ele->prev = tail;
-	 tail = ele;*/
-	  
-	  ele->prev = tail;
-	  tail = ele;	  
-	  ele->prev->next = ele;	  
-	  
-  }    
+    tail = ele;*/
+
+     ele->prev = tail;
+     tail = ele;
+     ele->prev->next = ele;
+
+  }
   return ele;
 }
 
 template<class S> ListElement<S>* List<S>::Prepend(S* v){
   if(!v){
-    cerr << "ERROR: cannot add null Body to list" << endl;
+    std::cerr << "ERROR: cannot add null Body to list" << std::endl;
     exit(0);
   }
 
@@ -170,9 +169,9 @@ template<class S> ListElement<S>* List<S>::Prepend(S* v){
   if(numelements==1)
     head = tail = ele;
   else{
-	  ele->next = head;
-	  head = ele;
-	  ele->next->prev = ele;
+     ele->next = head;
+     head = ele;
+     ele->next->prev = ele;
   }
   return ele;
 }
@@ -190,15 +189,15 @@ template<class S> S** List<S>::CreateArray(){
 
 template<class S> S* List<S>::operator()(int id){
   if(id >= numelements){
-    cerr << "ERROR: subscript out of bounds" << endl;
+    std::cerr << "ERROR: subscript out of bounds" << std::endl;
     exit(0);
   }
-  
+
   ListElement<S>* ele = head;
   for(int i=0;i<id;i++){
     ele = ele->next;
   }
-  
+
   return ele->value;
 }
 
@@ -214,16 +213,15 @@ template<class S> void List<S>::RemoveElementAndDeleteValue(ListElement<S>* ele)
 }
 
 template<class S> void List<S>::PrintList(){
-	cout<<"Printing List "<<endl;
-	ListElement<S>* ele = head;
-	cout<<*(ele->value)<<" ";
-	ele = ele->next;
-	for(int k =2; k<numelements; k++){
-		cout<<*(ele->value)<<" ";	
-		ele = ele->next;
-	}
-	cout<<*(ele->value)<<endl;
+   std::cout << "Printing List " << std::endl;
+   ListElement<S>* ele = head;
+   std::cout << *(ele->value) << " ";
+   ele = ele->next;
+   for(int k =2; k<numelements; k++){
+     std::cout << *(ele->value) << " ";
+     ele = ele->next;
+   }
+   std::cout << *(ele->value) << std::endl;
+}
 }
-
-
 #endif
diff --git a/lib/poems/poemsnodelib.h b/lib/poems/poemsnodelib.h
index 06cf4530d2..d323b701b0 100644
--- a/lib/poems/poemsnodelib.h
+++ b/lib/poems/poemsnodelib.h
@@ -20,9 +20,7 @@
 
 #include <iostream>
 
-using namespace std;
-
-
+namespace POEMS {
 TreeNode *GetTreeNode(int item,TreeNode *lptr = NULL,TreeNode *rptr =NULL); 
 
 void FreeTreeNode(TreeNode *p);
@@ -81,7 +79,7 @@ TreeNode *GetTreeNode(int item,TreeNode *lptr,TreeNode *rptr)
 	// if insufficient memory, terminatewith an error message
 	if (p == NULL)
 	{
-		cerr << "Memory allocation failure!\n";
+                std::cerr << "Memory allocation failure!\n";
 		exit(1);
 	}
 
@@ -140,7 +138,7 @@ void IndentBlanks(int num)
 //	const int indentblock = 6;
 
 	for(int i = 0; i < num; i++)
-		cout << " ";
+		std::cout << " ";
 }
 
 void PrintTree (TreeNode *t, int level)
@@ -153,10 +151,11 @@ void PrintTree (TreeNode *t, int level)
 		PrintTree(t->Right(),level + 1);
 		// indent to current level; output node data
 		IndentBlanks(indentUnit*level);
-		cout << t->GetData() << endl;
+		std::cout << t->GetData() << std::endl;
 		// print left branch of tree t
 		PrintTree(t->Left(),level + 1);
 	}
 }
+}
 #endif
 
diff --git a/lib/poems/poemsobject.cpp b/lib/poems/poemsobject.cpp
index 7c3f1ca872..dfdd95953e 100644
--- a/lib/poems/poemsobject.cpp
+++ b/lib/poems/poemsobject.cpp
@@ -19,6 +19,8 @@
 #include "poemsobject.h"
 #include <cstring>
 
+using namespace POEMS;
+
 POEMSObject::POEMSObject(){
   name = 0;
   ChangeName((const char*)"unnamed");
diff --git a/lib/poems/poemsobject.h b/lib/poems/poemsobject.h
index 63b2915638..173d9da7ab 100644
--- a/lib/poems/poemsobject.h
+++ b/lib/poems/poemsobject.h
@@ -19,7 +19,7 @@
 #ifndef POEMSOBJECT_H
 #define POEMSOBJECT_H
 
-
+namespace POEMS {
 class POEMSObject {
   char* name;
   int ID;
@@ -31,5 +31,5 @@ public:
   int GetID();
   void SetID(int id);
 };
-
+}
 #endif
diff --git a/lib/poems/poemstree.h b/lib/poems/poemstree.h
index 8f8e80ab66..d49a0d0a88 100644
--- a/lib/poems/poemstree.h
+++ b/lib/poems/poemstree.h
@@ -21,7 +21,7 @@
 #include "poemstreenode.h"
 #include "poemsnodelib.h"
 
-
+namespace POEMS {
 // constants to indicate the balance factor of a node
 const int leftheavy = -1;
 const int balanced = 0;
@@ -609,5 +609,5 @@ void Tree::ClearList(void)
 	delete current;
 	size = 0;
 }
-
+}
 #endif
diff --git a/lib/poems/poemstreenode.cpp b/lib/poems/poemstreenode.cpp
index 051e13d6e0..ca4307f69f 100644
--- a/lib/poems/poemstreenode.cpp
+++ b/lib/poems/poemstreenode.cpp
@@ -17,6 +17,8 @@
 
 #include "poemstreenode.h"
 
+using namespace POEMS;
+
 // constructor; initialize the data and pointer fields
 // The pointer value NULL assigns a empty subtree
 TreeNode::TreeNode (const int & item, TreeNode *lptr,TreeNode *rptr,
diff --git a/lib/poems/poemstreenode.h b/lib/poems/poemstreenode.h
index e22bd9b97d..69f3caf958 100644
--- a/lib/poems/poemstreenode.h
+++ b/lib/poems/poemstreenode.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: poemstreenode.h                                         *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,36 +11,37 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
 #ifndef TREENODE_H
 #define TREENODE_H
 
+namespace POEMS {
 // declares a tree node object for a binary tree
 class TreeNode{
 
 private:
 // points to the left and right children of the node
-	TreeNode *left;
-	TreeNode *right;
+        TreeNode *left;
+        TreeNode *right;
 
-	int balanceFactor;
-	int data;
-	void * aux_data;
-public:	
-	// make Tree a friend because it needs access to left and right pointer fields of a node
-	friend class Tree;
-	TreeNode * Left();
-	TreeNode * Right();
-	int GetData();
-	void * GetAuxData() {return aux_data;};
-	void SetAuxData(void * AuxData) {aux_data = AuxData;};
-	int GetBalanceFactor();
-	TreeNode(const int &item, TreeNode *lptr, TreeNode *rptr, int balfac = 0);
-	//friend class DCASolver;
+        int balanceFactor;
+        int data;
+        void * aux_data;
+public:
+        // make Tree a friend because it needs access to left and right pointer fields of a node
+        friend class Tree;
+        TreeNode * Left();
+        TreeNode * Right();
+        int GetData();
+        void * GetAuxData() {return aux_data;};
+        void SetAuxData(void * AuxData) {aux_data = AuxData;};
+        int GetBalanceFactor();
+        TreeNode(const int &item, TreeNode *lptr, TreeNode *rptr, int balfac = 0);
+        //friend class DCASolver;
 };
-
+}
 #endif
 
diff --git a/lib/poems/point.cpp b/lib/poems/point.cpp
index d7eed94efd..f648b58d3d 100644
--- a/lib/poems/point.cpp
+++ b/lib/poems/point.cpp
@@ -20,6 +20,8 @@
 #include "point.h"
 #include "vect3.h"
 
+using namespace POEMS;
+
 Point::Point(){
   position.Zeros();
 }
@@ -35,7 +37,7 @@ void Point::WriteOut(std::ostream& out){
   WriteOutPointData(out);
 }
 
-Point* NewPoint(int type){
+Point* POEMS::NewPoint(int type){
   switch( PointType(type) )
   {
     case FIXEDPOINT :  // A Fixed Point
diff --git a/lib/poems/point.h b/lib/poems/point.h
index 746a527cd7..c08e92db0d 100644
--- a/lib/poems/point.h
+++ b/lib/poems/point.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: point.h                                                  *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,7 +11,7 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
@@ -22,21 +22,21 @@
 #include "poemsobject.h"
 #include "vect3.h"
 
-
+namespace POEMS {
 // emumerated type
 enum PointType {
   FIXEDPOINT = 0
 };
 
 class Point : public POEMSObject {
-public: 
+public:
   Vect3 position;
 
-  Point();   
-  bool ReadIn(std::istream& in);  
+  Point();
+  bool ReadIn(std::istream& in);
   void WriteOut(std::ostream& out);
-  
-  virtual ~Point();  
+
+  virtual ~Point();
   virtual PointType GetType() = 0;
   virtual Vect3 GetPoint() = 0;
   virtual bool ReadInPointData(std::istream& in) = 0;
@@ -45,5 +45,5 @@ public:
 
 // global point functions
 Point* NewPoint(int type);
-
+}
 #endif
diff --git a/lib/poems/prismaticjoint.cpp b/lib/poems/prismaticjoint.cpp
index c7d4bed568..3e649f0d22 100644
--- a/lib/poems/prismaticjoint.cpp
+++ b/lib/poems/prismaticjoint.cpp
@@ -25,6 +25,8 @@
 #include "mat3x3.h"
 #include "virtualmatrix.h"
 
+using namespace POEMS;
+
 PrismaticJoint::PrismaticJoint(){
   q.Dim(1);
   qdot.Dim(1);
diff --git a/lib/poems/prismaticjoint.h b/lib/poems/prismaticjoint.h
index c17be350cb..df3ab3195c 100644
--- a/lib/poems/prismaticjoint.h
+++ b/lib/poems/prismaticjoint.h
@@ -24,8 +24,7 @@
 #include "vect3.h"
 #include "matrix.h"
 
-
-
+namespace POEMS {
 class PrismaticJoint : public Joint  {
   Vect3 axis_pk; // unit vector in body1 basis
   Vect3 axis_k;  // unit vector in body2 basis
@@ -45,5 +44,5 @@ public:
   void ForwardKinematics();
   void BackwardKinematics();
 };
-
+}
 #endif
diff --git a/lib/poems/revolutejoint.cpp b/lib/poems/revolutejoint.cpp
index 955a2f0dea..6b11acbf0f 100644
--- a/lib/poems/revolutejoint.cpp
+++ b/lib/poems/revolutejoint.cpp
@@ -25,6 +25,8 @@
 #include "mat3x3.h"
 #include "virtualmatrix.h"
 
+using namespace POEMS;
+
 RevoluteJoint::RevoluteJoint(){
   DimQandU(1);
   Vect3 axis;
diff --git a/lib/poems/revolutejoint.h b/lib/poems/revolutejoint.h
index ac1c06d5c5..10a38114df 100644
--- a/lib/poems/revolutejoint.h
+++ b/lib/poems/revolutejoint.h
@@ -23,6 +23,7 @@
 #include "vect3.h"
 #include "matrix.h"
 
+namespace POEMS {
 class VirtualMatrix;
 
 class RevoluteJoint : public Joint  {
@@ -44,5 +45,5 @@ public:
   void ForwardKinematics();
   void BackwardKinematics();
 };
-
+}
 #endif
diff --git a/lib/poems/rigidbody.cpp b/lib/poems/rigidbody.cpp
index 53acfa8bd4..a7d4f9f19f 100644
--- a/lib/poems/rigidbody.cpp
+++ b/lib/poems/rigidbody.cpp
@@ -20,6 +20,8 @@
 #include "virtualmatrix.h"
 
 using namespace std;
+using namespace POEMS;
+
 
 RigidBody::RigidBody(){
 }
diff --git a/lib/poems/rigidbody.h b/lib/poems/rigidbody.h
index 2f4430e264..d4be92ebcf 100644
--- a/lib/poems/rigidbody.h
+++ b/lib/poems/rigidbody.h
@@ -22,6 +22,7 @@
 
 #include "body.h"
 
+namespace POEMS {
 class RigidBody : public Body  {
 public: 
   RigidBody();
@@ -30,5 +31,5 @@ public:
   bool ReadInBodyData(std::istream& in);
   void WriteOutBodyData(std::ostream& out);
 };
-
+}
 #endif
diff --git a/lib/poems/rowmatrix.cpp b/lib/poems/rowmatrix.cpp
index 9fa57e773f..eac8fbef5d 100644
--- a/lib/poems/rowmatrix.cpp
+++ b/lib/poems/rowmatrix.cpp
@@ -21,6 +21,8 @@
 #include <cstdlib>
 
 using namespace std;
+using namespace POEMS;
+
 
 RowMatrix::RowMatrix(){
   numcols = 0;
diff --git a/lib/poems/rowmatrix.h b/lib/poems/rowmatrix.h
index 1dbb3d11e0..8e8d7b5f99 100644
--- a/lib/poems/rowmatrix.h
+++ b/lib/poems/rowmatrix.h
@@ -24,6 +24,7 @@
 #include "virtualrowmatrix.h"
 #include "virtualmatrix.h"
 
+namespace POEMS {
 class RowMatrix : public VirtualRowMatrix  {
   double* elements;
 public:
@@ -52,5 +53,5 @@ public:
   RowMatrix& operator=(const VirtualMatrix& A); // overloaded =
   RowMatrix& operator*=(double b);
 };
-
+}
 #endif
diff --git a/lib/poems/solver.cpp b/lib/poems/solver.cpp
index 9d2245149f..04afb27945 100644
--- a/lib/poems/solver.cpp
+++ b/lib/poems/solver.cpp
@@ -20,6 +20,8 @@
 #include "system.h"
 #include "onsolver.h"
 
+using namespace POEMS;
+
 Solver::Solver(){
 
 }
diff --git a/lib/poems/solver.h b/lib/poems/solver.h
index 6a2debf497..489e6d41f1 100644
--- a/lib/poems/solver.h
+++ b/lib/poems/solver.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: solver.h                                                *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,7 +11,7 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:       anderk5@rpi.edu                                     *
  *_________________________________________________________________________*/
 
@@ -21,6 +21,7 @@
 #include "colmatmap.h"
 #include "defines.h"
 
+namespace POEMS {
 class System;
 class Matrix;
 
@@ -28,11 +29,11 @@ class Solver{
 protected:
   System* system;
 
-  
+
   double time;
   ColMatMap state;
   ColMatMap statedot;
-  ColMatMap statedoubledot;  
+  ColMatMap statedoubledot;
 
   SolverType type;
 
@@ -47,13 +48,13 @@ public:
 
   virtual void DeleteModel() = 0;
   virtual void CreateModel() = 0;
-  virtual void Solve(double time, Matrix& FF) = 0;  
-  
-    
-    
+  virtual void Solve(double time, Matrix& FF) = 0;
+
+
+
   ColMatMap* GetState();
   ColMatMap* GetStateDerivative();
   ColMatMap* GetStateDerivativeDerivative();
 };
-
+}
 #endif
diff --git a/lib/poems/sphericaljoint.cpp b/lib/poems/sphericaljoint.cpp
index cba6d00a6e..996bb5a887 100644
--- a/lib/poems/sphericaljoint.cpp
+++ b/lib/poems/sphericaljoint.cpp
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: sphericaljoint.cpp                                      *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,10 +11,10 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
- 
+
 
 #include "sphericaljoint.h"
 
@@ -30,6 +30,8 @@
 #include "vect3.h"
 #include "virtualmatrix.h"
 
+using namespace std;
+using namespace POEMS;
 
 SphericalJoint::SphericalJoint(){
   DimQandU(4,3);
@@ -48,24 +50,24 @@ bool SphericalJoint::ReadInJointData(std::istream& in){
 void SphericalJoint::WriteOutJointData(std::ostream& out){
 }
 
-Matrix SphericalJoint::GetForward_sP(){  
+Matrix SphericalJoint::GetForward_sP(){
   Mat3x3 sPa,sPl;
   Matrix sP(6,3);
-  sPa.Identity();    
+  sPa.Identity();
   sPl.Zeros();
-  Vect3 temp = -(point2->position);  
-  
+  Vect3 temp = -(point2->position);
+
   sPl(1,2) = temp(3);
   sPl(1,3) = -temp(2);
-  
+
   sPl(2,1) = -temp(3);
   sPl(2,3) = temp(1);
-  
-  sPl(3,1) = temp(2);  
+
+  sPl(3,1) = temp(2);
   sPl(3,2) = -temp(1);
-         
-  sP=Stack(sPa,sPl);  
-  return sP;    
+
+  sP=Stack(sPa,sPl);
+  return sP;
 }
 
 void SphericalJoint::UpdateForward_sP( Matrix& sP){
@@ -79,15 +81,15 @@ Matrix SphericalJoint::GetBackward_sP(){
 	cout<<" -----------"<<endl;	
   Mat3x3 sPa,sPl;
   Matrix sP;
-  sPa.Identity();    
-  sPl.Zeros();          
-  sPl(3,2)=(point2->position(1));  
-  sPl(2,3)=-(point2->position(1));   
+  sPa.Identity();
+  sPl.Zeros();
+  sPl(3,2)=(point2->position(1));
+  sPl(2,3)=-(point2->position(1));
   sP=Stack(sPa,sPl);
-  return sP;    
+  return sP;
 }
 
-void SphericalJoint::UpdateBackward_sP( Matrix& sP){  
+void SphericalJoint::UpdateBackward_sP( Matrix& sP){
   // sP is constant, do nothing.
 }
 
@@ -96,106 +98,106 @@ void SphericalJoint::ComputeLocalTransform(){
 	// Obtain the transformation matrix from euler parameters
 	EP_Transformation(q, ko_C_k);
 	FastMult(pk_C_ko,ko_C_k,pk_C_k);
-	}  
+	}
 
 
 void SphericalJoint::ForwardKinematics(){
   Vect3 result1,result2,result3,result4,result5;
   Vect3 pk_w_k;
-  
+
   //cout<<"Check in spherical "<<q<<" "<<qdot<<endl;
   EP_Normalize(q);
 
 
   // orientations
-  ComputeForwardTransforms();    
+  ComputeForwardTransforms();
+
 
-  
   //----------------------------------//
-  // COMPUTE POSITION VECTOR R12 aka GAMMA    
-  
-  FastNegMult(pk_C_k,point2->position,result1); // parents basis      
-  FastAdd(result1,point1->position,r12);      
-  
-  // compute position vector r21 
-  FastNegMult(k_C_pk,r12,r21); 
+  // COMPUTE POSITION VECTOR R12 aka GAMMA
+
+  FastNegMult(pk_C_k,point2->position,result1); // parents basis
+  FastAdd(result1,point1->position,r12);
+
+  // compute position vector r21
+  FastNegMult(k_C_pk,r12,r21);
+
+
 
-  
-  
   //----------------------------------//
-  // COMPUTE GLOBAL LOCATION  
-  FastMult(body1->n_C_k,(body1->GetPoint(2))->position,result1);    
-  FastAdd(result1,body1->r,result1);  
-  FastNegMult(body2->n_C_k,(body2->GetPoint(1))->position,result2);  
-  FastAdd(result1,result2,body2->r);  
-  
-  qdot_to_u(q, u, qdot);  
-  
+  // COMPUTE GLOBAL LOCATION
+  FastMult(body1->n_C_k,(body1->GetPoint(2))->position,result1);
+  FastAdd(result1,body1->r,result1);
+  FastNegMult(body2->n_C_k,(body2->GetPoint(1))->position,result2);
+  FastAdd(result1,result2,body2->r);
+
+  qdot_to_u(q, u, qdot);
+
+
+  //-----------------------------------
+  // angular velocities
 
-  //-----------------------------------  
-  // angular velocities  
-  
   FastAssign(u,pk_w_k);
-  FastTMult(pk_C_k,body1->omega_k,result1);        
+  FastTMult(pk_C_k,body1->omega_k,result1);
   FastAdd(result1,pk_w_k,body2->omega_k);
-  FastMult(body2->n_C_k,body2->omega_k,body2->omega); 
+  FastMult(body2->n_C_k,body2->omega_k,body2->omega);
 
-  
-  
-  //-----------------------------------  
-  
-  // compute velocities  
-  FastCross(body1->omega_k,(body1->GetPoint(2))->position,result1);  
+
+
+  //-----------------------------------
+
+  // compute velocities
+  FastCross(body1->omega_k,(body1->GetPoint(2))->position,result1);
   FastAdd(body1->v_k,result1,result2);
-  FastTMult(pk_C_k,result2,result1); // In body basis  
-  
-  FastCross((body2->GetPoint(1))->position,body2->omega_k,result2);                 
-  FastAdd(result1,result2,body2->v_k);    // In body basis     
-  FastMult(body2->n_C_k,body2->v_k,body2->v); 
-  
-  
+  FastTMult(pk_C_k,result2,result1); // In body basis
+
+  FastCross((body2->GetPoint(1))->position,body2->omega_k,result2);
+  FastAdd(result1,result2,body2->v_k);    // In body basis
+  FastMult(body2->n_C_k,body2->v_k,body2->v);
+
+
   //------------------------------------------
-  //Compute the KE   
+  //Compute the KE
   Matrix tempke;
-  tempke = T(body2->v)*(body2->v);    
+  tempke = T(body2->v)*(body2->v);
   double ke = 0.0;
-  ke = body2->mass*tempke(1,1);  
+  ke = body2->mass*tempke(1,1);
   FastMult(body2->inertia,body2->omega_k,result1);
-  tempke= T(body2->omega_k)*result1;  
+  tempke= T(body2->omega_k)*result1;
   ke = 0.5*ke + 0.5*tempke(1,1);
   body2->KE=ke;
 
-  //-----------------------------------  
-  // compute state explicit angular acceleration  // Has to be in body basis    
-  FastTMult(pk_C_k,body1->alpha_t,result2);  
+  //-----------------------------------
+  // compute state explicit angular acceleration  // Has to be in body basis
+  FastTMult(pk_C_k,body1->alpha_t,result2);
   FastCross(body2->omega_k,pk_w_k,result1);
-  FastAdd(result1,result2,body2->alpha_t); 
-   
-  //-----------------------------------  
-  // compute state explicit acceleration   
-  // NEED TO DO THIS COMPLETELY IN BODY BASIS 
-  
+  FastAdd(result1,result2,body2->alpha_t);
+
+  //-----------------------------------
+  // compute state explicit acceleration
+  // NEED TO DO THIS COMPLETELY IN BODY BASIS
+
   FastCross(body1->omega_k,(body1->GetPoint(2))->position,result1);
-  FastCross(body1->omega_k,result1,result2); 
-  FastTMult(pk_C_k,result2,result1); 
-  
-  //FastCross(body2->omega_k,-(body2->GetPoint(1))->position,result3);  
-  FastCross((body2->GetPoint(1))->position,body2->omega_k,result3);  
-  FastCross(body2->omega_k,result3,result2);    
-  FastAdd(result1,result2,result3); //wxwxr in body basis  
-  
-  FastCross(body1->alpha_t,(body1->GetPoint(2))->position,result4); 
-  FastTMult(pk_C_k,result4,result5);  
+  FastCross(body1->omega_k,result1,result2);
+  FastTMult(pk_C_k,result2,result1);
+
+  //FastCross(body2->omega_k,-(body2->GetPoint(1))->position,result3);
+  FastCross((body2->GetPoint(1))->position,body2->omega_k,result3);
+  FastCross(body2->omega_k,result3,result2);
+  FastAdd(result1,result2,result3); //wxwxr in body basis
+
+  FastCross(body1->alpha_t,(body1->GetPoint(2))->position,result4);
+  FastTMult(pk_C_k,result4,result5);
   FastAssign(result5,result4);
-  
-  FastCross((body2->GetPoint(1))->position,body2->alpha_t,result2);   
-  FastAdd(result2,result4,result1); //alphaxr in body basis  
-  
+
+  FastCross((body2->GetPoint(1))->position,body2->alpha_t,result2);
+  FastAdd(result2,result4,result1); //alphaxr in body basis
+
   FastTMult(pk_C_k,body1->a_t,result2);
-  FastTripleSum(result3,result1,result2,body2->a_t);     // in body basis  
-  
-  
-  //-----------------------------------    
+  FastTripleSum(result3,result1,result2,body2->a_t);     // in body basis
+
+
+  //-----------------------------------
 }
 
 // NOTE: NOT USING BACKWARDKINEMATICS AT PRESENT
@@ -206,13 +208,13 @@ void SphericalJoint::BackwardKinematics(){
 
   // orientations
   ComputeBackwardTransforms();
-  
+
 
   // compute position vector r21
   //r21 = point2->position - k_C_pk * point1->position;
   FastMult(k_C_pk,point1->position,result1);
   FastSubt(point2->position,result1,r21);
-  
+
 
   // compute position vector r21
   FastNegMult(pk_C_k,r21,r12);
@@ -231,7 +233,7 @@ void SphericalJoint::BackwardKinematics(){
   EP_Derivatives(q,u,qdot);
 
   // angular velocities
- 
+
   FastMult(body2->n_C_k,u,result2);
   FastAdd(body2->omega,result2,body1->omega);
   FastAssign(u,k_w_pk);
@@ -251,7 +253,7 @@ void SphericalJoint::BackwardKinematics(){
   FastCross(body1->omega_k,k_w_pk,result1);
   FastMult(pk_C_k,body2->alpha_t,result2);
   FastAdd(result1,result2,body1->alpha_t);
-  
+
   // compute state explicit acceleration
   FastCross(body2->alpha_t,point2->position,result1);
   FastCross(body2->omega_k,point2->position,result2);
@@ -264,5 +266,5 @@ void SphericalJoint::BackwardKinematics(){
   FastCross(body1->omega_k,result2,result3);
 
   FastTripleSum(result5,result1,result3,body1->a_t);
-  
+
 }
diff --git a/lib/poems/sphericaljoint.h b/lib/poems/sphericaljoint.h
index 90f1547873..db048e5544 100644
--- a/lib/poems/sphericaljoint.h
+++ b/lib/poems/sphericaljoint.h
@@ -22,6 +22,7 @@
 #include "joint.h"
 #include "matrix.h"
 
+namespace POEMS {
 class SphericalJoint : public Joint  {
   Matrix const_sP;
 public:
@@ -38,5 +39,5 @@ public:
   void ForwardKinematics();
   void BackwardKinematics();
 };
-
+}
 #endif
diff --git a/lib/poems/system.cpp b/lib/poems/system.cpp
index 1ac9ce1236..8515588424 100644
--- a/lib/poems/system.cpp
+++ b/lib/poems/system.cpp
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: system.cpp                                              *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,7 +11,7 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
@@ -36,11 +36,12 @@
 #include "vect3.h"
 #include "virtualmatrix.h"
 
-class Point;
-
+namespace POEMS { class Point; }
+using namespace std;
+using namespace POEMS;
 
 System::System(){
-	mappings = NULL;
+        mappings = NULL;
 }
 
 System::~System(){
@@ -59,7 +60,7 @@ int System::GetNumBodies(){
 
 int * System::GetMappings()
 {
-	return mappings;
+        return mappings;
 }
 
 void System::AddBody(Body* body){
@@ -78,16 +79,16 @@ double System::GetTime(){
   return time;
 }
 
-void System::ComputeForces(){	
-	// NOT DOING ANYTHING AT THIS TIME
-  }  
-  
+void System::ComputeForces(){
+        // NOT DOING ANYTHING AT THIS TIME
+  }
+
 bool System::ReadIn(istream& in){
   int numbodies;
   Body* body;
   int bodytype;
   char bodyname[256];
-  int index;  
+  int index;
 
   // get number of bodies
   in >> numbodies;
@@ -165,12 +166,12 @@ bool System::ReadIn(istream& in){
       delete [] bodyarray;
       return false;
     }
-    
+
     joint->SetBodies(bodyarray[body1], bodyarray[body2]);
 
     bodyarray[body1]->AddJoint(joint);
     bodyarray[body2]->AddJoint(joint);
-    
+
     in >> point1 >> point2;
 
     joint->SetPoints(bodyarray[body1]->GetPoint(point1),bodyarray[body2]->GetPoint(point2));
@@ -203,7 +204,7 @@ void System::WriteOut(ostream& out){
 
     // set the body ID for later identification
     body->SetID(i);
-    
+
     // write out the data
     body->WriteOut(out);
 
@@ -211,7 +212,7 @@ void System::WriteOut(ostream& out){
   }
 
   // number of joints
-  out << joints.GetNumElements() << endl;  
+  out << joints.GetNumElements() << endl;
 
   // joints loop
   i = 0;
@@ -226,7 +227,7 @@ void System::WriteOut(ostream& out){
 
     // write out the data
     joint->WriteOut(out);
-    
+
     i++; j_ele = j_ele->next;
   }
 }
@@ -252,68 +253,68 @@ void System::ClearJointIDs(){
 
 void System::Create_DegenerateSystem(int& nfree, int*freelist, double *&masstotal, double **&inertia, double **&xcm, double **&vcm, double **&omega, double **&ex_space, double **&ey_space, double **&ez_space){
 
- //-------------------------------------------------------------------------//    
-  // Declaring Temporary Entities 
-  //-------------------------------------------------------------------------//   
+ //-------------------------------------------------------------------------//
+  // Declaring Temporary Entities
+  //-------------------------------------------------------------------------//
      Body* body = NULL;
      Body* prev;
      Body* Inertial;
      Point* origin;
      Joint* joint;
-     Point* point_CM; 
-     Point* point_p;  
-     Point* point_k;  
+     Point* point_CM;
+     Point* point_p;
+     Point* point_k;
      Point* point_ch = NULL;
-     Vect3 r1,r2,r3,v1,v2,v3; 
+     Vect3 r1,r2,r3,v1,v2,v3;
      Mat3x3 IM, N, PKCK,PKCN ;
      ColMatrix qo, uo, q, qdot,w;
-     
-	 mappings = new int[nfree];
-	 for(int i = 0; i < nfree; i++)
-	 {
-		 mappings[i] = freelist[i];		 
-	 }
+
+         mappings = new int[nfree];
+         for(int i = 0; i < nfree; i++)
+         {
+                 mappings[i] = freelist[i];
+         }
      qo.Dim(4);
      uo.Dim(3);
      q.Dim(4);
-     qdot.Dim(4);     
+     qdot.Dim(4);
      PKCN.Identity();
      PKCK.Identity();
      w.Dim(3);
-	
-//-------------------------------------------------------------------------//    
-  // Setting up Inertial Frame, gravity and Origin  
-  //-------------------------------------------------------------------------//
-     Inertial= new InertialFrame;  
-     AddBody(Inertial);
-     
-     Vect3 temp1;  
-     temp1.Zeros();  
-     ((InertialFrame*) Inertial)->SetGravity(temp1);    
-     origin= new FixedPoint(temp1);
-     Inertial->AddPoint(origin);    
+
+//-------------------------------------------------------------------------//
+  // Setting up Inertial Frame, gravity and Origin
+  //-------------------------------------------------------------------------//
+     Inertial= new InertialFrame;
+     AddBody(Inertial);
+
+     Vect3 temp1;
+     temp1.Zeros();
+     ((InertialFrame*) Inertial)->SetGravity(temp1);
+     origin= new FixedPoint(temp1);
+     Inertial->AddPoint(origin);
+//-------------------------------------------------------------------------//
+        double ** xh1 = new double*[nfree];
+        double ** xh2 = new double*[nfree];
+
+        for (int i=0; i<nfree; i++){
+                xh1[i] = new double[3];
+                xh2[i] = new double[3];
+        }
+        for (int i=0; i<nfree; i++){
+                for (int j=0; j<3; j++){
+                        xh1[i][j] = xcm[mappings[i]-1][j];
+                }
+        }
+
 //-------------------------------------------------------------------------//
-	double ** xh1 = new double*[nfree];
-	double ** xh2 = new double*[nfree];
-	
-	for (int i=0; i<nfree; i++){		
- 		xh1[i] = new double[3];
-		xh2[i] = new double[3];
-	}	
-	for (int i=0; i<nfree; i++){		
-		for (int j=0; j<3; j++){
-			xh1[i][j] = xcm[mappings[i]-1][j];
-		}
-	}
-	
-//-------------------------------------------------------------------------//    
 // Begin looping over each body for recursive kinematics
 //-------------------------------------------------------------------------//
     for(int i=0;i<nfree;i++){
-	    prev=Inertial;
-	    point_p=origin;               
-      
-     body = new RigidBody;   
+            prev=Inertial;
+            point_p=origin;
+
+     body = new RigidBody;
      body->mass=masstotal[mappings[i]-1];
      IM(1,1)=inertia[mappings[i]-1][0];
      IM(2,2)=inertia[mappings[i]-1][1];
@@ -324,49 +325,49 @@ void System::Create_DegenerateSystem(int& nfree, int*freelist, double *&masstota
      IM(2,1)=IM(1,2);
      IM(3,1)=IM(1,3);
      IM(3,2)=IM(2,3);
-     body->inertia = IM;   
+     body->inertia = IM;
+
+//-------------------------------------------------------//
+
+
+            for (int k=0;k<3;k++){
+          r1(k+1)=xh1[i][k]-xcm[mappings[i]-1][k];
+                         r3(k+1) = xcm[mappings[i]-1][k];
+                         r3(k+1)=xh2[i][k]-xcm[mappings[i]-1][k];
+            }
+
+     r2.Zeros();
+
+     for (int k=1;k<=3;k++){
+          N(k,1)=ex_space[mappings[i]-1][k-1];
+          N(k,2)=ey_space[mappings[i]-1][k-1];
+          N(k,3)=ez_space[mappings[i]-1][k-1];
+         }
 
-//-------------------------------------------------------//                  
-    
-    
-	    for (int k=0;k<3;k++){ 
-          r1(k+1)=xh1[i][k]-xcm[mappings[i]-1][k];     
-			 r3(k+1) = xcm[mappings[i]-1][k];
-			 r3(k+1)=xh2[i][k]-xcm[mappings[i]-1][k]; 	  
-	    }
-	  
-     r2.Zeros(); 
-     
-     for (int k=1;k<=3;k++){     
-          N(k,1)=ex_space[mappings[i]-1][k-1];              
-          N(k,2)=ey_space[mappings[i]-1][k-1];              
-          N(k,3)=ez_space[mappings[i]-1][k-1];              
-         }                  
-     
      PKCK=T(N);
      PKCN=T(N);
-    
-     q.Zeros();     
+
+     q.Zeros();
      EP_FromTransformation(q,N);
-                    
+
      r1=PKCN*r1;
      r3=PKCN*r3;
-     
+
      for (int k=1;k<=3;k++){
-     w(k)=omega[mappings[i]-1][k-1];                                   
-     }     
-     
+     w(k)=omega[mappings[i]-1][k-1];
+     }
+
      Vect3 cart_r, cart_v;
      for (int k=1;k<=3;k++){
-	     cart_r(k)=xcm[mappings[i]-1][k-1];	     
-	     cart_v(k)=vcm[mappings[i]-1][k-1];
-	     }
-	     	     
-	     w=PKCN*w;     
-	     EP_Derivatives(q,w,qdot);
-	     
-     
-//-------------------------------------------------------------------------//    
+             cart_r(k)=xcm[mappings[i]-1][k-1];
+             cart_v(k)=vcm[mappings[i]-1][k-1];
+             }
+
+             w=PKCN*w;
+             EP_Derivatives(q,w,qdot);
+
+
+//-------------------------------------------------------------------------//
 // Create bodies and joints with associated properties for POEMS
 //-------------------------------------------------------------------------//
 
@@ -376,12 +377,12 @@ void System::Create_DegenerateSystem(int& nfree, int*freelist, double *&masstota
      body->AddPoint(point_CM);
      body->AddPoint(point_k);
      body->AddPoint(point_ch);
-     AddBody(body);    
-   
+     AddBody(body);
+
      Mat3x3 One;
-     One.Identity();	
-	  ColMatrix qq=Stack(q,cart_r);
-          ColMatrix vv=Stack(qdot,cart_v);          
+     One.Identity();
+          ColMatrix qq=Stack(q,cart_r);
+          ColMatrix vv=Stack(qdot,cart_v);
           joint=new FreeBodyJoint;
           AddJoint(joint);
           joint->SetBodies(prev,body);
@@ -391,225 +392,225 @@ void System::Create_DegenerateSystem(int& nfree, int*freelist, double *&masstota
           joint->SetZeroOrientation(One);
           joint->DimQandU(7,6);
           joint->SetInitialState(qq,vv);
-          joint->ForwardKinematics();                             
+          joint->ForwardKinematics();
   }
   for(int i = 0; i < nfree; i++) {
-	  delete [] xh1[i];
-	  delete [] xh2[i];
+          delete [] xh1[i];
+          delete [] xh2[i];
   }
   delete [] xh1;
-  delete [] xh2;  
+  delete [] xh2;
 }
 
 
 void System::Create_System_LAMMPS(int numbodies, double *mass,double **inertia, double ** xcm, double ** xjoint,double **vcm,double **omega,double **ex_space, double **ey_space, double **ez_space, int b, int * mapping, int count){
 
-	//-------------------------------------------------------------------------//    
-  // Declaring Temporary Entities 
-	//-------------------------------------------------------------------------//   
-  	
-	Body* body = NULL;
-	Body* prev;
-	Body* Inertial;
-	Point* origin;
-	Joint* joint;
-	Point* point_CM;   
-	Point* point_p;    
-	Point* point_k;    
-	Point* point_ch = NULL;  
-	Vect3 r1,r2,r3,v1,v2,v3; 
-	Mat3x3 IM, N, PKCK,PKCN ;
-	ColMatrix qo, uo, q, qdot,w;
-	Vect3 cart_r, cart_v;
-	mappings = new int[b];
-	for(int i = 0; i < b; i++){
-		mappings[i] = mapping[i];		 
-	}	 
-	 
-     
-	qo.Dim(4);
-	uo.Dim(3);
-	q.Dim(4);
-	qdot.Dim(4);     
-	PKCN.Identity();
-	PKCK.Identity();
-	w.Dim(3);
-	
-	//-------------------------------------------------------------------------//    
-  // Setting up Inertial Frame, gravity and Origin  
-	//-------------------------------------------------------------------------//
-	Inertial= new InertialFrame;  
-	AddBody(Inertial);
-     
-	Vect3 temp1;  
-	temp1.Zeros();  
-	((InertialFrame*) Inertial)->SetGravity(temp1);    
-	origin= new FixedPoint(temp1);
-	Inertial->AddPoint(origin);    
-	//-------------------------------------------------------------------------//
+        //-------------------------------------------------------------------------//
+  // Declaring Temporary Entities
+        //-------------------------------------------------------------------------//
 
-	double ** xh1;
-	double ** xh2;
-		  
-	xh1 = new double*[b];		  
-	xh2 = new double*[b];
-		  
-	
-	for (int i=0; i<b; i++){			
-		xh1[i] = new double[3];
-		xh2[i] = new double[3];			
-	}	
-		  
-		  
-		  
-	for (int j=0; j<3; j++){
-		xh1[0][j] = xcm[mapping[0]-1][j];
-		xh2[b-1][j] = xcm[mapping[b-1]-1][j];
-	}
-	
-	for (int i=0; i<b-1; i++){
-		for (int j=0; j<3; j++){
-			xh1[i+1][j] = xjoint[mapping[i]-count-1][j];
-		}
-	}		  
-	
-	for (int i=0; i<b-1; i++){
-		for (int j=0; j<3; j++){
-			xh2[i][j] = xjoint[mapping[i]-count-1][j];
-		}
-	}
+        Body* body = NULL;
+        Body* prev;
+        Body* Inertial;
+        Point* origin;
+        Joint* joint;
+        Point* point_CM;
+        Point* point_p;
+        Point* point_k;
+        Point* point_ch = NULL;
+        Vect3 r1,r2,r3,v1,v2,v3;
+        Mat3x3 IM, N, PKCK,PKCN ;
+        ColMatrix qo, uo, q, qdot,w;
+        Vect3 cart_r, cart_v;
+        mappings = new int[b];
+        for(int i = 0; i < b; i++){
+                mappings[i] = mapping[i];
+        }
 
-	
-	//-------------------------------------------------------------------------//    
+
+        qo.Dim(4);
+        uo.Dim(3);
+        q.Dim(4);
+        qdot.Dim(4);
+        PKCN.Identity();
+        PKCK.Identity();
+        w.Dim(3);
+
+        //-------------------------------------------------------------------------//
+  // Setting up Inertial Frame, gravity and Origin
+        //-------------------------------------------------------------------------//
+        Inertial= new InertialFrame;
+        AddBody(Inertial);
+
+        Vect3 temp1;
+        temp1.Zeros();
+        ((InertialFrame*) Inertial)->SetGravity(temp1);
+        origin= new FixedPoint(temp1);
+        Inertial->AddPoint(origin);
+        //-------------------------------------------------------------------------//
+
+        double ** xh1;
+        double ** xh2;
+
+        xh1 = new double*[b];
+        xh2 = new double*[b];
+
+
+        for (int i=0; i<b; i++){
+                xh1[i] = new double[3];
+                xh2[i] = new double[3];
+        }
+
+
+
+        for (int j=0; j<3; j++){
+                xh1[0][j] = xcm[mapping[0]-1][j];
+                xh2[b-1][j] = xcm[mapping[b-1]-1][j];
+        }
+
+        for (int i=0; i<b-1; i++){
+                for (int j=0; j<3; j++){
+                        xh1[i+1][j] = xjoint[mapping[i]-count-1][j];
+                }
+        }
+
+        for (int i=0; i<b-1; i++){
+                for (int j=0; j<3; j++){
+                        xh2[i][j] = xjoint[mapping[i]-count-1][j];
+                }
+        }
+
+
+        //-------------------------------------------------------------------------//
 // Begin looping over each body for recursive kinematics
-	//-------------------------------------------------------------------------//
-	for(int i=0;i<b;i++){	    
-		if (i == 0){
-			prev=Inertial;
-			point_p=origin;
-		}    
-		else{
-			prev = body;            
-			point_p = point_ch;     
-		}   
-    
-  	  
-		body = new RigidBody;  
-		body->mass=mass[mapping[i]-1];
-		IM(1,1)=inertia[mapping[i]-1][0];
-		IM(2,2)=inertia[mapping[i]-1][1];
-		IM(3,3)=inertia[mapping[i]-1][2];
-		IM(1,2)=0.0;
-		IM(1,3)=0.0;
-		IM(2,3)=0.0;
-		IM(2,1)=IM(1,2);
-		IM(3,1)=IM(1,3);
-		IM(3,2)=IM(2,3);
-		body->inertia = IM;   
-     
-		//-------------------------------------------------------//
-             
-		for (int k=0;k<3;k++){
-			r1(k+1)=xh1[i][k]-xcm[mapping[i]-1][k]; 
-			r3(k+1)=xh2[i][k]-xcm[mapping[i]-1][k]; 
-		}
-		r2.Zeros();   
-     
-		for (int k=1;k<=3;k++){     
-			N(k,1)=ex_space[mapping[i]-1][k-1];              
-			N(k,2)=ey_space[mapping[i]-1][k-1];              
-			N(k,3)=ez_space[mapping[i]-1][k-1];              
-		}
-	 
-         
-		if (i==0){
-			PKCK=T(N);     
-			PKCN=T(N);     
-	     
-			q.Zeros();	     
-			EP_FromTransformation(q,N);
-	     
-			r1=PKCN*r1;
-			r3=PKCN*r3;
-	     
-			for (int k=1;k<=3;k++){     
-				w(k)=omega[mappings[i]-1][k-1];
-			}	     
-	     
-			for (int k=1;k<=3;k++){
-				cart_r(k)=xcm[mappings[i]-1][k-1];		     
-				cart_v(k)=vcm[mappings[i]-1][k-1];
-			}
-			w=PKCN*w;     
-			EP_Derivatives(q,w,qdot);    
-		 
-		}                 
-		else{	     
-			PKCK=PKCN*N;   
-			PKCN=T(N);     
-    
-			q.Zeros();
-			EP_FromTransformation(q,PKCK);
-               
-			r1=PKCN*r1;
-			r3=PKCN*r3;
-     
-			for (int k=1;k<=3;k++){     
-				w(k)=omega[mapping[i]-1][k-1]-omega[mapping[i-1]-1][k-1];                         
-			}
-    	
-			w=PKCN*w; 
-			EP_Derivatives(q, w, qdot);        
-		}  
-        
-	
-		//-------------------------------------------------------------------------//    
-// Create bodies and joints with associated properties for POEMS
-		//-------------------------------------------------------------------------//
+        //-------------------------------------------------------------------------//
+        for(int i=0;i<b;i++){
+                if (i == 0){
+                        prev=Inertial;
+                        point_p=origin;
+                }
+                else{
+                        prev = body;
+                        point_p = point_ch;
+                }
+
+
+                body = new RigidBody;
+                body->mass=mass[mapping[i]-1];
+                IM(1,1)=inertia[mapping[i]-1][0];
+                IM(2,2)=inertia[mapping[i]-1][1];
+                IM(3,3)=inertia[mapping[i]-1][2];
+                IM(1,2)=0.0;
+                IM(1,3)=0.0;
+                IM(2,3)=0.0;
+                IM(2,1)=IM(1,2);
+                IM(3,1)=IM(1,3);
+                IM(3,2)=IM(2,3);
+                body->inertia = IM;
+
+                //-------------------------------------------------------//
+
+                for (int k=0;k<3;k++){
+                        r1(k+1)=xh1[i][k]-xcm[mapping[i]-1][k];
+                        r3(k+1)=xh2[i][k]-xcm[mapping[i]-1][k];
+                }
+                r2.Zeros();
+
+                for (int k=1;k<=3;k++){
+                        N(k,1)=ex_space[mapping[i]-1][k-1];
+                        N(k,2)=ey_space[mapping[i]-1][k-1];
+                        N(k,3)=ez_space[mapping[i]-1][k-1];
+                }
+
+
+                if (i==0){
+                        PKCK=T(N);
+                        PKCN=T(N);
+
+                        q.Zeros();
+                        EP_FromTransformation(q,N);
+
+                        r1=PKCN*r1;
+                        r3=PKCN*r3;
+
+                        for (int k=1;k<=3;k++){
+                                w(k)=omega[mappings[i]-1][k-1];
+                        }
+
+                        for (int k=1;k<=3;k++){
+                                cart_r(k)=xcm[mappings[i]-1][k-1];
+                                cart_v(k)=vcm[mappings[i]-1][k-1];
+                        }
+                        w=PKCN*w;
+                        EP_Derivatives(q,w,qdot);
+
+                }
+                else{
+                        PKCK=PKCN*N;
+                        PKCN=T(N);
+
+                        q.Zeros();
+                        EP_FromTransformation(q,PKCK);
+
+                        r1=PKCN*r1;
+                        r3=PKCN*r3;
+
+                        for (int k=1;k<=3;k++){
+                                w(k)=omega[mapping[i]-1][k-1]-omega[mapping[i-1]-1][k-1];
+                        }
+
+                        w=PKCN*w;
+                        EP_Derivatives(q, w, qdot);
+                }
+
+
+                //-------------------------------------------------------------------------//
+// Create bodies and joints with associated properties for POEMS
+                //-------------------------------------------------------------------------//
+
+                point_CM = new FixedPoint(r2);
+                point_k = new FixedPoint(r1);
+                point_ch = new FixedPoint(r3);
+                body->AddPoint(point_CM);
+                body->AddPoint(point_k);
+                body->AddPoint(point_ch);
+                AddBody(body);
+
+                Mat3x3 One;
+                One.Identity();
+                if (i==0){
+                        ColMatrix qq=Stack(q,cart_r);
+                        ColMatrix vv=Stack(qdot,cart_v);
+                        joint=new FreeBodyJoint;
+                        AddJoint(joint);
+                        joint->SetBodies(prev,body);
+                        body->AddJoint(joint);
+                        prev->AddJoint(joint);
+                        joint->SetPoints(point_p,point_k);
+                        joint->SetZeroOrientation(One);
+                        joint->DimQandU(7,6);
+                        joint->SetInitialState(qq,vv);
+                        joint->ForwardKinematics();
+                }
+                else{
+                        joint= new SphericalJoint;
+                        AddJoint(joint);
+                        joint->SetBodies(prev,body);
+                        body->AddJoint(joint);
+                        prev->AddJoint(joint);
+                        joint->SetPoints(point_p,point_k);
+                        joint->SetZeroOrientation(One);
+                        joint->DimQandU(4,3);
+                        joint->SetInitialState(q,qdot);
+                        joint->ForwardKinematics();
+                }
+        }
+        for(int i = 0; i < b; i++)
+        {
+                delete [] xh1[i];
+                delete [] xh2[i];
+        }
+        delete [] xh1;
+        delete [] xh2;
 
-		point_CM = new FixedPoint(r2);
-		point_k = new FixedPoint(r1);
-		point_ch = new FixedPoint(r3);
-		body->AddPoint(point_CM);
-		body->AddPoint(point_k);
-		body->AddPoint(point_ch);
-		AddBody(body);    
-   
-		Mat3x3 One;
-		One.Identity();	
-		if (i==0){	  
-			ColMatrix qq=Stack(q,cart_r);
-			ColMatrix vv=Stack(qdot,cart_v);          
-			joint=new FreeBodyJoint;
-			AddJoint(joint);
-			joint->SetBodies(prev,body);
-			body->AddJoint(joint);
-			prev->AddJoint(joint);
-			joint->SetPoints(point_p,point_k);
-			joint->SetZeroOrientation(One);
-			joint->DimQandU(7,6);
-			joint->SetInitialState(qq,vv);  	  
-			joint->ForwardKinematics();  
-		}    
-		else{	     
-			joint= new SphericalJoint;
-			AddJoint(joint);
-			joint->SetBodies(prev,body);
-			body->AddJoint(joint);
-			prev->AddJoint(joint);
-			joint->SetPoints(point_p,point_k);          
-			joint->SetZeroOrientation(One);
-			joint->DimQandU(4,3);
-			joint->SetInitialState(q,qdot);
-			joint->ForwardKinematics();  
-		}
-	}
-	for(int i = 0; i < b; i++)
-	{
-		delete [] xh1[i];
-		delete [] xh2[i];
-	}
-	delete [] xh1;
-	delete [] xh2;
-  
 }
diff --git a/lib/poems/system.h b/lib/poems/system.h
index e5ce27194b..2402e5852f 100644
--- a/lib/poems/system.h
+++ b/lib/poems/system.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: system.h                                                *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,7 +11,7 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
@@ -22,6 +22,7 @@
 #include <iostream>
 #include "poemslist.h"
 
+namespace POEMS {
      class Body;
      class Joint;
 
@@ -33,37 +34,37 @@
        double time;
        List<Body> bodies;
        List<Joint> joints;
-  
+
        System();
        ~System();
        void Delete();
-       
+
        int GetNumBodies();
-       
+
 	   int * GetMappings();
 
        void AddBody(Body* body);
-       
+
        void AddJoint(Joint* joint);
-       
+
        void SetTime(double t);
-       
+
        double GetTime();
-       
+
        void ComputeForces();
-       
+
        bool ReadIn(std::istream& in);
-       
+
        void WriteOut(std::ostream& out);
-       
+
        void ClearBodyIDs();
-       
-       void ClearJointIDs();              
+
+       void ClearJointIDs();
 
        void Create_System_LAMMPS(int numbodies, double *mass,double **inertia, double ** xcm, double ** xjoint,double **vh1,double **omega,double **ex_space, double **ey_space, double **ez_space, int b, int * mapping, int count);
-       
+
        void Create_DegenerateSystem(int& nfree, int*freelist, double *&masstotal, double **&inertia, double **&xcm, double **&vcm, double **&omega, double **&ex_space, double **&ey_space, double **&ez_space);
 
 };
-
+}
 #endif
diff --git a/lib/poems/vect3.cpp b/lib/poems/vect3.cpp
index 7ce3230238..66c2f902d5 100644
--- a/lib/poems/vect3.cpp
+++ b/lib/poems/vect3.cpp
@@ -21,6 +21,8 @@
 #include <iostream>
 
 using namespace std;
+using namespace POEMS;
+
 
 Vect3::Vect3(){
   numrows = 3; numcols = 1;
diff --git a/lib/poems/vect3.h b/lib/poems/vect3.h
index a28832db6c..6b0efc371a 100644
--- a/lib/poems/vect3.h
+++ b/lib/poems/vect3.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: vect3.h                                              *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,7 +11,7 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
@@ -23,6 +23,7 @@
 #include "virtualcolmatrix.h"
 #include "virtualmatrix.h"
 
+namespace POEMS {
 class Matrix;
 class Mat3x3;
 class Mat6x6;
@@ -82,7 +83,6 @@ public:
   friend void FastMult(Mat3x3& A, ColMatrix& B, Vect3& C);
   friend void FastAssign(ColMatrix&A, Vect3& C);
   friend void FastMult(Mat3x3& A, Vect3& B, ColMatrix& C);
-  
 };
-
+}
 #endif
diff --git a/lib/poems/vect4.cpp b/lib/poems/vect4.cpp
index b045998017..8e903d57bd 100644
--- a/lib/poems/vect4.cpp
+++ b/lib/poems/vect4.cpp
@@ -21,6 +21,8 @@
 #include <iostream>
 
 using namespace std;
+using namespace POEMS;
+
 
 Vect4::Vect4(){
   numrows = 4; numcols = 1;
diff --git a/lib/poems/vect4.h b/lib/poems/vect4.h
index 0b28994043..754706f10a 100644
--- a/lib/poems/vect4.h
+++ b/lib/poems/vect4.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: vect4.h                                                 *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,7 +11,7 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
@@ -23,6 +23,7 @@
 #include "virtualcolmatrix.h"
 #include "virtualmatrix.h"
 
+namespace POEMS {
 class Matrix;
 class Mat4x4;
 
@@ -41,7 +42,7 @@ public:
   void BasicSet_1int(int i, double value);
   void BasicIncrement_1int(int i, double value);
 
-  
+
   void Const(double value);
   MatrixType GetType() const;
   std::ostream& WriteData(std::ostream& c) const;
@@ -69,5 +70,5 @@ public:
   friend void FastAdd(Vect4& A, Vect4& B, Vect4& C);
   friend void FastSubt(Vect4& A, Vect4& B, Vect4& C);
 };
-
+}
 #endif
diff --git a/lib/poems/vect6.cpp b/lib/poems/vect6.cpp
index 5baadafc34..5aabe520e1 100644
--- a/lib/poems/vect6.cpp
+++ b/lib/poems/vect6.cpp
@@ -21,6 +21,8 @@
 #include <iostream>
 
 using namespace std;
+using namespace POEMS;
+
 
 Vect6::Vect6(){
   numrows = 6; numcols = 1;
diff --git a/lib/poems/vect6.h b/lib/poems/vect6.h
index 8047c1bbf8..5def4c6b1b 100644
--- a/lib/poems/vect6.h
+++ b/lib/poems/vect6.h
@@ -23,6 +23,7 @@
 #include "virtualcolmatrix.h"
 #include "virtualmatrix.h"
 
+namespace POEMS {
 class Matrix;
 class Mat6x6;
 class ColMatrix;
@@ -69,5 +70,5 @@ public:
 
   friend void OnPopulateSVect(Vect3& angular, Vect3& linear, Vect6& sV);
 };
-
+}
 #endif
diff --git a/lib/poems/virtualcolmatrix.cpp b/lib/poems/virtualcolmatrix.cpp
index 3a6c6e22cd..fa430b81f7 100644
--- a/lib/poems/virtualcolmatrix.cpp
+++ b/lib/poems/virtualcolmatrix.cpp
@@ -20,6 +20,8 @@
 #include <cstdlib>
 
 using namespace std;
+using namespace POEMS;
+
 
 VirtualColMatrix::VirtualColMatrix(){
   numcols = 1;
diff --git a/lib/poems/virtualcolmatrix.h b/lib/poems/virtualcolmatrix.h
index 1548794ca0..040fce051e 100644
--- a/lib/poems/virtualcolmatrix.h
+++ b/lib/poems/virtualcolmatrix.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: virtualcolmatrix.h                                      *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,7 +11,7 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
@@ -21,25 +21,25 @@
 
 #include "virtualmatrix.h"
 
-
+namespace POEMS {
 class VirtualColMatrix : public VirtualMatrix  {
-public: 
-	VirtualColMatrix();
-	~VirtualColMatrix();
-	double& operator_2int (int i, int j); // array access
-	double Get_2int (int i, int j) const;
-	void Set_2int (int i, int j, double value);
-	double BasicGet_2int(int i, int j) const;
-	void BasicSet_2int(int i, int j, double value);
-	void BasicIncrement_2int(int i, int j, double value);
+public:
+        VirtualColMatrix();
+        ~VirtualColMatrix();
+        double& operator_2int (int i, int j); // array access
+        double Get_2int (int i, int j) const;
+        void Set_2int (int i, int j, double value);
+        double BasicGet_2int(int i, int j) const;
+        void BasicSet_2int(int i, int j, double value);
+        void BasicIncrement_2int(int i, int j, double value);
 
-	virtual double& operator_1int (int i) = 0; // array access
-	virtual double Get_1int(int i) const = 0;
-	virtual void Set_1int(int i, double value) = 0;
-	virtual double BasicGet_1int(int i) const = 0;
-	virtual void BasicSet_1int(int i, double value) = 0;
-	virtual void BasicIncrement_1int(int i, double value) = 0;
+        virtual double& operator_1int (int i) = 0; // array access
+        virtual double Get_1int(int i) const = 0;
+        virtual void Set_1int(int i, double value) = 0;
+        virtual double BasicGet_1int(int i) const = 0;
+        virtual void BasicSet_1int(int i, double value) = 0;
+        virtual void BasicIncrement_1int(int i, double value) = 0;
 
 };
-
+}
 #endif
diff --git a/lib/poems/virtualmatrix.cpp b/lib/poems/virtualmatrix.cpp
index 19d5946961..bfffeb3d74 100644
--- a/lib/poems/virtualmatrix.cpp
+++ b/lib/poems/virtualmatrix.cpp
@@ -21,6 +21,8 @@
 #include <cstdlib>
 
 using namespace std;
+using namespace POEMS;
+
 
 VirtualMatrix::VirtualMatrix(){
   numrows = numcols = 0;
@@ -140,14 +142,14 @@ istream& VirtualMatrix::ReadData(istream& c){
 // operators and functions
 //
 
-ostream& operator<< (ostream& c, const VirtualMatrix& A){  //output
+ostream& POEMS::operator<< (ostream& c, const VirtualMatrix& A){  //output
   c << A.GetType() << ' ';
   A.WriteData(c);
   c << endl;
   return c;
 }
 
-istream& operator>> (istream& c, VirtualMatrix& A){  //input
+istream& POEMS::operator>> (istream& c, VirtualMatrix& A){  //input
   VirtualMatrix* vm;
   int matrixtype;
   c >> matrixtype;
diff --git a/lib/poems/virtualmatrix.h b/lib/poems/virtualmatrix.h
index a27c7de047..706b943d66 100644
--- a/lib/poems/virtualmatrix.h
+++ b/lib/poems/virtualmatrix.h
@@ -20,6 +20,7 @@
 #define VIRTUALMATRIX_H
 #include <iostream>
 
+namespace POEMS {
 enum MatrixType {
 	MATRIX = 0,
 	COLMATRIX = 1,
@@ -83,5 +84,5 @@ protected:
 // overloaded operators
 std::ostream& operator<< (std::ostream& c, const VirtualMatrix& A); // output
 std::istream& operator>> (std::istream& c, VirtualMatrix& A); // input
-
+}
 #endif
diff --git a/lib/poems/virtualrowmatrix.cpp b/lib/poems/virtualrowmatrix.cpp
index 1a12cfbfd9..740c747407 100644
--- a/lib/poems/virtualrowmatrix.cpp
+++ b/lib/poems/virtualrowmatrix.cpp
@@ -21,6 +21,8 @@
 #include <cstdlib>
 
 using namespace std;
+using namespace POEMS;
+
 
 VirtualRowMatrix::VirtualRowMatrix(){
   numrows = 1;
diff --git a/lib/poems/virtualrowmatrix.h b/lib/poems/virtualrowmatrix.h
index 68b39f137d..0a3c877f64 100644
--- a/lib/poems/virtualrowmatrix.h
+++ b/lib/poems/virtualrowmatrix.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: virtualrowmatrix.h                                      *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,7 +11,7 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
@@ -21,8 +21,9 @@
 
 #include "virtualmatrix.h"
 
+namespace POEMS {
 class VirtualRowMatrix : public VirtualMatrix  {
-public: 
+public:
   VirtualRowMatrix();
   ~VirtualRowMatrix();
   double& operator_2int (int i, int j); // array access
@@ -31,7 +32,7 @@ public:
   double BasicGet_2int(int i, int j) const;
   void BasicSet_2int(int i, int j, double value);
   void BasicIncrement_2int(int i, int j, double value);
-  
+
   virtual double& operator_1int (int i) = 0; // array access
   virtual double Get_1int(int i) const = 0;
   virtual void Set_1int(int i, double value) = 0;
@@ -39,5 +40,5 @@ public:
   virtual void BasicSet_1int(int i, double value) = 0;
   virtual void BasicIncrement_1int(int i, double value) = 0;
 };
-
+}
 #endif
diff --git a/lib/poems/workspace.cpp b/lib/poems/workspace.cpp
index 29a1c95832..03f18a7895 100644
--- a/lib/poems/workspace.cpp
+++ b/lib/poems/workspace.cpp
@@ -38,6 +38,8 @@
 
 
 using namespace std;
+using namespace POEMS;
+
 
 void Workspace::allocateNewSystem() {
   currentIndex++;
diff --git a/lib/poems/workspace.h b/lib/poems/workspace.h
index 6213815d81..2510b2fc3a 100644
--- a/lib/poems/workspace.h
+++ b/lib/poems/workspace.h
@@ -3,7 +3,7 @@
  *      POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE     *
  *      DESCRIPTION: SEE READ-ME                                           *
  *      FILE NAME: workspace.h                                             *
- *      AUTHORS: See Author List                                           * 
+ *      AUTHORS: See Author List                                           *
  *      GRANTS: See Grants List                                            *
  *      COPYRIGHT: (C) 2005 by Authors as listed in Author's List          *
  *      LICENSE: Please see License Agreement                              *
@@ -11,7 +11,7 @@
  *      ADMINISTRATOR: Prof. Kurt Anderson                                 *
  *                     Computational Dynamics Lab                          *
  *                     Rensselaer Polytechnic Institute                    *
- *                     110 8th St. Troy NY 12180                           * 
+ *                     110 8th St. Troy NY 12180                           *
  *      CONTACT:        anderk5@rpi.edu                                    *
  *_________________________________________________________________________*/
 
@@ -19,61 +19,62 @@
 #ifndef WORKSPACE_H
 #define WORKSPACE_H
 
+namespace POEMS {
 class System;
 class Solver;
 
 struct SysData{
-	System * system;
-	int solver;
-	int integrator;
+        System * system;
+        int solver;
+        int integrator;
 };
 
 class Workspace {
-	SysData * system; // the multibody systems data
-	int currentIndex;
-	int maxAlloc;
-	
+        SysData * system; // the multibody systems data
+        int currentIndex;
+        int maxAlloc;
+
 public:
      Workspace();
      ~Workspace();
-     
+
      double Thalf;
      double Tfull;
      double ConFac;
      double KE_val;
-	  int FirstTime;
-     
+          int FirstTime;
+
      bool LoadFile(char* filename);
-     
+
      bool SaveFile(char* filename, int index = -1);
 
      System* GetSystem(int index = -1);
-     
-     void AddSolver(Solver* s, int index = -1);        
-    
 
-     void LobattoOne(double **&xcm, double **&vcm,double **&omega,double **&torque, double **&fcm, double **&ex_space, double **&ey_space, double **&ez_space);	
-     
-     void LobattoTwo(double **&vcm,double **&omega,double **&torque, double **&fcm);	 
-     
- 
+     void AddSolver(Solver* s, int index = -1);
+
+
+     void LobattoOne(double **&xcm, double **&vcm,double **&omega,double **&torque, double **&fcm, double **&ex_space, double **&ey_space, double **&ez_space);
+
+     void LobattoTwo(double **&vcm,double **&omega,double **&torque, double **&fcm);
+
+
      bool MakeSystem(int& nbody, double *&masstotal, double **&inertia, double **&xcm, double **&vcm, double **&omega, double **&ex_space, double **&ey_space, double **&ez_space, int &njoint, int **&jointbody, double **&xjoint, int& nfree, int*freelist, double dthalf, double dtv, double tempcon, double KE);
-     																																							
-	  
-	  bool SaveSystem(int& nbody, double *&masstotal, double **&inertia, double **&xcm, double **&xjoint, double **&vcm, double **&omega, double **&ex_space, double **&ey_space, double **&ez_space, double **&acm, double **&alpha, double **&torque, double **&fcm, int **&jointbody, int &njoint);
-	  																		
-	 bool MakeDegenerateSystem(int& nfree, int*freelist, double *&masstotal, double **&inertia, double **&xcm, double **&vcm, double **&omega, double **&ex_space, double **&ey_space, double **&ez_space);
+
+
+          bool SaveSystem(int& nbody, double *&masstotal, double **&inertia, double **&xcm, double **&xjoint, double **&vcm, double **&omega, double **&ex_space, double **&ey_space, double **&ez_space, double **&acm, double **&alpha, double **&torque, double **&fcm, int **&jointbody, int &njoint);
+
+         bool MakeDegenerateSystem(int& nfree, int*freelist, double *&masstotal, double **&inertia, double **&xcm, double **&vcm, double **&omega, double **&ex_space, double **&ey_space, double **&ez_space);
      int getNumberOfSystems();
-     
+
      void SetLammpsValues(double dtv, double dthalf, double tempcon);
      void SetKE(int temp, double SysKE);
-	  
-	  void RKStep(double **&xcm, double **&vcm,double **&omega,double **&torque, double **&fcm, double **&ex_space, double **&ey_space, double **&ez_space);	
-	  
-	  void WriteFile(char* filename);
+
+          void RKStep(double **&xcm, double **&vcm,double **&omega,double **&torque, double **&fcm, double **&ex_space, double **&ey_space, double **&ez_space);
+
+          void WriteFile(char* filename);
 
 private:
-	void allocateNewSystem(); //helper function to handle vector resizing and such for the array of system pointers
+        void allocateNewSystem(); //helper function to handle vector resizing and such for the array of system pointers
 };
-
+}
 #endif
diff --git a/src/POEMS/fix_poems.cpp b/src/POEMS/fix_poems.cpp
index 6d943a79db..21c12f700b 100644
--- a/src/POEMS/fix_poems.cpp
+++ b/src/POEMS/fix_poems.cpp
@@ -20,7 +20,6 @@
 #include "fix_poems.h"
 #include <mpi.h>
 #include <cmath>
-#include <cstdio>
 #include <cstring>
 #include <cstdlib>
 #include "workspace.h"
@@ -30,7 +29,6 @@
 #include "respa.h"
 #include "modify.h"
 #include "force.h"
-#include "output.h"
 #include "group.h"
 #include "comm.h"
 #include "citeme.h"
@@ -265,7 +263,7 @@ FixPOEMS::FixPOEMS(LAMMPS *lmp, int narg, char **arg) :
 
   // create POEMS instance
 
-  poems = new Workspace;
+  poems = new POEMS::Workspace;
 
   // compute per body forces and torques inside final_integrate() by default
 
diff --git a/src/POEMS/fix_poems.h b/src/POEMS/fix_poems.h
index 7b2c77a1bd..87f9557ac0 100644
--- a/src/POEMS/fix_poems.h
+++ b/src/POEMS/fix_poems.h
@@ -22,6 +22,7 @@ FixStyle(poems,FixPOEMS)
 
 #include "fix.h"
 
+namespace POEMS { class Workspace; }
 namespace LAMMPS_NS {
 
 class FixPOEMS : public Fix  {
@@ -94,7 +95,7 @@ class FixPOEMS : public Fix  {
 
   // POEMS object
 
-  class Workspace *poems;
+  POEMS::Workspace *poems;
 
   // internal class functions