Skip to content
Snippets Groups Projects
Commit b68817e8 authored by Daniel Brown's avatar Daniel Brown
Browse files

adding romhom greedypoint file parameter, changing SIfloat to convert numpy...

adding romhom greedypoint file parameter, changing SIfloat to convert numpy arrays, beam_param can accept numpy arrays for computing multiple points
parent ac57bb98
Branches
No related tags found
No related merge requests found
import os import os
import re import re
import pykat.exceptions as pkex import pykat.exceptions as pkex
import numpy as np
__suffix = {'y': 'e-24', # yocto __suffix = {'y': 'e-24', # yocto
'z': 'e-21', # zepto 'z': 'e-21', # zepto
...@@ -23,10 +24,19 @@ def SIfloat(value): ...@@ -23,10 +24,19 @@ def SIfloat(value):
if value==None: if value==None:
return value return value
if type(value)==list: value = np.array(value)
return [convertToFloat(s) for s in value]
v = np.vectorize(convertToFloat)
if value.size == 1:
return float(v(value))
a = v(value)
if len(a) == 1:
return a[0]
else: else:
return convertToFloat(value) return a
def convertToFloat(value): def convertToFloat(value):
......
...@@ -128,7 +128,7 @@ class surfacemap(object): ...@@ -128,7 +128,7 @@ class surfacemap(object):
raise BasePyKatException("Map type needs handling") 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: if interpolate:
from scipy.interpolate import interp2d from scipy.interpolate import interp2d
...@@ -169,12 +169,12 @@ class surfacemap(object): ...@@ -169,12 +169,12 @@ class surfacemap(object):
EI = {} 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: if symm:
EI["ym"] = EI["xm"] EI["ym"] = EI["xm"]
else: 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 EI["limits"] = EI["xm"].limits
......
...@@ -48,7 +48,7 @@ class gauss_param(object): ...@@ -48,7 +48,7 @@ class gauss_param(object):
elif len(kwargs) == 2: elif len(kwargs) == 2:
if "w0" in kwargs and "z" in kwargs: 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: 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: elif "rc" in kwargs and "w" in kwargs:
...@@ -80,7 +80,7 @@ class gauss_param(object): ...@@ -80,7 +80,7 @@ class gauss_param(object):
@property @property
def w(self): 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): def beamsize(self, z=None, wavelength=None, nr=None, w0=None):
...@@ -135,15 +135,20 @@ class gauss_param(object): ...@@ -135,15 +135,20 @@ class gauss_param(object):
@property @property
def w0(self): 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 @property
def Rc(self): def Rc(self):
if self.__q.real != 0: def __rc(z, zr):
return self.z * (1 + (self.zr/self.z)**2) if z != 0:
return z * (1 + (zr/z)**2)
else: else:
return float("inf") return float("inf")
v = np.vectorize(__rc)
return v(self.z, self.zr)
def curvature(self, z=None, wavelength=None, nr=None, w0=None): def curvature(self, z=None, wavelength=None, nr=None, w0=None):
if z == None: if z == None:
z = self.z z = self.z
...@@ -314,13 +319,13 @@ class HG_beam(object): ...@@ -314,13 +319,13 @@ class HG_beam(object):
def _calc_constants(self): def _calc_constants(self):
self.__xpre_const = math.pow(2.0/math.pi, 0.25) 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 *= np.sqrt(1.0/(self._qx.w0 * 2**(self._n) * np.factorial(self._n)))
self.__xpre_const *= cmath.sqrt(1j*self._qx.imag / self._qx.q) 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.__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.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 *= np.sqrt(1.0/(self._qy.w0 * 2**(self._m) * np.factorial(self._m)))
self.__ypre_const *= cmath.sqrt(1j*self._qy.imag / self._qy.q) 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.__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 self.__sqrt2_wxz = math.sqrt(2) / self._qx.w
......
...@@ -115,8 +115,11 @@ def u_star_u(re_q1, re_q2, w0_1, w0_2, n1, n2, x, x2=None): ...@@ -115,8 +115,11 @@ 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() 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 greedyfile != None:
greedypts = str(greedyfile)
else:
if isModeMatched: if isModeMatched:
greedypts = 'matched20.txt' greedypts = 'matched20.txt'
else: else:
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment