diff --git a/src/include/xfor/xfor.h b/src/include/xfor/xfor.h
index d4fb92f..59ec7a6 100644
--- a/src/include/xfor/xfor.h
+++ b/src/include/xfor/xfor.h
@@ -36,9 +36,10 @@ namespace MultiArrayHelper
         //virtual size_t rootSteps() const = 0;
         virtual std::shared_ptr<ExtBase> operator+(const ExtBase& in) const = 0;
         virtual std::shared_ptr<ExtBase> operator*(size_t in) const = 0;
-
+	
 	template <class ExtType>
-	const ExtType& expl() const { return dynamic_cast<const ExtT<ExtType>*>(this)->ext(); }
+	const ExtType& expl() const;
+
     };
 
     typedef std::shared_ptr<ExtBase> DExt;
diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt
index 3a338da..42c9eac 100644
--- a/src/tests/CMakeLists.txt
+++ b/src/tests/CMakeLists.txt
@@ -21,6 +21,16 @@ add_dependencies(oputest multiarray)
 target_link_libraries(oputest ${GTEST_BOTH_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT} multiarray)
 add_test(NAME oputest COMMAND oputest)
 
+add_executable(op2utest op2_unit_test.cc)
+add_dependencies(op2utest multiarray)
+target_link_libraries(op2utest ${GTEST_BOTH_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT} multiarray)
+add_test(NAME op2utest COMMAND op2utest)
+
+add_executable(op3utest op3_unit_test.cc)
+add_dependencies(op3utest multiarray)
+target_link_libraries(op3utest ${GTEST_BOTH_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT} multiarray)
+add_test(NAME op3utest COMMAND op3utest)
+
 add_executable(opptest op_perf_test.cc)
 add_dependencies(opptest multiarray)
 target_link_libraries(opptest multiarray)
diff --git a/src/tests/op2_unit_test.cc b/src/tests/op2_unit_test.cc
new file mode 100644
index 0000000..3c044d4
--- /dev/null
+++ b/src/tests/op2_unit_test.cc
@@ -0,0 +1,250 @@
+
+#include "test_header.h"
+
+namespace
+{
+    class OpTest_Sub : public ::testing::Test
+    {
+    protected:
+
+	typedef SingleRangeFactory<char,SpaceType::ANY> SRF;
+	typedef SRF::oType SRange;
+
+        OpTest_Sub()
+	{
+	    swapFactory<SRF>(rfbptr, {'x', 'l'} );
+	    sr1ptr = std::dynamic_pointer_cast<SRange>( rfbptr->create() );
+
+	    swapFactory<SRF>(rfbptr, {'1', '2', '3'} );
+	    sr2ptr = std::dynamic_pointer_cast<SRange>( rfbptr->create() );
+
+	    swapFactory<SRF>(rfbptr, {'a', 'b'} );
+	    sr3ptr = std::dynamic_pointer_cast<SRange>( rfbptr->create() );
+
+        }
+
+	std::shared_ptr<RangeFactoryBase> rfbptr;
+        std::shared_ptr<SRange> sr1ptr;
+	std::shared_ptr<SRange> sr2ptr;
+	std::shared_ptr<SRange> sr3ptr;
+
+        vector<double> v1 = { 2.917, 9.436, 0.373, 0.353, 4.005, 1.070,
+                                   -14.364, -1.868, -25.703, 13.836, 23.563, 41.339 };
+
+        vector<double> v2 = { 0.353, 4.005, 1.070, 2.310, 9.243, 2.911 };
+    };
+
+    
+    class MapTest : public ::testing::Test
+    {
+    protected:
+
+	typedef SingleRangeFactory<size_t,SpaceType::ANY> SRF;
+	typedef SRF::oType SRange;
+        typedef typename SRange::IndexType SIndex;
+
+        std::shared_ptr<SRange> sr1ptr;
+	std::shared_ptr<SRange> sr2ptr;
+        std::shared_ptr<SIndex> si1ptr;
+        std::shared_ptr<SIndex> si2ptr;
+        
+	vector<double> v1 = { -31.71, -77.16, -18.81,
+                              -67.06, 72.31, -54.48,
+                              -50.91, -11.62, -59.57,
+                              -42.53, 80.41, 6.35 };
+
+        typedef std::remove_reference<decltype(*mkMapR( mkMapOp(std::make_shared<plus<size_t>>(),getIndex(sr1ptr),getIndex(sr2ptr)) , sr1ptr, sr2ptr ))>::type MpRange;
+        std::shared_ptr<MpRange> mpr1ptr;
+        typedef MpRange::ORType TRange;
+        
+	MapTest()
+	{
+	    SRF srf1( { 1,  3,  7,  10 } );
+	    SRF srf2( { 2,  6,  8 } );
+
+	    sr1ptr = std::dynamic_pointer_cast<SRange>( srf1.create() );
+	    sr2ptr = std::dynamic_pointer_cast<SRange>( srf2.create() );
+            si1ptr = getIndex(sr1ptr);
+            si2ptr = getIndex(sr2ptr);
+            
+            mpr1ptr = mkMapR( mkMapOp(std::make_shared<plus<size_t>>(),si1ptr,si2ptr) , sr1ptr, sr2ptr );
+	}
+
+    };
+    
+
+    class MetaOp_Test : public ::testing::Test
+    {
+    protected:
+
+	typedef SingleRangeFactory<double,SpaceType::ANY> SRF;
+	typedef SRF::oType SR;
+
+	const size_t s1 = 10;
+	const size_t s2 = 17;
+	
+	MetaOp_Test()
+	{
+	    mv1 = { 2.476, 9.665, 1.289, 2.89, 77.04, -11.09, 100.4, 2.0, -26.5, -0.001 };
+	    mv2 = { 44.56, 23.097, -117.3, -0.0765, 3.445, 0.02389, -4.0112, 10.567, 8.99, -177.2, 475.3,
+		    11.111, 13.108, -35.6, 64.32, 2.44, -12.};
+
+	    assert(mv1.size() == s1); // just to prevent typos...
+	    assert(mv2.size() == s2);
+	    
+	    swapFactory<SRF>(rfbptr, mv1);
+	    sr1ptr = std::dynamic_pointer_cast<SR>(rfbptr->create());
+
+	    swapFactory<SRF>(rfbptr, mv2);
+	    sr2ptr = std::dynamic_pointer_cast<SR>(rfbptr->create());
+	}
+
+	std::shared_ptr<RangeFactoryBase> rfbptr;
+	std::shared_ptr<SR> sr1ptr;
+	std::shared_ptr<SR> sr2ptr;
+
+	vector<double> mv1;
+	vector<double> mv2;
+
+	double vcontract(size_t i)
+	{
+	    double res = 0.;
+	    for(auto& x: mv2){
+		res += pow(mv1[i],x);
+	    }
+	    return res;
+	}
+    };
+
+    TEST_F(MapTest, Exec1)
+    {
+	MultiArray<double,SRange,SRange> ma1(sr1ptr,sr2ptr,v1);
+
+
+	auto ii1 = getIndex( rptr<0>( ma1 ) );
+	auto ii2 = getIndex( rptr<1>( ma1 ) );
+        
+        auto mr = mkMapR( mkMapOp(std::make_shared<plus<size_t>>(),ii1,ii2) , sr1ptr, sr2ptr );
+	//auto mr = mkMapR(std::make_shared<plus<size_t>>(),sr1ptr,sr2ptr);
+	MultiArray<double,MpRange> res(mr);
+	MultiArray<double,MpRange> res2(mr);
+	auto jj = getIndex( mr );
+	(*jj)(ii1,ii2);
+	///auto jj = mkMapI( std::make_shared<plus<size_t> >(), ii1, ii1 );
+		
+	res(jj) += ma1(ii1,ii2);
+	auto mult = mr->explMapMultiplicity();
+        res2(jj) += ma1(ii1,ii2) / staticcast<double>( mult(jj) );
+
+	MultiArray<double,TRange> form = res.format(mpr1ptr->outRange());
+	MultiArray<double,TRange> form2 = res2.format(mpr1ptr->outRange());
+	
+	EXPECT_EQ( jj->range()->outRange()->size(), static_cast<size_t>( 10 ) );
+ 	EXPECT_EQ( jj->range()->mapMultiplicity().at(9), static_cast<size_t>( 3 ) );
+	EXPECT_EQ( jj->range()->mapMultiplicity().at(3), static_cast<size_t>( 1 ) );
+        
+	EXPECT_EQ( form.at(3), -31.71 );
+	EXPECT_EQ( form.at(7), -77.16 );
+	EXPECT_EQ( form.at(9), -18.81 + 72.31 -50.91 );
+	EXPECT_EQ( form.at(5), -67.06 );
+	EXPECT_EQ( form.at(11), -54.48 );
+	EXPECT_EQ( form.at(13), -11.62 );
+	EXPECT_EQ( form.at(15), -59.57 );
+	EXPECT_EQ( form.at(12), -42.53 );
+	EXPECT_EQ( form.at(16), 80.41 );
+	EXPECT_EQ( form.at(18), 6.35 );
+        
+	EXPECT_EQ( form2.at(3), -31.71 );
+	EXPECT_EQ( form2.at(7), -77.16 );
+	EXPECT_EQ( xround( form2.at(9) ), xround( (-18.81 + 72.31 -50.91) / 3. ) );
+	EXPECT_EQ( form2.at(5), -67.06 );
+	EXPECT_EQ( form2.at(11), -54.48 );
+	EXPECT_EQ( form2.at(13), -11.62 );
+	EXPECT_EQ( form2.at(15), -59.57 );
+	EXPECT_EQ( form2.at(12), -42.53 );
+	EXPECT_EQ( form2.at(16), 80.41 );
+	EXPECT_EQ( form2.at(18), 6.35 );
+    }
+    
+    TEST_F(MetaOp_Test, SimpleCall)
+    {
+	FunctionalMultiArray<double,Pow<double>,SR,SR> fma(sr1ptr, sr2ptr);
+
+	auto i = fma.begin();
+	
+	EXPECT_EQ( xround( fma[ i.at( mkt(9.665, -0.0765) ) ] ), xround( pow(9.665, -0.0765) ) );
+    }
+
+    TEST_F(MetaOp_Test, Operation)
+    {
+	FunctionalMultiArray<double,Pow<double>,SR,SR> fma(sr1ptr, sr2ptr);
+	MultiArray<double,SR> res( sr1ptr );
+
+	auto i1 = MAT::getIndex(sr1ptr);
+	auto i2 = MAT::getIndex(sr2ptr);
+	
+	res(i1) += fma(i1,i2).c(i2);
+
+	auto i = res.begin();
+	
+	EXPECT_EQ( xround( res[ i.at( 2.476 ) ] ), xround( vcontract(0) ) );
+	EXPECT_EQ( xround( res[ i.at( 9.665 ) ] ), xround( vcontract(1) ) );
+	EXPECT_EQ( xround( res[ i.at( 1.289 ) ] ), xround( vcontract(2) ) );
+	EXPECT_EQ( xround( res[ i.at( 2.89 ) ] ), xround( vcontract(3) ) );
+	EXPECT_EQ( xround( res[ i.at( 77.04 ) ] ), xround( vcontract(4) ) );
+	
+	EXPECT_EQ( xround( res[ i.at( -11.09 ) ] ), xround( vcontract(5) ) );
+	EXPECT_EQ( xround( res[ i.at( 100.4 ) ] ), xround( vcontract(6) ) );
+	EXPECT_EQ( xround( res[ i.at( 2.0 ) ] ), xround( vcontract(7) ) );
+	EXPECT_EQ( xround( res[ i.at( -26.5 ) ] ), xround( vcontract(8) ) );
+	EXPECT_EQ( xround( res[ i.at( -0.001 ) ] ), xround( vcontract(9) ) );
+    }
+
+    TEST_F(OpTest_Sub, Exec)
+    {
+        MultiArray<double,SRange,SRange,SRange> ma1(sr1ptr, sr2ptr, sr3ptr, v1);
+        MultiArray<double,SRange,SRange> ma2(sr3ptr, sr2ptr, v2);
+
+        SubRangeFactory<SRange> subf(sr2ptr, vector<size_t>({0,2}));
+        auto subptr = MAT::createExplicit(subf);
+        
+        MultiArray<double,SRange,SRange> res(sr3ptr,sr1ptr,0.);
+	MultiArray<double,SRange,SubRange<SRange>,SRange> res2(sr3ptr,subptr,sr1ptr,0.);
+
+        auto i1 = MAT::getIndex( sr1ptr );
+        auto i2 = MAT::getIndex( sr2ptr );
+        auto i3 = MAT::getIndex( sr3ptr );
+        auto si = MAT::getIndex( subptr );
+        (*si)(i2);
+
+        res(i3,i1) += (ma2(i3,i2) - ma1(i1,i2,i3)).c(si);
+	res2(i3,si,i1) = ma2(i3,i2) - ma1(i1,i2,i3);
+
+        EXPECT_EQ( res2.size(), static_cast<size_t>(8) );
+        EXPECT_EQ( res2.vdata().size(), static_cast<size_t>(8) );
+        EXPECT_EQ( MAT::rptr<1>( res2 )->size(), static_cast<size_t>(2) );
+
+        EXPECT_EQ( MAT::rptr<1>( res2 )->isMeta('1'), true );
+        EXPECT_EQ( MAT::rptr<1>( res2 )->isMeta('3'), true );
+        
+        EXPECT_EQ( xround( res.at(mkt('a','x')) ), xround((ma2.at(mkt('a','1')) - ma1.at(mkt('x','1','a'))) + (ma2.at(mkt('a','3')) - ma1.at(mkt('x','3','a')) ) ));
+	EXPECT_EQ( xround( res.at(mkt('a','l')) ), xround((ma2.at(mkt('a','1')) - ma1.at(mkt('l','1','a'))) + (ma2.at(mkt('a','3')) - ma1.at(mkt('l','3','a')) ) ));
+
+        EXPECT_EQ( xround( res.at(mkt('b','x')) ), xround((ma2.at(mkt('b','1')) - ma1.at(mkt('x','1','b'))) + (ma2.at(mkt('b','3')) - ma1.at(mkt('x','3','b')) ) ) );
+        EXPECT_EQ( xround( res.at(mkt('b','l')) ), xround((ma2.at(mkt('b','1')) - ma1.at(mkt('l','1','b'))) + (ma2.at(mkt('b','3')) - ma1.at(mkt('l','3','b')) ) ) );
+
+	EXPECT_EQ( xround( res2.at(mkt('a','1','x')) ), xround(ma2.at(mkt('a','1')) - ma1.at(mkt('x','1','a'))) );
+        EXPECT_EQ( xround( res2.at(mkt('a','1','l')) ), xround(ma2.at(mkt('a','1')) - ma1.at(mkt('l','1','a'))) );
+
+        EXPECT_EQ( xround( res2.at(mkt('b','1','x')) ), xround(ma2.at(mkt('b','1')) - ma1.at(mkt('x','1','b'))) );
+        EXPECT_EQ( xround( res2.at(mkt('b','1','l')) ), xround(ma2.at(mkt('b','1')) - ma1.at(mkt('l','1','b'))) );
+
+        EXPECT_EQ( xround( res2.at(mkt('a','3','x')) ), xround(ma2.at(mkt('a','3')) - ma1.at(mkt('x','3','a'))) );
+        EXPECT_EQ( xround( res2.at(mkt('a','3','l')) ), xround(ma2.at(mkt('a','3')) - ma1.at(mkt('l','3','a'))) );
+
+        EXPECT_EQ( xround( res2.at(mkt('b','3','x')) ), xround(ma2.at(mkt('b','3')) - ma1.at(mkt('x','3','b'))) );
+        EXPECT_EQ( xround( res2.at(mkt('b','3','l')) ), xround(ma2.at(mkt('b','3')) - ma1.at(mkt('l','3','b'))) );
+
+    }
+
+}
diff --git a/src/tests/op3_unit_test.cc b/src/tests/op3_unit_test.cc
new file mode 100644
index 0000000..876f9fd
--- /dev/null
+++ b/src/tests/op3_unit_test.cc
@@ -0,0 +1,118 @@
+
+#include "test_header.h"
+
+namespace
+{
+    class OpTest_Spin : public ::testing::Test
+    {
+    protected:
+
+	typedef SpinRF SRF;
+	typedef SpinRange SR;
+	typedef MultiRangeFactory<SR,SR,SR,SR,SR,SR,SR,SR> SR8F;
+	typedef SR8F::oType SR8;
+
+	static const size_t s = 65536;
+	
+	OpTest_Spin()
+	{
+	    data.resize(s);
+	    for(size_t i = 0; i != s; ++i){
+		const double arg = static_cast<double>( s - i ) - 0.1; 
+		data[i] = sin(arg)*10;
+	    }
+	    SRF f;
+	    sr = std::dynamic_pointer_cast<SR>(f.create());
+	}
+
+	vector<double> data;
+	std::shared_ptr<SR> sr;
+    };
+
+    TEST_F(OpTest_Spin, Contract)
+    {
+	MultiArray<double,SR,SR,SR,SR,SR,SR,SR,SR> ma(sr, sr, sr, sr, sr, sr, sr, sr, data);
+	MultiArray<double,SR,SR> res1( sr, sr );
+	MultiArray<double,SR,SR> res2( sr, sr );
+	
+	auto alpha = MAT::getIndex<SR>();
+	auto beta = MAT::getIndex<SR>();
+	auto gamma = MAT::getIndex<SR>();
+	auto delta = MAT::getIndex<SR>();
+	auto deltap = MAT::getIndex<SR>();
+	
+	auto mix = MAT::mkMIndex( alpha, beta, gamma );
+
+	std::clock_t begin = std::clock();
+	res1(delta, deltap) += ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
+	std::clock_t end = std::clock();
+	std::cout << "MultiArray time: " << static_cast<double>( end - begin ) / CLOCKS_PER_SEC
+		  << std::endl;
+
+	res2(delta, deltap).par() += ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(alpha).c(beta).c(gamma);
+	
+	vector<double> vres(4*4);
+
+	std::clock_t begin2 = std::clock();
+	for(size_t d = 0; d != 4; ++d){
+	    for(size_t p = 0; p != 4; ++p){
+		const size_t tidx = d*4 + p;
+		vres[tidx] = 0.;
+		for(size_t a = 0; a != 4; ++a){
+		    for(size_t b = 0; b != 4; ++b){
+			for(size_t c = 0; c != 4; ++c){
+			    const size_t sidx = d*4*4*4*4*4*4*4 + a*5*4*4*4*4*4 + b*5*4*4*4 + c*5*4 + p;
+			    vres[tidx] += data[sidx];
+			}
+		    }
+		}
+	    }
+	}
+	std::clock_t end2 = std::clock();
+
+	EXPECT_EQ( xround(res1.at(mkts(0,0))), xround(vres[0]) );
+	EXPECT_EQ( xround(res1.at(mkts(0,1))), xround(vres[1]) );
+	EXPECT_EQ( xround(res1.at(mkts(0,2))), xround(vres[2]) );
+	EXPECT_EQ( xround(res1.at(mkts(0,3))), xround(vres[3]) );
+	
+	EXPECT_EQ( xround(res1.at(mkts(1,0))), xround(vres[4]) );
+	EXPECT_EQ( xround(res1.at(mkts(1,1))), xround(vres[5]) );
+	EXPECT_EQ( xround(res1.at(mkts(1,2))), xround(vres[6]) );
+	EXPECT_EQ( xround(res1.at(mkts(1,3))), xround(vres[7]) );
+
+	EXPECT_EQ( xround(res1.at(mkts(2,0))), xround(vres[8]) );
+	EXPECT_EQ( xround(res1.at(mkts(2,1))), xround(vres[9]) );
+	EXPECT_EQ( xround(res1.at(mkts(2,2))), xround(vres[10]) );
+	EXPECT_EQ( xround(res1.at(mkts(2,3))), xround(vres[11]) );
+
+	EXPECT_EQ( xround(res1.at(mkts(3,0))), xround(vres[12]) );
+	EXPECT_EQ( xround(res1.at(mkts(3,1))), xround(vres[13]) );
+	EXPECT_EQ( xround(res1.at(mkts(3,2))), xround(vres[14]) );
+	EXPECT_EQ( xround(res1.at(mkts(3,3))), xround(vres[15]) );
+
+	EXPECT_EQ( xround(res2.at(mkts(0,0))), xround(vres[0]) );
+	EXPECT_EQ( xround(res2.at(mkts(0,1))), xround(vres[1]) );
+	EXPECT_EQ( xround(res2.at(mkts(0,2))), xround(vres[2]) );
+	EXPECT_EQ( xround(res2.at(mkts(0,3))), xround(vres[3]) );
+	
+	EXPECT_EQ( xround(res2.at(mkts(1,0))), xround(vres[4]) );
+	EXPECT_EQ( xround(res2.at(mkts(1,1))), xround(vres[5]) );
+	EXPECT_EQ( xround(res2.at(mkts(1,2))), xround(vres[6]) );
+	EXPECT_EQ( xround(res2.at(mkts(1,3))), xround(vres[7]) );
+
+	EXPECT_EQ( xround(res2.at(mkts(2,0))), xround(vres[8]) );
+	EXPECT_EQ( xround(res2.at(mkts(2,1))), xround(vres[9]) );
+	EXPECT_EQ( xround(res2.at(mkts(2,2))), xround(vres[10]) );
+	EXPECT_EQ( xround(res2.at(mkts(2,3))), xround(vres[11]) );
+
+	EXPECT_EQ( xround(res2.at(mkts(3,0))), xround(vres[12]) );
+	EXPECT_EQ( xround(res2.at(mkts(3,1))), xround(vres[13]) );
+	EXPECT_EQ( xround(res2.at(mkts(3,2))), xround(vres[14]) );
+	EXPECT_EQ( xround(res2.at(mkts(3,3))), xround(vres[15]) );
+
+	std::cout << "vector - for loop time: " << static_cast<double>( end2 - begin2 ) / CLOCKS_PER_SEC
+		  << std::endl;
+	std::cout << "ratio: " << static_cast<double>( end - begin ) / static_cast<double>( end2 - begin2 ) << std::endl;
+    }    
+
+}
diff --git a/src/tests/op_unit_test.cc b/src/tests/op_unit_test.cc
index 4440366..c45891b 100644
--- a/src/tests/op_unit_test.cc
+++ b/src/tests/op_unit_test.cc
@@ -1,95 +1,9 @@
 // -*- C++ -*-
 
-#include <cstdlib>
-#include "gtest/gtest.h"
-#include <iostream>
-
-#include "multi_array_header.h"
-
-#include <ctime>
-#include <cmath>
-
-namespace MAT = MultiArrayTools;
-
-namespace {
-
-    double xround(double arg)
-    {
-	if(std::isnan(arg)) { return 0.; }
-	return roundf(arg * 100000.) / 100000.;
-    }
-    
-    using namespace MAT;
-
-    template <class Factory>
-    void swapFactory(std::shared_ptr<RangeFactoryBase>& fptr)
-    {
-	auto nptr = std::make_shared<Factory>();
-	fptr = nptr;
-    }
-    
-    template <class Factory, typename T>
-    void swapFactory(std::shared_ptr<RangeFactoryBase>& fptr, std::initializer_list<T> ilist)
-    {
-	vector<T> tmp = ilist;
-	auto nptr = std::make_shared<Factory>( tmp );
-	fptr = nptr;
-    }
-
-    template <class Factory, typename T>
-    void swapFactory(std::shared_ptr<RangeFactoryBase>& fptr, vector<T>& ilist)
-    {
-	vector<T> tmp = ilist;
-	auto nptr = std::make_shared<Factory>( tmp );
-	fptr = nptr;
-    }
-
-    
-    template <class Factory, class... Rs>
-    void swapMFactory(std::shared_ptr<RangeFactoryBase>& fptr, const Rs&... rs)
-    {
-	auto nptr = std::make_shared<Factory>( rs... );
-	fptr = nptr;
-    }
-
-    template <typename... Ts>
-    auto mkt(Ts&&... ts) -> decltype(std::make_tuple(ts...))
-    {
-	return std::make_tuple(ts...);
-    }
-
-    template <typename... Ts>
-    auto mkts(Ts&&... ts) -> decltype(std::make_tuple(ts...))
-    {
-	return std::make_tuple(static_cast<size_t>( ts )...);
-    }
-
-    //typedef Expressions1 EC1;
-    
-    template <typename T>
-    struct Pow
-    {
-	static constexpr bool FISSTATIC = true;
-	typedef T value_type;
-	
-	static inline T apply(T b, T e)
-	{
-	    return pow(b, e);
-	}
-
-	template <class... Ops>
-	static auto mk(const Ops&... ops)
-	    -> Operation<T,Pow<T>,Ops...>
-	{
-	    return Operation<T,Pow<T>,Ops...>(ops...);
-	}
-
-	static inline T apply(const std::tuple<double, double>& arg)
-	{
-	    return pow(std::get<0>(arg), std::get<1>(arg));
-	}
-    };
+#include "test_header.h"
 
+namespace
+{
     
     class OpTest_1Dim : public ::testing::Test
     {
@@ -159,376 +73,7 @@ namespace {
 				   -14.364, -1.868, -25.703, 13.836, 23.563, 41.339 };
     };
 
-    class OpTest_Sub : public ::testing::Test
-    {
-    protected:
 
-	typedef SingleRangeFactory<char,SpaceType::ANY> SRF;
-	typedef SRF::oType SRange;
-
-        OpTest_Sub()
-	{
-	    swapFactory<SRF>(rfbptr, {'x', 'l'} );
-	    sr1ptr = std::dynamic_pointer_cast<SRange>( rfbptr->create() );
-
-	    swapFactory<SRF>(rfbptr, {'1', '2', '3'} );
-	    sr2ptr = std::dynamic_pointer_cast<SRange>( rfbptr->create() );
-
-	    swapFactory<SRF>(rfbptr, {'a', 'b'} );
-	    sr3ptr = std::dynamic_pointer_cast<SRange>( rfbptr->create() );
-
-        }
-
-	std::shared_ptr<RangeFactoryBase> rfbptr;
-        std::shared_ptr<SRange> sr1ptr;
-	std::shared_ptr<SRange> sr2ptr;
-	std::shared_ptr<SRange> sr3ptr;
-
-        vector<double> v1 = { 2.917, 9.436, 0.373, 0.353, 4.005, 1.070,
-                                   -14.364, -1.868, -25.703, 13.836, 23.563, 41.339 };
-
-        vector<double> v2 = { 0.353, 4.005, 1.070, 2.310, 9.243, 2.911 };
-    };
-
-    
-    class MapTest : public ::testing::Test
-    {
-    protected:
-
-	typedef SingleRangeFactory<size_t,SpaceType::ANY> SRF;
-	typedef SRF::oType SRange;
-        typedef typename SRange::IndexType SIndex;
-        /*
-	typedef FunctionalMultiArray<size_t,plus<size_t>,SRange,SRange> MapF;
-	
-	typedef MapRangeFactory<MapF,SRange,SRange> MpRF;
-	typedef MpRF::oType MpRange;
-        */
-	
-
-        std::shared_ptr<SRange> sr1ptr;
-	std::shared_ptr<SRange> sr2ptr;
-        std::shared_ptr<SIndex> si1ptr;
-        std::shared_ptr<SIndex> si2ptr;
-        
-	//std::shared_ptr<MpRange> mpr1ptr;
-	vector<double> v1 = { -31.71, -77.16, -18.81,
-                              -67.06, 72.31, -54.48,
-                              -50.91, -11.62, -59.57,
-                              -42.53, 80.41, 6.35 };
-
-        typedef std::remove_reference<decltype(*mkMapR( mkMapOp(std::make_shared<plus<size_t>>(),getIndex(sr1ptr),getIndex(sr2ptr)) , sr1ptr, sr2ptr ))>::type MpRange;
-        std::shared_ptr<MpRange> mpr1ptr;
-        typedef MpRange::ORType TRange;
-        
-	MapTest()
-	{
-	    SRF srf1( { 1,  3,  7,  10 } );
-	    SRF srf2( { 2,  6,  8 } );
-
-	    sr1ptr = std::dynamic_pointer_cast<SRange>( srf1.create() );
-	    sr2ptr = std::dynamic_pointer_cast<SRange>( srf2.create() );
-            si1ptr = getIndex(sr1ptr);
-            si2ptr = getIndex(sr2ptr);
-            
-            mpr1ptr = mkMapR( mkMapOp(std::make_shared<plus<size_t>>(),si1ptr,si2ptr) , sr1ptr, sr2ptr );
-	}
-
-    };
-    
-    class OpTest_Performance : public ::testing::Test
-    {
-    protected:
-
-	typedef SingleRangeFactory<size_t,SpaceType::ANY> SRF;
-	typedef SRF::oType SRange;
-
-	typedef MultiRangeFactory<SRange,SRange> MRF;
-	typedef MRF::oType MRange;
-
-	OpTest_Performance()
-	{
-	
-	    vector<size_t> initvec1(vs1);
-	    cv1.resize(vs1);
-	    for(size_t i = 0; i != vs1; ++i){
-		initvec1[i] = i;
-		cv1[i] = sqrt( static_cast<double>(i)*0.53 );
-	    }
-
-	    vector<size_t> initvec2(vs2);
-	    cv2.resize(vs2*vs1);
-	    for(size_t i = 0; i != vs2; ++i){
-		initvec2[i] = i;
-		for(size_t j = 0; j != vs1; ++j){
-		    cv2[i*vs1 + j] = static_cast<double>(i) * sin(static_cast<double>(j)*0.4);
-		}
-	    }
-
-	    swapFactory<SRF>(rfbptr, initvec1);
-	    sr1ptr = std::dynamic_pointer_cast<SRange>(rfbptr->create());
-
-	    swapFactory<SRF>(rfbptr, initvec2);
-	    sr2ptr = std::dynamic_pointer_cast<SRange>(rfbptr->create());
-
-	    swapMFactory<MRF>(rfbptr, sr2ptr, sr1ptr);
-	    mrptr = std::dynamic_pointer_cast<MRange>(rfbptr->create());
-	}
-
-	//const size_t vs1 = 10000;
-	//const size_t vs2 = 1000;
-	const size_t vs1 = 4000;
-	const size_t vs2 = 2500;
-
-	std::shared_ptr<RangeFactoryBase> rfbptr;
-	std::shared_ptr<SRange> sr1ptr;
-	std::shared_ptr<SRange> sr2ptr;
-	std::shared_ptr<MRange> mrptr;
-
-	vector<double> cv1;
-	vector<double> cv2;
-    };
-
-    class OpTest_Spin : public ::testing::Test
-    {
-    protected:
-
-	typedef SpinRF SRF;
-	typedef SpinRange SR;
-	typedef MultiRangeFactory<SR,SR,SR,SR,SR,SR,SR,SR> SR8F;
-	typedef SR8F::oType SR8;
-
-	static const size_t s = 65536;
-	
-	OpTest_Spin()
-	{
-	    data.resize(s);
-	    for(size_t i = 0; i != s; ++i){
-		const double arg = static_cast<double>( s - i ) - 0.1; 
-		data[i] = sin(arg)*10;
-	    }
-	    SRF f;
-	    sr = std::dynamic_pointer_cast<SR>(f.create());
-	}
-
-	vector<double> data;
-	std::shared_ptr<SR> sr;
-    };
-
-    class MetaOp_Test : public ::testing::Test
-    {
-    protected:
-
-	typedef SingleRangeFactory<double,SpaceType::ANY> SRF;
-	typedef SRF::oType SR;
-
-	const size_t s1 = 10;
-	const size_t s2 = 17;
-	
-	MetaOp_Test()
-	{
-	    mv1 = { 2.476, 9.665, 1.289, 2.89, 77.04, -11.09, 100.4, 2.0, -26.5, -0.001 };
-	    mv2 = { 44.56, 23.097, -117.3, -0.0765, 3.445, 0.02389, -4.0112, 10.567, 8.99, -177.2, 475.3,
-		    11.111, 13.108, -35.6, 64.32, 2.44, -12.};
-
-	    assert(mv1.size() == s1); // just to prevent typos...
-	    assert(mv2.size() == s2);
-	    
-	    swapFactory<SRF>(rfbptr, mv1);
-	    sr1ptr = std::dynamic_pointer_cast<SR>(rfbptr->create());
-
-	    swapFactory<SRF>(rfbptr, mv2);
-	    sr2ptr = std::dynamic_pointer_cast<SR>(rfbptr->create());
-	}
-
-	std::shared_ptr<RangeFactoryBase> rfbptr;
-	std::shared_ptr<SR> sr1ptr;
-	std::shared_ptr<SR> sr2ptr;
-
-	vector<double> mv1;
-	vector<double> mv2;
-
-	double vcontract(size_t i)
-	{
-	    double res = 0.;
-	    for(auto& x: mv2){
-		res += pow(mv1[i],x);
-	    }
-	    return res;
-	}
-    };
-
-
-    TEST_F(MapTest, Exec1)
-    {
-	MultiArray<double,SRange,SRange> ma1(sr1ptr,sr2ptr,v1);
-
-
-	auto ii1 = getIndex( rptr<0>( ma1 ) );
-	auto ii2 = getIndex( rptr<1>( ma1 ) );
-        
-        auto mr = mkMapR( mkMapOp(std::make_shared<plus<size_t>>(),ii1,ii2) , sr1ptr, sr2ptr );
-	//auto mr = mkMapR(std::make_shared<plus<size_t>>(),sr1ptr,sr2ptr);
-	MultiArray<double,MpRange> res(mr);
-	MultiArray<double,MpRange> res2(mr);
-	auto jj = getIndex( mr );
-	(*jj)(ii1,ii2);
-	///auto jj = mkMapI( std::make_shared<plus<size_t> >(), ii1, ii1 );
-		
-	res(jj) += ma1(ii1,ii2);
-	auto mult = mr->explMapMultiplicity();
-        res2(jj) += ma1(ii1,ii2) / staticcast<double>( mult(jj) );
-
-	MultiArray<double,TRange> form = res.format(mpr1ptr->outRange());
-	MultiArray<double,TRange> form2 = res2.format(mpr1ptr->outRange());
-	
-	EXPECT_EQ( jj->range()->outRange()->size(), static_cast<size_t>( 10 ) );
- 	EXPECT_EQ( jj->range()->mapMultiplicity().at(9), static_cast<size_t>( 3 ) );
-	EXPECT_EQ( jj->range()->mapMultiplicity().at(3), static_cast<size_t>( 1 ) );
-        
-	EXPECT_EQ( form.at(3), -31.71 );
-	EXPECT_EQ( form.at(7), -77.16 );
-	EXPECT_EQ( form.at(9), -18.81 + 72.31 -50.91 );
-	EXPECT_EQ( form.at(5), -67.06 );
-	EXPECT_EQ( form.at(11), -54.48 );
-	EXPECT_EQ( form.at(13), -11.62 );
-	EXPECT_EQ( form.at(15), -59.57 );
-	EXPECT_EQ( form.at(12), -42.53 );
-	EXPECT_EQ( form.at(16), 80.41 );
-	EXPECT_EQ( form.at(18), 6.35 );
-        
-	EXPECT_EQ( form2.at(3), -31.71 );
-	EXPECT_EQ( form2.at(7), -77.16 );
-	EXPECT_EQ( xround( form2.at(9) ), xround( (-18.81 + 72.31 -50.91) / 3. ) );
-	EXPECT_EQ( form2.at(5), -67.06 );
-	EXPECT_EQ( form2.at(11), -54.48 );
-	EXPECT_EQ( form2.at(13), -11.62 );
-	EXPECT_EQ( form2.at(15), -59.57 );
-	EXPECT_EQ( form2.at(12), -42.53 );
-	EXPECT_EQ( form2.at(16), 80.41 );
-	EXPECT_EQ( form2.at(18), 6.35 );
-    }
-    
-    TEST_F(MetaOp_Test, SimpleCall)
-    {
-	FunctionalMultiArray<double,Pow<double>,SR,SR> fma(sr1ptr, sr2ptr);
-
-	auto i = fma.begin();
-	
-	EXPECT_EQ( xround( fma[ i.at( mkt(9.665, -0.0765) ) ] ), xround( pow(9.665, -0.0765) ) );
-    }
-
-    TEST_F(MetaOp_Test, Operation)
-    {
-	FunctionalMultiArray<double,Pow<double>,SR,SR> fma(sr1ptr, sr2ptr);
-	MultiArray<double,SR> res( sr1ptr );
-
-	auto i1 = MAT::getIndex(sr1ptr);
-	auto i2 = MAT::getIndex(sr2ptr);
-	
-	res(i1) += fma(i1,i2).c(i2);
-
-	auto i = res.begin();
-	
-	EXPECT_EQ( xround( res[ i.at( 2.476 ) ] ), xround( vcontract(0) ) );
-	EXPECT_EQ( xround( res[ i.at( 9.665 ) ] ), xround( vcontract(1) ) );
-	EXPECT_EQ( xround( res[ i.at( 1.289 ) ] ), xround( vcontract(2) ) );
-	EXPECT_EQ( xround( res[ i.at( 2.89 ) ] ), xround( vcontract(3) ) );
-	EXPECT_EQ( xround( res[ i.at( 77.04 ) ] ), xround( vcontract(4) ) );
-	
-	EXPECT_EQ( xround( res[ i.at( -11.09 ) ] ), xround( vcontract(5) ) );
-	EXPECT_EQ( xround( res[ i.at( 100.4 ) ] ), xround( vcontract(6) ) );
-	EXPECT_EQ( xround( res[ i.at( 2.0 ) ] ), xround( vcontract(7) ) );
-	EXPECT_EQ( xround( res[ i.at( -26.5 ) ] ), xround( vcontract(8) ) );
-	EXPECT_EQ( xround( res[ i.at( -0.001 ) ] ), xround( vcontract(9) ) );
-    }
-
-    
-    TEST_F(OpTest_Spin, Contract)
-    {
-	MultiArray<double,SR,SR,SR,SR,SR,SR,SR,SR> ma(sr, sr, sr, sr, sr, sr, sr, sr, data);
-	MultiArray<double,SR,SR> res1( sr, sr );
-	MultiArray<double,SR,SR> res2( sr, sr );
-	
-	auto alpha = MAT::getIndex<SR>();
-	auto beta = MAT::getIndex<SR>();
-	auto gamma = MAT::getIndex<SR>();
-	auto delta = MAT::getIndex<SR>();
-	auto deltap = MAT::getIndex<SR>();
-	
-	auto mix = MAT::mkMIndex( alpha, beta, gamma );
-
-	std::clock_t begin = std::clock();
-	res1(delta, deltap) += ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
-	std::clock_t end = std::clock();
-	std::cout << "MultiArray time: " << static_cast<double>( end - begin ) / CLOCKS_PER_SEC
-		  << std::endl;
-
-	res2(delta, deltap).par() += ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(alpha).c(beta).c(gamma);
-	
-	vector<double> vres(4*4);
-
-	std::clock_t begin2 = std::clock();
-	for(size_t d = 0; d != 4; ++d){
-	    for(size_t p = 0; p != 4; ++p){
-		const size_t tidx = d*4 + p;
-		vres[tidx] = 0.;
-		for(size_t a = 0; a != 4; ++a){
-		    for(size_t b = 0; b != 4; ++b){
-			for(size_t c = 0; c != 4; ++c){
-			    const size_t sidx = d*4*4*4*4*4*4*4 + a*5*4*4*4*4*4 + b*5*4*4*4 + c*5*4 + p;
-			    vres[tidx] += data[sidx];
-			}
-		    }
-		}
-	    }
-	}
-	std::clock_t end2 = std::clock();
-
-	EXPECT_EQ( xround(res1.at(mkts(0,0))), xround(vres[0]) );
-	EXPECT_EQ( xround(res1.at(mkts(0,1))), xround(vres[1]) );
-	EXPECT_EQ( xround(res1.at(mkts(0,2))), xround(vres[2]) );
-	EXPECT_EQ( xround(res1.at(mkts(0,3))), xround(vres[3]) );
-	
-	EXPECT_EQ( xround(res1.at(mkts(1,0))), xround(vres[4]) );
-	EXPECT_EQ( xround(res1.at(mkts(1,1))), xround(vres[5]) );
-	EXPECT_EQ( xround(res1.at(mkts(1,2))), xround(vres[6]) );
-	EXPECT_EQ( xround(res1.at(mkts(1,3))), xround(vres[7]) );
-
-	EXPECT_EQ( xround(res1.at(mkts(2,0))), xround(vres[8]) );
-	EXPECT_EQ( xround(res1.at(mkts(2,1))), xround(vres[9]) );
-	EXPECT_EQ( xround(res1.at(mkts(2,2))), xround(vres[10]) );
-	EXPECT_EQ( xround(res1.at(mkts(2,3))), xround(vres[11]) );
-
-	EXPECT_EQ( xround(res1.at(mkts(3,0))), xround(vres[12]) );
-	EXPECT_EQ( xround(res1.at(mkts(3,1))), xround(vres[13]) );
-	EXPECT_EQ( xround(res1.at(mkts(3,2))), xround(vres[14]) );
-	EXPECT_EQ( xround(res1.at(mkts(3,3))), xround(vres[15]) );
-
-	EXPECT_EQ( xround(res2.at(mkts(0,0))), xround(vres[0]) );
-	EXPECT_EQ( xround(res2.at(mkts(0,1))), xround(vres[1]) );
-	EXPECT_EQ( xround(res2.at(mkts(0,2))), xround(vres[2]) );
-	EXPECT_EQ( xround(res2.at(mkts(0,3))), xround(vres[3]) );
-	
-	EXPECT_EQ( xround(res2.at(mkts(1,0))), xround(vres[4]) );
-	EXPECT_EQ( xround(res2.at(mkts(1,1))), xround(vres[5]) );
-	EXPECT_EQ( xround(res2.at(mkts(1,2))), xround(vres[6]) );
-	EXPECT_EQ( xround(res2.at(mkts(1,3))), xround(vres[7]) );
-
-	EXPECT_EQ( xround(res2.at(mkts(2,0))), xround(vres[8]) );
-	EXPECT_EQ( xround(res2.at(mkts(2,1))), xround(vres[9]) );
-	EXPECT_EQ( xround(res2.at(mkts(2,2))), xround(vres[10]) );
-	EXPECT_EQ( xround(res2.at(mkts(2,3))), xround(vres[11]) );
-
-	EXPECT_EQ( xround(res2.at(mkts(3,0))), xround(vres[12]) );
-	EXPECT_EQ( xround(res2.at(mkts(3,1))), xround(vres[13]) );
-	EXPECT_EQ( xround(res2.at(mkts(3,2))), xround(vres[14]) );
-	EXPECT_EQ( xround(res2.at(mkts(3,3))), xround(vres[15]) );
-
-	std::cout << "vector - for loop time: " << static_cast<double>( end2 - begin2 ) / CLOCKS_PER_SEC
-		  << std::endl;
-	std::cout << "ratio: " << static_cast<double>( end - begin ) / static_cast<double>( end2 - begin2 ) << std::endl;
-    }    
 
     TEST_F(OpTest_1Dim, ExecOp)
     {
@@ -546,52 +91,6 @@ namespace {
 	EXPECT_EQ( fabs( res.at('g') - (7.192+5.063) ) < 0.0001, true );
     }
 
-    TEST_F(OpTest_Sub, Exec)
-    {
-        MultiArray<double,SRange,SRange,SRange> ma1(sr1ptr, sr2ptr, sr3ptr, v1);
-        MultiArray<double,SRange,SRange> ma2(sr3ptr, sr2ptr, v2);
-
-        SubRangeFactory<SRange> subf(sr2ptr, vector<size_t>({0,2}));
-        auto subptr = MAT::createExplicit(subf);
-        
-        MultiArray<double,SRange,SRange> res(sr3ptr,sr1ptr,0.);
-	MultiArray<double,SRange,SubRange<SRange>,SRange> res2(sr3ptr,subptr,sr1ptr,0.);
-
-        auto i1 = MAT::getIndex( sr1ptr );
-        auto i2 = MAT::getIndex( sr2ptr );
-        auto i3 = MAT::getIndex( sr3ptr );
-        auto si = MAT::getIndex( subptr );
-        (*si)(i2);
-
-        res(i3,i1) += (ma2(i3,i2) - ma1(i1,i2,i3)).c(si);
-	res2(i3,si,i1) = ma2(i3,i2) - ma1(i1,i2,i3);
-
-        EXPECT_EQ( res2.size(), static_cast<size_t>(8) );
-        EXPECT_EQ( res2.vdata().size(), static_cast<size_t>(8) );
-        EXPECT_EQ( MAT::rptr<1>( res2 )->size(), static_cast<size_t>(2) );
-
-        EXPECT_EQ( MAT::rptr<1>( res2 )->isMeta('1'), true );
-        EXPECT_EQ( MAT::rptr<1>( res2 )->isMeta('3'), true );
-        
-        EXPECT_EQ( xround( res.at(mkt('a','x')) ), xround((ma2.at(mkt('a','1')) - ma1.at(mkt('x','1','a'))) + (ma2.at(mkt('a','3')) - ma1.at(mkt('x','3','a')) ) ));
-	EXPECT_EQ( xround( res.at(mkt('a','l')) ), xround((ma2.at(mkt('a','1')) - ma1.at(mkt('l','1','a'))) + (ma2.at(mkt('a','3')) - ma1.at(mkt('l','3','a')) ) ));
-
-        EXPECT_EQ( xround( res.at(mkt('b','x')) ), xround((ma2.at(mkt('b','1')) - ma1.at(mkt('x','1','b'))) + (ma2.at(mkt('b','3')) - ma1.at(mkt('x','3','b')) ) ) );
-        EXPECT_EQ( xround( res.at(mkt('b','l')) ), xround((ma2.at(mkt('b','1')) - ma1.at(mkt('l','1','b'))) + (ma2.at(mkt('b','3')) - ma1.at(mkt('l','3','b')) ) ) );
-
-	EXPECT_EQ( xround( res2.at(mkt('a','1','x')) ), xround(ma2.at(mkt('a','1')) - ma1.at(mkt('x','1','a'))) );
-        EXPECT_EQ( xround( res2.at(mkt('a','1','l')) ), xround(ma2.at(mkt('a','1')) - ma1.at(mkt('l','1','a'))) );
-
-        EXPECT_EQ( xround( res2.at(mkt('b','1','x')) ), xround(ma2.at(mkt('b','1')) - ma1.at(mkt('x','1','b'))) );
-        EXPECT_EQ( xround( res2.at(mkt('b','1','l')) ), xround(ma2.at(mkt('b','1')) - ma1.at(mkt('l','1','b'))) );
-
-        EXPECT_EQ( xround( res2.at(mkt('a','3','x')) ), xround(ma2.at(mkt('a','3')) - ma1.at(mkt('x','3','a'))) );
-        EXPECT_EQ( xround( res2.at(mkt('a','3','l')) ), xround(ma2.at(mkt('a','3')) - ma1.at(mkt('l','3','a'))) );
-
-        EXPECT_EQ( xround( res2.at(mkt('b','3','x')) ), xround(ma2.at(mkt('b','3')) - ma1.at(mkt('x','3','b'))) );
-        EXPECT_EQ( xround( res2.at(mkt('b','3','l')) ), xround(ma2.at(mkt('b','3')) - ma1.at(mkt('l','3','b'))) );
-
-    }
     
     TEST_F(OpTest_MDim, ExecOp1)
     {
@@ -686,24 +185,6 @@ namespace {
 	    EXPECT_EQ( xround( maa2.at(i1->meta()) ), xround( maa.at(i1->meta()) * maa.at(i1->meta()) ) );
 	}
     }
-
-    template <typename T>
-    struct Monopole
-    {
-	static constexpr bool FISSTATIC = true;
-
-	static inline T apply(T f0, T x, T n)
-	{
-	    return f0 / ( 1 + x / n );
-	}
-
-	template <class... Ops>
-	static auto mk(const Ops&... ops)
-	    -> Operation<T,Monopole<T>,Ops...>
-	{
-	    return Operation<T,Monopole<T>,Ops...>(ops...);
-	}
-    };
     
     TEST_F(OpTest_MDim, ExecFOp)
     {
@@ -882,40 +363,7 @@ namespace {
 		   xround( ma4.at(mkt('l','2','b','B')) ) );
 	EXPECT_EQ( xround( ma7.at(mkt('3','B')).at(mkt('l','b')) ),
 		   xround( ma4.at(mkt('l','3','b','B')) ) );
-}
-
-    /*
-    TEST_F(OpTest_MDim, ExecAnonOp1)
-    {
-	MultiArray<double,MRange,SRange> res(mr1ptr,sr4ptr);
-	MultiArray<double,MRange> ma1(mr1ptr, v3);
-	MultiArray<double,SRange> ma2(sr2ptr, v1);
-	MultiArray<double,SRange> ma3(sr4ptr, v4);
-
-	auto si1 = MAT::getIndex( sr2ptr );
-	auto si2 = MAT::getIndex( sr3ptr );
-	auto si3 = MAT::getIndex( sr4ptr );
-	auto mi = MAT::getIndex( mr1ptr );
-	mi->operator()(si1,si2);
-	
-	res(mi,si3) = ma1(mi) + ma2(si1) + ma3(si3);
-
-	EXPECT_EQ( xround( res.at(mkt(mkt('1','a'),'A')) ), xround(0.353 + 2.917 + 1.470) );
-	EXPECT_EQ( xround( res.at(mkt(mkt('1','a'),'B')) ), xround(0.353 + 2.917 + 2.210) );
-	EXPECT_EQ( xround( res.at(mkt(mkt('1','b'),'A')) ), xround(4.005 + 2.917 + 1.470) );
-	EXPECT_EQ( xround( res.at(mkt(mkt('1','b'),'B')) ), xround(4.005 + 2.917 + 2.210) );
-
-	EXPECT_EQ( xround( res.at(mkt(mkt('2','a'),'A')) ), xround(1.070 + 9.436 + 1.470) );
-	EXPECT_EQ( xround( res.at(mkt(mkt('2','a'),'B')) ), xround(1.070 + 9.436 + 2.210) );
-	EXPECT_EQ( xround( res.at(mkt(mkt('2','b'),'A')) ), xround(2.310 + 9.436 + 1.470) );
-	EXPECT_EQ( xround( res.at(mkt(mkt('2','b'),'B')) ), xround(2.310 + 9.436 + 2.210) );
-
-	EXPECT_EQ( xround( res.at(mkt(mkt('3','a'),'A')) ), xround(9.243 + 0.373 + 1.470) );
-	EXPECT_EQ( xround( res.at(mkt(mkt('3','a'),'B')) ), xround(9.243 + 0.373 + 2.210) );
-	EXPECT_EQ( xround( res.at(mkt(mkt('3','b'),'A')) ), xround(2.911 + 0.373 + 1.470) );
-	EXPECT_EQ( xround( res.at(mkt(mkt('3','b'),'B')) ), xround(2.911 + 0.373 + 2.210) );
     }
-    */
     
 } // anonymous namspace
 
diff --git a/src/tests/test_header.h b/src/tests/test_header.h
new file mode 100644
index 0000000..663ed8f
--- /dev/null
+++ b/src/tests/test_header.h
@@ -0,0 +1,109 @@
+
+
+#include <cstdlib>
+#include "gtest/gtest.h"
+#include <iostream>
+
+#include "multi_array_header.h"
+
+#include <ctime>
+#include <cmath>
+
+namespace MAT = MultiArrayTools;
+
+double xround(double arg)
+{
+    if(std::isnan(arg)) { return 0.; }
+    return roundf(arg * 100000.) / 100000.;
+}
+    
+using namespace MAT;
+
+template <class Factory>
+void swapFactory(std::shared_ptr<RangeFactoryBase>& fptr)
+{
+    auto nptr = std::make_shared<Factory>();
+    fptr = nptr;
+}
+    
+template <class Factory, typename T>
+void swapFactory(std::shared_ptr<RangeFactoryBase>& fptr, std::initializer_list<T> ilist)
+{
+    vector<T> tmp = ilist;
+    auto nptr = std::make_shared<Factory>( tmp );
+    fptr = nptr;
+}
+
+template <class Factory, typename T>
+void swapFactory(std::shared_ptr<RangeFactoryBase>& fptr, vector<T>& ilist)
+{
+    vector<T> tmp = ilist;
+    auto nptr = std::make_shared<Factory>( tmp );
+    fptr = nptr;
+}
+
+    
+template <class Factory, class... Rs>
+void swapMFactory(std::shared_ptr<RangeFactoryBase>& fptr, const Rs&... rs)
+{
+    auto nptr = std::make_shared<Factory>( rs... );
+    fptr = nptr;
+}
+
+template <typename... Ts>
+auto mkt(Ts&&... ts) -> decltype(std::make_tuple(ts...))
+{
+    return std::make_tuple(ts...);
+}
+
+template <typename... Ts>
+auto mkts(Ts&&... ts) -> decltype(std::make_tuple(ts...))
+{
+    return std::make_tuple(static_cast<size_t>( ts )...);
+}
+
+//typedef Expressions1 EC1;
+    
+template <typename T>
+struct Pow
+{
+    static constexpr bool FISSTATIC = true;
+    typedef T value_type;
+	
+    static inline T apply(T b, T e)
+    {
+	return pow(b, e);
+    }
+
+    template <class... Ops>
+    static auto mk(const Ops&... ops)
+	-> Operation<T,Pow<T>,Ops...>
+    {
+	return Operation<T,Pow<T>,Ops...>(ops...);
+    }
+
+    static inline T apply(const std::tuple<double, double>& arg)
+    {
+	return pow(std::get<0>(arg), std::get<1>(arg));
+    }
+};
+
+
+template <typename T>
+struct Monopole
+{
+    static constexpr bool FISSTATIC = true;
+
+    static inline T apply(T f0, T x, T n)
+    {
+	return f0 / ( 1 + x / n );
+    }
+
+    template <class... Ops>
+    static auto mk(const Ops&... ops)
+	-> Operation<T,Monopole<T>,Ops...>
+    {
+	return Operation<T,Monopole<T>,Ops...>(ops...);
+    }
+};
+