From b68817e85ef9189a883c8267666af4fc4103e80a Mon Sep 17 00:00:00 2001
From: Daniel Brown <ddb@star.sr.bham.ac.uk>
Date: Wed, 17 Sep 2014 12:48:07 +0100
Subject: [PATCH] adding romhom greedypoint file parameter, changing SIfloat to
 convert numpy arrays, beam_param can accept numpy arrays for computing
 multiple points

---
 pykat/SIfloat.py                         | 16 ++++++++++---
 pykat/utilities/maps.py                  |  6 ++---
 pykat/utilities/optics/gaussian_beams.py | 29 ++++++++++++++----------
 pykat/utilities/romhom.py                | 11 +++++----
 4 files changed, 40 insertions(+), 22 deletions(-)

diff --git a/pykat/SIfloat.py b/pykat/SIfloat.py
index 5a3556b..cac016f 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 6939b82..a9e190c 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 b97a3bb..ac62b38 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 350fb6a..af8ad94 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)
     
-- 
GitLab