diff --git a/.gitignore b/.gitignore index 94487b9516e57164cde23bec81c7aef250fd2a8f..20c5256ec7e862e466a66e2db1adacd9a82b0efb 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,9 @@ *.pyc +.DS_store +*.dat +*.txt +*.log +*.kat +*.pdf +*.out +*.bak \ No newline at end of file diff --git a/bin/test_knm.py b/bin/test_knm.py new file mode 100644 index 0000000000000000000000000000000000000000..07d1d7c64a0beb822eafc7a62efda9404a2dc40a --- /dev/null +++ b/bin/test_knm.py @@ -0,0 +1,39 @@ +from pykat import * +from pykat.utilities.knm import knmHG, makeCouplingMatrix, plot_knm_matrix +from pykat.utilities.maps import aperturemap +import numpy as np +import time + +q1 = beam_param(w0=5e-2, z=0) +q2 = beam_param(w0=5e-2, z=0) + +aperture = 0.1 +s = 1000 +size = np.array([s, s]) +stepsize = 0.3/(size-1) +smap = aperturemap("tilt", size, stepsize, aperture) + +couplings = makeCouplingMatrix(1) + +params = {"usepolar":True, "aperture":aperture, "epsabs": 1e-3, "epsrel": 1e-3} + +t0 = time.time() +kbh = knmHG(couplings, q1, q2, method="adaptive", verbose=True, params=params) +print time.time() - t0 + +t0 = time.time() +kr = knmHG(couplings, q1, q2, surface_map=smap, method="riemann", verbose=True) +print time.time() - t0 + +smap.generateROMWeights(isModeMatched=True, verbose=True) +t0 = time.time() +krm = knmHG(couplings, q1, q2, surface_map=smap, method="romhom") +print time.time() - t0 + +print kbh +print kr +print krm + +plot_knm_matrix(couplings, np.log10(np.abs(kbh - kr))) + +print params \ No newline at end of file diff --git a/pykat/utilities/knm.py b/pykat/utilities/knm.py index 1607bc26d990bff8813f08f01536d50637f48bbd..6374fd16f3010cd02ca82098ef195f746568bbf2 100644 --- a/pykat/utilities/knm.py +++ b/pykat/utilities/knm.py @@ -1,6 +1,10 @@ from itertools import combinations_with_replacement as combinations from pykat.utilities.optics.gaussian_beams import beam_param, HG_beam from pykat.exceptions import BasePyKatException +from romhom import u_star_u +from progressbar import ProgressBar, ETA, Percentage, Bar +from scipy.interpolate import interp2d +from scipy.integrate import dblquad import time import maps @@ -11,10 +15,7 @@ import collections import math import cmath -from romhom import u_star_u -from progressbar import ProgressBar, ETA, Percentage, Bar - -def makeCouplingMatrix(max_order): +def makeCouplingMatrix(max_order, Neven=True, Nodd=True, Meven=True, Modd=True): max_order = int(max_order) c = [] for n in range(0, max_order+1): @@ -30,13 +31,119 @@ def makeCouplingMatrix(max_order): for j in c: e = list(i) e.extend(j) + + if not Neven and (e[0]-e[2]) % 2 == 0: continue + if not Nodd and (e[0]-e[2]) % 2 == 1: continue + if not Meven and (e[1]-e[3]) % 2 == 0: continue + if not Modd and (e[1]-e[3]) % 2 == 1: continue + row.append(e) - M.append(row) + + M.append(np.array(row).squeeze()) return np.array(M) -def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None, cache=None, delta=(0,0)): +def adaptive_knm(mode_in, mode_out, q1, q2, q1y=None, q2y=None, smap=None, delta=(0,0), params={}): + + if q1y == None: + q1y = q1 + + if q2y == None: + q2y = q2 + + if "epsabs" not in params: params["epsabs"] = 1e-6 + if "epsrel" not in params: params["epsrel"] = 1e-6 + if "usepolar" not in params: params["usepolar"] = False + + if len(mode_in) != 2 or len(mode_out) != 2: + raise BasePyKatException("Both mode in and out should be a container with modes [n m]") + + Hg_in = HG_beam(qx=q1, qy=q1y, n=mode_in[0], m=mode_in[1]) + Hg_out = HG_beam(qx=q2, qy=q2y, n=mode_out[0], m=mode_out[1]) + + Nfuncs = [] + Nfuncs.append(0) + + if smap != None: + + if not params["usepolar"]: + xlims = (min(smap.x), max(smap.x)) + ylims = (min(smap.y), max(smap.y)) + + def Rfunc(y,x): + Nfuncs[-1] += len(x) + return (Hg_in.Unm(x+delta[0], y+delta[1]) * smap.z_xy(x=x,y=y) * Hg_out.Unm(x, y).conjugate()).real + + def Ifunc(y,x): + Nfuncs[-1] += len(x) + return (Hg_in.Unm(x+delta[0], y+delta[1]) * smap.z_xy(x=x,y=y) * Hg_out.Unm(x, y).conjugate()).imag + + else: + xlims = (0, 2*math.pi) + ylims = (0, params["aperture"]) + + def Rfunc(r, phi): + Nfuncs[-1] += len(x) + x = r*np.cos(phi) + y = r*np.sin(phi) + return (r * Hg_in.Unm(x, y) * smap.z_xy(x=x,y=y) * Hg_out.Unm(x, y).conjugate()).real + + def Ifunc(r, phi): + Nfuncs[-1] += len(x) + x = r*np.cos(phi) + y = r*np.sin(phi) + return (r * Hg_in.Unm(x, y) * smap.z_xy(x=x,y=y) * Hg_out.Unm(x, y).conjugate()).imag + + else: + if not params["usepolar"]: + _x = 4 * math.sqrt(1+max(mode_in[0],mode_in[1])) * q1.w + _y = 4 * math.sqrt(1+max(mode_in[0],mode_in[1])) * q1y.w + + xlims = (-_x, _x) + ylims = (-_y, _y) + + def Rfunc(y, x): + Nfuncs[-1] += len(r) + return (Hg_in.Unm(x+delta[0], y+delta[1]) * Hg_out.Unm(x, y).conjugate()).real + + def Ifunc(y,x): + Nfuncs[-1] += len(r) + return (Hg_in.Unm(x+delta[0], y+delta[1]) * Hg_out.Unm(x, y).conjugate()).imag + else: + xlims = (0, 2*math.pi) + ylims = (0, params["aperture"]) + + def Rfunc(r, phi): + + if hasattr(r, "__len__"): + Nfuncs[-1] += len(r) + else: + Nfuncs[-1] += 1 + + x = r*np.cos(phi) + y = r*np.sin(phi) + return (r * Hg_in.Unm(x, y) * Hg_out.Unm(x, y).conjugate()).real + + def Ifunc(r, phi): + if hasattr(r, "__len__"): + Nfuncs[-1] += len(r) + else: + Nfuncs[-1] += 1 + + x = r*np.cos(phi) + y = r*np.sin(phi) + return (r * Hg_in.Unm(x, y) * Hg_out.Unm(x, y).conjugate()).imag + + R, errR = dblquad(Rfunc, xlims[0], xlims[1], lambda y: ylims[0], lambda y: ylims[1], epsabs=params["epsabs"], epsrel=params["epsrel"]) + I, errI = dblquad(Ifunc, xlims[0], xlims[1], lambda y: ylims[0], lambda y: ylims[1], epsabs=params["epsabs"], epsrel=params["epsrel"]) + + params["Nfuncs"] = Nfuncs[0] + params["errors"] = (errR, errI) + + return R + 1j * I + +def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None, cache=None, delta=(0,0), params={}): if Axy == None: Axy == np.ones((len(x), len(y))) @@ -71,7 +178,7 @@ def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None -def __gen_riemann_knm_cache(x, y, couplings, q1, q2, q1y=None, q2y=None, delta=(0,0)): +def __gen_riemann_knm_cache(x, y, couplings, q1, q2, q1y=None, q2y=None, delta=(0,0), params={}): if q1y == None: q1y = q1 @@ -311,7 +418,9 @@ def bayerhelms_HG_knm(mode_in, mode_out, q1, q2, q1y=None, q2y=None, gamma=(0,0) return __bayerhelms_kn(n,_n, q1, q2, 2*gamma[0]) * __bayerhelms_kn(m, _m, q1y, q2y, 2*gamma[1]) -def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riemann", verbose=False, profile=False, gamma=(0,0), delta=(0,0)): + + +def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riemann", verbose=False, profile=False, gamma=(0,0), delta=(0,0), params={}): if q1y == None: q1y = q1 @@ -339,7 +448,7 @@ def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riema __fac_cache.append(math.factorial(n)) if surface_map != None: - Axy = surface_map.z_xy(q1.wavelength) + Axy = surface_map.z_xy(wavelength=q1.wavelength) x = surface_map.x y = surface_map.y @@ -388,12 +497,15 @@ def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riema mode_in = [int(it.next()), int(it.next())] mode_out = [int(it.next()), int(it.next())] + if method == "riemann": K[i] = riemann_HG_knm(x, y, mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, Axy=Axy, cache=cache, delta=delta) elif method == "romhom": K[i] = ROM_HG_knm(weights, mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, cache=cache) elif method == "bayerhelms": K[i] = bayerhelms_HG_knm(mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, gamma=gamma) + elif method == "adaptive": + K[i] = adaptive_knm(mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, smap=surface_map, delta=delta, params=params) else: raise BasePyKatException("method value '%s' not accepted" % method) diff --git a/pykat/utilities/maps.py b/pykat/utilities/maps.py index 397cfe0772b5494a9427bd0fadaef7b7a207ba0e..6939b82c587d9faea1838bf51a324350632b47a5 100644 --- a/pykat/utilities/maps.py +++ b/pykat/utilities/maps.py @@ -1,4 +1,5 @@ from pykat.utilities.romhom import makeReducedBasis, makeEmpiricalInterpolant, makeWeights +from scipy.interpolate import interp2d import numpy as np import math @@ -11,7 +12,8 @@ class surfacemap(object): self.center = center self.step_size = step_size self.scaling = scaling - + self.__interp = None + if data == None: self.data = np.zeros(size) else: @@ -36,6 +38,42 @@ class surfacemap(object): mapfile.write("%.15g " % self.data[i,j]) mapfile.write("\n") + @property + def data(self): + return self.__data + + @data.setter + def data(self, value): + self.__data = value + self.__interp = None + + @property + def center(self): + return self.__center + + @center.setter + def center(self, value): + self.__center = value + self.__interp = None + + @property + def step_size(self): + return self.__step_size + + @step_size.setter + def step_size(self, value): + self.__step_size = value + self.__interp = None + + @property + def scaling(self): + return self.__scaling + + @scaling.setter + def scaling(self, value): + self.__scaling = value + self.__interp = None + @property def x(self): return self.step_size[0] * (np.array(range(1, self.data.shape[0]+1)) - self.center[0]) @@ -56,22 +94,34 @@ class surfacemap(object): def ROMWeights(self): return self._rom_weights - def z_xy(self, wavelength=1064e-9, direction="reflection", nr1=1.0, nr2=1.0): + def z_xy(self, x=None, y=None, wavelength=1064e-9, direction="reflection", nr1=1.0, nr2=1.0): + if x == None and y == None: + data = self.scaling * self.data + else: + if self.__interp == None: + self.__interp = interp2d(self.x, self.y, self.data * self.scaling) + + data = self.__interp(x, y) + if direction == "reflection": if "phase" in self.type: k = math.pi * 2 / wavelength - return np.exp(2j * k * self.scaling * self.data) + return np.exp(2j * k * data) + elif "absorption" in self.type: - return np.sqrt(1.0 - self.scaling * self.data) + return np.sqrt(1.0 - data) + else: raise BasePyKatException("Map type needs handling") elif direction == "transmission": if "phase" in self.type: k = math.pi * 2 / wavelength - return np.exp((nr1-nr2)*k * self.scaling * self.data) + return np.exp((nr1-nr2)*k * data) + elif "absorption" in self.type: - return np.sqrt(1.0 - self.scaling * self.data) + return np.sqrt(1.0 - data) + else: raise BasePyKatException("Map type needs handling") else: diff --git a/pykat/utilities/optics/gaussian_beams.py b/pykat/utilities/optics/gaussian_beams.py index 0ff65548a2db0b5be2eff7b40c2575b895e3ebf6..b97a3bbe59c8287bd337b687961d8127949c49e3 100644 --- a/pykat/utilities/optics/gaussian_beams.py +++ b/pykat/utilities/optics/gaussian_beams.py @@ -314,12 +314,12 @@ 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.w*2**self._n * math.factorial(self._n))) + 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 *= ((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.w*2**self._m * math.factorial(self._m))) + 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 *= ((1j*self._qy.imag * self._qy.q.conjugate())/(-1j*self._qy.imag * self._qy.q)) **(self._m/2.0)