diff --git a/build_linux64/libamsculib3.linux64.a b/build_linux64/libamsculib3.linux64.a
index 1e90a39..1c8e36b 100644
Binary files a/build_linux64/libamsculib3.linux64.a and b/build_linux64/libamsculib3.linux64.a differ
diff --git a/build_linux64/objstore/amscu_comp128.o b/build_linux64/objstore/amscu_comp128.o
index 1ef784f..2a47d18 100644
Binary files a/build_linux64/objstore/amscu_comp128.o and b/build_linux64/objstore/amscu_comp128.o differ
diff --git a/build_linux64/objstore/amscu_comp64.o b/build_linux64/objstore/amscu_comp64.o
index ea550fa..a217bae 100644
Binary files a/build_linux64/objstore/amscu_comp64.o and b/build_linux64/objstore/amscu_comp64.o differ
diff --git a/build_linux64/objstore/amscu_cudafunctions.o b/build_linux64/objstore/amscu_cudafunctions.o
index 9a319ce..db83c5b 100644
Binary files a/build_linux64/objstore/amscu_cudafunctions.o and b/build_linux64/objstore/amscu_cudafunctions.o differ
diff --git a/build_linux64/objstore/amscu_random.o b/build_linux64/objstore/amscu_random.o
index f47b889..0e407a7 100644
Binary files a/build_linux64/objstore/amscu_random.o and b/build_linux64/objstore/amscu_random.o differ
diff --git a/build_linux64/objstore/amscuarray.o b/build_linux64/objstore/amscuarray.o
index 86f9af8..365e905 100644
Binary files a/build_linux64/objstore/amscuarray.o and b/build_linux64/objstore/amscuarray.o differ
diff --git a/build_linux64/objstore/amscuarray_dops.o b/build_linux64/objstore/amscuarray_dops.o
index 524d91c..6e79477 100644
Binary files a/build_linux64/objstore/amscuarray_dops.o and b/build_linux64/objstore/amscuarray_dops.o differ
diff --git a/build_linux64/objstore/amscugeom.o b/build_linux64/objstore/amscugeom.o
index 21bb4e8..8896ad7 100644
Binary files a/build_linux64/objstore/amscugeom.o and b/build_linux64/objstore/amscugeom.o differ
diff --git a/build_linux64/objstore/amsculib2.o b/build_linux64/objstore/amsculib2.o
deleted file mode 100644
index f02f1dc..0000000
Binary files a/build_linux64/objstore/amsculib2.o and /dev/null differ
diff --git a/build_linux64/objstore/amsculib3.o b/build_linux64/objstore/amsculib3.o
index 3a8daef..a0141c8 100644
Binary files a/build_linux64/objstore/amsculib3.o and b/build_linux64/objstore/amsculib3.o differ
diff --git a/build_linux64/objstore/amscumath.o b/build_linux64/objstore/amscumath.o
index 4605968..8dc15b4 100644
Binary files a/build_linux64/objstore/amscumath.o and b/build_linux64/objstore/amscumath.o differ
diff --git a/build_linux64/objstore/amscurarray.o b/build_linux64/objstore/amscurarray.o
index d3ed6c5..519248e 100644
Binary files a/build_linux64/objstore/amscurarray.o and b/build_linux64/objstore/amscurarray.o differ
diff --git a/build_linux64/objstore/cuvec2f.o b/build_linux64/objstore/cuvec2f.o
new file mode 100644
index 0000000..17771d0
Binary files /dev/null and b/build_linux64/objstore/cuvec2f.o differ
diff --git a/build_linux64/objstore/cuvect3f.o b/build_linux64/objstore/cuvec3f.o
similarity index 53%
rename from build_linux64/objstore/cuvect3f.o
rename to build_linux64/objstore/cuvec3f.o
index dc38945..b3cdf72 100644
Binary files a/build_linux64/objstore/cuvect3f.o and b/build_linux64/objstore/cuvec3f.o differ
diff --git a/build_linux64/objstore/cuvect4f.o b/build_linux64/objstore/cuvec4f.o
similarity index 57%
rename from build_linux64/objstore/cuvect4f.o
rename to build_linux64/objstore/cuvec4f.o
index 4b50b03..5bde412 100644
Binary files a/build_linux64/objstore/cuvect4f.o and b/build_linux64/objstore/cuvec4f.o differ
diff --git a/build_linux64/objstore/cuvect2.o b/build_linux64/objstore/cuvect2.o
deleted file mode 100644
index 7a4f9fc..0000000
Binary files a/build_linux64/objstore/cuvect2.o and /dev/null differ
diff --git a/build_linux64/objstore/cuvect2f.o b/build_linux64/objstore/cuvect2f.o
deleted file mode 100644
index 5d0683f..0000000
Binary files a/build_linux64/objstore/cuvect2f.o and /dev/null differ
diff --git a/build_linux64/objstore/cuvect3.o b/build_linux64/objstore/cuvect3.o
deleted file mode 100644
index b1b1237..0000000
Binary files a/build_linux64/objstore/cuvect3.o and /dev/null differ
diff --git a/build_linux64/objstore/cuvect4.o b/build_linux64/objstore/cuvect4.o
deleted file mode 100644
index ec09b88..0000000
Binary files a/build_linux64/objstore/cuvect4.o and /dev/null differ
diff --git a/build_linux64/test b/build_linux64/test
index ab4136b..a3941b7 100644
Binary files a/build_linux64/test and b/build_linux64/test differ
diff --git a/include/amsculib3/amsculib3.hpp b/include/amsculib3/amsculib3.hpp
index 067b25a..ec2ffaa 100644
--- a/include/amsculib3/amsculib3.hpp
+++ b/include/amsculib3/amsculib3.hpp
@@ -19,9 +19,9 @@
 class cuvect2;
 class cuvect3;
 class cuvect4;
-class cuvect2f;
-class cuvect3f;
-class cuvect4f;
+class cuvec2f;
+class cuvec3f;
+class cuvec4f;
 
 //Need a way to define the same symbols using both host and device code
 //A solution was found here: https://stackoverflow.com/questions/9457572/cuda-host-and-device-using-same-constant-memory
@@ -46,8 +46,8 @@ namespace amscuda
 //Components
 #include <amsculib3/amscu_cudafunctions.hpp>
 #include <amsculib3/math/amscumath.hpp>
+#include <amsculib3/geom/amscugeom.hpp>
 
-#include <amsculib3/amscugeom.hpp>
 #include <amsculib3/amscuarray.hpp>
 #include <amsculib3/amscuda_binarrrw.hpp>
 #include <amsculib3/amscu_random.hpp>
diff --git a/include/amsculib3/amscugeom.hpp b/include/amsculib3/geom/amscugeom.hpp
similarity index 100%
rename from include/amsculib3/amscugeom.hpp
rename to include/amsculib3/geom/amscugeom.hpp
diff --git a/include/amsculib3/math/amscu_comp128.hpp b/include/amsculib3/math/amscu_comp128.hpp
index 05bb2c6..a100d1c 100644
--- a/include/amsculib3/math/amscu_comp128.hpp
+++ b/include/amsculib3/math/amscu_comp128.hpp
@@ -15,7 +15,8 @@ namespace cmp
         __host__ __device__ cucomp128();
         __host__ __device__ ~cucomp128();
         __host__ __device__ cucomp128(const cucomp128 &other);
-        __host__ __device__ cucomp128(const double &other);
+        __host__ __device__ explicit cucomp128(const double &other);
+        __host__ __device__ explicit cucomp128(const double &_real, const double &_imag);
 
         __host__ __device__ cucomp128& operator=(cucomp128& other);
         __host__ __device__ const cucomp128& operator=(const cucomp128& other);
@@ -58,7 +59,6 @@ namespace cmp
 
     __host__ __device__ double arg(cucomp128 z);
 
-    __host__ __device__ cucomp128 dtocomp(double _r, double _i);
     __host__ __device__ double real(cucomp128 z);
     __host__ __device__ double imag(cucomp128 z);
     __host__ __device__ cucomp128 sin(cucomp128 z);
diff --git a/include/amsculib3/math/amscu_comp64.hpp b/include/amsculib3/math/amscu_comp64.hpp
index 53c86c0..7b8cf29 100644
--- a/include/amsculib3/math/amscu_comp64.hpp
+++ b/include/amsculib3/math/amscu_comp64.hpp
@@ -15,7 +15,8 @@ namespace cmp
         __host__ __device__ cucomp64();
         __host__ __device__ ~cucomp64();
         __host__ __device__ cucomp64(const cucomp64 &other);
-        __host__ __device__ cucomp64(const float &other);
+        __host__ __device__ explicit cucomp64(const float &other);
+        __host__ __device__ explicit cucomp64(const float &_real, const float &_imag);
 
         __host__ __device__ cucomp64& operator=(cucomp64& other);
         __host__ __device__ const cucomp64& operator=(const cucomp64& other);
@@ -53,31 +54,42 @@ namespace cmp
         __host__ __device__ bool iszero() const;
         __host__ __device__ float arg() const;
         __host__ __device__ float mag() const;
+        __host__ __device__ float abs() const;
         __host__ __device__ cucomp64 conj() const;
+
+        //accumulation operators
+        __host__ __device__ cucomp64& operator+=(const cucomp64& z);
+        __host__ __device__ cucomp64& operator-=(const cucomp64& z);
+        __host__ __device__ cucomp64& operator*=(const cucomp64& z);
+        __host__ __device__ cucomp64& operator/=(const cucomp64& z);
+
+        __host__ __device__ cucomp64& operator+=(const float& z);
+        __host__ __device__ cucomp64& operator-=(const float& z);
+        __host__ __device__ cucomp64& operator*=(const float& z);
+        __host__ __device__ cucomp64& operator/=(const float& z);
+        
     };
 
-    __host__ __device__ float arg(cucomp64 z);
-
-    __host__ __device__ cucomp64 dtocomp64(float _r, float _i);
-    __host__ __device__ float real(cucomp64 z);
-    __host__ __device__ float imag(cucomp64 z);
-    __host__ __device__ cucomp64 sin(cucomp64 z);
-    __host__ __device__ cucomp64 cos(cucomp64 z);
-    __host__ __device__ cucomp64 tan(cucomp64 z);
-    __host__ __device__ cucomp64 exp(cucomp64 z);
-    __host__ __device__ cucomp64 log(cucomp64 z);
-    __host__ __device__ float abs(cucomp64 z);
-    __host__ __device__ cucomp64 conj(cucomp64 z);
+    __host__ __device__ float arg(const cucomp64 &z);
+    __host__ __device__ float real(const cucomp64 &z);
+    __host__ __device__ float imag(const cucomp64 &z);
+    __host__ __device__ cucomp64 sin(const cucomp64 &z);
+    __host__ __device__ cucomp64 cos(const cucomp64 &z);
+    __host__ __device__ cucomp64 tan(const cucomp64 &z);
+    __host__ __device__ cucomp64 exp(const cucomp64 &z);
+    __host__ __device__ cucomp64 log(const cucomp64 &z);
+    __host__ __device__ float abs(const cucomp64 &z);
+    __host__ __device__ cucomp64 conj(const cucomp64 &z);
 
     // //need hyperbolic trig Functions
-    __host__ __device__ cucomp64 cosh(cucomp64 z);
-    __host__ __device__ cucomp64 sinh(cucomp64 z);
-    __host__ __device__ cucomp64 tanh(cucomp64 z);
+    __host__ __device__ cucomp64 cosh(const cucomp64 &z);
+    __host__ __device__ cucomp64 sinh(const cucomp64 &z);
+    __host__ __device__ cucomp64 tanh(const cucomp64 &z);
 
-    __host__ __device__ cucomp64 pow(cucomp64 z1, cucomp64 z2);
+    __host__ __device__ cucomp64 pow(const cucomp64 &z1, const cucomp64 &z2);
 
     // //returns "complex sign" of complex number - 0, or a unit number with same argument
-    __host__ __device__ cucomp64 csgn(cucomp64 z);
+    __host__ __device__ cucomp64 csgn(const cucomp64 &z);
 
 void test_cucomp64_1();
 
diff --git a/include/amsculib3/math/amscumath.hpp b/include/amsculib3/math/amscumath.hpp
index 6452455..df934b8 100644
--- a/include/amsculib3/math/amscumath.hpp
+++ b/include/amsculib3/math/amscumath.hpp
@@ -59,12 +59,12 @@ namespace amscuda
 //component headers
 #include <amsculib3/math/amscu_comp64.hpp>
 #include <amsculib3/math/amscu_comp128.hpp>
-#include <amsculib3/math/cuvect2.hpp>
-#include <amsculib3/math/cuvect3.hpp>
-#include <amsculib3/math/cuvect4.hpp>
-#include <amsculib3/math/cuvect2f.hpp>
-#include <amsculib3/math/cuvect3f.hpp>
-#include <amsculib3/math/cuvect4f.hpp>
+//#include <amsculib3/math/cuvect2.hpp>
+//#include <amsculib3/math/cuvect3.hpp>
+//#include <amsculib3/math/cuvect4.hpp>
+#include <amsculib3/math/cuvec2f.hpp>
+#include <amsculib3/math/cuvec3f.hpp>
+#include <amsculib3/math/cuvec4f.hpp>
 
 #endif
 
diff --git a/include/amsculib3/math/cuvect2f.hpp b/include/amsculib3/math/cuvec2f.hpp
similarity index 52%
rename from include/amsculib3/math/cuvect2f.hpp
rename to include/amsculib3/math/cuvec2f.hpp
index ed6f58c..00a88e5 100644
--- a/include/amsculib3/math/cuvect2f.hpp
+++ b/include/amsculib3/math/cuvec2f.hpp
@@ -1,33 +1,33 @@
-#ifndef __CUVECT2F_HPP__
-#define __CUVECT2F_HPP__
+#ifndef __CUVEC2F_HPP__
+#define __CUVEC2F_HPP__
 
 namespace amscuda
 {
 
-    class cuvect2f
+    class cuvec2f
     {
         public:
         float x;
         float y;
 
-        __host__ __device__ cuvect2f();
-        __host__ __device__ ~cuvect2f();
-        __host__ __device__ cuvect2f(const float &_x, const float &_y);
+        __host__ __device__ cuvec2f();
+        __host__ __device__ ~cuvec2f();
+        __host__ __device__ cuvec2f(const float &_x, const float &_y);
         
         
         __host__ __device__ float& operator[](const int &I);
         __host__ __device__ const float& operator[](const int &I) const;
 
-        __host__ __device__ cuvect2f operator+(const cuvect2f &rhs);
-        __host__ __device__ cuvect2f operator-(const cuvect2f &rhs);
-        __host__ __device__ cuvect2f operator*(const float &rhs);
-        __host__ __device__ cuvect2f operator/(const float &rhs);
-        __host__ __device__ friend cuvect2f operator-(const cuvect2f &rhs);
+        __host__ __device__ cuvec2f operator+(const cuvec2f &rhs);
+        __host__ __device__ cuvec2f operator-(const cuvec2f &rhs);
+        __host__ __device__ cuvec2f operator*(const float &rhs);
+        __host__ __device__ cuvec2f operator/(const float &rhs);
+        __host__ __device__ friend cuvec2f operator-(const cuvec2f &rhs);
 
-        __host__ __device__ cuvect2f& operator+=(const cuvect2f &rhs);
-        __host__ __device__ cuvect2f& operator-=(const cuvect2f &rhs);
-        __host__ __device__ cuvect2f& operator/=(const float &rhs);
-        __host__ __device__ cuvect2f& operator*=(const float &rhs);
+        __host__ __device__ cuvec2f& operator+=(const cuvec2f &rhs);
+        __host__ __device__ cuvec2f& operator-=(const cuvec2f &rhs);
+        __host__ __device__ cuvec2f& operator/=(const float &rhs);
+        __host__ __device__ cuvec2f& operator*=(const float &rhs);
 
     };
 
@@ -59,7 +59,7 @@ namespace amscuda
         __host__ __device__ cumat2f operator-(const cumat2f &rhs);
         __host__ __device__ cumat2f operator*(const float &rhs);
         __host__ __device__ cumat2f operator/(const float &rhs);
-        __host__ __device__ cuvect2f operator*(const cuvect2f &rhs);
+        __host__ __device__ cuvec2f operator*(const cuvec2f &rhs);
         __host__ __device__ cumat2f operator*(const cumat2f &rhs);
         __host__ __device__ friend cumat2f operator-(const cumat2f &rhs);
         
@@ -78,45 +78,18 @@ namespace amscuda
         __host__ __device__ cumat2f& operator*=(const cumat2f &rhs);
     };
 
-    __host__ __device__ float cuvect2f_dot(const cuvect2f &a, const cuvect2f &b);
-    __host__ __device__ float cuvect2f_cross(const cuvect2f &a, const cuvect2f &b);
-    __host__ __device__ float cuvect2f_norm(const cuvect2f &a);
-    __host__ __device__ cuvect2f cuvect2f_normalize(const cuvect2f &a);
-    __host__ __device__ cuvect2f cuvect2f_proj(const cuvect2f &a, const cuvect2f &b);
-
+    __host__ __device__ float cuvec2f_dot(const cuvec2f &a, const cuvec2f &b);
+    __host__ __device__ float cuvec2f_cross(const cuvec2f &a, const cuvec2f &b);
+    __host__ __device__ float cuvec2f_norm(const cuvec2f &a);
+    __host__ __device__ cuvec2f cuvec2f_normalize(const cuvec2f &a);
+    __host__ __device__ cuvec2f cuvec2f_proj(const cuvec2f &a, const cuvec2f &b);
     __host__ __device__ cumat2f cumat2f_rot_from_angle(const float &angle);
 
 
     ///////////
     // Tests //
     ///////////
-    void test_cuvect2f_1();
-
-
-    ///////////////////////////
-    //legacy array operations//
-    ///////////////////////////
-
-    //2x2 matrix operations
-    //matrix order is assumed to be mat[I,J] = mat[I+2*J]
-
-    //transpose a 2x2 matrix in place
-    __host__ __device__ void mat2f_transpose(float *mat2inout);
-    
-    //copies src to dest
-    __host__ __device__ void mat2f_copy(float *mat2f_dest, const float *mat2f_src);
-    
-    //inverts mat?inout[4]
-    __host__ __device__ void mat2f_inverse(float *mat2inout);
-    
-    //rotatin matrix from angle
-    __host__ __device__ void mat2f_rot_from_angle(float angle, float *mat2);
-
-    //multiplies c = a*b
-    __host__ __device__ void mat2f_mult(float *mat2a, float *mat2b, float *mat2c);
-
-    // ret = a*b
-    __host__ __device__ cuvect2f mat2f_mult(float *mat2a, const cuvect2f &b);
+    void test_cuvec2f_1();
 
 
 
diff --git a/include/amsculib3/math/cuvect3f.hpp b/include/amsculib3/math/cuvec3f.hpp
similarity index 61%
rename from include/amsculib3/math/cuvect3f.hpp
rename to include/amsculib3/math/cuvec3f.hpp
index 8a3008e..08962b4 100644
--- a/include/amsculib3/math/cuvect3f.hpp
+++ b/include/amsculib3/math/cuvec3f.hpp
@@ -1,34 +1,34 @@
-#ifndef __CUVECT3F_HPP__
-#define __CUVECT3F_HPP__
+#ifndef __CUVEC3F_HPP__
+#define __CUVEC3F_HPP__
 
 namespace amscuda
 {
 
-    class cuvect3f
+    class cuvec3f
     {
         public:
         float x;
         float y;
         float z;
 
-        __host__ __device__ cuvect3f();
-        __host__ __device__ ~cuvect3f();
-        __host__ __device__ cuvect3f(const float &_x, const float &_y, const float &_z);
+        __host__ __device__ cuvec3f();
+        __host__ __device__ ~cuvec3f();
+        __host__ __device__ cuvec3f(const float &_x, const float &_y, const float &_z);
         
         
         __host__ __device__ float& operator[](const int &I);
         __host__ __device__ const float& operator[](const int &I) const;
 
-        __host__ __device__ cuvect3f operator+(const cuvect3f &rhs);
-        __host__ __device__ cuvect3f operator-(const cuvect3f &rhs);
-        __host__ __device__ cuvect3f operator*(const float &rhs);
-        __host__ __device__ cuvect3f operator/(const float &rhs);
-        __host__ __device__ friend cuvect3f operator-(const cuvect3f &rhs);
+        __host__ __device__ cuvec3f operator+(const cuvec3f &rhs);
+        __host__ __device__ cuvec3f operator-(const cuvec3f &rhs);
+        __host__ __device__ cuvec3f operator*(const float &rhs);
+        __host__ __device__ cuvec3f operator/(const float &rhs);
+        __host__ __device__ friend cuvec3f operator-(const cuvec3f &rhs);
 
-        __host__ __device__ cuvect3f& operator+=(const cuvect3f &rhs);
-        __host__ __device__ cuvect3f& operator-=(const cuvect3f &rhs);
-        __host__ __device__ cuvect3f& operator/=(const float &rhs);
-        __host__ __device__ cuvect3f& operator*=(const float &rhs);
+        __host__ __device__ cuvec3f& operator+=(const cuvec3f &rhs);
+        __host__ __device__ cuvec3f& operator-=(const cuvec3f &rhs);
+        __host__ __device__ cuvec3f& operator/=(const float &rhs);
+        __host__ __device__ cuvec3f& operator*=(const float &rhs);
 
     };
 
@@ -62,7 +62,7 @@ namespace amscuda
         __host__ __device__ cumat3f operator-(const cumat3f &rhs);
         __host__ __device__ cumat3f operator*(const float &rhs);
         __host__ __device__ cumat3f operator/(const float &rhs);
-        __host__ __device__ cuvect3f operator*(const cuvect3f &rhs);
+        __host__ __device__ cuvec3f operator*(const cuvec3f &rhs);
         __host__ __device__ cumat3f operator*(const cumat3f &rhs);
         __host__ __device__ friend cumat3f operator-(const cumat3f &rhs);
         
@@ -82,21 +82,19 @@ namespace amscuda
         
     };
 
-    __host__ __device__ float cuvect3f_dot(const cuvect3f &a,const cuvect3f &b);
-    __host__ __device__ cuvect3f cuvect3f_cross(const cuvect3f &a, const cuvect3f &b);
-    __host__ __device__ float cuvect3f_norm(const cuvect3f &a);
-    __host__ __device__ cuvect3f cuvect3f_normalize(const cuvect3f &a);
-    __host__ __device__ cuvect3f cuvect3f_proj(const cuvect3f &a, const cuvect3f &b);
+    __host__ __device__ float cuvec3f_dot(const cuvec3f &a,const cuvec3f &b);
+    __host__ __device__ cuvec3f cuvec3f_cross(const cuvec3f &a, const cuvec3f &b);
+    __host__ __device__ float cuvec3f_norm(const cuvec3f &a);
+    __host__ __device__ cuvec3f cuvec3f_normalize(const cuvec3f &a);
+    __host__ __device__ cuvec3f cuvec3f_proj(const cuvec3f &a, const cuvec3f &b);
 
-    __host__ __device__ cumat3f hodge_dual(const cuvect3f &vin);
-    __host__ __device__ cuvect3f hodge_dual(const cumat3f &min);
-    __host__ __device__ cumat3f rotmat_from_axisangle(const cuvect3f &axis, const float &angle);
+    __host__ __device__ cumat3f hodge_dual(const cuvec3f &vin);
+    __host__ __device__ cuvec3f hodge_dual(const cumat3f &min);
+    __host__ __device__ cumat3f rotmat_from_axisangle(const cuvec3f &axis, const float &angle);
 
     __host__ void test_cudavectf_logic1();
 
 };
 
-#include <amsculib2/cuvect3f_legacy.hpp>
-
 #endif
 
diff --git a/include/amsculib3/math/cuvect4f.hpp b/include/amsculib3/math/cuvec4f.hpp
similarity index 68%
rename from include/amsculib3/math/cuvect4f.hpp
rename to include/amsculib3/math/cuvec4f.hpp
index 4139caf..d79f298 100644
--- a/include/amsculib3/math/cuvect4f.hpp
+++ b/include/amsculib3/math/cuvec4f.hpp
@@ -1,10 +1,10 @@
-#ifndef __CUVECT4F_HPP__
-#define __CUVECT4F_HPP__
+#ifndef __CUVEC4F_HPP__
+#define __CUVEC4F_HPP__
 
 namespace amscuda
 {
 
-    class cuvect4f
+    class cuvec4f
     {
         public:
         float x;
@@ -12,24 +12,24 @@ namespace amscuda
         float z;
         float w;
 
-        __host__ __device__ cuvect4f();
-        __host__ __device__ ~cuvect4f();
-        __host__ __device__ cuvect4f(const float &_x, const float &_y, const float &_z, const float &_w);
+        __host__ __device__ cuvec4f();
+        __host__ __device__ ~cuvec4f();
+        __host__ __device__ cuvec4f(const float &_x, const float &_y, const float &_z, const float &_w);
         
         
         __host__ __device__ float& operator[](const int &I);
         __host__ __device__ const float& operator[](const int &I) const;
 
-        __host__ __device__ cuvect4f operator+(const cuvect4f &rhs);
-        __host__ __device__ cuvect4f operator-(const cuvect4f &rhs);
-        __host__ __device__ cuvect4f operator*(const float &rhs);
-        __host__ __device__ cuvect4f operator/(const float &rhs);
-        __host__ __device__ friend cuvect4f operator-(const cuvect4f &rhs);
+        __host__ __device__ cuvec4f operator+(const cuvec4f &rhs);
+        __host__ __device__ cuvec4f operator-(const cuvec4f &rhs);
+        __host__ __device__ cuvec4f operator*(const float &rhs);
+        __host__ __device__ cuvec4f operator/(const float &rhs);
+        __host__ __device__ friend cuvec4f operator-(const cuvec4f &rhs);
 
-        __host__ __device__ cuvect4f& operator+=(const cuvect4f &rhs);
-        __host__ __device__ cuvect4f& operator-=(const cuvect4f &rhs);
-        __host__ __device__ cuvect4f& operator/=(const float &rhs);
-        __host__ __device__ cuvect4f& operator*=(const float &rhs);
+        __host__ __device__ cuvec4f& operator+=(const cuvec4f &rhs);
+        __host__ __device__ cuvec4f& operator-=(const cuvec4f &rhs);
+        __host__ __device__ cuvec4f& operator/=(const float &rhs);
+        __host__ __device__ cuvec4f& operator*=(const float &rhs);
     };
 
     class cumat4f
@@ -66,7 +66,7 @@ namespace amscuda
         __host__ __device__ cumat4f operator-(const cumat4f &rhs);
         __host__ __device__ cumat4f operator*(const float &rhs);
         __host__ __device__ cumat4f operator/(const float &rhs);
-        __host__ __device__ cuvect4f operator*(const cuvect4f &rhs);
+        __host__ __device__ cuvec4f operator*(const cuvec4f &rhs);
         __host__ __device__ cumat4f operator*(const cumat4f &rhs);
         __host__ __device__ friend cumat4f operator-(const cumat4f &rhs);
         
@@ -85,10 +85,10 @@ namespace amscuda
         __host__ __device__ cumat4f& operator*=(const cumat4f &rhs);
     };
 
-    __host__ __device__ float cuvect4f_dot(cuvect4f &a, cuvect4f &b);
-    __host__ __device__ float cuvect4f_norm(cuvect4f &a);
-    __host__ __device__ cuvect4f cuvect4f_normalize(cuvect4f &a);
-    __host__ __device__ cuvect4f cuvect4f_proj(cuvect4f &a, cuvect4f &b);
+    __host__ __device__ float cuvec4f_dot(cuvec4f &a, cuvec4f &b);
+    __host__ __device__ float cuvec4f_norm(cuvec4f &a);
+    __host__ __device__ cuvec4f cuvec4f_normalize(cuvec4f &a);
+    __host__ __device__ cuvec4f cuvec4f_proj(cuvec4f &a, cuvec4f &b);
 
 
 };
diff --git a/include/amsculib3/math/cuvect2.hpp b/include/amsculib3/math/cuvect2.hpp
deleted file mode 100644
index 170897e..0000000
--- a/include/amsculib3/math/cuvect2.hpp
+++ /dev/null
@@ -1,87 +0,0 @@
-#ifndef __CUVECT2_HPP__
-#define __CUVECT2_HPP__
-
-namespace amscuda
-{
-
-    class cuvect2
-    {
-        public:
-        double x;
-        double y;
-
-
-        __host__ __device__ cuvect2();
-        __host__ __device__ ~cuvect2();
-        __host__ __device__ cuvect2(double _x, double _y);
-
-        __host__ __device__ double& operator[](const int I);
-        __host__ __device__ const double& operator[](const int I) const;
-
-        __host__ __device__ cuvect2 operator+(cuvect2 lhs);
-        __host__ __device__ cuvect2 operator-(cuvect2 lhs);
-        __host__ __device__ cuvect2 operator*(double lhs);
-        __host__ __device__ cuvect2 operator/(double lhs);
-        __host__ __device__ friend cuvect2 operator-(cuvect2 rhs);
-    };
-
-    class cumat2
-    {
-        public:
-        double dat[4];
-
-        __host__ __device__ cumat2();
-        __host__ __device__ ~cumat2();
-        __host__ __device__ double& operator[](const int I);
-        __host__ __device__ double& operator()(const int I, const int J);
-        __host__ __device__ double& at(const int I, const int J);
-
-        __host__ __device__ cumat2 operator+(cumat2 lhs);
-        __host__ __device__ cumat2 operator-(cumat2 lhs);
-        __host__ __device__ cumat2 operator*(double lhs);
-        __host__ __device__ cumat2 operator/(double lhs);
-        __host__ __device__ cuvect2 operator*(cuvect2 lhs);
-        __host__ __device__ cumat2 operator*(cumat2 lhs);
-        __host__ __device__ friend cumat2 operator-(cumat2 rhs);
-
-        __host__ __device__ double det();
-        __host__ __device__ cumat2 transpose();
-        __host__ __device__ cumat2 inverse();
-    };
-
-    __host__ __device__ double cuvect2_dot(cuvect2 a, cuvect2 b);
-    __host__ __device__ double cuvect2_cross(cuvect2 a, cuvect2 b);
-    __host__ __device__ double cuvect2_norm(cuvect2 a);
-    __host__ __device__ cuvect2 cuvect2_normalize(cuvect2 a);
-    __host__ __device__ cuvect2 cuvect2_proj(cuvect2 a, cuvect2 b);
-
-    //2x2 matrix operations
-    //matrix order is assumed to be mat[I,J] = mat[I+3*J]
-
-    //transpose a 2x2 matrix in place
-    __host__ __device__ void mat2_transpose(double *mat2inout);
-    
-    //copies src to dest
-    __host__ __device__ void mat2_copy(double *mat2_dest, const double *mat2_src);
-    
-    //inverts mat?inout[4]
-    __host__ __device__ void mat2_inverse(double *mat2inout);
-    
-    //rotatin matrix from angle
-    __host__ __device__ void mat2_rot_from_angle(double angle, double *mat2);
-
-    //multiplies c = a*b
-    __host__ __device__ void mat2_mult(double *mat2a, double *mat2b, double *mat2c);
-
-    // ret = a*b
-    __host__ __device__ cuvect2 mat2_mult(double *mat2a, cuvect2 b);
-
-
-
-    void test_cuvect2_1();
-    
-
-}; //end namespace amscuda
-
-#endif
-
diff --git a/include/amsculib3/math/cuvect3.hpp b/include/amsculib3/math/cuvect3.hpp
deleted file mode 100644
index 00f3834..0000000
--- a/include/amsculib3/math/cuvect3.hpp
+++ /dev/null
@@ -1,122 +0,0 @@
-#ifndef __CUVECT3_HPP__
-#define __CUVECT3_HPP__
-
-namespace amscuda
-{
-
-    class cuvect3
-    {
-        public:
-        double x;
-        double y;
-        double z;
-
-        __host__ __device__ cuvect3();
-        __host__ __device__ ~cuvect3();
-        __host__ __device__ cuvect3(const double &_x, const double &_y, const double &_z);
-        
-        
-        __host__ __device__ double& operator[](const int &I);
-        __host__ __device__ const double& operator[](const int &I) const;
-
-        __host__ __device__ cuvect3 operator+(const cuvect3 &rhs);
-        __host__ __device__ cuvect3 operator-(const cuvect3 &rhs);
-        __host__ __device__ cuvect3 operator*(const double &rhs);
-        __host__ __device__ cuvect3 operator/(const double &rhs);
-        __host__ __device__ friend cuvect3 operator-(const cuvect3 &rhs);
-
-        __host__ __device__ cuvect3& operator+=(const cuvect3 &rhs);
-        __host__ __device__ cuvect3& operator-=(const cuvect3 &rhs);
-        __host__ __device__ cuvect3& operator/=(const double &rhs);
-        __host__ __device__ cuvect3& operator*=(const double &rhs);
-
-    };
-
-    class cumat3
-    {
-        public:
-        double m00,m10,m20; //named references to force register use?
-        double m01,m11,m21; //switched to column-major-order to match GLSL/lapack
-        double m02,m12,m22;
-
-        __host__ __device__ cumat3();
-        __host__ __device__ ~cumat3();
-        __host__ __device__ cumat3(
-            const double & _m00, const double & _m01, const double & _m02,
-            const double & _m10, const double & _m11, const double & _m12,
-            const double & _m20, const double & _m21, const double & _m22
-        );
-
-        __host__ __device__ double& operator[](const int &I);
-        __host__ __device__ double& operator()(const int &I, const int &J);
-        __host__ __device__ double& at(const int &I, const int &J);
-
-        __host__ __device__ const double& operator[](const int &I) const;
-        __host__ __device__ const double& operator()(const int &I, const int &J) const;
-        __host__ __device__ const double& at(const int &I, const int &J) const;
-
-        __host__ __device__ cumat3 operator+(const cumat3 &rhs);
-        __host__ __device__ cumat3 operator-(const cumat3 &rhs);
-        __host__ __device__ cumat3 operator*(const double &rhs);
-        __host__ __device__ cumat3 operator/(const double &rhs);
-        __host__ __device__ cuvect3 operator*(const cuvect3 &rhs);
-        __host__ __device__ cumat3 operator*(const cumat3 &rhs);
-        __host__ __device__ friend cumat3 operator-(const cumat3 &rhs);
-        
-        __host__ __device__ double det();
-        __host__ __device__ cumat3 transpose();
-        __host__ __device__ cumat3 inverse();
-
-        __host__ __device__ double* data(); //pointer to double[9] representation of matrix
-        __host__ __device__ const double* data() const; //pointer to double[9] representation of matrix
-
-        //In place operations (to save GPU register use)
-        __host__ __device__ cumat3& operator+=(const cumat3 &rhs);
-        __host__ __device__ cumat3& operator-=(const cumat3 &rhs);
-        __host__ __device__ cumat3& operator/=(const double &rhs);
-        __host__ __device__ cumat3& operator*=(const double &rhs);
-        __host__ __device__ cumat3& operator*=(const cumat3 &rhs);
-        
-    };
-
-    __host__ __device__ double cuvect3_dot(const cuvect3 &a,const cuvect3 &b);
-    __host__ __device__ cuvect3 cuvect3_cross(const cuvect3 &a, const cuvect3 &b);
-    __host__ __device__ double cuvect3_norm(const cuvect3 &a);
-    __host__ __device__ cuvect3 cuvect3_normalize(const cuvect3 &a);
-    __host__ __device__ cuvect3 cuvect3_proj(const cuvect3 &a, const cuvect3 &b);
-
-    __host__ __device__ cumat3 hodge_dual(const cuvect3 &vin);
-    __host__ __device__ cuvect3 hodge_dual(const cumat3 &min);
-    __host__ __device__ cumat3 rotmat_from_axisangle(const cuvect3 &axis, const double &angle);
-
-
-    //3x3 matrix operations
-    //matrix order is assumed to be mat[I,J] = mat[I+3*J]
-    
-    //transposes a 3x3 (9 element) matrix
-    __host__ __device__ void mat3f_transpose(double *mat3inout);
-
-    //copies src to dest
-    __host__ __device__ void mat3f_copy(double *mat3f_dest, const double *mat3f_src);
-
-    //returns determinant of 3x3 matrix
-    __host__ __device__ double mat3f_det(double *mat3in);
-
-    //inverts a 3x3 (9 element) matrix
-    __host__ __device__ void mat3f_inverse(double *mat3inout);
-
-    __host__ __device__ cuvect3 mat3f_mult(double *mat3in, const cuvect3 &cvin);
-    __host__ __device__ void mat3f_mult(double *matina, double *matinb, double *matout);
-
-    __host__ __device__ void mat3f_hodgedual(const cuvect3 &vecin, double *matout);
-    __host__ __device__ void mat3f_hodgedual(double *matin, cuvect3 &vecout);
-
-    //returns direction cosine rotation matrix from axis and angle
-    __host__ __device__ void mat3f_rot_from_axisangle(cuvect3 axis, double angle, double *matout);
-
-    __host__ void test_cudavect_logic1();
-
-}; //end namespace amscuda
-
-#endif
-
diff --git a/include/amsculib3/math/cuvect4.hpp b/include/amsculib3/math/cuvect4.hpp
deleted file mode 100644
index 2125852..0000000
--- a/include/amsculib3/math/cuvect4.hpp
+++ /dev/null
@@ -1,61 +0,0 @@
-#ifndef __CUVECT4_HPP__
-#define __CUVECT4_HPP__
-
-namespace amscuda
-{
-
-    class cuvect4
-    {
-        public:
-        double x;
-        double y;
-        double z;
-        double w;
-
-        __host__ __device__ cuvect4();
-        __host__ __device__ ~cuvect4();
-        __host__ __device__ cuvect4(double _x, double _y, double _z, double _w);
-                
-        __host__ __device__ double& operator[](const int I);
-        __host__ __device__ const double& operator[](const int I) const;
-
-        __host__ __device__ cuvect4 operator+(cuvect4 lhs);
-        __host__ __device__ cuvect4 operator-(cuvect4 lhs);
-        __host__ __device__ cuvect4 operator*(double lhs);
-        __host__ __device__ cuvect4 operator/(double lhs);
-        __host__ __device__ friend cuvect4 operator-(cuvect4 rhs);
-    };
-
-    class cumat4
-    {
-        public:
-        double dat[16];
-
-        __host__ __device__ cumat4();
-        __host__ __device__ ~cumat4();
-        __host__ __device__ double& operator[](const int I);
-        __host__ __device__ double& operator()(const int I, const int J);
-        __host__ __device__ double& at(const int I, const int J);
-
-        __host__ __device__ cumat4 operator+(cumat4 lhs);
-        __host__ __device__ cumat4 operator-(cumat4 lhs);
-        __host__ __device__ cumat4 operator*(double lhs);
-        __host__ __device__ cumat4 operator/(double lhs);
-        __host__ __device__ cuvect4 operator*(cuvect4 lhs);
-        __host__ __device__ cumat4 operator*(cumat4 lhs);
-        __host__ __device__ friend cumat4 operator-(cumat4 rhs);
-
-        __host__ __device__ double det();
-        __host__ __device__ cumat4 transpose();
-        __host__ __device__ cumat4 inverse();
-    };
-
-    __host__ __device__ double cuvect4_dot(cuvect4 a, cuvect4 b);
-    __host__ __device__ double cuvect4_norm(cuvect4 a);
-    __host__ __device__ cuvect4 cuvect4_normalize(cuvect4 a);
-    __host__ __device__ cuvect4 cuvect4_proj(cuvect4 a, cuvect4 b);
-
-}; //end namespace amscuda
-
-#endif
-
diff --git a/old/9apr26_prerefactor/include/amsculib2/amsculib2.hpp b/old/9apr26_prerefactor/include/amsculib2/amsculib2.hpp
index 2c026c4..0796683 100644
--- a/old/9apr26_prerefactor/include/amsculib2/amsculib2.hpp
+++ b/old/9apr26_prerefactor/include/amsculib2/amsculib2.hpp
@@ -19,9 +19,9 @@
 class cuvect2;
 class cuvect3;
 class cuvect4;
-class cuvect2f;
-class cuvect3f;
-class cuvect4f;
+class cuvec2f;
+class cuvec3f;
+class cuvec4f;
 
 //Need a way to define the same symbols using both host and device code
 //A solution was found here: https://stackoverflow.com/questions/9457572/cuda-host-and-device-using-same-constant-memory
@@ -51,9 +51,9 @@ namespace amscuda
 #include <amsculib2/cuvect2.hpp>
 #include <amsculib2/cuvect3.hpp>
 #include <amsculib2/cuvect4.hpp>
-#include <amsculib2/cuvect2f.hpp>
-#include <amsculib2/cuvect3f.hpp>
-#include <amsculib2/cuvect4f.hpp>
+#include <amsculib2/cuvec2f.hpp>
+#include <amsculib2/cuvec3f.hpp>
+#include <amsculib2/cuvec4f.hpp>
 #include <amsculib2/amscugeom.hpp>
 #include <amsculib2/amscuarray.hpp>
 #include <amsculib2/amscuda_binarrrw.hpp>
diff --git a/old/9apr26_prerefactor/include/amsculib2/cuvect2f.hpp b/old/9apr26_prerefactor/include/amsculib2/cuvect2f.hpp
index a0ba6ce..88a0a40 100644
--- a/old/9apr26_prerefactor/include/amsculib2/cuvect2f.hpp
+++ b/old/9apr26_prerefactor/include/amsculib2/cuvect2f.hpp
@@ -1,28 +1,28 @@
-#ifndef __CUVECT2F_HPP__
-#define __CUVECT2F_HPP__
+#ifndef __cuvec2f_HPP__
+#define __cuvec2f_HPP__
 
 namespace amscuda
 {
 
-    class cuvect2f
+    class cuvec2f
     {
         public:
         float x;
         float y;
 
 
-        __host__ __device__ cuvect2f();
-        __host__ __device__ ~cuvect2f();
-        __host__ __device__ cuvect2f(const float &_x, const float &_y);
+        __host__ __device__ cuvec2f();
+        __host__ __device__ ~cuvec2f();
+        __host__ __device__ cuvec2f(const float &_x, const float &_y);
 
         __host__ __device__ float& operator[](const int &I);
         __host__ __device__ const float& operator[](const int &I) const;
 
-        __host__ __device__ cuvect2f operator+(const cuvect2f &lhs);
-        __host__ __device__ cuvect2f operator-(const cuvect2f &lhs);
-        __host__ __device__ cuvect2f operator*(const float &lhs);
-        __host__ __device__ cuvect2f operator/(const float &lhs);
-        __host__ __device__ friend cuvect2f operator-(const cuvect2f &rhs);
+        __host__ __device__ cuvec2f operator+(const cuvec2f &lhs);
+        __host__ __device__ cuvec2f operator-(const cuvec2f &lhs);
+        __host__ __device__ cuvec2f operator*(const float &lhs);
+        __host__ __device__ cuvec2f operator/(const float &lhs);
+        __host__ __device__ friend cuvec2f operator-(const cuvec2f &rhs);
     };
 
     class cumat2f
@@ -40,7 +40,7 @@ namespace amscuda
         __host__ __device__ cumat2f operator-(cumat2f lhs);
         __host__ __device__ cumat2f operator*(float lhs);
         __host__ __device__ cumat2f operator/(float lhs);
-        __host__ __device__ cuvect2f operator*(cuvect2f lhs);
+        __host__ __device__ cuvec2f operator*(cuvec2f lhs);
         __host__ __device__ cumat2f operator*(cumat2f lhs);
         __host__ __device__ friend cumat2f operator-(cumat2f rhs);
 
@@ -49,11 +49,11 @@ namespace amscuda
         __host__ __device__ cumat2f inverse();
     };
 
-    __host__ __device__ float cuvect2f_dot(const cuvect2f &a, const cuvect2f &b);
-    __host__ __device__ float cuvect2f_cross(const cuvect2f &a, const cuvect2f &b);
-    __host__ __device__ float cuvect2f_norm(const cuvect2f &a);
-    __host__ __device__ cuvect2f cuvect2f_normalize(const cuvect2f &a);
-    __host__ __device__ cuvect2f cuvect2f_proj(const cuvect2f &a, const cuvect2f &b);
+    __host__ __device__ float cuvec2f_dot(const cuvec2f &a, const cuvec2f &b);
+    __host__ __device__ float cuvec2f_cross(const cuvec2f &a, const cuvec2f &b);
+    __host__ __device__ float cuvec2f_norm(const cuvec2f &a);
+    __host__ __device__ cuvec2f cuvec2f_normalize(const cuvec2f &a);
+    __host__ __device__ cuvec2f cuvec2f_proj(const cuvec2f &a, const cuvec2f &b);
 
     //2x2 matrix operations
     //matrix order is assumed to be mat[I,J] = mat[I+3*J]
@@ -74,10 +74,10 @@ namespace amscuda
     __host__ __device__ void mat2f_mult(float *mat2a, float *mat2b, float *mat2c);
 
     // ret = a*b
-    __host__ __device__ cuvect2f mat2f_mult(float *mat2a, const cuvect2f &b);
+    __host__ __device__ cuvec2f mat2f_mult(float *mat2a, const cuvec2f &b);
 
 
-    void test_cuvect2f_1();
+    void test_cuvec2f_1();
 
 
 };
diff --git a/old/9apr26_prerefactor/include/amsculib2/cuvect3f.hpp b/old/9apr26_prerefactor/include/amsculib2/cuvect3f.hpp
index 0e06f45..b133293 100644
--- a/old/9apr26_prerefactor/include/amsculib2/cuvect3f.hpp
+++ b/old/9apr26_prerefactor/include/amsculib2/cuvect3f.hpp
@@ -1,34 +1,34 @@
-#ifndef __CUVECT3F_HPP__
-#define __CUVECT3F_HPP__
+#ifndef __cuvec3f_HPP__
+#define __cuvec3f_HPP__
 
 namespace amscuda
 {
 
-    class cuvect3f
+    class cuvec3f
     {
         public:
         float x;
         float y;
         float z;
 
-        __host__ __device__ cuvect3f();
-        __host__ __device__ ~cuvect3f();
-        __host__ __device__ cuvect3f(const float &_x, const float &_y, const float &_z);
+        __host__ __device__ cuvec3f();
+        __host__ __device__ ~cuvec3f();
+        __host__ __device__ cuvec3f(const float &_x, const float &_y, const float &_z);
         
         
         __host__ __device__ float& operator[](const int &I);
         __host__ __device__ const float& operator[](const int &I) const;
 
-        __host__ __device__ cuvect3f operator+(const cuvect3f &rhs);
-        __host__ __device__ cuvect3f operator-(const cuvect3f &rhs);
-        __host__ __device__ cuvect3f operator*(const float &rhs);
-        __host__ __device__ cuvect3f operator/(const float &rhs);
-        __host__ __device__ friend cuvect3f operator-(const cuvect3f &rhs);
+        __host__ __device__ cuvec3f operator+(const cuvec3f &rhs);
+        __host__ __device__ cuvec3f operator-(const cuvec3f &rhs);
+        __host__ __device__ cuvec3f operator*(const float &rhs);
+        __host__ __device__ cuvec3f operator/(const float &rhs);
+        __host__ __device__ friend cuvec3f operator-(const cuvec3f &rhs);
 
-        __host__ __device__ cuvect3f& operator+=(const cuvect3f &rhs);
-        __host__ __device__ cuvect3f& operator-=(const cuvect3f &rhs);
-        __host__ __device__ cuvect3f& operator/=(const float &rhs);
-        __host__ __device__ cuvect3f& operator*=(const float &rhs);
+        __host__ __device__ cuvec3f& operator+=(const cuvec3f &rhs);
+        __host__ __device__ cuvec3f& operator-=(const cuvec3f &rhs);
+        __host__ __device__ cuvec3f& operator/=(const float &rhs);
+        __host__ __device__ cuvec3f& operator*=(const float &rhs);
 
     };
 
@@ -59,7 +59,7 @@ namespace amscuda
         __host__ __device__ cumat3f operator-(const cumat3f &rhs);
         __host__ __device__ cumat3f operator*(const float &rhs);
         __host__ __device__ cumat3f operator/(const float &rhs);
-        __host__ __device__ cuvect3f operator*(const cuvect3f &rhs);
+        __host__ __device__ cuvec3f operator*(const cuvec3f &rhs);
         __host__ __device__ cumat3f operator*(const cumat3f &rhs);
         __host__ __device__ friend cumat3f operator-(const cumat3f &rhs);
         
@@ -79,15 +79,15 @@ namespace amscuda
         
     };
 
-    __host__ __device__ float cuvect3f_dot(const cuvect3f &a,const cuvect3f &b);
-    __host__ __device__ cuvect3f cuvect3f_cross(const cuvect3f &a, const cuvect3f &b);
-    __host__ __device__ float cuvect3f_norm(const cuvect3f &a);
-    __host__ __device__ cuvect3f cuvect3f_normalize(const cuvect3f &a);
-    __host__ __device__ cuvect3f cuvect3f_proj(const cuvect3f &a, const cuvect3f &b);
+    __host__ __device__ float cuvec3f_dot(const cuvec3f &a,const cuvec3f &b);
+    __host__ __device__ cuvec3f cuvec3f_cross(const cuvec3f &a, const cuvec3f &b);
+    __host__ __device__ float cuvec3f_norm(const cuvec3f &a);
+    __host__ __device__ cuvec3f cuvec3f_normalize(const cuvec3f &a);
+    __host__ __device__ cuvec3f cuvec3f_proj(const cuvec3f &a, const cuvec3f &b);
 
-    __host__ __device__ cumat3f hodge_dual(const cuvect3f &vin);
-    __host__ __device__ cuvect3f hodge_dual(const cumat3f &min);
-    __host__ __device__ cumat3f rotmat_from_axisangle(const cuvect3f &axis, const float &angle);
+    __host__ __device__ cumat3f hodge_dual(const cuvec3f &vin);
+    __host__ __device__ cuvec3f hodge_dual(const cumat3f &min);
+    __host__ __device__ cumat3f rotmat_from_axisangle(const cuvec3f &axis, const float &angle);
 
 
     //3x3 matrix operations
@@ -105,14 +105,14 @@ namespace amscuda
     //inverts a 3x3 (9 element) matrix
     __host__ __device__ void mat3f_inverse(float *mat3inout);
 
-    __host__ __device__ cuvect3f mat3f_mult(float *mat3in, const cuvect3f &cvin);
+    __host__ __device__ cuvec3f mat3f_mult(float *mat3in, const cuvec3f &cvin);
     __host__ __device__ void mat3f_mult(float *matina, float *matinb, float *matout);
 
-    __host__ __device__ void mat3f_hodgedual(const cuvect3f &vecin, float *matout);
-    __host__ __device__ void mat3f_hodgedual(float *matin, cuvect3f &vecout);
+    __host__ __device__ void mat3f_hodgedual(const cuvec3f &vecin, float *matout);
+    __host__ __device__ void mat3f_hodgedual(float *matin, cuvec3f &vecout);
 
     //returns direction cosine rotation matrix from axis and angle
-    __host__ __device__ void mat3f_rot_from_axisangle(cuvect3f axis, float angle, float *matout);
+    __host__ __device__ void mat3f_rot_from_axisangle(cuvec3f axis, float angle, float *matout);
 
     __host__ void test_cudavectf_logic1();
 
diff --git a/old/9apr26_prerefactor/include/amsculib2/cuvect4f.hpp b/old/9apr26_prerefactor/include/amsculib2/cuvect4f.hpp
index e70b2f9..cb06685 100644
--- a/old/9apr26_prerefactor/include/amsculib2/cuvect4f.hpp
+++ b/old/9apr26_prerefactor/include/amsculib2/cuvect4f.hpp
@@ -1,10 +1,10 @@
-#ifndef __CUVECT4F_HPP__
-#define __CUVECT4F_HPP__
+#ifndef __cuvec4f_HPP__
+#define __cuvec4f_HPP__
 
 namespace amscuda
 {
 
-    class cuvect4f
+    class cuvec4f
     {
         public:
         float x;
@@ -12,18 +12,18 @@ namespace amscuda
         float z;
         float w;
 
-        __host__ __device__ cuvect4f();
-        __host__ __device__ ~cuvect4f();
-        __host__ __device__ cuvect4f(float _x, float _y, float _z, float _w);
+        __host__ __device__ cuvec4f();
+        __host__ __device__ ~cuvec4f();
+        __host__ __device__ cuvec4f(float _x, float _y, float _z, float _w);
                 
         __host__ __device__ float& operator[](const int I);
         __host__ __device__ const float& operator[](const int I) const;
 
-        __host__ __device__ cuvect4f operator+(cuvect4f lhs);
-        __host__ __device__ cuvect4f operator-(cuvect4f lhs);
-        __host__ __device__ cuvect4f operator*(float lhs);
-        __host__ __device__ cuvect4f operator/(float lhs);
-        __host__ __device__ friend cuvect4f operator-(cuvect4f rhs);
+        __host__ __device__ cuvec4f operator+(cuvec4f lhs);
+        __host__ __device__ cuvec4f operator-(cuvec4f lhs);
+        __host__ __device__ cuvec4f operator*(float lhs);
+        __host__ __device__ cuvec4f operator/(float lhs);
+        __host__ __device__ friend cuvec4f operator-(cuvec4f rhs);
     };
 
     class cumat4f
@@ -41,7 +41,7 @@ namespace amscuda
         __host__ __device__ cumat4f operator-(cumat4f lhs);
         __host__ __device__ cumat4f operator*(float lhs);
         __host__ __device__ cumat4f operator/(float lhs);
-        __host__ __device__ cuvect4f operator*(cuvect4f lhs);
+        __host__ __device__ cuvec4f operator*(cuvec4f lhs);
         __host__ __device__ cumat4f operator*(cumat4f lhs);
         __host__ __device__ friend cumat4f operator-(cumat4f rhs);
 
@@ -50,10 +50,10 @@ namespace amscuda
         __host__ __device__ cumat4f inverse();
     };
 
-    __host__ __device__ float cuvect4f_dot(cuvect4f a, cuvect4f b);
-    __host__ __device__ float cuvect4f_norm(cuvect4f a);
-    __host__ __device__ cuvect4f cuvect4f_normalize(cuvect4f a);
-    __host__ __device__ cuvect4f cuvect4f_proj(cuvect4f a, cuvect4f b);
+    __host__ __device__ float cuvec4f_dot(cuvec4f a, cuvec4f b);
+    __host__ __device__ float cuvec4f_norm(cuvec4f a);
+    __host__ __device__ cuvec4f cuvec4f_normalize(cuvec4f a);
+    __host__ __device__ cuvec4f cuvec4f_proj(cuvec4f a, cuvec4f b);
 
 
 };
diff --git a/old/9apr26_prerefactor/src/amsculib2/cuvect2f.cu b/old/9apr26_prerefactor/src/amsculib2/cuvect2f.cu
index 9b53d79..d18e1b0 100644
--- a/old/9apr26_prerefactor/src/amsculib2/cuvect2f.cu
+++ b/old/9apr26_prerefactor/src/amsculib2/cuvect2f.cu
@@ -3,91 +3,91 @@
 namespace amscuda
 {
 
-    __host__ __device__ cuvect2f::cuvect2f()
+    __host__ __device__ cuvec2f::cuvec2f()
     {
         x = 0.0; y = 0.0;
         return;
     }
 
-    __host__ __device__ cuvect2f::~cuvect2f()
+    __host__ __device__ cuvec2f::~cuvec2f()
     {
         x = 0.0; y = 0.0;
         return;
     }
 
-    __host__ __device__ cuvect2f::cuvect2f(const float &_x, const float &_y)
+    __host__ __device__ cuvec2f::cuvec2f(const float &_x, const float &_y)
     {
         x = _x; y = _y;
         return;
     }
 
-    __host__ __device__ float& cuvect2f::operator[](const int &I)
+    __host__ __device__ float& cuvec2f::operator[](const int &I)
     {
         if(I==0) return x;
         if(I==1) return y;
         return x;
     }
 
-    __host__ __device__ const float& cuvect2f::operator[](const int &I) const
+    __host__ __device__ const float& cuvec2f::operator[](const int &I) const
     {
         if(I==0) return x;
         if(I==1) return y;
         return x;
     }
 
-    __host__ __device__ cuvect2f cuvect2f::operator+(const cuvect2f &lhs)
+    __host__ __device__ cuvec2f cuvec2f::operator+(const cuvec2f &lhs)
     {
-        cuvect2f ret;
+        cuvec2f ret;
         ret.x = x + lhs.x;
         ret.y = y + lhs.y;
         return ret;
     }
 
-    __host__ __device__ cuvect2f cuvect2f::operator-(const cuvect2f &lhs)
+    __host__ __device__ cuvec2f cuvec2f::operator-(const cuvec2f &lhs)
     {
-        cuvect2f ret;
+        cuvec2f ret;
         ret.x = x - lhs.x;
         ret.y = y - lhs.y;
         return ret;
     }
-    __host__ __device__ cuvect2f cuvect2f::operator*(const float &lhs)
+    __host__ __device__ cuvec2f cuvec2f::operator*(const float &lhs)
     {
-        cuvect2f ret;
+        cuvec2f ret;
         ret.x = x*lhs;
         ret.y = y*lhs;
         return ret;
     }
 
-    __host__ __device__ cuvect2f cuvect2f::operator/(const float &lhs)
+    __host__ __device__ cuvec2f cuvec2f::operator/(const float &lhs)
     {
-        cuvect2f ret;
+        cuvec2f ret;
         ret.x = x/lhs;
         ret.y = y/lhs;
         return ret;
     }
 
-    __host__ __device__ float cuvect2f_dot(const cuvect2f &a, const cuvect2f &b)
+    __host__ __device__ float cuvec2f_dot(const cuvec2f &a, const cuvec2f &b)
     {
         float ret = a.x*b.x+a.y*b.y;
         return ret;
     }
 
-    __host__ __device__ float cuvect2f_cross(const cuvect2f &a, const cuvect2f &b)
+    __host__ __device__ float cuvec2f_cross(const cuvec2f &a, const cuvec2f &b)
     {
         float ret = a.x*b.y-a.y*b.x;
         return ret;
     }
 
-    __host__ __device__ float cuvect2f_norm(const cuvect2f &a)
+    __host__ __device__ float cuvec2f_norm(const cuvec2f &a)
     {
         float ret = ::sqrtf(a.x*a.x+a.y*a.y);
         return ret;
     }
 
-    __host__ __device__ cuvect2f cuvect2f_normalize(const cuvect2f &a)
+    __host__ __device__ cuvec2f cuvec2f_normalize(const cuvec2f &a)
     {
-        cuvect2f ret;
-        float m = cuvect2f_norm(a);
+        cuvec2f ret;
+        float m = cuvec2f_norm(a);
         if(m>0.0)
         {
             ret.x = a.x/m; ret.y = a.y/m;
@@ -99,11 +99,11 @@ namespace amscuda
         return ret;
     }
 
-    __host__ __device__ cuvect2f cuvect2f_proj(const cuvect2f &a, const cuvect2f &b)
+    __host__ __device__ cuvec2f cuvec2f_proj(const cuvec2f &a, const cuvec2f &b)
     {
-        cuvect2f ret;
-        cuvect2f bn = cuvect2f_normalize(b);
-        float m = cuvect2f_dot(a,bn);
+        cuvec2f ret;
+        cuvec2f bn = cuvec2f_normalize(b);
+        float m = cuvec2f_dot(a,bn);
         ret = bn*m;
         return ret;
     }
@@ -189,9 +189,9 @@ __host__ __device__ float& cumat2f::at(const int I, const int J)
 }
         
 
-__host__ __device__ cuvect2f cumat2f::operator*(cuvect2f lhs)
+__host__ __device__ cuvec2f cumat2f::operator*(cuvec2f lhs)
 {
-    cuvect2f ret = cuvect2f(0.0,0.0);
+    cuvec2f ret = cuvec2f(0.0,0.0);
     int I,J;
     for(I=0;I<2;I++)
     {
@@ -343,17 +343,17 @@ __host__ __device__ void mat2f_mult(float *mat2a, float *mat2b, float *mat2c)
 }
 
 // ret = a*b
-__host__ __device__ cuvect2f mat2f_mult(float *mat2a, const cuvect2f &b)
+__host__ __device__ cuvec2f mat2f_mult(float *mat2a, const cuvec2f &b)
 {
-    cuvect2f ret;
+    cuvec2f ret;
     ret.x = b.x*mat2a[0+0*2] + b.y*mat2a[1+0*2];
     ret.y = b.x*mat2a[0+1*2] + b.y*mat2a[1+1*2];
     return ret;
 }
 
-__host__ __device__ cuvect2f operator-(cuvect2f rhs)
+__host__ __device__ cuvec2f operator-(cuvec2f rhs)
 {
-    cuvect2f ret;
+    cuvec2f ret;
     ret[0] = -rhs[0];
     ret[1] = -rhs[1];
     return ret;
@@ -367,7 +367,7 @@ __host__ __device__ cumat2f operator-(cumat2f rhs)
     return ret;
 }
 
-void test_cuvect2f_1()
+void test_cuvec2f_1()
 {
 
 
diff --git a/old/9apr26_prerefactor/src/amsculib2/cuvect3f.cu b/old/9apr26_prerefactor/src/amsculib2/cuvect3f.cu
index 309d6bb..a92c3d5 100644
--- a/old/9apr26_prerefactor/src/amsculib2/cuvect3f.cu
+++ b/old/9apr26_prerefactor/src/amsculib2/cuvect3f.cu
@@ -3,19 +3,19 @@
 namespace amscuda
 {
 
-    __host__ __device__ cuvect3f::cuvect3f()
+    __host__ __device__ cuvec3f::cuvec3f()
     {
         x = 0.0; y = 0.0; z = 0.0;
         return;
     }
 
-    __host__ __device__ cuvect3f::~cuvect3f()
+    __host__ __device__ cuvec3f::~cuvec3f()
     {
         x = 0.0; y = 0.0; z = 0.0;
         return;
     }
 
-    __host__ __device__ float& cuvect3f::operator[](const int &I)
+    __host__ __device__ float& cuvec3f::operator[](const int &I)
     {
         if(I==0) return x;
         if(I==1) return y;
@@ -23,7 +23,7 @@ namespace amscuda
         return x;
     }
 
-    __host__ __device__ const float& cuvect3f::operator[](const int &I) const
+    __host__ __device__ const float& cuvec3f::operator[](const int &I) const
     {
         if(I==0) return x;
         if(I==1) return y;
@@ -31,9 +31,9 @@ namespace amscuda
         return x;
     }
 
-    __host__ __device__ cuvect3f cuvect3f::operator+(const cuvect3f &rhs)
+    __host__ __device__ cuvec3f cuvec3f::operator+(const cuvec3f &rhs)
     {
-        cuvect3f ret;
+        cuvec3f ret;
         ret.x = x+rhs.x;
         ret.y = y+rhs.y;
         ret.z = z+rhs.z;
@@ -41,9 +41,9 @@ namespace amscuda
         return ret;
     }
 
-    __host__ __device__ cuvect3f cuvect3f::operator-(const cuvect3f &rhs)
+    __host__ __device__ cuvec3f cuvec3f::operator-(const cuvec3f &rhs)
     {
-        cuvect3f ret;
+        cuvec3f ret;
         ret.x = x-rhs.x;
         ret.y = y-rhs.y;
         ret.z = z-rhs.z;
@@ -51,25 +51,25 @@ namespace amscuda
         return ret;
     }
 
-    __host__ __device__ cuvect3f cuvect3f::operator*(const float &rhs)
+    __host__ __device__ cuvec3f cuvec3f::operator*(const float &rhs)
     {
-        cuvect3f ret;
+        cuvec3f ret;
         ret.x = x*rhs;
         ret.y = y*rhs;
         ret.z = z*rhs;
         return ret;
     }
 
-    __host__ __device__ cuvect3f cuvect3f::operator/(const float &rhs)
+    __host__ __device__ cuvec3f cuvec3f::operator/(const float &rhs)
     {
-        cuvect3f ret;
+        cuvec3f ret;
         ret.x = x/rhs;
         ret.y = y/rhs;
         ret.z = z/rhs;
         return ret;
     }
 
-    __host__ __device__ cuvect3f& cuvect3f::operator+=(const cuvect3f &rhs)
+    __host__ __device__ cuvec3f& cuvec3f::operator+=(const cuvec3f &rhs)
     {
         x = x + rhs.x;
         y = y + rhs.y;
@@ -77,7 +77,7 @@ namespace amscuda
         return *this;
     }
 
-    __host__ __device__ cuvect3f& cuvect3f::operator-=(const cuvect3f &rhs)
+    __host__ __device__ cuvec3f& cuvec3f::operator-=(const cuvec3f &rhs)
     {
         x = x - rhs.x;
         y = y - rhs.y;
@@ -85,7 +85,7 @@ namespace amscuda
         return *this;
     }
 
-    __host__ __device__ cuvect3f& cuvect3f::operator*=(const float &rhs)
+    __host__ __device__ cuvec3f& cuvec3f::operator*=(const float &rhs)
     {
         x = x * rhs;
         y = y * rhs;
@@ -93,7 +93,7 @@ namespace amscuda
         return *this;
     }
 
-    __host__ __device__ cuvect3f& cuvect3f::operator/=(const float &rhs)
+    __host__ __device__ cuvec3f& cuvec3f::operator/=(const float &rhs)
     {
         x = x / rhs;
         y = y / rhs;
@@ -102,23 +102,23 @@ namespace amscuda
     }
 
 
-    __host__ __device__ cuvect3f::cuvect3f(const float &_x, const float &_y, const float &_z)
+    __host__ __device__ cuvec3f::cuvec3f(const float &_x, const float &_y, const float &_z)
     {
         x = _x; y = _y; z = _z;
         return;
     }
 
 
-    __host__ __device__ float cuvect3f_dot(const cuvect3f &a, const cuvect3f &b)
+    __host__ __device__ float cuvec3f_dot(const cuvec3f &a, const cuvec3f &b)
     {
         float ret = a.x*b.x+a.y*b.y+a.z*b.z;
 
         return ret;
     }
 
-    __host__ __device__ cuvect3f cuvect3f_cross(const cuvect3f &a, const cuvect3f &b)
+    __host__ __device__ cuvec3f cuvec3f_cross(const cuvec3f &a, const cuvec3f &b)
     {
-        cuvect3f ret;
+        cuvec3f ret;
         ret[0] = a[1]*b[2]-a[2]*b[1];
         ret[1] = a[2]*b[0]-a[0]*b[2];
         ret[2] = a[0]*b[1]-a[1]*b[0];
@@ -126,16 +126,16 @@ namespace amscuda
         return ret;
     }
 
-    __host__ __device__ float cuvect3f_norm(const cuvect3f &a)
+    __host__ __device__ float cuvec3f_norm(const cuvec3f &a)
     {
         float ret;
         ret = ::sqrtf(a.x*a.x+a.y*a.y+a.z*a.z);
         return ret;
     }
 
-    __host__ __device__ cuvect3f cuvect3f_normalize(const cuvect3f &a)
+    __host__ __device__ cuvec3f cuvec3f_normalize(const cuvec3f &a)
     {
-        cuvect3f ret;
+        cuvec3f ret;
         float m;
         m = ::sqrtf(a.x*a.x+a.y*a.y+a.z*a.z);
         if(m>0.0)
@@ -150,11 +150,11 @@ namespace amscuda
         return ret;
     }
 
-    __host__ __device__ cuvect3f cuvect3f_proj(const cuvect3f &a, const cuvect3f &b)
+    __host__ __device__ cuvec3f cuvec3f_proj(const cuvec3f &a, const cuvec3f &b)
     {
-        cuvect3f ret;
-        cuvect3f bn = cuvect3f_normalize(b);
-        float m = cuvect3f_dot(a,bn);
+        cuvec3f ret;
+        cuvec3f bn = cuvec3f_normalize(b);
+        float m = cuvec3f_dot(a,bn);
         ret = bn*m;
         return ret;
     }
@@ -346,9 +346,9 @@ __host__ __device__ cumat3f cumat3f::operator/(const float &rhs)
     return ret;
 }
         
-__host__ __device__ cuvect3f cumat3f::operator*(const cuvect3f &rhs)
+__host__ __device__ cuvec3f cumat3f::operator*(const cuvec3f &rhs)
 {
-    cuvect3f ret;
+    cuvec3f ret;
     
     ret.x = m00*rhs.x + m01*rhs.y + m02*rhs.z;
     ret.y = m10*rhs.x + m11*rhs.y + m12*rhs.z;
@@ -649,10 +649,10 @@ __host__ __device__ void mat3f_inverse(float *mat3inout)
     return;
 }
 
-__host__ __device__ cuvect3f mat3f_mult(float *mat3in, const cuvect3f &cvin)
+__host__ __device__ cuvec3f mat3f_mult(float *mat3in, const cuvec3f &cvin)
 {
     int I,J;
-    cuvect3f ret;
+    cuvec3f ret;
     for(I=0;I<3;I++)
     {
         ret[I] = 0.0;
@@ -700,7 +700,7 @@ __host__ __device__ void mat3f_mult(float *matina, float *matinb, float *matout)
     return;
 }
 
-__host__ __device__ void mat3f_hodgedual(const cuvect3f &vecin, float *matout)
+__host__ __device__ void mat3f_hodgedual(const cuvec3f &vecin, float *matout)
 {
     matout[0 + 0*3] = 0.0f;
     matout[1 + 0*3] = -vecin[2];
@@ -716,7 +716,7 @@ __host__ __device__ void mat3f_hodgedual(const cuvect3f &vecin, float *matout)
     return;
 }
 
-__host__ __device__ void mat3f_hodgedual(float *matin, cuvect3f &vecout)
+__host__ __device__ void mat3f_hodgedual(float *matin, cuvec3f &vecout)
 {
     vecout[0] = 0.5*(matin[1 + 2*3] - matin[2 + 1*3]);
     vecout[1] = 0.5*(matin[2 + 0*3] - matin[0 + 2*3]);
@@ -726,7 +726,7 @@ __host__ __device__ void mat3f_hodgedual(float *matin, cuvect3f &vecout)
 }
 
 //returns direction cosine rotation matrix from axis and angle
-__host__ __device__ void mat3f_rot_from_axisangle(cuvect3f axis, float angle, float *matout)
+__host__ __device__ void mat3f_rot_from_axisangle(cuvec3f axis, float angle, float *matout)
 {
     int I;
     float H[9];
@@ -738,7 +738,7 @@ __host__ __device__ void mat3f_rot_from_axisangle(cuvect3f axis, float angle, fl
     II[1+1*3] = 1.0;
     II[2+2*3] = 1.0;
 
-    axis = cuvect3f_normalize(axis);
+    axis = cuvec3f_normalize(axis);
     
     mat3f_hodgedual(axis,H);
     mat3f_mult(H,H,Hsq);
@@ -759,7 +759,7 @@ __host__ void test_cudavectf_logic1()
 
     printf("3 dim vector and matrix functional tests on host side\n");
     
-    cuvect3f a,b,c;
+    cuvec3f a,b,c;
     float ma[9],mb[9],mc[9];
 
     int I,J;
@@ -798,7 +798,7 @@ __host__ void test_cudavectf_logic1()
         }
     }
 
-    a = cuvect3f(1,1,1);
+    a = cuvec3f(1,1,1);
     b = mat3f_mult(ma,a);
     b = mat3f_mult(mb,b);
     
@@ -807,8 +807,8 @@ __host__ void test_cudavectf_logic1()
         printf("a[%d] = %1.3f, b[%d] = %1.3f\n",I,a[I],I,b[I]);
     }
 
-    a = cuvect3f(1,0,1);
-    b = cuvect3f(0,1,-1);
+    a = cuvec3f(1,0,1);
+    b = cuvec3f(0,1,-1);
     c = a+b;
 
     for(I=0;I<3;I++)
@@ -823,31 +823,31 @@ __host__ void test_cudavectf_logic1()
         printf("a[%d] = %1.3f, b[%d] = %1.3f, c[%d] = %1.3f\n",I,a[I],I,b[I],I,c[I]);
     }
 
-    c = cuvect3f_cross(a,b);
+    c = cuvec3f_cross(a,b);
 
     for(I=0;I<3;I++)
     {
         printf("a[%d] = %1.3f, b[%d] = %1.3f, c[%d] = %1.3f\n",I,a[I],I,b[I],I,c[I]);
     }
 
-    printf("c dot a = %1.3f, c dot b = %1.3f\n",cuvect3f_dot(c,a),cuvect3f_dot(c,b));
+    printf("c dot a = %1.3f, c dot b = %1.3f\n",cuvec3f_dot(c,a),cuvec3f_dot(c,b));
 
-    printf("norm(a)=%1.3f, norm(b)=%1.3f, norm(c)=%1.3f\n",cuvect3f_norm(a),cuvect3f_norm(b),cuvect3f_norm(c));
-    a = cuvect3f_normalize(a);
-    b = cuvect3f_normalize(b);
-    c = cuvect3f_normalize(c);
+    printf("norm(a)=%1.3f, norm(b)=%1.3f, norm(c)=%1.3f\n",cuvec3f_norm(a),cuvec3f_norm(b),cuvec3f_norm(c));
+    a = cuvec3f_normalize(a);
+    b = cuvec3f_normalize(b);
+    c = cuvec3f_normalize(c);
 
     for(I=0;I<3;I++)
     {
         printf("a[%d] = %1.3f, b[%d] = %1.3f, c[%d] = %1.3f\n",I,a[I],I,b[I],I,c[I]);
     }
-    printf("c dot a = %1.3f, c dot b = %1.3f\n",cuvect3f_dot(c,a),cuvect3f_dot(c,b));
-    printf("norm(a)=%1.3f, norm(b)=%1.3f, norm(c)=%1.3f\n",cuvect3f_norm(a),cuvect3f_norm(b),cuvect3f_norm(c));
+    printf("c dot a = %1.3f, c dot b = %1.3f\n",cuvec3f_dot(c,a),cuvec3f_dot(c,b));
+    printf("norm(a)=%1.3f, norm(b)=%1.3f, norm(c)=%1.3f\n",cuvec3f_norm(a),cuvec3f_norm(b),cuvec3f_norm(c));
 
     return;
 }
 
-__host__ __device__ cumat3f hodge_dual(const cuvect3f &vin)
+__host__ __device__ cumat3f hodge_dual(const cuvec3f &vin)
 {
     cumat3f ret;
 
@@ -866,9 +866,9 @@ __host__ __device__ cumat3f hodge_dual(const cuvect3f &vin)
     return ret;
 }
 
-__host__ __device__ cuvect3f hodge_dual(const cumat3f &min)
+__host__ __device__ cuvec3f hodge_dual(const cumat3f &min)
 {
-    cuvect3f ret;
+    cuvec3f ret;
 
     ret.x = 0.5f*(min.m12 - min.m21);
     ret.y = 0.5f*(min.m20 - min.m02);
@@ -893,13 +893,13 @@ __host__ __device__ const cumat3f cumat3f_zeros()
     return ret;
 }
 
-__host__ __device__ cumat3f rotmat_from_axisangle(const cuvect3f &axis, const float &angle)
+__host__ __device__ cumat3f rotmat_from_axisangle(const cuvec3f &axis, const float &angle)
 {
     cumat3f ret = cumat3f_zeros();
     cumat3f H;
-    cuvect3f _naxis;
+    cuvec3f _naxis;
 
-    _naxis = cuvect3f_normalize(axis);
+    _naxis = cuvec3f_normalize(axis);
     H = hodge_dual(_naxis);
 
     ret += H*sinf(angle);
diff --git a/old/9apr26_prerefactor/src/amsculib2/cuvect4f.cu b/old/9apr26_prerefactor/src/amsculib2/cuvect4f.cu
index 3718341..9d3a182 100644
--- a/old/9apr26_prerefactor/src/amsculib2/cuvect4f.cu
+++ b/old/9apr26_prerefactor/src/amsculib2/cuvect4f.cu
@@ -4,28 +4,28 @@ namespace amscuda
 {
 
 ////////////
-//cuvect4ff//
+//cuvec4ff//
 ////////////
 
-__host__ __device__ cuvect4f::cuvect4f()
+__host__ __device__ cuvec4f::cuvec4f()
 {
     x = 0.0; y = 0.0; z = 0.0; w = 0.0;
     return;
 }
 
-__host__ __device__ cuvect4f::~cuvect4f()
+__host__ __device__ cuvec4f::~cuvec4f()
 {
     x = 0.0; y = 0.0; z = 0.0; w = 0.0;
     return;
 }
 
-__host__ __device__ cuvect4f::cuvect4f(float _x, float _y, float _z, float _w)
+__host__ __device__ cuvec4f::cuvec4f(float _x, float _y, float _z, float _w)
 {
     x = _x; y = _y; z = _z; w = _w;
     return;
 }
         
-__host__ __device__ float& cuvect4f::operator[](const int I)
+__host__ __device__ float& cuvec4f::operator[](const int I)
 {
     if(I==0) return x;
     else if(I==1) return y;
@@ -34,7 +34,7 @@ __host__ __device__ float& cuvect4f::operator[](const int I)
     return x;
 }
 
-__host__ __device__ const float& cuvect4f::operator[](const int I) const
+__host__ __device__ const float& cuvec4f::operator[](const int I) const
 {
     if(I==0) return x;
     else if(I==1) return y;
@@ -43,9 +43,9 @@ __host__ __device__ const float& cuvect4f::operator[](const int I) const
     return x;
 }
 
-__host__ __device__ cuvect4f cuvect4f::operator+(cuvect4f lhs)
+__host__ __device__ cuvec4f cuvec4f::operator+(cuvec4f lhs)
 {
-    cuvect4f ret;
+    cuvec4f ret;
     ret.x = this->x + lhs.x;
     ret.y = this->y + lhs.y;
     ret.z = this->z + lhs.z;
@@ -53,9 +53,9 @@ __host__ __device__ cuvect4f cuvect4f::operator+(cuvect4f lhs)
     return ret;
 }
 
-__host__ __device__ cuvect4f cuvect4f::operator-(cuvect4f lhs)
+__host__ __device__ cuvec4f cuvec4f::operator-(cuvec4f lhs)
 {
-    cuvect4f ret;
+    cuvec4f ret;
     ret.x = this->x - lhs.x;
     ret.y = this->y - lhs.y;
     ret.z = this->z - lhs.z;
@@ -63,9 +63,9 @@ __host__ __device__ cuvect4f cuvect4f::operator-(cuvect4f lhs)
     return ret;
 }
 
-__host__ __device__ cuvect4f cuvect4f::operator*(float lhs)
+__host__ __device__ cuvec4f cuvec4f::operator*(float lhs)
 {
-    cuvect4f ret;
+    cuvec4f ret;
     ret.x = this->x*lhs;
     ret.y = this->y*lhs;
     ret.z = this->z*lhs;
@@ -73,9 +73,9 @@ __host__ __device__ cuvect4f cuvect4f::operator*(float lhs)
     return ret;
 }
 
-__host__ __device__ cuvect4f cuvect4f::operator/(float lhs)
+__host__ __device__ cuvec4f cuvec4f::operator/(float lhs)
 {
-    cuvect4f ret;
+    cuvec4f ret;
     ret.x = this->x/lhs;
     ret.y = this->y/lhs;
     ret.z = this->z/lhs;
@@ -162,9 +162,9 @@ __host__ __device__ cumat4f cumat4f::operator/(float lhs)
     return ret;
 }
 
-__host__ __device__ cuvect4f cumat4f::operator*(cuvect4f lhs)
+__host__ __device__ cuvec4f cumat4f::operator*(cuvec4f lhs)
 {
-    cuvect4f ret = cuvect4f(0.0,0.0,0.0,0.0);
+    cuvec4f ret = cuvec4f(0.0,0.0,0.0,0.0);
     int I,J;
 
     for(I=0;I<4;I++)
@@ -382,41 +382,41 @@ __host__ __device__ cumat4f cumat4f::inverse()
 }
 
 
-__host__ __device__ float cuvect4f_dot(cuvect4f a, cuvect4f b)
+__host__ __device__ float cuvec4f_dot(cuvec4f a, cuvec4f b)
 {
     float ret = 0.0;
     ret = a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w;
     return ret;
 }
 
-__host__ __device__ float cuvect4f_norm(cuvect4f a)
+__host__ __device__ float cuvec4f_norm(cuvec4f a)
 {
     float ret = 0.0;
-    ret = ::sqrtf(cuvect4f_dot(a,a));
+    ret = ::sqrtf(cuvec4f_dot(a,a));
     return ret;
 }
 
-__host__ __device__ cuvect4f cuvect4f_normalize(cuvect4f a)
+__host__ __device__ cuvec4f cuvec4f_normalize(cuvec4f a)
 {
-    cuvect4f ret = cuvect4f(0.0f,0.0f,0.0f,0.0f);
-    float nrm = cuvect4f_norm(a);
+    cuvec4f ret = cuvec4f(0.0f,0.0f,0.0f,0.0f);
+    float nrm = cuvec4f_norm(a);
     if(nrm>0.0)
         ret = a/nrm;
     return ret;
 }
 
-__host__ __device__ cuvect4f cuvect4f_proj(cuvect4f a, cuvect4f b)
+__host__ __device__ cuvec4f cuvec4f_proj(cuvec4f a, cuvec4f b)
 {
-    cuvect4f ret;
-    cuvect4f bn = cuvect4f_normalize(b);
-    float d = cuvect4f_dot(a,bn);
+    cuvec4f ret;
+    cuvec4f bn = cuvec4f_normalize(b);
+    float d = cuvec4f_dot(a,bn);
     ret = bn*d;
     return ret;
 }
 
-__host__ __device__ cuvect4f operator-(cuvect4f rhs)
+__host__ __device__ cuvec4f operator-(cuvec4f rhs)
 {
-    cuvect4f ret;
+    cuvec4f ret;
     ret[0] = -rhs[0];
     ret[1] = -rhs[1];
     ret[2] = -rhs[2];
diff --git a/src/amsculib3/amscu_random.cu b/src/amsculib3/amscu_random.cu
index 5e363ad..8fabb1f 100644
--- a/src/amsculib3/amscu_random.cu
+++ b/src/amsculib3/amscu_random.cu
@@ -163,59 +163,59 @@ __host__ __device__ float dpr64_randf(int64_t *seedinout)
 
 void test_dprg64()
 {
-    printf("Tests for dprg:\n");
-    long I;
-    int64_t seed = 133LL;
-    double d;
-    float f;
-    cuvect3 qv;
+    // printf("Tests for dprg:\n");
+    // long I;
+    // int64_t seed = 133LL;
+    // double d;
+    // float f;
+    // cuvec3 qv;
 
-    printf("dpr64_randd test\n");
-    seed = 133LL;
-    for(I=0;I<10;I++)
-    {
-        d = dpr64_randd(&seed);
-        printf("seed: %lld rand: %1.4f\n",(long long)seed,d);
-    }
+    // printf("dpr64_randd test\n");
+    // seed = 133LL;
+    // for(I=0;I<10;I++)
+    // {
+    //     d = dpr64_randd(&seed);
+    //     printf("seed: %lld rand: %1.4f\n",(long long)seed,d);
+    // }
 
-    printf("\n\n");
-    printf("dpr64_randf test\n");
-    seed = 133LL;
-    for(I=0;I<10;I++)
-    {
-        f = dpr64_randf(&seed);
-        printf("seed: %lld rand: %1.4f\n",(long long)seed,f);
-    }
+    // printf("\n\n");
+    // printf("dpr64_randf test\n");
+    // seed = 133LL;
+    // for(I=0;I<10;I++)
+    // {
+    //     f = dpr64_randf(&seed);
+    //     printf("seed: %lld rand: %1.4f\n",(long long)seed,f);
+    // }
 
-    return;
+    // return;
 }
 
 void test_dprg32()
 {
-    printf("Tests for dprg:\n");
-    long I;
-    int32_t seed = 133;
-    double d;
-    float f;
-    cuvect3 qv;
+    // printf("Tests for dprg:\n");
+    // long I;
+    // int32_t seed = 133;
+    // double d;
+    // float f;
+    // cuvec3 qv;
 
-    printf("dpr32_randf test\n");
-    seed = 133;
-    for(I=0;I<10;I++)
-    {
-        f = dpr32_randf(&seed);
-        printf("seed: %lld rand: %1.4f\n",(long long)seed,f);
-    }
+    // printf("dpr32_randf test\n");
+    // seed = 133;
+    // for(I=0;I<10;I++)
+    // {
+    //     f = dpr32_randf(&seed);
+    //     printf("seed: %lld rand: %1.4f\n",(long long)seed,f);
+    // }
 
-    printf("\n\ndpr32_randnf test\n");
-    seed = 133;
-    for(I=0;I<10;I++)
-    {
-        f = dpr32_randnf(&seed);
-        printf("seed: %lld rand: %1.4f\n",(long long)seed,f);
-    }
+    // printf("\n\ndpr32_randnf test\n");
+    // seed = 133;
+    // for(I=0;I<10;I++)
+    // {
+    //     f = dpr32_randnf(&seed);
+    //     printf("seed: %lld rand: %1.4f\n",(long long)seed,f);
+    // }
 
-    return;
+    // return;
 }
 
 
diff --git a/src/amsculib3/amscugeom.cu b/src/amsculib3/geom/amscugeom.cu
similarity index 100%
rename from src/amsculib3/amscugeom.cu
rename to src/amsculib3/geom/amscugeom.cu
diff --git a/src/amsculib3/math/amscu_comp64.cu b/src/amsculib3/math/amscu_comp64.cu
index 1fa673c..e5789f5 100644
--- a/src/amsculib3/math/amscu_comp64.cu
+++ b/src/amsculib3/math/amscu_comp64.cu
@@ -7,15 +7,15 @@ namespace cmp
 
 __host__ __device__ cucomp64::cucomp64()
 {
-    real = 0.0;
-    imag = 0.0;
+    real = 0.0f;
+    imag = 0.0f;
     return;
 }
 
 __host__ __device__ cucomp64::~cucomp64()
 {
-    real = 0.0;
-    imag = 0.0;
+    real = 0.0f;
+    imag = 0.0f;
     return;
 }
 
@@ -30,7 +30,14 @@ __host__ __device__ cucomp64::cucomp64(const cucomp64 &other)
 __host__ __device__ cucomp64::cucomp64(const float &other)
 {
     real = other;
-    imag = 0.0;
+    imag = 0.0f;
+    return;
+}
+
+__host__ __device__ cucomp64::cucomp64(const float &_real, const float &_imag)
+{
+    real = _real;
+    imag = _imag;
     return;
 }
 
@@ -51,14 +58,14 @@ __host__ __device__ const cucomp64& cucomp64::operator=(const cucomp64& other)
 __host__ __device__ cucomp64& cucomp64::operator=(float& other)
 {
     real = other;
-    imag = 0.0;
+    imag = 0.0f;
     return *this;
 }
 
 __host__ __device__ const cucomp64& cucomp64::operator=(const float& other) 
 {
     this->real = other;
-    this->imag = 0.0;
+    this->imag = 0.0f;
     return *this;
 }
 
@@ -307,14 +314,21 @@ __host__ __device__ bool cucomp64::iszero() const
 
 __host__ __device__ float cucomp64::arg() const
 {
-    float ret = 0.0;
+    float ret = 0.0f;
     ret = (float) amscuda::arg((double)real,(double)imag);
     return ret;
 }
 
 __host__ __device__ float cucomp64::mag() const
 {
-    float ret = 0.0;
+    float ret = 0.0f;
+    ret = ::sqrt(real*real+imag*imag);
+    return ret;
+}
+
+__host__ __device__ float cucomp64::abs() const
+{
+    float ret = 0.0f;
     ret = ::sqrt(real*real+imag*imag);
     return ret;
 }
@@ -327,62 +341,54 @@ __host__ __device__ cucomp64 cucomp64::conj() const
     return ret;
 }
 
-__host__ __device__ float arg(cucomp64 z)
+__host__ __device__ float arg(const cucomp64 &z)
 {
     return z.arg();
 }
 
-__host__ __device__ float abs(cucomp64 z)
+__host__ __device__ float abs(const cucomp64 &z)
 {
     return z.mag();
 }
 
-__host__ __device__ cucomp64 dtocomp64(float _r, float _i)
-{
-    cucomp64 ret;
-    ret.real = _r;
-    ret.imag = _i;
-    return ret;
-}
-
-__host__ __device__ float real(cucomp64 z)
+__host__ __device__ float real(const cucomp64 &z)
 {
     return z.real;
 }
 
-__host__ __device__ float imag(cucomp64 z)
+__host__ __device__ float imag(const cucomp64 &z)
 {
     return z.imag;
 }
 
-__host__ __device__ cucomp64 sin(cucomp64 z)
+__host__ __device__ cucomp64 sin(const cucomp64 &z)
 {
     cucomp64 ret;
-    cucomp64 im1 = dtocomp64(0.0f,1.0f);
-    cucomp64 div = dtocomp64(0.0f,2.0f);
+    cucomp64 im1 = cucomp64(0.0f,1.0f);
+    cucomp64 div = cucomp64(0.0f,2.0f);
 
     ret = (exp(im1*z)-exp(-im1*z))/div;
 
     return ret;
 }
 
-__host__ __device__ cucomp64 cos(cucomp64 z)
+__host__ __device__ cucomp64 cos(const cucomp64 &z)
 {
     cucomp64 ret;
-    cucomp64 im1 = dtocomp64(0.0f,1.0f);
-    cucomp64 div = dtocomp64(2.0f,0.0f);
+    cucomp64 im1 = cucomp64(0.0f,1.0f);
+    cucomp64 div = cucomp64(2.0f,0.0f);
 
     ret = (exp(im1*z)+exp(-im1*z))/div;
 
     return ret;
 }
 
-__host__ __device__ cucomp64 tan(cucomp64 z)
+__host__ __device__ cucomp64 tan(const cucomp64 &z)
 {
     return sin(z)/cos(z);
 }
 
-__host__ __device__ cucomp64 exp(cucomp64 z)
+__host__ __device__ cucomp64 exp(const cucomp64 &z)
 {
     cucomp64 ret;
     ret.real = ::exp(z.real)*::cos(z.imag);
@@ -390,7 +396,7 @@ __host__ __device__ cucomp64 exp(cucomp64 z)
     return ret;
 }
 
-__host__ __device__ cucomp64 log(cucomp64 z)
+__host__ __device__ cucomp64 log(const cucomp64 &z)
 {
     cucomp64 ret;
     ret.real = ::log(::sqrt(z.real*z.real+z.imag*z.imag));
@@ -398,46 +404,103 @@ __host__ __device__ cucomp64 log(cucomp64 z)
     return ret;
 }
 
-__host__ __device__ cucomp64 conj(cucomp64 z)
+__host__ __device__ cucomp64 conj(const cucomp64 &z)
 {
     return z.conj();
 }
 
-__host__ __device__ cucomp64 cosh(cucomp64 z)
+__host__ __device__ cucomp64 cosh(const cucomp64 &z)
 {
     cucomp64 ret;
-    cucomp64 div = dtocomp64(2.0f,0.0f);
+    cucomp64 div = cucomp64(2.0f,0.0f);
 
     ret = (exp(z)+exp(-z))/div;
 
     return ret;
 }
 
-__host__ __device__ cucomp64 sinh(cucomp64 z)
+__host__ __device__ cucomp64 sinh(const cucomp64 &z)
 {
     cucomp64 ret;
-    cucomp64 div = dtocomp64(2.0f,0.0f);
+    cucomp64 div = cucomp64(2.0f,0.0f);
 
     ret = (exp(z)-exp(-z))/div;
 
     return ret;
 }
 
-__host__ __device__ cucomp64 tanh(cucomp64 z)
+__host__ __device__ cucomp64 tanh(const cucomp64 &z)
 {
     return sinh(z)/cosh(z);
 }
 
-__host__ __device__ cucomp64 pow(cucomp64 z1, cucomp64 z2)
+__host__ __device__ cucomp64 pow(const cucomp64 &z1, const cucomp64 &z2)
 {
 	cucomp64 ret;
-	if(z1.mag()>0.0)
-		ret = exp(z2*log(z1));
+	if(z1.mag()>0.0f)
+		ret = amscuda::cmp::exp(amscuda::cmp::log(z1)*z2);
 	else
-		ret = dtocomp64(0.0f,0.0f);
+		ret = cucomp64(0.0f,0.0f);
 	return ret;
 }
 
+// //returns "complex sign" of complex number - 0, or a unit number with same argument
+__host__ __device__ cucomp64 csgn(const cucomp64 &z)
+{
+    cucomp64 ret;
+    float mag = z.mag();
+    if(mag>0.0f)
+    {
+        ret = z;
+        ret = ret/mag;
+    }
+    return ret;
+}
+
+//////////////////////////
+//accumulation operators//
+//////////////////////////
+
+__host__ __device__ cucomp64& cucomp64::operator+=(const cucomp64& z)
+{
+
+    return (*this);
+}
+
+__host__ __device__ cucomp64& cucomp64::operator-=(const cucomp64& z)
+{
+    return (*this);
+}
+
+__host__ __device__ cucomp64& cucomp64::operator*=(const cucomp64& z)
+{
+    return (*this);
+}
+
+__host__ __device__ cucomp64& cucomp64::operator/=(const cucomp64& z)
+{
+    return (*this);
+}
+
+__host__ __device__ cucomp64& cucomp64::operator+=(const float& z)
+{
+    return (*this);
+}
+
+__host__ __device__ cucomp64& cucomp64::operator-=(const float& z)
+{
+    return (*this);
+}
+
+__host__ __device__ cucomp64& cucomp64::operator*=(const float& z)
+{
+    return (*this);
+}
+
+__host__ __device__ cucomp64& cucomp64::operator/=(const float& z)
+{
+    return (*this);
+}
 
 void test_cucomp64_1()
 {
@@ -448,11 +511,10 @@ void test_cucomp64_1()
 
     printf("sizeof double=%ld\n",(long)(8*sizeof(f1)));
     printf("sizeof double=%ld\n",(long)(8*sizeof(d1)));
-    printf("sizeof complex=%ld\n",(long)(8*sizeof(z1)));
-    printf("sizeof cucomp128=%ld\n",(long)(8*sizeof(a)));
+    printf("sizeof complex64=%ld\n",(long)(8*sizeof(z1)));
 
-    a = dtocomp64(1.0,1.0);
-    b = dtocomp64(1.0,-1.0);
+    a = cucomp64(1.0,1.0);
+    b = cucomp64(1.0,-1.0);
 
     printf("a=%1.4f + %1.4fi\n",a[0],a[1]);
     printf("b=%1.4f + %1.4fi\n",b[0],b[1]);
diff --git a/src/amsculib3/amscumath.cu b/src/amsculib3/math/amscumath.cu
similarity index 100%
rename from src/amsculib3/amscumath.cu
rename to src/amsculib3/math/amscumath.cu
diff --git a/src/amsculib3/math/cuvect2f.cu b/src/amsculib3/math/cuvec2f.cu
similarity index 68%
rename from src/amsculib3/math/cuvect2f.cu
rename to src/amsculib3/math/cuvec2f.cu
index b1185de..2f1df96 100644
--- a/src/amsculib3/math/cuvect2f.cu
+++ b/src/amsculib3/math/cuvec2f.cu
@@ -3,88 +3,88 @@
 namespace amscuda
 {
 
-    __host__ __device__ cuvect2f::cuvect2f()
+    __host__ __device__ cuvec2f::cuvec2f()
     {
         x = 0.0f; y = 0.0f;
         return;
     }
 
-    __host__ __device__ cuvect2f::~cuvect2f()
+    __host__ __device__ cuvec2f::~cuvec2f()
     {
         x = 0.0f; y = 0.0f;
         return;
     }
 
-    __host__ __device__ float& cuvect2f::operator[](const int &I)
+    __host__ __device__ float& cuvec2f::operator[](const int &I)
     {
         if(I==0) return x;
         if(I==1) return y;
         return x;
     }
 
-    __host__ __device__ const float& cuvect2f::operator[](const int &I) const
+    __host__ __device__ const float& cuvec2f::operator[](const int &I) const
     {
         if(I==0) return x;
         if(I==1) return y;
         return x;
     }
 
-    __host__ __device__ cuvect2f cuvect2f::operator+(const cuvect2f &rhs)
+    __host__ __device__ cuvec2f cuvec2f::operator+(const cuvec2f &rhs)
     {
-        cuvect2f ret;
+        cuvec2f ret;
         ret.x = x+rhs.x;
         ret.y = y+rhs.y;
 
         return ret;
     }
 
-    __host__ __device__ cuvect2f cuvect2f::operator-(const cuvect2f &rhs)
+    __host__ __device__ cuvec2f cuvec2f::operator-(const cuvec2f &rhs)
     {
-        cuvect2f ret;
+        cuvec2f ret;
         ret.x = x-rhs.x;
         ret.y = y-rhs.y;
 
         return ret;
     }
 
-    __host__ __device__ cuvect2f cuvect2f::operator*(const float &rhs)
+    __host__ __device__ cuvec2f cuvec2f::operator*(const float &rhs)
     {
-        cuvect2f ret;
+        cuvec2f ret;
         ret.x = x*rhs;
         ret.y = y*rhs;
         return ret;
     }
 
-    __host__ __device__ cuvect2f cuvect2f::operator/(const float &rhs)
+    __host__ __device__ cuvec2f cuvec2f::operator/(const float &rhs)
     {
-        cuvect2f ret;
+        cuvec2f ret;
         ret.x = x/rhs;
         ret.y = y/rhs;
         return ret;
     }
 
-    __host__ __device__ cuvect2f& cuvect2f::operator+=(const cuvect2f &rhs)
+    __host__ __device__ cuvec2f& cuvec2f::operator+=(const cuvec2f &rhs)
     {
         x = x + rhs.x;
         y = y + rhs.y;
         return *this;
     }
 
-    __host__ __device__ cuvect2f& cuvect2f::operator-=(const cuvect2f &rhs)
+    __host__ __device__ cuvec2f& cuvec2f::operator-=(const cuvec2f &rhs)
     {
         x = x - rhs.x;
         y = y - rhs.y;
         return *this;
     }
 
-    __host__ __device__ cuvect2f& cuvect2f::operator*=(const float &rhs)
+    __host__ __device__ cuvec2f& cuvec2f::operator*=(const float &rhs)
     {
         x = x * rhs;
         y = y * rhs;
         return *this;
     }
 
-    __host__ __device__ cuvect2f& cuvect2f::operator/=(const float &rhs)
+    __host__ __device__ cuvec2f& cuvec2f::operator/=(const float &rhs)
     {
         x = x / rhs;
         y = y / rhs;
@@ -92,34 +92,34 @@ namespace amscuda
     }
 
 
-    __host__ __device__ cuvect2f::cuvect2f(const float &_x, const float &_y)
+    __host__ __device__ cuvec2f::cuvec2f(const float &_x, const float &_y)
     {
         x = _x; y = _y;
         return;
     }
 
-    __host__ __device__ float cuvect2f_dot(const cuvect2f &a, const cuvect2f &b)
+    __host__ __device__ float cuvec2f_dot(const cuvec2f &a, const cuvec2f &b)
     {
         float ret = a.x*b.x+a.y*b.y;
         return ret;
     }
 
-    __host__ __device__ float cuvect2f_cross(const cuvect2f &a, const cuvect2f &b)
+    __host__ __device__ float cuvec2f_cross(const cuvec2f &a, const cuvec2f &b)
     {
         float ret = a.x*b.y-a.y*b.x;
         return ret;
     }
 
-    __host__ __device__ float cuvect2f_norm(const cuvect2f &a)
+    __host__ __device__ float cuvec2f_norm(const cuvec2f &a)
     {
         float ret = ::sqrtf(a.x*a.x+a.y*a.y);
         return ret;
     }
 
-    __host__ __device__ cuvect2f cuvect2f_normalize(const cuvect2f &a)
+    __host__ __device__ cuvec2f cuvec2f_normalize(const cuvec2f &a)
     {
-        cuvect2f ret;
-        float m = cuvect2f_norm(a);
+        cuvec2f ret;
+        float m = cuvec2f_norm(a);
         if(m>0.0)
         {
             ret.x = a.x/m; ret.y = a.y/m;
@@ -131,11 +131,11 @@ namespace amscuda
         return ret;
     }
 
-    __host__ __device__ cuvect2f cuvect2f_proj(const cuvect2f &a, const cuvect2f &b)
+    __host__ __device__ cuvec2f cuvec2f_proj(const cuvec2f &a, const cuvec2f &b)
     {
-        cuvect2f ret;
-        cuvect2f bn = cuvect2f_normalize(b);
-        float m = cuvect2f_dot(a,bn);
+        cuvec2f ret;
+        cuvec2f bn = cuvec2f_normalize(b);
+        float m = cuvec2f_dot(a,bn);
         ret = bn*m;
         return ret;
     }
@@ -264,9 +264,9 @@ namespace amscuda
         return ret;
     }
             
-    __host__ __device__ cuvect2f cumat2f::operator*(const cuvect2f &rhs)
+    __host__ __device__ cuvec2f cumat2f::operator*(const cuvec2f &rhs)
     {
-        cuvect2f ret;
+        cuvec2f ret;
         
         ret.x = m00*rhs.x + m01*rhs.y;
         ret.y = m10*rhs.x + m11*rhs.y;
@@ -440,94 +440,16 @@ namespace amscuda
         return R;
     }
 
-///////////////////////////
-//legacy array operations//
-///////////////////////////
-
-//copies src to dest
-__host__ __device__ void mat2f_copy(float *mat2f_dest, const float *mat2f_src)
-{
-    mat2f_dest[0+0*2] = mat2f_src[0+0*2];
-    mat2f_dest[1+0*2] = mat2f_src[1+0*2];
-    mat2f_dest[0+1*2] = mat2f_src[0+1*2];
-    mat2f_dest[1+1*2] = mat2f_src[1+1*2];
-    
-    return;
-}
-
-//inverts mat?inout[4]
-__host__ __device__ void mat2f_inverse(float *mat2inout)
-{
-    float det = mat2inout[0+0*2]*mat2inout[1+1*2]-mat2inout[0+1*2]*mat2inout[1+0*2];
-    float mat2in[4];
-
-    mat2f_copy(mat2in,mat2inout);
-    mat2inout[0+0*2] = mat2inout[1+1*2]/det;
-    mat2inout[1+0*2] = -mat2inout[1+0*2]/det;
-    mat2inout[0+1*2] = -mat2inout[0+1*2]/det;
-    mat2inout[1+1*2] = mat2inout[0+0*2]/det;
-    
-    return;
-}
-
-//rotatin matrix from angle
-__host__ __device__ void mat2f_rot_from_angle(float angle, float *mat2)
-{
-    mat2[0+0*2] = ::cosf(angle);
-    mat2[1+0*2] = ::sinf(angle);
-    mat2[0+1*2] = -::sinf(angle);
-    mat2[1+1*2] = ::cosf(angle);
-    return;
-}
-
-//multiplies c = a*b
-__host__ __device__ void mat2f_mult(float *mat2a, float *mat2b, float *mat2c)
-{
-    float mat2a_in[4];
-    float mat2b_in[4];
-
-    if(mat2a==NULL || mat2b==NULL || mat2c==NULL)
+    __host__ __device__ cuvec2f operator-(const cuvec2f &rhs)
     {
-        return;
+        cuvec2f ret;
+        ret[0] = -rhs[0];
+        ret[1] = -rhs[1];
+        return ret;
     }
 
-    mat2f_copy(mat2a_in,mat2a);
-    mat2f_copy(mat2b_in,mat2b);
-    
-    mat2c[0+0*2] = mat2a_in[0+0*2]*mat2b_in[0+0*2] + mat2a_in[1+0*2]*mat2b_in[0+1*2];
-    mat2c[1+0*2] = mat2a_in[0+0*2]*mat2b_in[1+0*2] + mat2a_in[1+0*2]*mat2b_in[1+1*2];
-    mat2c[0+1*2] = mat2a_in[0+1*2]*mat2b_in[0+0*2] + mat2a_in[1+1*2]*mat2b_in[0+1*2];
-    mat2c[1+1*2] = mat2a_in[0+1*2]*mat2b_in[1+0*2] + mat2a_in[1+1*2]*mat2b_in[1+1*2];
 
-    return;
-}
-
-// ret = a*b
-__host__ __device__ cuvect2f mat2f_mult(float *mat2a, const cuvect2f &b)
-{
-    cuvect2f ret;
-    ret.x = b.x*mat2a[0+0*2] + b.y*mat2a[1+0*2];
-    ret.y = b.x*mat2a[0+1*2] + b.y*mat2a[1+1*2];
-    return ret;
-}
-
-__host__ __device__ cuvect2f operator-(cuvect2f rhs)
-{
-    cuvect2f ret;
-    ret[0] = -rhs[0];
-    ret[1] = -rhs[1];
-    return ret;
-}
-
-__host__ __device__ cumat2f operator-(cumat2f rhs)
-{
-    cumat2f ret;
-    int I;
-    for(I=0;I<4;I++) ret[I] = -rhs[I];
-    return ret;
-}
-
-void test_cuvect2f_1()
+void test_cuvec2f_1()
 {
 
 
diff --git a/src/amsculib3/math/cuvect3f.cu b/src/amsculib3/math/cuvec3f.cu
similarity index 85%
rename from src/amsculib3/math/cuvect3f.cu
rename to src/amsculib3/math/cuvec3f.cu
index b296dbd..fb69594 100644
--- a/src/amsculib3/math/cuvect3f.cu
+++ b/src/amsculib3/math/cuvec3f.cu
@@ -3,19 +3,19 @@
 namespace amscuda
 {
 
-    __host__ __device__ cuvect3f::cuvect3f()
+    __host__ __device__ cuvec3f::cuvec3f()
     {
         x = 0.0f; y = 0.0f; z = 0.0f;
         return;
     }
 
-    __host__ __device__ cuvect3f::~cuvect3f()
+    __host__ __device__ cuvec3f::~cuvec3f()
     {
         x = 0.0f; y = 0.0f; z = 0.0f;
         return;
     }
 
-    __host__ __device__ float& cuvect3f::operator[](const int &I)
+    __host__ __device__ float& cuvec3f::operator[](const int &I)
     {
         if(I==0) return x;
         if(I==1) return y;
@@ -23,7 +23,7 @@ namespace amscuda
         return x;
     }
 
-    __host__ __device__ const float& cuvect3f::operator[](const int &I) const
+    __host__ __device__ const float& cuvec3f::operator[](const int &I) const
     {
         if(I==0) return x;
         if(I==1) return y;
@@ -31,9 +31,9 @@ namespace amscuda
         return x;
     }
 
-    __host__ __device__ cuvect3f cuvect3f::operator+(const cuvect3f &rhs)
+    __host__ __device__ cuvec3f cuvec3f::operator+(const cuvec3f &rhs)
     {
-        cuvect3f ret;
+        cuvec3f ret;
         ret.x = x+rhs.x;
         ret.y = y+rhs.y;
         ret.z = z+rhs.z;
@@ -41,9 +41,9 @@ namespace amscuda
         return ret;
     }
 
-    __host__ __device__ cuvect3f cuvect3f::operator-(const cuvect3f &rhs)
+    __host__ __device__ cuvec3f cuvec3f::operator-(const cuvec3f &rhs)
     {
-        cuvect3f ret;
+        cuvec3f ret;
         ret.x = x-rhs.x;
         ret.y = y-rhs.y;
         ret.z = z-rhs.z;
@@ -51,25 +51,25 @@ namespace amscuda
         return ret;
     }
 
-    __host__ __device__ cuvect3f cuvect3f::operator*(const float &rhs)
+    __host__ __device__ cuvec3f cuvec3f::operator*(const float &rhs)
     {
-        cuvect3f ret;
+        cuvec3f ret;
         ret.x = x*rhs;
         ret.y = y*rhs;
         ret.z = z*rhs;
         return ret;
     }
 
-    __host__ __device__ cuvect3f cuvect3f::operator/(const float &rhs)
+    __host__ __device__ cuvec3f cuvec3f::operator/(const float &rhs)
     {
-        cuvect3f ret;
+        cuvec3f ret;
         ret.x = x/rhs;
         ret.y = y/rhs;
         ret.z = z/rhs;
         return ret;
     }
 
-    __host__ __device__ cuvect3f& cuvect3f::operator+=(const cuvect3f &rhs)
+    __host__ __device__ cuvec3f& cuvec3f::operator+=(const cuvec3f &rhs)
     {
         x = x + rhs.x;
         y = y + rhs.y;
@@ -77,7 +77,7 @@ namespace amscuda
         return *this;
     }
 
-    __host__ __device__ cuvect3f& cuvect3f::operator-=(const cuvect3f &rhs)
+    __host__ __device__ cuvec3f& cuvec3f::operator-=(const cuvec3f &rhs)
     {
         x = x - rhs.x;
         y = y - rhs.y;
@@ -85,7 +85,7 @@ namespace amscuda
         return *this;
     }
 
-    __host__ __device__ cuvect3f& cuvect3f::operator*=(const float &rhs)
+    __host__ __device__ cuvec3f& cuvec3f::operator*=(const float &rhs)
     {
         x = x * rhs;
         y = y * rhs;
@@ -93,7 +93,7 @@ namespace amscuda
         return *this;
     }
 
-    __host__ __device__ cuvect3f& cuvect3f::operator/=(const float &rhs)
+    __host__ __device__ cuvec3f& cuvec3f::operator/=(const float &rhs)
     {
         x = x / rhs;
         y = y / rhs;
@@ -102,23 +102,23 @@ namespace amscuda
     }
 
 
-    __host__ __device__ cuvect3f::cuvect3f(const float &_x, const float &_y, const float &_z)
+    __host__ __device__ cuvec3f::cuvec3f(const float &_x, const float &_y, const float &_z)
     {
         x = _x; y = _y; z = _z;
         return;
     }
 
 
-    __host__ __device__ float cuvect3f_dot(const cuvect3f &a, const cuvect3f &b)
+    __host__ __device__ float cuvec3f_dot(const cuvec3f &a, const cuvec3f &b)
     {
         float ret = a.x*b.x+a.y*b.y+a.z*b.z;
 
         return ret;
     }
 
-    __host__ __device__ cuvect3f cuvect3f_cross(const cuvect3f &a, const cuvect3f &b)
+    __host__ __device__ cuvec3f cuvec3f_cross(const cuvec3f &a, const cuvec3f &b)
     {
-        cuvect3f ret;
+        cuvec3f ret;
         ret[0] = a[1]*b[2]-a[2]*b[1];
         ret[1] = a[2]*b[0]-a[0]*b[2];
         ret[2] = a[0]*b[1]-a[1]*b[0];
@@ -126,16 +126,16 @@ namespace amscuda
         return ret;
     }
 
-    __host__ __device__ float cuvect3f_norm(const cuvect3f &a)
+    __host__ __device__ float cuvec3f_norm(const cuvec3f &a)
     {
         float ret;
         ret = ::sqrtf(a.x*a.x+a.y*a.y+a.z*a.z);
         return ret;
     }
 
-    __host__ __device__ cuvect3f cuvect3f_normalize(const cuvect3f &a)
+    __host__ __device__ cuvec3f cuvec3f_normalize(const cuvec3f &a)
     {
-        cuvect3f ret;
+        cuvec3f ret;
         float m;
         m = ::sqrtf(a.x*a.x+a.y*a.y+a.z*a.z);
         if(m>0.0)
@@ -150,11 +150,11 @@ namespace amscuda
         return ret;
     }
 
-    __host__ __device__ cuvect3f cuvect3f_proj(const cuvect3f &a, const cuvect3f &b)
+    __host__ __device__ cuvec3f cuvec3f_proj(const cuvec3f &a, const cuvec3f &b)
     {
-        cuvect3f ret;
-        cuvect3f bn = cuvect3f_normalize(b);
-        float m = cuvect3f_dot(a,bn);
+        cuvec3f ret;
+        cuvec3f bn = cuvec3f_normalize(b);
+        float m = cuvec3f_dot(a,bn);
         ret = bn*m;
         return ret;
     }
@@ -361,9 +361,9 @@ __host__ __device__ cumat3f cumat3f::operator/(const float &rhs)
     return ret;
 }
         
-__host__ __device__ cuvect3f cumat3f::operator*(const cuvect3f &rhs)
+__host__ __device__ cuvec3f cumat3f::operator*(const cuvec3f &rhs)
 {
-    cuvect3f ret;
+    cuvec3f ret;
     
     ret.x = m00*rhs.x + m01*rhs.y + m02*rhs.z;
     ret.y = m10*rhs.x + m11*rhs.y + m12*rhs.z;
@@ -597,7 +597,7 @@ __host__ void test_cudavectf_logic1()
 
     // printf("3 dim vector and matrix functional tests on host side\n");
     
-    // cuvect3f a,b,c;
+    // cuvec3f a,b,c;
     // float ma[9],mb[9],mc[9];
 
     // int I,J;
@@ -636,7 +636,7 @@ __host__ void test_cudavectf_logic1()
     //     }
     // }
 
-    // a = cuvect3f(1,1,1);
+    // a = cuvec3f(1,1,1);
     // b = mat3f_mult(ma,a);
     // b = mat3f_mult(mb,b);
     
@@ -645,8 +645,8 @@ __host__ void test_cudavectf_logic1()
     //     printf("a[%d] = %1.3f, b[%d] = %1.3f\n",I,a[I],I,b[I]);
     // }
 
-    // a = cuvect3f(1,0,1);
-    // b = cuvect3f(0,1,-1);
+    // a = cuvec3f(1,0,1);
+    // b = cuvec3f(0,1,-1);
     // c = a+b;
 
     // for(I=0;I<3;I++)
@@ -661,31 +661,31 @@ __host__ void test_cudavectf_logic1()
     //     printf("a[%d] = %1.3f, b[%d] = %1.3f, c[%d] = %1.3f\n",I,a[I],I,b[I],I,c[I]);
     // }
 
-    // c = cuvect3f_cross(a,b);
+    // c = cuvec3f_cross(a,b);
 
     // for(I=0;I<3;I++)
     // {
     //     printf("a[%d] = %1.3f, b[%d] = %1.3f, c[%d] = %1.3f\n",I,a[I],I,b[I],I,c[I]);
     // }
 
-    // printf("c dot a = %1.3f, c dot b = %1.3f\n",cuvect3f_dot(c,a),cuvect3f_dot(c,b));
+    // printf("c dot a = %1.3f, c dot b = %1.3f\n",cuvec3f_dot(c,a),cuvec3f_dot(c,b));
 
-    // printf("norm(a)=%1.3f, norm(b)=%1.3f, norm(c)=%1.3f\n",cuvect3f_norm(a),cuvect3f_norm(b),cuvect3f_norm(c));
-    // a = cuvect3f_normalize(a);
-    // b = cuvect3f_normalize(b);
-    // c = cuvect3f_normalize(c);
+    // printf("norm(a)=%1.3f, norm(b)=%1.3f, norm(c)=%1.3f\n",cuvec3f_norm(a),cuvec3f_norm(b),cuvec3f_norm(c));
+    // a = cuvec3f_normalize(a);
+    // b = cuvec3f_normalize(b);
+    // c = cuvec3f_normalize(c);
 
     // for(I=0;I<3;I++)
     // {
     //     printf("a[%d] = %1.3f, b[%d] = %1.3f, c[%d] = %1.3f\n",I,a[I],I,b[I],I,c[I]);
     // }
-    // printf("c dot a = %1.3f, c dot b = %1.3f\n",cuvect3f_dot(c,a),cuvect3f_dot(c,b));
-    // printf("norm(a)=%1.3f, norm(b)=%1.3f, norm(c)=%1.3f\n",cuvect3f_norm(a),cuvect3f_norm(b),cuvect3f_norm(c));
+    // printf("c dot a = %1.3f, c dot b = %1.3f\n",cuvec3f_dot(c,a),cuvec3f_dot(c,b));
+    // printf("norm(a)=%1.3f, norm(b)=%1.3f, norm(c)=%1.3f\n",cuvec3f_norm(a),cuvec3f_norm(b),cuvec3f_norm(c));
 
     // return;
 }
 
-__host__ __device__ cumat3f hodge_dual(const cuvect3f &vin)
+__host__ __device__ cumat3f hodge_dual(const cuvec3f &vin)
 {
     cumat3f ret;
 
@@ -704,9 +704,9 @@ __host__ __device__ cumat3f hodge_dual(const cuvect3f &vin)
     return ret;
 }
 
-__host__ __device__ cuvect3f hodge_dual(const cumat3f &min)
+__host__ __device__ cuvec3f hodge_dual(const cumat3f &min)
 {
-    cuvect3f ret;
+    cuvec3f ret;
 
     ret.x = 0.5f*(min.m12 - min.m21);
     ret.y = 0.5f*(min.m20 - min.m02);
@@ -731,13 +731,13 @@ __host__ __device__ const cumat3f cumat3f_zeros()
     return ret;
 }
 
-__host__ __device__ cumat3f rotmat_from_axisangle(const cuvect3f &axis, const float &angle)
+__host__ __device__ cumat3f rotmat_from_axisangle(const cuvec3f &axis, const float &angle)
 {
     cumat3f ret = cumat3f_zeros();
     cumat3f H;
-    cuvect3f _naxis;
+    cuvec3f _naxis;
 
-    _naxis = cuvect3f_normalize(axis);
+    _naxis = cuvec3f_normalize(axis);
     H = hodge_dual(_naxis);
 
     ret += H*sinf(angle);
diff --git a/src/amsculib3/math/cuvect4f.cu b/src/amsculib3/math/cuvec4f.cu
similarity index 90%
rename from src/amsculib3/math/cuvect4f.cu
rename to src/amsculib3/math/cuvec4f.cu
index 666d539..588a03e 100644
--- a/src/amsculib3/math/cuvect4f.cu
+++ b/src/amsculib3/math/cuvec4f.cu
@@ -4,22 +4,22 @@ namespace amscuda
 {
 
 ////////////
-//cuvect4ff//
+//cuvec4ff//
 ////////////
 
-__host__ __device__ cuvect4f::cuvect4f()
+__host__ __device__ cuvec4f::cuvec4f()
 {
     x = 0.0f; y = 0.0f; z = 0.0f; w = 0.0f;
     return;
 }
 
-__host__ __device__ cuvect4f::~cuvect4f()
+__host__ __device__ cuvec4f::~cuvec4f()
 {
     x = 0.0f; y = 0.0f; z = 0.0f; w = 0.0f;
     return;
 }
 
-__host__ __device__ float& cuvect4f::operator[](const int &I)
+__host__ __device__ float& cuvec4f::operator[](const int &I)
 {
     if(I==0) return x;
     if(I==1) return y;
@@ -28,7 +28,7 @@ __host__ __device__ float& cuvect4f::operator[](const int &I)
     return x;
 }
 
-__host__ __device__ const float& cuvect4f::operator[](const int &I) const
+__host__ __device__ const float& cuvec4f::operator[](const int &I) const
 {
     if(I==0) return x;
     if(I==1) return y;
@@ -37,9 +37,9 @@ __host__ __device__ const float& cuvect4f::operator[](const int &I) const
     return x;
 }
 
-__host__ __device__ cuvect4f cuvect4f::operator+(const cuvect4f &rhs)
+__host__ __device__ cuvec4f cuvec4f::operator+(const cuvec4f &rhs)
 {
-    cuvect4f ret;
+    cuvec4f ret;
     ret.x = x+rhs.x;
     ret.y = y+rhs.y;
     ret.z = z+rhs.z;
@@ -48,9 +48,9 @@ __host__ __device__ cuvect4f cuvect4f::operator+(const cuvect4f &rhs)
     return ret;
 }
 
-__host__ __device__ cuvect4f cuvect4f::operator-(const cuvect4f &rhs)
+__host__ __device__ cuvec4f cuvec4f::operator-(const cuvec4f &rhs)
 {
-    cuvect4f ret;
+    cuvec4f ret;
     ret.x = x-rhs.x;
     ret.y = y-rhs.y;
     ret.z = z-rhs.z;
@@ -59,9 +59,9 @@ __host__ __device__ cuvect4f cuvect4f::operator-(const cuvect4f &rhs)
     return ret;
 }
 
-__host__ __device__ cuvect4f cuvect4f::operator*(const float &rhs)
+__host__ __device__ cuvec4f cuvec4f::operator*(const float &rhs)
 {
-    cuvect4f ret;
+    cuvec4f ret;
     ret.x = x*rhs;
     ret.y = y*rhs;
     ret.z = z*rhs;
@@ -69,9 +69,9 @@ __host__ __device__ cuvect4f cuvect4f::operator*(const float &rhs)
     return ret;
 }
 
-__host__ __device__ cuvect4f cuvect4f::operator/(const float &rhs)
+__host__ __device__ cuvec4f cuvec4f::operator/(const float &rhs)
 {
-    cuvect4f ret;
+    cuvec4f ret;
     ret.x = x/rhs;
     ret.y = y/rhs;
     ret.z = z/rhs;
@@ -79,7 +79,7 @@ __host__ __device__ cuvect4f cuvect4f::operator/(const float &rhs)
     return ret;
 }
 
-__host__ __device__ cuvect4f& cuvect4f::operator+=(const cuvect4f &rhs)
+__host__ __device__ cuvec4f& cuvec4f::operator+=(const cuvec4f &rhs)
 {
     x = x + rhs.x;
     y = y + rhs.y;
@@ -88,7 +88,7 @@ __host__ __device__ cuvect4f& cuvect4f::operator+=(const cuvect4f &rhs)
     return *this;
 }
 
-__host__ __device__ cuvect4f& cuvect4f::operator-=(const cuvect4f &rhs)
+__host__ __device__ cuvec4f& cuvec4f::operator-=(const cuvec4f &rhs)
 {
     x = x - rhs.x;
     y = y - rhs.y;
@@ -97,7 +97,7 @@ __host__ __device__ cuvect4f& cuvect4f::operator-=(const cuvect4f &rhs)
     return *this;
 }
 
-__host__ __device__ cuvect4f& cuvect4f::operator*=(const float &rhs)
+__host__ __device__ cuvec4f& cuvec4f::operator*=(const float &rhs)
 {
     x = x * rhs;
     y = y * rhs;
@@ -106,7 +106,7 @@ __host__ __device__ cuvect4f& cuvect4f::operator*=(const float &rhs)
     return *this;
 }
 
-__host__ __device__ cuvect4f& cuvect4f::operator/=(const float &rhs)
+__host__ __device__ cuvec4f& cuvec4f::operator/=(const float &rhs)
 {
     x = x / rhs;
     y = y / rhs;
@@ -115,9 +115,9 @@ __host__ __device__ cuvect4f& cuvect4f::operator/=(const float &rhs)
     return *this;
 }
 
-__host__ __device__ cuvect4f operator-(const cuvect4f &rhs)
+__host__ __device__ cuvec4f operator-(const cuvec4f &rhs)
 {
-    cuvect4f ret;
+    cuvec4f ret;
     ret[0] = -rhs[0];
     ret[1] = -rhs[1];
     ret[2] = -rhs[2];
@@ -126,7 +126,7 @@ __host__ __device__ cuvect4f operator-(const cuvect4f &rhs)
 }
 
 
-__host__ __device__ cuvect4f::cuvect4f(const float &_x, const float &_y, const float &_z, const float &_w)
+__host__ __device__ cuvec4f::cuvec4f(const float &_x, const float &_y, const float &_z, const float &_w)
 {
     x = _x; y = _y; z = _z; w = _w;
     return;
@@ -404,9 +404,9 @@ __host__ __device__ cumat4f cumat4f::operator/(const float &rhs)
     return ret;
 }
         
-__host__ __device__ cuvect4f cumat4f::operator*(const cuvect4f &rhs)
+__host__ __device__ cuvec4f cumat4f::operator*(const cuvec4f &rhs)
 {
-    cuvect4f ret;
+    cuvec4f ret;
     ret.x = m00*rhs.x + m01*rhs.y + m02*rhs.z + m03*rhs.w;
     ret.y = m10*rhs.x + m11*rhs.y + m12*rhs.z + m13*rhs.w;
     ret.z = m20*rhs.x + m21*rhs.y + m22*rhs.z + m23*rhs.w;
@@ -674,34 +674,34 @@ __host__ __device__ cumat4f& cumat4f::operator*=(const cumat4f &rhs)
 }
 
 
-__host__ __device__ float cuvect4f_dot(cuvect4f &a, cuvect4f &b)
+__host__ __device__ float cuvec4f_dot(cuvec4f &a, cuvec4f &b)
 {
     float ret = 0.0f;
     ret = a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w;
     return ret;
 }
 
-__host__ __device__ float cuvect4f_norm(cuvect4f &a)
+__host__ __device__ float cuvec4f_norm(cuvec4f &a)
 {
     float ret = 0.0f;
-    ret = ::sqrtf(cuvect4f_dot(a,a));
+    ret = ::sqrtf(cuvec4f_dot(a,a));
     return ret;
 }
 
-__host__ __device__ cuvect4f cuvect4f_normalize(cuvect4f &a)
+__host__ __device__ cuvec4f cuvec4f_normalize(cuvec4f &a)
 {
-    cuvect4f ret = cuvect4f(0.0f,0.0f,0.0f,0.0f);
-    float nrm = cuvect4f_norm(a);
+    cuvec4f ret = cuvec4f(0.0f,0.0f,0.0f,0.0f);
+    float nrm = cuvec4f_norm(a);
     if(nrm>0.0f)
         ret = a/nrm;
     return ret;
 }
 
-__host__ __device__ cuvect4f cuvect4f_proj(cuvect4f &a, cuvect4f &b)
+__host__ __device__ cuvec4f cuvec4f_proj(cuvec4f &a, cuvec4f &b)
 {
-    cuvect4f ret;
-    cuvect4f bn = cuvect4f_normalize(b);
-    float d = cuvect4f_dot(a,bn);
+    cuvec4f ret;
+    cuvec4f bn = cuvec4f_normalize(b);
+    float d = cuvec4f_dot(a,bn);
     ret = bn*d;
     return ret;
 }
diff --git a/src/amsculib3/math/cuvect2.cu b/src/amsculib3/math/cuvect2.cu
deleted file mode 100644
index e9b5169..0000000
--- a/src/amsculib3/math/cuvect2.cu
+++ /dev/null
@@ -1,377 +0,0 @@
-#include <amsculib3/amsculib3.hpp>
-
-namespace amscuda
-{
-
-    __host__ __device__ cuvect2::cuvect2()
-    {
-        x = 0.0; y = 0.0;
-        return;
-    }
-
-    __host__ __device__ cuvect2::~cuvect2()
-    {
-        x = 0.0; y = 0.0;
-        return;
-    }
-
-    __host__ __device__ cuvect2::cuvect2(double _x, double _y)
-    {
-        x = _x; y = _y;
-        return;
-    }
-
-    __host__ __device__ double& cuvect2::operator[](const int I)
-    {
-        if(I==0) return x;
-        if(I==1) return y;
-        return x;
-    }
-
-    __host__ __device__ const double& cuvect2::operator[](const int I) const
-    {
-        if(I==0) return x;
-        if(I==1) return y;
-        return x;
-    }
-
-    __host__ __device__ cuvect2 cuvect2::operator+(cuvect2 lhs)
-    {
-        cuvect2 ret;
-        ret.x = x + lhs.x;
-        ret.y = y + lhs.y;
-        return ret;
-    }
-
-    __host__ __device__ cuvect2 cuvect2::operator-(cuvect2 lhs)
-    {
-        cuvect2 ret;
-        ret.x = x - lhs.x;
-        ret.y = y - lhs.y;
-        return ret;
-    }
-    __host__ __device__ cuvect2 cuvect2::operator*(double lhs)
-    {
-        cuvect2 ret;
-        ret.x = x*lhs;
-        ret.y = y*lhs;
-        return ret;
-    }
-
-    __host__ __device__ cuvect2 cuvect2::operator/(double lhs)
-    {
-        cuvect2 ret;
-        ret.x = x/lhs;
-        ret.y = y/lhs;
-        return ret;
-    }
-
-    __host__ __device__ double cuvect2_dot(cuvect2 a, cuvect2 b)
-    {
-        double ret = a.x*b.x+a.y*b.y;
-        return ret;
-    }
-
-    __host__ __device__ double cuvect2_cross(cuvect2 a, cuvect2 b)
-    {
-        double ret = a.x*b.y-a.y*b.x;
-        return ret;
-    }
-
-    __host__ __device__ double cuvect2_norm(cuvect2 a)
-    {
-        double ret = ::sqrt(a.x*a.x+a.y*a.y);
-        return ret;
-    }
-
-    __host__ __device__ cuvect2 cuvect2_normalize(cuvect2 a)
-    {
-        cuvect2 ret;
-        double m = cuvect2_norm(a);
-        if(m>0.0)
-        {
-            ret.x = a.x/m; ret.y = a.y/m;
-        }
-        else
-        {
-            ret.x = 0.0; ret.y = 0.0;
-        }
-        return ret;
-    }
-
-    __host__ __device__ cuvect2 cuvect2_proj(cuvect2 a, cuvect2 b)
-    {
-        cuvect2 ret;
-        cuvect2 bn = cuvect2_normalize(b);
-        double m = cuvect2_dot(a,bn);
-        ret = bn*m;
-        return ret;
-    }
-
-
-__host__ __device__ cumat2::cumat2()
-{
-    int I;
-    for(I=0;I<4;I++)
-    {
-        dat[I] = 0.0;
-    }
-    return;
-}
-
-__host__ __device__ cumat2::~cumat2()
-{
-    int I;
-    for(I=0;I<4;I++)
-    {
-        dat[I] = 0.0;
-    }
-    return;
-}
-
-__host__ __device__ double& cumat2::operator[](const int I)
-{
-    return dat[I];
-}
-
-__host__ __device__ double& cumat2::operator()(const int I, const int J)
-{
-    return dat[I+2*J];
-}
-
-__host__ __device__ cumat2 cumat2::operator+(cumat2 lhs)
-{
-    int I;
-    cumat2 ret;
-    for(I=0;I<4;I++)
-    {
-        ret.dat[I] = this->dat[I] + lhs.dat[I];
-    }
-    return ret;
-}
-
-__host__ __device__ cumat2 cumat2::operator-(cumat2 lhs)
-{
-    int I;
-    cumat2 ret;
-    for(I=0;I<4;I++)
-    {
-        ret.dat[I] = this->dat[I] - lhs.dat[I];
-    }
-    return ret;
-}
-
-__host__ __device__ cumat2 cumat2::operator*(double lhs)
-{
-    cumat2 ret;
-    int I;
-    for(I=0;I<4;I++)
-    {
-        ret.dat[I] = this->dat[I]*lhs;
-    }
-    return ret;
-}
-
-__host__ __device__ cumat2 cumat2::operator/(double lhs)
-{
-    cumat2 ret;
-    int I;
-    for(I=0;I<4;I++)
-    {
-        ret.dat[I] = this->dat[I]/lhs;
-    }
-    return ret;
-}
-
-__host__ __device__ double& cumat2::at(const int I, const int J)
-{
-    return dat[I+2*J];
-}
-        
-
-__host__ __device__ cuvect2 cumat2::operator*(cuvect2 lhs)
-{
-    cuvect2 ret = cuvect2(0.0,0.0);
-    int I,J;
-    for(I=0;I<2;I++)
-    {
-        for(J=0;J<2;J++)
-        {
-            ret[I] = ret[I] + this->at(I,J)*lhs[J];
-        }
-    }
-    return ret;
-}
-
-__host__ __device__ cumat2 cumat2::operator*(cumat2 lhs)
-{
-    cumat2 ret;
-    int I,J,K;
-
-    for(I=0;I<2;I++)
-    {
-        for(J=0;J<2;J++)
-        {
-            ret(I,J) = 0.0;
-            for(K=0;K<2;K++)
-            {
-                ret(I,J) = ret(I,J) + this->at(I,K)*lhs(K,J);
-            }
-        }
-    }
-
-    return ret;
-}
-
-__host__ __device__ double cumat2::det()
-{
-    double ret = 0.0;
-    ret = ret + this->at(0,0)*this->at(1,1);
-    ret = ret - this->at(1,0)*this->at(0,1);
-    return ret;
-}
-
-__host__ __device__ cumat2 cumat2::transpose()
-{
-    cumat2 q;
-    int I,J;
-    for(I=0;I<2;I++)
-    {
-        for(J=0;J<2;J++)
-        {
-            q.at(I,J) = this->at(J,I);
-        }
-    }
-    return q;
-}
-
-__host__ __device__ cumat2 cumat2::inverse()
-{
-    cumat2 q;
-    double dt = q.det();
-    if(dt!=0.0)
-    {
-        q(0,0) = this->at(1,1)/dt;
-        q(0,1) = -this->at(0,1)/dt;
-        q(1,0) = -this->at(1,0)/dt;
-        q(1,1) = this->at(0,0)/dt;
-    }
-    else
-    {
-        q(0,0) = inf;
-        q(0,1) = inf;
-        q(1,0) = inf;
-        q(1,1) = inf;
-    }
-
-    return q;
-}
-
-//2x2 matrix operations
-//matrix order is assumed to be mat[I,J] = mat[I+3*J]
-
-//transpose a 2x2 matrix in place
-__host__ __device__ void mat2_transpose(double *mat2inout)
-{
-    double mat2_in[4];
-
-    mat2_copy(mat2_in,mat2inout);
-    mat2inout[0+0*2] = mat2_in[0+0*2];
-    mat2inout[1+0*2] = mat2_in[0+1*2];
-    mat2inout[0+1*2] = mat2_in[1+0*2];
-    mat2inout[1+1*2] = mat2_in[1+1*2];
-    
-    return;
-}
-
-//copies src to dest
-__host__ __device__ void mat2_copy(double *mat2_dest, const double *mat2_src)
-{
-    mat2_dest[0+0*2] = mat2_src[0+0*2];
-    mat2_dest[1+0*2] = mat2_src[1+0*2];
-    mat2_dest[0+1*2] = mat2_src[0+1*2];
-    mat2_dest[1+1*2] = mat2_src[1+1*2];
-    
-    return;
-}
-
-//inverts mat?inout[4]
-__host__ __device__ void mat2_inverse(double *mat2inout)
-{
-    double det = mat2inout[0+0*2]*mat2inout[1+1*2]-mat2inout[0+1*2]*mat2inout[1+0*2];
-    double mat2in[4];
-
-    mat2_copy(mat2in,mat2inout);
-    mat2inout[0+0*2] = mat2inout[1+1*2]/det;
-    mat2inout[1+0*2] = -mat2inout[1+0*2]/det;
-    mat2inout[0+1*2] = -mat2inout[0+1*2]/det;
-    mat2inout[1+1*2] = mat2inout[0+0*2]/det;
-    
-    return;
-}
-
-//rotatin matrix from angle
-__host__ __device__ void mat2_rot_from_angle(double angle, double *mat2)
-{
-    mat2[0+0*2] = ::cos(angle);
-    mat2[1+0*2] = ::sin(angle);
-    mat2[0+1*2] = -::sin(angle);
-    mat2[1+1*2] = ::cos(angle);
-    return;
-}
-
-//multiplies c = a*b
-__host__ __device__ void mat2_mult(double *mat2a, double *mat2b, double *mat2c)
-{
-    double mat2a_in[4];
-    double mat2b_in[4];
-
-    if(mat2a==NULL || mat2b==NULL || mat2c==NULL)
-    {
-        return;
-    }
-
-    mat2_copy(mat2a_in,mat2a);
-    mat2_copy(mat2b_in,mat2b);
-    
-    mat2c[0+0*2] = mat2a_in[0+0*2]*mat2b_in[0+0*2] + mat2a_in[1+0*2]*mat2b_in[0+1*2];
-    mat2c[1+0*2] = mat2a_in[0+0*2]*mat2b_in[1+0*2] + mat2a_in[1+0*2]*mat2b_in[1+1*2];
-    mat2c[0+1*2] = mat2a_in[0+1*2]*mat2b_in[0+0*2] + mat2a_in[1+1*2]*mat2b_in[0+1*2];
-    mat2c[1+1*2] = mat2a_in[0+1*2]*mat2b_in[1+0*2] + mat2a_in[1+1*2]*mat2b_in[1+1*2];
-
-    return;
-}
-
-// ret = a*b
-__host__ __device__ cuvect2 mat2_mult(double *mat2a, cuvect2 b)
-{
-    cuvect2 ret;
-    ret.x = b.x*mat2a[0+0*2] + b.y*mat2a[1+0*2];
-    ret.y = b.x*mat2a[0+1*2] + b.y*mat2a[1+1*2];
-    return ret;
-}
-
-__host__ __device__ cuvect2 operator-(cuvect2 rhs)
-{
-    cuvect2 ret;
-    ret[0] = -rhs[0];
-    ret[1] = -rhs[1];
-    return ret;
-}
-
-__host__ __device__ cumat2 operator-(cumat2 rhs)
-{
-    cumat2 ret;
-    int I;
-    for(I=0;I<4;I++) ret[I] = -rhs[I];
-    return ret;
-}
-
-void test_cuvect2_1()
-{
-
-
-    return;
-}
-
-};
\ No newline at end of file
diff --git a/src/amsculib3/math/cuvect3.cu b/src/amsculib3/math/cuvect3.cu
deleted file mode 100644
index 34159c5..0000000
--- a/src/amsculib3/math/cuvect3.cu
+++ /dev/null
@@ -1,860 +0,0 @@
-#include <amsculib3/amsculib3.hpp>
-
-namespace amscuda
-{
-
-    __host__ __device__ cuvect3::cuvect3()
-    {
-        x = 0.0; y = 0.0; z = 0.0;
-        return;
-    }
-
-    __host__ __device__ cuvect3::~cuvect3()
-    {
-        x = 0.0; y = 0.0; z = 0.0;
-        return;
-    }
-
-    __host__ __device__ double& cuvect3::operator[](const int &I)
-    {
-        if(I==0) return x;
-        if(I==1) return y;
-        if(I==2) return z;
-        return x;
-    }
-
-    __host__ __device__ const double& cuvect3::operator[](const int &I) const
-    {
-        if(I==0) return x;
-        if(I==1) return y;
-        if(I==2) return z;
-        return x;
-    }
-
-    __host__ __device__ cuvect3 cuvect3::operator+(const cuvect3 &rhs)
-    {
-        cuvect3 ret;
-        ret.x = x+rhs.x;
-        ret.y = y+rhs.y;
-        ret.z = z+rhs.z;
-
-        return ret;
-    }
-
-    __host__ __device__ cuvect3 cuvect3::operator-(const cuvect3 &rhs)
-    {
-        cuvect3 ret;
-        ret.x = x-rhs.x;
-        ret.y = y-rhs.y;
-        ret.z = z-rhs.z;
-
-        return ret;
-    }
-
-    __host__ __device__ cuvect3 cuvect3::operator*(const double &rhs)
-    {
-        cuvect3 ret;
-        ret.x = x*rhs;
-        ret.y = y*rhs;
-        ret.z = z*rhs;
-        return ret;
-    }
-
-    __host__ __device__ cuvect3 cuvect3::operator/(const double &rhs)
-    {
-        cuvect3 ret;
-        ret.x = x/rhs;
-        ret.y = y/rhs;
-        ret.z = z/rhs;
-        return ret;
-    }
-
-    __host__ __device__ cuvect3& cuvect3::operator+=(const cuvect3 &rhs)
-    {
-        x = x + rhs.x;
-        y = y + rhs.y;
-        z = z + rhs.z;
-        return *this;
-    }
-
-    __host__ __device__ cuvect3& cuvect3::operator-=(const cuvect3 &rhs)
-    {
-        x = x - rhs.x;
-        y = y - rhs.y;
-        z = z - rhs.z;
-        return *this;
-    }
-
-    __host__ __device__ cuvect3& cuvect3::operator*=(const double &rhs)
-    {
-        x = x * rhs;
-        y = y * rhs;
-        z = z * rhs;
-        return *this;
-    }
-
-    __host__ __device__ cuvect3& cuvect3::operator/=(const double &rhs)
-    {
-        x = x / rhs;
-        y = y / rhs;
-        z = z / rhs;
-        return *this;
-    }
-
-
-    __host__ __device__ cuvect3::cuvect3(const double &_x, const double &_y, const double &_z)
-    {
-        x = _x; y = _y; z = _z;
-        return;
-    }
-
-
-    __host__ __device__ double cuvect3_dot(const cuvect3 &a, const cuvect3 &b)
-    {
-        double ret = a.x*b.x+a.y*b.y+a.z*b.z;
-
-        return ret;
-    }
-
-    __host__ __device__ cuvect3 cuvect3_cross(const cuvect3 &a, const cuvect3 &b)
-    {
-        cuvect3 ret;
-        ret[0] = a[1]*b[2]-a[2]*b[1];
-        ret[1] = a[2]*b[0]-a[0]*b[2];
-        ret[2] = a[0]*b[1]-a[1]*b[0];
-
-        return ret;
-    }
-
-    __host__ __device__ double cuvect3_norm(const cuvect3 &a)
-    {
-        double ret;
-        ret = ::sqrtf(a.x*a.x+a.y*a.y+a.z*a.z);
-        return ret;
-    }
-
-    __host__ __device__ cuvect3 cuvect3_normalize(const cuvect3 &a)
-    {
-        cuvect3 ret;
-        double m;
-        m = ::sqrtf(a.x*a.x+a.y*a.y+a.z*a.z);
-        if(m>0.0)
-        {
-            ret.x = a.x/m; ret.y = a.y/m; ret.z = a.z/m;
-        }
-        else
-        {
-            ret.x = 0.0; ret.y = 0.0; ret.z = 0.0;
-        }
-
-        return ret;
-    }
-
-    __host__ __device__ cuvect3 cuvect3_proj(const cuvect3 &a, const cuvect3 &b)
-    {
-        cuvect3 ret;
-        cuvect3 bn = cuvect3_normalize(b);
-        double m = cuvect3_dot(a,bn);
-        ret = bn*m;
-        return ret;
-    }
-
-
-
-
-__host__ __device__ cumat3::cumat3()
-{
-    m00 = 0.0;
-    m01 = 0.0;
-    m02 = 0.0;
-    m10 = 0.0;
-    m11 = 0.0;
-    m12 = 0.0;
-    m20 = 0.0;
-    m21 = 0.0;
-    m22 = 0.0;
-
-    return;
-}
-
-__host__ __device__ cumat3::~cumat3()
-{
-    m00 = 0.0;
-    m01 = 0.0;
-    m02 = 0.0;
-    m10 = 0.0;
-    m11 = 0.0;
-    m12 = 0.0;
-    m20 = 0.0;
-    m21 = 0.0;
-    m22 = 0.0;
-    return;
-}
-
-__host__ __device__ double& cumat3::operator[](const int &I)
-{
-    if(I==0) return m00;
-    if(I==1) return m10;
-    if(I==2) return m20;
-    if(I==3) return m01;
-    if(I==4) return m11;
-    if(I==5) return m21;
-    if(I==6) return m02;
-    if(I==7) return m12;
-    if(I==8) return m22;
-    
-    return m00;
-}
-
-__host__ __device__ double& cumat3::operator()(const int &I, const int &J)
-{
-    if(I==0 && J==0) return m00;
-    if(I==1 && J==0) return m10;
-    if(I==2 && J==0) return m20;
-    if(I==0 && J==1) return m01;
-    if(I==1 && J==1) return m11;
-    if(I==2 && J==1) return m21;
-    if(I==0 && J==2) return m02;
-    if(I==1 && J==2) return m12;
-    if(I==2 && J==2) return m22;
-
-    return m00;
-}
-
-
-__host__ __device__ double& cumat3::at(const int &I, const int &J)
-{
-    if(I==0 && J==0) return m00;
-    if(I==1 && J==0) return m10;
-    if(I==2 && J==0) return m20;
-    if(I==0 && J==1) return m01;
-    if(I==1 && J==1) return m11;
-    if(I==2 && J==1) return m21;
-    if(I==0 && J==2) return m02;
-    if(I==1 && J==2) return m12;
-    if(I==2 && J==2) return m22;
-
-    return m00;
-}
-
-__host__ __device__ const double& cumat3::operator[](const int &I) const 
-{
-    if(I==0) return m00;
-    if(I==1) return m10;
-    if(I==2) return m20;
-    if(I==3) return m01;
-    if(I==4) return m11;
-    if(I==5) return m21;
-    if(I==6) return m02;
-    if(I==7) return m12;
-    if(I==8) return m22;
-    
-    return m00;
-}
-
-__host__ __device__ const double& cumat3::operator()(const int &I, const int &J) const 
-{
-    if(I==0 && J==0) return m00;
-    if(I==1 && J==0) return m10;
-    if(I==2 && J==0) return m20;
-    if(I==0 && J==1) return m01;
-    if(I==1 && J==1) return m11;
-    if(I==2 && J==1) return m21;
-    if(I==0 && J==2) return m02;
-    if(I==1 && J==2) return m12;
-    if(I==2 && J==2) return m22;
-
-    return m00;
-}
-
-
-__host__ __device__ const double& cumat3::at(const int &I, const int &J) const 
-{
-    if(I==0 && J==0) return m00;
-    if(I==1 && J==0) return m10;
-    if(I==2 && J==0) return m20;
-    if(I==0 && J==1) return m01;
-    if(I==1 && J==1) return m11;
-    if(I==2 && J==1) return m21;
-    if(I==0 && J==2) return m02;
-    if(I==1 && J==2) return m12;
-    if(I==2 && J==2) return m22;
-
-    return m00;
-}
-
-
-__host__ __device__ cumat3 cumat3::operator+(const cumat3 &rhs)
-{
-    cumat3 ret;
-    ret.m00 = m00 + rhs.m00;
-    ret.m10 = m10 + rhs.m10;
-    ret.m20 = m20 + rhs.m20;
-    ret.m01 = m01 + rhs.m01;
-    ret.m11 = m11 + rhs.m11;
-    ret.m21 = m21 + rhs.m21;
-    ret.m02 = m02 + rhs.m02;
-    ret.m12 = m12 + rhs.m12;
-    ret.m22 = m22 + rhs.m22;
-    
-
-    return ret;
-}
-
-__host__ __device__ cumat3 cumat3::operator-(const cumat3 &rhs)
-{
-    cumat3 ret;
-    ret.m00 = m00 - rhs.m00;
-    ret.m10 = m10 - rhs.m10;
-    ret.m20 = m20 - rhs.m20;
-    ret.m01 = m01 - rhs.m01;
-    ret.m11 = m11 - rhs.m11;
-    ret.m21 = m21 - rhs.m21;
-    ret.m02 = m02 - rhs.m02;
-    ret.m12 = m12 - rhs.m12;
-    ret.m22 = m22 - rhs.m22;
-    return ret;
-}
-
-__host__ __device__ cumat3 cumat3::operator*(const double &rhs)
-{
-    cumat3 ret;
-    ret.m00 = m00 * rhs;
-    ret.m10 = m10 * rhs;
-    ret.m20 = m20 * rhs;
-    ret.m01 = m01 * rhs;
-    ret.m11 = m11 * rhs;
-    ret.m21 = m21 * rhs;
-    ret.m02 = m02 * rhs;
-    ret.m12 = m12 * rhs;
-    ret.m22 = m22 * rhs;
-    return ret;
-}
-
-__host__ __device__ cumat3 cumat3::operator/(const double &rhs)
-{
-    cumat3 ret;
-    ret.m00 = m00 / rhs;
-    ret.m10 = m10 / rhs;
-    ret.m20 = m20 / rhs;
-    ret.m01 = m01 / rhs;
-    ret.m11 = m11 / rhs;
-    ret.m21 = m21 / rhs;
-    ret.m02 = m02 / rhs;
-    ret.m12 = m12 / rhs;
-    ret.m22 = m22 / rhs;
-    return ret;
-}
-        
-__host__ __device__ cuvect3 cumat3::operator*(const cuvect3 &rhs)
-{
-    cuvect3 ret;
-    
-    ret.x = m00*rhs.x + m01*rhs.y + m02*rhs.z;
-    ret.y = m10*rhs.x + m11*rhs.y + m12*rhs.z;
-    ret.z = m20*rhs.x + m21*rhs.y + m22*rhs.z;
-
-    return ret;
-}
-
-__host__ __device__ cumat3 cumat3::operator*(const cumat3 &rhs)
-{
-    cumat3 ret;
-
-    ret.m00 = m00*rhs.m00 + m01*rhs.m10 + m02*rhs.m20;
-    ret.m01 = m00*rhs.m01 + m01*rhs.m11 + m02*rhs.m21;
-    ret.m02 = m00*rhs.m02 + m01*rhs.m12 + m02*rhs.m22;
-    ret.m10 = m10*rhs.m00 + m11*rhs.m10 + m12*rhs.m20;
-    ret.m11 = m10*rhs.m01 + m11*rhs.m11 + m12*rhs.m21;
-    ret.m12 = m10*rhs.m02 + m11*rhs.m12 + m12*rhs.m22;
-    ret.m20 = m20*rhs.m00 + m21*rhs.m10 + m22*rhs.m20;
-    ret.m21 = m20*rhs.m01 + m21*rhs.m11 + m22*rhs.m21;
-    ret.m22 = m20*rhs.m02 + m21*rhs.m12 + m22*rhs.m22;
-
-    return ret;
-}
-
-__host__ __device__ double cumat3::det()
-{
-    double ret = 0.0;
-
-    ret += m00*m11*m22;
-    ret += m01*m12*m20;
-    ret += m02*m10*m21;
-    ret -= m20*m11*m02;
-    ret -= m21*m12*m00;
-    ret -= m22*m10*m01;
-
-    return ret;
-}
-
-__host__ __device__ cumat3 cumat3::transpose()
-{
-    cumat3 ret;
-
-    ret.m00 = m00;
-    ret.m01 = m10;
-    ret.m02 = m20;
-    ret.m10 = m01;
-    ret.m11 = m11;
-    ret.m12 = m21;
-    ret.m20 = m02;
-    ret.m21 = m12;
-    ret.m22 = m22;
-
-    return ret;
-}
-
-__host__ __device__ cumat3 cumat3::inverse()
-{
-    cumat3 q;
-    double dt = det();
-    if(dt!=0.0)
-    {
-        q(0,0) =  (at(1,1)*at(2,2)-at(1,2)*at(2,1))/dt;
-        q(1,0) = -(at(1,0)*at(2,2)-at(1,2)*at(2,0))/dt;
-        q(2,0) =  (at(1,0)*at(2,1)-at(1,1)*at(2,0))/dt;
-        q(0,1) = -(at(0,1)*at(2,2)-at(0,2)*at(2,1))/dt;
-        q(1,1) =  (at(0,0)*at(2,2)-at(0,2)*at(2,0))/dt;
-        q(2,1) = -(at(0,0)*at(2,1)-at(0,1)*at(2,0))/dt;
-        q(0,2) =  (at(0,1)*at(1,2)-at(0,2)*at(1,1))/dt;
-        q(1,2) = -(at(0,0)*at(1,2)-at(0,2)*at(1,0))/dt;
-        q(2,2) =  (at(0,0)*at(1,1)-at(0,1)*at(1,0))/dt;
-        // q(0,0) =  (at(1,1)*at(2,2)-at(1,2)*at(2,1))/dt;
-        // q(0,1) = -(at(1,0)*at(2,2)-at(1,2)*at(2,0))/dt;
-        // q(0,2) =  (at(1,0)*at(2,1)-at(1,1)*at(2,0))/dt;
-        // q(1,0) = -(at(0,1)*at(2,2)-at(0,2)*at(2,1))/dt;
-        // q(1,1) =  (at(0,0)*at(2,2)-at(0,2)*at(2,0))/dt;
-        // q(1,2) = -(at(0,0)*at(2,1)-at(0,1)*at(2,0))/dt;
-        // q(2,0) =  (at(0,1)*at(1,2)-at(0,2)*at(1,1))/dt;
-        // q(2,1) = -(at(0,0)*at(1,2)-at(0,2)*at(1,0))/dt;
-        // q(2,2) =  (at(0,0)*at(1,1)-at(0,1)*at(1,0))/dt;
-        // q = q.transpose();
-    }
-    else
-    {
-        q(0,0) = inf;
-        q(0,1) = inf;
-        q(0,2) = inf;
-        q(1,0) = inf;
-        q(1,1) = inf;
-        q(1,2) = inf;
-        q(2,0) = inf;
-        q(2,1) = inf;
-        q(2,2) = inf;
-    }
-
-    return q;
-}
-
-__host__ __device__ cumat3 operator-(const cumat3 &rhs)
-{
-    cumat3 ret;
-    ret.m00 = -rhs.m00;
-    ret.m10 = -rhs.m10;
-    ret.m20 = -rhs.m20;
-    ret.m01 = -rhs.m01;
-    ret.m11 = -rhs.m11;
-    ret.m21 = -rhs.m21;
-    ret.m02 = -rhs.m02;
-    ret.m12 = -rhs.m12;
-    ret.m22 = -rhs.m22;
-
-    return ret;
-}
-
-__host__ __device__ cumat3& cumat3::operator+=(const cumat3 &rhs)
-{
-    m00 += rhs.m00;
-    m10 += rhs.m10;
-    m20 += rhs.m20;
-    m01 += rhs.m01;
-    m11 += rhs.m11;
-    m21 += rhs.m21;
-    m02 += rhs.m02;
-    m12 += rhs.m12;
-    m22 += rhs.m22;
-
-    return *this;
-}
-
-__host__ __device__ cumat3& cumat3::operator-=(const cumat3 &rhs)
-{
-    m00 -= rhs.m00;
-    m10 -= rhs.m10;
-    m20 -= rhs.m20;
-    m01 -= rhs.m01;
-    m11 -= rhs.m11;
-    m21 -= rhs.m21;
-    m02 -= rhs.m02;
-    m12 -= rhs.m12;
-    m22 -= rhs.m22;
-
-    return *this;
-}
-
-
-
-__host__ __device__ cumat3& cumat3::operator/=(const double &rhs)
-{
-    m00 /= rhs;
-    m10 /= rhs;
-    m20 /= rhs;
-    m01 /= rhs;
-    m11 /= rhs;
-    m21 /= rhs;
-    m02 /= rhs;
-    m12 /= rhs;
-    m22 /= rhs;
-
-    return *this;
-}
-
-__host__ __device__ cumat3& cumat3::operator*=(const double &rhs)
-{
-    m00 *= rhs;
-    m10 *= rhs;
-    m20 *= rhs;
-    m01 *= rhs;
-    m11 *= rhs;
-    m21 *= rhs;
-    m02 *= rhs;
-    m12 *= rhs;
-    m22 *= rhs;
-
-    return *this;
-}
-
-__host__ __device__ cumat3& cumat3::operator*=(const cumat3 &rhs)
-{
-    cumat3 tmp;
-
-    tmp.m00 = m00*rhs.m00 + m01*rhs.m10 + m02*rhs.m20;
-    tmp.m01 = m00*rhs.m01 + m01*rhs.m11 + m02*rhs.m21;
-    tmp.m02 = m00*rhs.m02 + m01*rhs.m12 + m02*rhs.m22;
-    tmp.m10 = m10*rhs.m00 + m11*rhs.m10 + m12*rhs.m20;
-    tmp.m11 = m10*rhs.m01 + m11*rhs.m11 + m12*rhs.m21;
-    tmp.m12 = m10*rhs.m02 + m11*rhs.m12 + m12*rhs.m22;
-    tmp.m20 = m20*rhs.m00 + m21*rhs.m10 + m22*rhs.m20;
-    tmp.m21 = m20*rhs.m01 + m21*rhs.m11 + m22*rhs.m21;
-    tmp.m22 = m20*rhs.m02 + m21*rhs.m12 + m22*rhs.m22;
-
-    (*this) = tmp;
-
-    return *this;
-}
-
-__host__ __device__ cumat3::cumat3(
-    const double & _m00, const double & _m01, const double & _m02,
-    const double & _m10, const double & _m11, const double & _m12,
-    const double & _m20, const double & _m21, const double & _m22
-)
-{
-    m00 = _m00;
-    m01 = _m01;
-    m02 = _m02;
-    m10 = _m10;
-    m11 = _m11;
-    m12 = _m12;
-    m20 = _m20;
-    m21 = _m21;
-    m22 = _m22;
-}
-
-__host__ __device__ double* cumat3::data()
-{
-     //pointer to double[9] representation of matrix
-    return (double*) this;
-}
-
-__host__ __device__ const double* cumat3::data() const
-{
-     //pointer to double[9] representation of matrix
-    return (const double*) this;
-}
-
-//transposes a 3x3 (9 element) matrix
-__host__ __device__ void mat3f_transpose(double *mat3inout)
-{
-    int I,J;
-    double matint[9];
-    for(I=0;I<9;I++)
-    {
-        matint[I] = mat3inout[I];
-    }
-
-    for(I=0;I<3;I++)
-    {
-        for(J=0;J<3;J++)
-        {
-            mat3inout[I+J*3] = matint[J+I*3];
-        }
-    }
-
-    return;
-}
-
-//copies src to dest
-__host__ __device__ void mat3f_copy(double *mat3f_dest, const double *mat3f_src)
-{
-    int I;
-    if(mat3f_dest==NULL || mat3f_src==NULL)
-        return;
-
-    for(I=0;I<9;I++)
-        mat3f_dest[I] = mat3f_src[I];
-
-    return;
-}
-
-
-__host__ __device__ double mat3f_det(double *mat3in)
-{
-    double ret = 0.0;
-
-    ret = ret + mat3in[0+0*3]*mat3in[1+1*3]*mat3in[2+2*3];
-    ret = ret + mat3in[0+1*3]*mat3in[1+2*3]*mat3in[2+0*3];
-    ret = ret + mat3in[0+2*3]*mat3in[1+0*3]*mat3in[2+1*3];
-    ret = ret - mat3in[0+0*3]*mat3in[1+2*3]*mat3in[2+1*3];
-    ret = ret - mat3in[0+1*3]*mat3in[1+0*3]*mat3in[2+2*3];
-    ret = ret - mat3in[0+2*3]*mat3in[1+1*3]*mat3in[2+0*3];
-
-    return ret;
-}
-
-//inverts a 3x3 (9 element) matrix
-__host__ __device__ void mat3f_inverse(double *mat3inout)
-{
-    int I;
-    double matint[9];
-    double det = mat3f_det(mat3inout);
-
-    for(I=0;I<9;I++)
-    {
-        matint[I] = mat3inout[I];
-    }
-
-    mat3inout[0+0*3] =  (matint[1+1*3]*matint[2+2*3]-matint[1+2*3]*matint[2+1*3])/det;
-    mat3inout[0+1*3] = -(matint[1+0*3]*matint[2+2*3]-matint[1+2*3]*matint[2+0*3])/det;
-    mat3inout[0+2*3] =  (matint[1+0*3]*matint[2+1*3]-matint[1+1*3]*matint[2+0*3])/det;
-    mat3inout[1+0*3] = -(matint[0+1*3]*matint[2+2*3]-matint[0+2*3]*matint[2+1*3])/det;
-    mat3inout[1+1*3] =  (matint[0+0*3]*matint[2+2*3]-matint[0+2*3]*matint[2+0*3])/det;
-    mat3inout[1+2*3] = -(matint[0+0*3]*matint[2+1*3]-matint[0+1*3]*matint[2+0*3])/det;
-    mat3inout[2+0*3] =  (matint[0+1*3]*matint[1+2*3]-matint[0+2*3]*matint[1+1*3])/det;
-    mat3inout[2+1*3] = -(matint[0+0*3]*matint[1+2*3]-matint[0+2*3]*matint[1+0*3])/det;
-    mat3inout[2+2*3] =  (matint[0+0*3]*matint[1+1*3]-matint[0+1*3]*matint[1+0*3])/det;
-
-    mat3f_transpose(mat3inout);
-
-    return;
-}
-
-__host__ __device__ cuvect3 mat3f_mult(double *mat3in, const cuvect3 &cvin)
-{
-    int I,J;
-    cuvect3 ret;
-    for(I=0;I<3;I++)
-    {
-        ret[I] = 0.0;
-        for(J=0;J<3;J++)
-        {
-            ret[I] = ret[I] + mat3in[I+3*J]*cvin[J];
-        }
-    }
-
-    return ret;
-}
-
-__host__ __device__ void mat3f_mult(double *matina, double *matinb, double *matout)
-{
-    double wrk[9];
-    int I,J,K;
-
-    for(I=0;I<3;I++)
-    {
-        for(J=0;J<3;J++)
-        {
-            wrk[I+3*J] = 0.0;
-        }
-    }
-
-    for(I=0;I<3;I++)
-    {
-        for(J=0;J<3;J++)
-        {
-            for(K=0;K<3;K++)
-            {
-                wrk[I+3*K] = wrk[I+3*K] + matina[I+3*J]*matinb[J+3*K];
-            }
-        }
-    }
-
-    for(I=0;I<3;I++)
-    {
-        for(J=0;J<3;J++)
-        {
-            matout[I+3*J] = wrk[I+3*J];
-        }
-    }
-
-    return;
-}
-
-__host__ void test_cudavect_logic1()
-{
-    //3 dim vector and matrix functional tests on host side
-
-    // printf("3 dim vector and matrix functional tests on host side\n");
-    
-    // cuvect3 a,b,c;
-    // double ma[9],mb[9],mc[9];
-
-    // int I,J;
-
-    // for(I=0;I<3;I++)
-    // {
-    //     for(J=0;J<3;J++)
-    //     {
-    //         ma[I+3*J] = ((double) rand())/((double) RAND_MAX);
-    //         mb[I+3*J] = ma[I+3*J];
-    //     }
-    // }
-
-    // mat3f_inverse(mb);
-    // mat3f_mult(ma,mb,mc);
-
-    // for(I=0;I<3;I++)
-    // {
-    //     for(J=0;J<3;J++)
-    //     {
-    //         printf("ma[%d,%d] = %1.3f\n",I,J,ma[I+3*J]);
-    //     }
-    // }
-    // for(I=0;I<3;I++)
-    // {
-    //     for(J=0;J<3;J++)
-    //     {
-    //         printf("mb[%d,%d] = %1.3f\n",I,J,mb[I+3*J]);
-    //     }
-    // }
-    // for(I=0;I<3;I++)
-    // {
-    //     for(J=0;J<3;J++)
-    //     {
-    //         printf("mc[%d,%d] = %1.3f\n",I,J,mc[I+3*J]);
-    //     }
-    // }
-
-    // a = cuvect3(1,1,1);
-    // b = mat3f_mult(ma,a);
-    // b = mat3f_mult(mb,b);
-    
-    // for(I=0;I<3;I++)
-    // {
-    //     printf("a[%d] = %1.3f, b[%d] = %1.3f\n",I,a[I],I,b[I]);
-    // }
-
-    // a = cuvect3(1,0,1);
-    // b = cuvect3(0,1,-1);
-    // c = a+b;
-
-    // for(I=0;I<3;I++)
-    // {
-    //     printf("a[%d] = %1.3f, b[%d] = %1.3f, c[%d] = %1.3f\n",I,a[I],I,b[I],I,c[I]);
-    // }
-
-    // c = c/2.0;
-
-    // for(I=0;I<3;I++)
-    // {
-    //     printf("a[%d] = %1.3f, b[%d] = %1.3f, c[%d] = %1.3f\n",I,a[I],I,b[I],I,c[I]);
-    // }
-
-    // c = cuvect3_cross(a,b);
-
-    // for(I=0;I<3;I++)
-    // {
-    //     printf("a[%d] = %1.3f, b[%d] = %1.3f, c[%d] = %1.3f\n",I,a[I],I,b[I],I,c[I]);
-    // }
-
-    // printf("c dot a = %1.3f, c dot b = %1.3f\n",cuvect3_dot(c,a),cuvect3_dot(c,b));
-
-    // printf("norm(a)=%1.3f, norm(b)=%1.3f, norm(c)=%1.3f\n",cuvect3_norm(a),cuvect3_norm(b),cuvect3_norm(c));
-    // a = cuvect3_normalize(a);
-    // b = cuvect3_normalize(b);
-    // c = cuvect3_normalize(c);
-
-    // for(I=0;I<3;I++)
-    // {
-    //     printf("a[%d] = %1.3f, b[%d] = %1.3f, c[%d] = %1.3f\n",I,a[I],I,b[I],I,c[I]);
-    // }
-    // printf("c dot a = %1.3f, c dot b = %1.3f\n",cuvect3_dot(c,a),cuvect3_dot(c,b));
-    // printf("norm(a)=%1.3f, norm(b)=%1.3f, norm(c)=%1.3f\n",cuvect3_norm(a),cuvect3_norm(b),cuvect3_norm(c));
-
-    return;
-}
-
-__host__ __device__ cumat3 hodge_dual(const cuvect3 &vin)
-{
-    cumat3 ret;
-
-    ret.m00 = 0.0;
-    ret.m11 = 0.0;
-    ret.m22 = 0.0;
-
-    ret.m01 = vin.z;
-    ret.m12 = vin.x;
-    ret.m20 = vin.y;
-
-    ret.m10 = -vin.z;
-    ret.m21 = -vin.x;
-    ret.m02 = -vin.y;
-
-    return ret;
-}
-
-__host__ __device__ cuvect3 hodge_dual(const cumat3 &min)
-{
-    cuvect3 ret;
-
-    ret.x = 0.5f*(min.m12 - min.m21);
-    ret.y = 0.5f*(min.m20 - min.m02);
-    ret.z = 0.5f*(min.m01 - min.m10);
-
-    return ret;
-}
-
-__host__ __device__ const cumat3 cumat3_eye()
-{
-    cumat3 ret;
-    ret.m00 = 1.0f;
-    ret.m11 = 1.0f;
-    ret.m22 = 1.0f;
-    
-    return ret;
-}
-
-__host__ __device__ const cumat3 cumat3_zeros()
-{
-    cumat3 ret;
-    return ret;
-}
-
-__host__ __device__ cumat3 rotmat_from_axisangle(const cuvect3 &axis, const double &angle)
-{
-    cumat3 ret = cumat3_zeros();
-    cumat3 H;
-    cuvect3 _naxis;
-
-    _naxis = cuvect3_normalize(axis);
-    H = hodge_dual(_naxis);
-
-    ret += H*sinf(angle);
-    H = H*H;
-    ret -= H*cosf(angle);
-    ret += (H + cumat3_eye());
-
-    return ret;
-}
-
-};
\ No newline at end of file
diff --git a/src/amsculib3/math/cuvect4.cu b/src/amsculib3/math/cuvect4.cu
deleted file mode 100644
index 859de8d..0000000
--- a/src/amsculib3/math/cuvect4.cu
+++ /dev/null
@@ -1,432 +0,0 @@
-#include <amsculib3/amsculib3.hpp>
-
-namespace amscuda
-{
-
-__host__ __device__ cuvect4::cuvect4()
-{
-    x = 0.0; y = 0.0; z = 0.0; w = 0.0;
-    return;
-}
-
-__host__ __device__ cuvect4::~cuvect4()
-{
-    x = 0.0; y = 0.0; z = 0.0; w = 0.0;
-    return;
-}
-
-__host__ __device__ cuvect4::cuvect4(double _x, double _y, double _z, double _w)
-{
-    x = _x; y = _y; z = _z; w = _w;
-    return;
-}
-        
-__host__ __device__ double& cuvect4::operator[](const int I)
-{
-    if(I==0) return x;
-    else if(I==1) return y;
-    else if(I==2) return z;
-    else if(I==3) return w;
-    return x;
-}
-
-__host__ __device__ const double& cuvect4::operator[](const int I) const
-{
-    if(I==0) return x;
-    else if(I==1) return y;
-    else if(I==2) return z;
-    else if(I==3) return w;
-    return x;
-}
-
-__host__ __device__ cuvect4 cuvect4::operator+(cuvect4 lhs)
-{
-    cuvect4 ret;
-    ret.x = this->x + lhs.x;
-    ret.y = this->y + lhs.y;
-    ret.z = this->z + lhs.z;
-    ret.w = this->w + lhs.w;
-    return ret;
-}
-
-__host__ __device__ cuvect4 cuvect4::operator-(cuvect4 lhs)
-{
-    cuvect4 ret;
-    ret.x = this->x - lhs.x;
-    ret.y = this->y - lhs.y;
-    ret.z = this->z - lhs.z;
-    ret.w = this->w - lhs.w;
-    return ret;
-}
-
-__host__ __device__ cuvect4 cuvect4::operator*(double lhs)
-{
-    cuvect4 ret;
-    ret.x = this->x*lhs;
-    ret.y = this->y*lhs;
-    ret.z = this->z*lhs;
-    ret.w = this->w*lhs;
-    return ret;
-}
-
-__host__ __device__ cuvect4 cuvect4::operator/(double lhs)
-{
-    cuvect4 ret;
-    ret.x = this->x/lhs;
-    ret.y = this->y/lhs;
-    ret.z = this->z/lhs;
-    ret.w = this->w/lhs;
-    return ret;
-}
-
-__host__ __device__ cumat4::cumat4()
-{
-    int I;
-    for(I=0;I<16;I++)
-    {
-        dat[I] = 0.0;
-    }
-    return;
-}
-
-__host__ __device__ cumat4::~cumat4()
-{
-    int I;
-    for(I=0;I<16;I++)
-    {
-        dat[I] = 0.0;
-    }
-    return;
-}
-
-__host__ __device__ double& cumat4::operator[](const int I)
-{
-    return dat[I];
-}
-
-__host__ __device__ double& cumat4::operator()(const int I, const int J)
-{
-    return dat[I+4*J];
-}
-
-__host__ __device__ double& cumat4::at(const int I, const int J)
-{
-    return dat[I+4*J];
-}
-
-__host__ __device__ cumat4 cumat4::operator+(cumat4 lhs)
-{
-    cumat4 ret;
-    int I;
-    for(I=0;I<16;I++)
-    {
-        ret.dat[I] = this->dat[I] + lhs.dat[I];
-    }
-    return ret;
-}
-
-__host__ __device__ cumat4 cumat4::operator-(cumat4 lhs)
-{
-    cumat4 ret;
-    int I;
-    for(I=0;I<16;I++)
-    {
-        ret.dat[I] = this->dat[I] - lhs.dat[I];
-    }
-    return ret;
-}
-
-__host__ __device__ cumat4 cumat4::operator*(double lhs)
-{
-    cumat4 ret;
-    int I;
-    for(I=0;I<16;I++)
-    {
-        ret.dat[I] = this->dat[I]*lhs;
-    }
-    return ret;
-}
-
-__host__ __device__ cumat4 cumat4::operator/(double lhs)
-{
-    cumat4 ret;
-    int I;
-    for(I=0;I<16;I++)
-    {
-        ret.dat[I] = this->dat[I]/lhs;
-    }
-    return ret;
-}
-
-__host__ __device__ cuvect4 cumat4::operator*(cuvect4 lhs)
-{
-    cuvect4 ret = cuvect4(0.0,0.0,0.0,0.0);
-    int I,J;
-
-    for(I=0;I<4;I++)
-    {
-        for(J=0;J<4;J++)
-        {
-            ret[I] = ret[I] + this->at(I,J)*lhs[J];
-        }
-    }
-
-    return ret;
-}
-
-__host__ __device__ cumat4 cumat4::operator*(cumat4 lhs)
-{
-    cumat4 ret;
-    int I,J,K;
-    for(I=0;I<4;I++)
-    {
-        for(J=0;J<4;J++)
-        {
-            ret(I,J) = 0;
-            for(K=0;K<4;K++)
-            {
-                ret(I,J) = ret(I,J) + this->at(I,K) * lhs(K,J);
-            }
-        }
-    }
-    return ret;
-}
-
-__host__ __device__ cumat4 cumat4::transpose()
-{
-    cumat4 q;
-    int I,J;
-    for(I=0;I<4;I++)
-    {
-        for(J=0;J<4;J++)
-        {
-            q(I,J) = this->at(J,I);
-        }
-    }
-    return q;
-}
-
-__host__ __device__ double cumat4::det()
-{
-    double a00,a01,a02,a03;
-    double a10,a11,a12,a13;
-    double a20,a21,a22,a23;
-    double a30,a31,a32,a33;
-    double det;
-
-    a00 = this->at(0,0);
-    a01 = this->at(0,1);
-    a02 = this->at(0,2);
-    a03 = this->at(0,3);
-    a10 = this->at(1,0);
-    a11 = this->at(1,1);
-    a12 = this->at(1,2);
-    a13 = this->at(1,3);
-    a20 = this->at(2,0);
-    a21 = this->at(2,1);
-    a22 = this->at(2,2);
-    a23 = this->at(2,3);
-    a30 = this->at(3,0);
-    a31 = this->at(3,1);
-    a32 = this->at(3,2);
-    a33 = this->at(3,3);
-
-    det = a03*a12*a21*a30 -
-        a02*a13*a21*a30 -
-        a03*a11*a22*a30 +
-        a01*a13*a22*a30 +
-        a02*a11*a23*a30 -
-        a01*a12*a23*a30 -
-        a03*a12*a20*a31 +
-        a02*a13*a20*a31 +
-        a03*a10*a22*a31 -
-        a00*a13*a22*a31 -
-        a02*a10*a23*a31 +
-        a00*a12*a23*a31 +
-        a03*a11*a20*a32 -
-        a01*a13*a20*a32 -
-        a03*a10*a21*a32 +
-        a00*a13*a21*a32 +
-        a01*a10*a23*a32 -
-        a00*a11*a23*a32 -
-        a02*a11*a20*a33 +
-        a01*a12*a20*a33 +
-        a02*a10*a21*a33 -
-        a00*a12*a21*a33 -
-        a01*a10*a22*a33 +
-        a00*a11*a22*a33;
-
-    return det;
-}
-
-__host__ __device__ cumat4 minverse(cumat4 ma)
-{
-  cumat4 mb;
-
-  double a00,a01,a02,a03;
-  double a10,a11,a12,a13;
-  double a20,a21,a22,a23;
-  double a30,a31,a32,a33;
-
-  double b00,b01,b02,b03;
-  double b10,b11,b12,b13;
-  double b20,b21,b22,b23;
-  double b30,b31,b32,b33;
-
-  double det = 0.0;
-
-  a00 = ma.at(0,0);
-  a01 = ma.at(0,1);
-  a02 = ma.at(0,2);
-  a03 = ma.at(0,3);
-  a10 = ma.at(1,0);
-  a11 = ma.at(1,1);
-  a12 = ma.at(1,2);
-  a13 = ma.at(1,3);
-  a20 = ma.at(2,0);
-  a21 = ma.at(2,1);
-  a22 = ma.at(2,2);
-  a23 = ma.at(2,3);
-  a30 = ma.at(3,0);
-  a31 = ma.at(3,1);
-  a32 = ma.at(3,2);
-  a33 = ma.at(3,3);
-
-  det = a03*a12*a21*a30 -
-        a02*a13*a21*a30 -
-        a03*a11*a22*a30 +
-        a01*a13*a22*a30 +
-        a02*a11*a23*a30 -
-        a01*a12*a23*a30 -
-        a03*a12*a20*a31 +
-        a02*a13*a20*a31 +
-        a03*a10*a22*a31 -
-        a00*a13*a22*a31 -
-        a02*a10*a23*a31 +
-        a00*a12*a23*a31 +
-        a03*a11*a20*a32 -
-        a01*a13*a20*a32 -
-        a03*a10*a21*a32 +
-        a00*a13*a21*a32 +
-        a01*a10*a23*a32 -
-        a00*a11*a23*a32 -
-        a02*a11*a20*a33 +
-        a01*a12*a20*a33 +
-        a02*a10*a21*a33 -
-        a00*a12*a21*a33 -
-        a01*a10*a22*a33 +
-        a00*a11*a22*a33;
-
-  if(det*det>1.0E-30)
-  {
-    b00 = -a13*a22*a31 + a12*a23*a31 + a13*a21*a32 - a11*a23*a32 - a12*a21*a33 + a11*a22*a33;
-    b01 =  a03*a22*a31 - a02*a23*a31 - a03*a21*a32 + a01*a23*a32 + a02*a21*a33 - a01*a22*a33;
-    b02 = -a03*a12*a31 + a02*a13*a31 + a03*a11*a32 - a01*a13*a32 - a02*a11*a33 + a01*a12*a33;
-    b03 =  a03*a12*a21 - a02*a13*a21 - a03*a11*a22 + a01*a13*a22 + a02*a11*a23 - a01*a12*a23;
-    b10 =  a13*a22*a30 - a12*a23*a30 - a13*a20*a32 + a10*a23*a32 + a12*a20*a33 - a10*a22*a33;
-    b11 = -a03*a22*a30 + a02*a23*a30 + a03*a20*a32 - a00*a23*a32 - a02*a20*a33 + a00*a22*a33;
-    b12 =  a03*a12*a30 - a02*a13*a30 - a03*a10*a32 + a00*a13*a32 + a02*a10*a33 - a00*a12*a33;
-    b13 = -a03*a12*a20 + a02*a13*a20 + a03*a10*a22 - a00*a13*a22 - a02*a10*a23 + a00*a12*a23;
-    b20 = -a13*a21*a30 + a11*a23*a30 + a13*a20*a31 - a10*a23*a31 - a11*a20*a33 + a10*a21*a33;
-    b21 =  a03*a21*a30 - a01*a23*a30 - a03*a20*a31 + a00*a23*a31 + a01*a20*a33 - a00*a21*a33;
-    b22 = -a03*a11*a30 + a01*a13*a30 + a03*a10*a31 - a00*a13*a31 - a01*a10*a33 + a00*a11*a33;
-    b23 =  a03*a11*a20 - a01*a13*a20 - a03*a10*a21 + a00*a13*a21 + a01*a10*a23 - a00*a11*a23;
-    b30 =  a12*a21*a30 - a11*a22*a30 - a12*a20*a31 + a10*a22*a31 + a11*a20*a32 - a10*a21*a32;
-    b31 = -a02*a21*a30 + a01*a22*a30 + a02*a20*a31 - a00*a22*a31 - a01*a20*a32 + a00*a21*a32;
-    b32 =  a02*a11*a30 - a01*a12*a30 - a02*a10*a31 + a00*a12*a31 + a01*a10*a32 - a00*a11*a32;
-    b33 = -a02*a11*a20 + a01*a12*a20 + a02*a10*a21 - a00*a12*a21 - a01*a10*a22 + a00*a11*a22;
-    b00 = b00/det;
-    b01 = b01/det;
-    b02 = b02/det;
-    b03 = b03/det;
-    b10 = b10/det;
-    b11 = b11/det;
-    b12 = b12/det;
-    b13 = b13/det;
-    b20 = b20/det;
-    b21 = b21/det;
-    b22 = b22/det;
-    b23 = b23/det;
-    b30 = b30/det;
-    b31 = b31/det;
-    b32 = b32/det;
-    b33 = b33/det;
-    mb.at(0,0) = b00;
-    mb.at(0,1) = b01;
-    mb.at(0,2) = b02;
-    mb.at(0,3) = b03;
-    mb.at(1,0) = b10;
-    mb.at(1,1) = b11;
-    mb.at(1,2) = b12;
-    mb.at(1,3) = b13;
-    mb.at(2,0) = b20;
-    mb.at(2,1) = b21;
-    mb.at(2,2) = b22;
-    mb.at(2,3) = b23;
-    mb.at(3,0) = b30;
-    mb.at(3,1) = b31;
-    mb.at(3,2) = b32;
-    mb.at(3,3) = b33;
-  }
-  //this was STUPID. Gaah. Computer algebra system saves the day? I'd be surprised if this didn't end up *slower* than gaussian elimination. Don't do this again!
-  return mb;
-}
-
-__host__ __device__ cumat4 cumat4::inverse()
-{
-    return minverse(*this);
-}
-
-
-__host__ __device__ double cuvect4_dot(cuvect4 a, cuvect4 b)
-{
-    double ret = 0.0;
-    ret = a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w;
-    return ret;
-}
-
-__host__ __device__ double cuvect4_norm(cuvect4 a)
-{
-    double ret = 0.0;
-    ret = ::sqrt(cuvect4_dot(a,a));
-    return ret;
-}
-
-__host__ __device__ cuvect4 cuvect4_normalize(cuvect4 a)
-{
-    cuvect4 ret = cuvect4(0.0f,0.0f,0.0f,0.0f);
-    double nrm = cuvect4_norm(a);
-    if(nrm>0.0)
-        ret = a/nrm;
-    return ret;
-}
-
-__host__ __device__ cuvect4 cuvect4_proj(cuvect4 a, cuvect4 b)
-{
-    cuvect4 ret;
-    cuvect4 bn = cuvect4_normalize(b);
-    double d = cuvect4_dot(a,bn);
-    ret = bn*d;
-    return ret;
-}
-
-__host__ __device__ cuvect4 operator-(cuvect4 rhs)
-{
-    cuvect4 ret;
-    ret[0] = -rhs[0];
-    ret[1] = -rhs[1];
-    ret[2] = -rhs[2];
-    ret[3] = -rhs[3];
-    return ret;
-}
-
-__host__ __device__ cumat4 operator-(cumat4 rhs)
-{
-    cumat4 ret;
-    int I;
-    for(I=0;I<16;I++) ret[I] = -rhs[I];
-    return ret;
-}
-
-
-};
\ No newline at end of file