diff --git a/pykat/SIfloat.py b/pykat/SIfloat.py
index 5a3556b8b548aa89b7d53dca0a2157e457fe6c33..cac016f0623c82601cb7b6c072cb0e3c5b56b010 100644
--- a/pykat/SIfloat.py
+++ b/pykat/SIfloat.py
@@ -1,6 +1,7 @@
 import os
 import re
 import pykat.exceptions as pkex
+import numpy as np
 
 __suffix = {'y': 'e-24',  # yocto
             'z': 'e-21',  # zepto
@@ -22,11 +23,20 @@ __suffix = {'y': 'e-24',  # yocto
 def SIfloat(value):
     if value==None: 
         return value
+        
+    value = np.array(value)
+    
+    v = np.vectorize(convertToFloat)
+    
+    if value.size == 1:
+        return float(v(value))
+        
+    a = v(value)
     
-    if type(value)==list:
-        return [convertToFloat(s) for s in value]
+    if len(a) == 1:
+        return a[0]
     else:
-        return convertToFloat(value)
+        return a
     
 def convertToFloat(value):
     
diff --git a/pykat/utilities/maps.py b/pykat/utilities/maps.py
index 6939b82c587d9faea1838bf51a324350632b47a5..a9e190cadc048a3ca50187814dc2fed6dd14b642 100644
--- a/pykat/utilities/maps.py
+++ b/pykat/utilities/maps.py
@@ -128,7 +128,7 @@ class surfacemap(object):
             raise BasePyKatException("Map type needs handling")
         
 
-    def generateROMWeights(self, isModeMatched=True, verbose=False, interpolate=False, tolerance = 1e-12, sigma = 1, sort=False):
+    def generateROMWeights(self, isModeMatched=True, verbose=False, interpolate=False, tolerance = 1e-12, sigma = 1, sort=False, greedyfile=None):
         
         if interpolate:
             from scipy.interpolate import interp2d
@@ -169,12 +169,12 @@ class surfacemap(object):
         
         EI = {}
         
-        EI["xm"] = makeEmpiricalInterpolant(makeReducedBasis(xm, isModeMatched=isModeMatched, tolerance = tolerance, sigma = sigma), sort=sort)
+        EI["xm"] = makeEmpiricalInterpolant(makeReducedBasis(xm, isModeMatched=isModeMatched, tolerance = tolerance, sigma = sigma, greedyfile=greedyfile), sort=sort)
         
         if symm:
             EI["ym"] = EI["xm"]
         else:
-            EI["ym"] = makeEmpiricalInterpolant(makeReducedBasis(ym, isModeMatched=isModeMatched, tolerance = tolerance, sigma = sigma), sort=sort)
+            EI["ym"] = makeEmpiricalInterpolant(makeReducedBasis(ym, isModeMatched=isModeMatched, tolerance = tolerance, sigma = sigma, greedyfile=greedyfile), sort=sort)
         
         EI["limits"] = EI["xm"].limits
         
diff --git a/pykat/utilities/optics/gaussian_beams.py b/pykat/utilities/optics/gaussian_beams.py
index b97a3bbe59c8287bd337b687961d8127949c49e3..ac62b382bd3b1015a8ae407f14a39c2d00824d10 100644
--- a/pykat/utilities/optics/gaussian_beams.py
+++ b/pykat/utilities/optics/gaussian_beams.py
@@ -48,9 +48,9 @@ class gauss_param(object):
         elif len(kwargs) == 2:        
             
             if "w0" in kwargs and "z" in kwargs:
-                q = SIfloat(kwargs["z"]) + 1j *float(math.pi*SIfloat(kwargs["w0"])**2/(self.__lambda/self.__nr) )
+                q = SIfloat(kwargs["z"]) + 1j * math.pi*SIfloat(kwargs["w0"])**2/(self.__lambda/self.__nr)
             elif "z" in kwargs and "zr" in kwargs:
-                q = SIfloat(kwargs["z"]) + 1j *SIfloat(kwargs["zr"]) 
+                q = SIfloat(kwargs["z"]) + 1j * SIfloat(kwargs["zr"]) 
             elif "rc" in kwargs and "w" in kwargs:
                 one_q = 1 / SIfloat(kwargs["rc"]) - 1j * self.__lamda / (math.pi * self.__nr * SIfloat(kwargs["w"])**2)
                 q = 1/one_q
@@ -80,7 +80,7 @@ class gauss_param(object):
     
     @property
     def w(self):
-        return abs(self.__q)*math.sqrt(self.__lambda / (self.__nr * math.pi * self.__q.imag))
+        return np.abs(self.__q)* np.sqrt(self.__lambda / (self.__nr * math.pi * self.__q.imag))
     
     def beamsize(self, z=None, wavelength=None, nr=None, w0=None):
 
@@ -135,14 +135,19 @@ class gauss_param(object):
         
     @property
     def w0(self):
-        return math.sqrt(self.__q.imag * self.__lambda / (self.__nr * math.pi))    
+        return np.sqrt(self.__q.imag * self.__lambda / (self.__nr * math.pi))    
 
     @property
     def Rc(self):
-        if self.__q.real != 0:
-            return self.z * (1 + (self.zr/self.z)**2)
-        else:
-            return float("inf")
+        def __rc(z, zr):
+            if z != 0:
+                return z * (1 + (zr/z)**2)
+            else:
+                return float("inf")
+                
+        v = np.vectorize(__rc)
+        
+        return v(self.z, self.zr)
     
     def curvature(self, z=None, wavelength=None, nr=None, w0=None):
         if z == None:
@@ -314,13 +319,13 @@ class HG_beam(object):
             
     def _calc_constants(self):
         self.__xpre_const = math.pow(2.0/math.pi, 0.25)
-        self.__xpre_const *= math.sqrt(1.0/(self._qx.w0 * 2**(self._n) * math.factorial(self._n)))
-        self.__xpre_const *= cmath.sqrt(1j*self._qx.imag / self._qx.q)
+        self.__xpre_const *= np.sqrt(1.0/(self._qx.w0 * 2**(self._n) * np.factorial(self._n)))
+        self.__xpre_const *= np.sqrt(1j*self._qx.imag / self._qx.q)
         self.__xpre_const *= ((1j*self._qx.imag * self._qx.q.conjugate())/(-1j*self._qx.imag * self._qx.q)) ** ( self._n/2.0)
         
         self.__ypre_const = math.pow(2.0/math.pi, 0.25)
-        self.__ypre_const *= math.sqrt(1.0/(self._qy.w0 * 2**(self._m) * math.factorial(self._m)))
-        self.__ypre_const *= cmath.sqrt(1j*self._qy.imag / self._qy.q)
+        self.__ypre_const *= np.sqrt(1.0/(self._qy.w0 * 2**(self._m) * np.factorial(self._m)))
+        self.__ypre_const *= np.sqrt(1j*self._qy.imag / self._qy.q)
         self.__ypre_const *= ((1j*self._qy.imag * self._qy.q.conjugate())/(-1j*self._qy.imag * self._qy.q)) **(self._m/2.0)
     
         self.__sqrt2_wxz = math.sqrt(2) / self._qx.w
diff --git a/pykat/utilities/romhom.py b/pykat/utilities/romhom.py
index 350fb6a6e061e4d364a8f6f0a24295e0388a1f49..af8ad94aaa9e215cc4021ed4a923a8bf35f9454c 100644
--- a/pykat/utilities/romhom.py
+++ b/pykat/utilities/romhom.py
@@ -115,12 +115,15 @@ def u_star_u(re_q1, re_q2, w0_1, w0_2, n1, n2, x, x2=None):
     return u(re_q1, w0_1, n1, x) * u(re_q2, w0_2, n2, x2).conjugate()
 
     
-def makeReducedBasis(x, isModeMatched=True, tolerance = 1e-12, sigma = 1):
+def makeReducedBasis(x, isModeMatched=True, tolerance = 1e-12, sigma = 1, greedyfile=None):
     
-    if isModeMatched:
-        greedypts = 'matched20.txt'
+    if greedyfile != None:
+        greedypts = str(greedyfile)
     else:
-        greedypts = 'mismatched20.txt'
+        if isModeMatched:
+            greedypts = 'matched20.txt'
+        else:
+            greedypts = 'mismatched20.txt'
     
     greedyfile = os.path.join(pykat.__path__[0],'utilities','greedypoints',greedypts)