diff --git a/bin/test_romhom.py b/bin/test_romhom.py
index 0b19cab2525613a5d7b80f5abd3c453e14c173e9..6275e6fbaf06d884fb808e15de2bf1ea61ea9592 100644
--- a/bin/test_romhom.py
+++ b/bin/test_romhom.py
@@ -1,4 +1,4 @@
-from pykat.utilities.maps import read_map
+from pykat.utilities.maps import read_map, tiltmap, surfacemap
from pykat.utilities.knm import *
from pykat.utilities.optics.gaussian_beams import beam_param
import time
@@ -7,20 +7,23 @@ couplings = np.array([[[0,0,0,0], [0,0,1,0], [0,0,0,1]],
[[1,0,0,0], [1,0,1,0], [1,0,0,1]],
[[0,1,0,0], [0,1,1,0], [0,1,0,1]]])
-q1 = beam_param(w0=4e-2, z=0)
-q2 = beam_param(w0=4e-2, z=0)
+q1 = beam_param(w0=3e-2, z=0)
+q2 = beam_param(w0=3e-2, z=0)
-m = read_map('etm08_virtual.txt')
+size = np.array([200,200])
+stepsize = 0.2/(size-1)
-m.data = np.zeros((100, 100))
-m.x = np.linspace(-0.15, 0.15, m.size[0])
-m.y = np.linspace(-0.15, 0.15, m.size[1])
-m.center = np.array(m.size)/2
-m.step_size = (m.x[1]-m.x[0], m.y[1]-m.y[0])
+# This map has points evenly spaced about the axes
+m = tiltmap("tilt", size, stepsize, (1e-6, 0))
+# Shifting the central point changes the results slightly
+# but this is about a 1e-7 change, probably due to extra
+# clipping
+m.center += 20
-print "generating weights..."
+print "Generating weights..."
t0 = time.time()
w, EI = m.generateROMWeights(isModeMatched=True)
+#w.writeToFile("testWeights.rom")
print "Completed in ", time.time()-t0
print "computing Riemann knm..."
@@ -39,7 +42,7 @@ print np.abs(K2)
print "Speed up", tr/tt
-#w.writeToFile("testWeights.rom")
+
diff --git a/pykat/utilities/knm.py b/pykat/utilities/knm.py
index 2303d6babf6c75be8fa3da8c3eab225ba2e44e61..d58d0180324dd4fc66eac259cbc15d8c707fa06f 100644
--- a/pykat/utilities/knm.py
+++ b/pykat/utilities/knm.py
@@ -94,47 +94,64 @@ def ROM_HG_knm(weights, mode_in, mode_out, q1, q2, q1y=None, q2y=None, cache=Non
npr = mode_in[1]
mpr = mode_out[1]
- # if cache == None:
- u_xm_nodes = u_star_u(q1.z, q2.z, q1.w0, q2.w0, n, m, weights.EI["xm"].nodes-0.1e-2)
- u_xp_nodes = u_star_u(q1.z, q2.z, q1.w0, q2.w0, n, m, weights.EI["xp"].nodes-0.1e-2)
- u_ym_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EI["ym"].nodes-0.1e-2)
- u_yp_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EI["yp"].nodes-0.1e-2)
-
- # else:
- # strx = "x[%i,%i]" % (mode_in[0], mode_out[0])
- # stry = "y[%i,%i]" % (mode_in[1], mode_out[1])
- #
- # u_x_nodes = cache[strx]
- # u_y_nodes = cache[stry]
-
- u_xy_Q1 = np.outer(u_xm_nodes, u_yp_nodes)
- u_xy_Q2 = np.outer(u_xp_nodes, u_yp_nodes)
- u_xy_Q3 = np.outer(u_xp_nodes, u_ym_nodes)
- u_xy_Q4 = np.outer(u_xm_nodes, u_ym_nodes)
-
- # n_mod_2 = n % 2
- # m_mod_2 = m % 2
- # npr_mod_2 = npr % 2
- # mpr_mod_2 = mpr % 2
- #
- # if n_mod_2 == m_mod_2 and npr_mod_2 == mpr_mod_2:
- # k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4)
- #
- # elif n_mod_2 != m_mod_2:
- # if npr_mod_2 == mpr_mod_2:
- # k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q4) + np.einsum('ij,ij', -u_xy_nodes, w_ij_Q2Q3)
- # else:
- # k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q2Q4) + np.einsum('ij,ij', -u_xy_nodes, w_ij_Q1Q3)
- #
- # elif npr_mod_2 != mpr_mod_2:
- # if n_mod_2 == m_mod_2:
- # k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q3Q4) + np.einsum('ij,ij', -u_xy_nodes, w_ij_Q1Q2)
- # else:
- # k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q2Q4) + np.einsum('ij,ij', -u_xy_nodes, w_ij_Q1Q3)
- k_ROQ = np.einsum('ij,ij', u_xy_Q1, weights.w_ij_Q1)
- k_ROQ += np.einsum('ij,ij', u_xy_Q2, weights.w_ij_Q2)
- k_ROQ += np.einsum('ij,ij', u_xy_Q3, weights.w_ij_Q3)
- k_ROQ += np.einsum('ij,ij', u_xy_Q4, weights.w_ij_Q4)
+ foundSymmetry = np.all(weights.EI["xp"].nodes == -weights.EI["xm"].nodes) and np.all(weights.EI["yp"].nodes == -weights.EI["ym"].nodes)
+
+ if foundSymmetry:
+ if cache == None:
+ u_x_nodes = u_star_u(q1.z, q2.z, q1.w0, q2.w0, n, m, weights.EI["xm"].nodes)
+ u_y_nodes = u_star_u(q1.z, q2.z, q1.w0, q2.w0, n, m, weights.EI["ym"].nodes)
+ w_ij_Q1Q3 = weights.w_ij_Q1 + weights.w_ij_Q3
+ w_ij_Q2Q4 = weights.w_ij_Q2 + weights.w_ij_Q4
+ w_ij_Q1Q2 = weights.w_ij_Q1 + weights.w_ij_Q2
+ w_ij_Q1Q4 = weights.w_ij_Q1 + weights.w_ij_Q4
+ w_ij_Q2Q3 = weights.w_ij_Q2 + weights.w_ij_Q3
+ w_ij_Q3Q4 = weights.w_ij_Q3 + weights.w_ij_Q4
+ w_ij_Q1Q2Q3Q4 = weights.w_ij_Q1 + weights.w_ij_Q2 + weights.w_ij_Q3 + weights.w_ij_Q4
+ else:
+ strx = "x[%i,%i]" % (mode_in[0], mode_out[0])
+ stry = "y[%i,%i]" % (mode_in[1], mode_out[1])
+
+ u_x_nodes = cache[strx]
+ u_y_nodes = cache[stry]
+
+ u_xy_nodes = np.outer(u_x_nodes, u_y_nodes)
+
+ n_mod_2 = n % 2
+ m_mod_2 = m % 2
+ npr_mod_2 = npr % 2
+ mpr_mod_2 = mpr % 2
+
+ if n_mod_2 == m_mod_2 and npr_mod_2 == mpr_mod_2:
+ k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4)
+
+ elif n_mod_2 != m_mod_2:
+ if npr_mod_2 == mpr_mod_2:
+ k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q4) + np.einsum('ij,ij', -u_xy_nodes, w_ij_Q2Q3)
+ else:
+ k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q2Q4) + np.einsum('ij,ij', -u_xy_nodes, w_ij_Q1Q3)
+
+ elif npr_mod_2 != mpr_mod_2:
+ if n_mod_2 == m_mod_2:
+ k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q3Q4) + np.einsum('ij,ij', -u_xy_nodes, w_ij_Q1Q2)
+ else:
+ k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q2Q4) + np.einsum('ij,ij', -u_xy_nodes, w_ij_Q1Q3)
+
+ else:
+ # If there is no symmetry about the axes then each quadrant is different
+ u_xm_nodes = u_star_u(q1.z, q2.z, q1.w0, q2.w0, n, m, weights.EI["xm"].nodes)
+ u_xp_nodes = u_star_u(q1.z, q2.z, q1.w0, q2.w0, n, m, weights.EI["xp"].nodes)
+ u_ym_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EI["ym"].nodes)
+ u_yp_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EI["yp"].nodes)
+
+ u_xy_Q1 = np.outer(u_xm_nodes, u_yp_nodes)
+ u_xy_Q2 = np.outer(u_xp_nodes, u_yp_nodes)
+ u_xy_Q3 = np.outer(u_xp_nodes, u_ym_nodes)
+ u_xy_Q4 = np.outer(u_xm_nodes, u_ym_nodes)
+
+ k_ROQ = np.einsum('ij,ij', u_xy_Q1, weights.w_ij_Q1)
+ k_ROQ += np.einsum('ij,ij', u_xy_Q2, weights.w_ij_Q2)
+ k_ROQ += np.einsum('ij,ij', u_xy_Q3, weights.w_ij_Q3)
+ k_ROQ += np.einsum('ij,ij', u_xy_Q4, weights.w_ij_Q4)
return k_ROQ
diff --git a/pykat/utilities/maps.py b/pykat/utilities/maps.py
index d19b032788b842971bf1d91e4979c988d13ca6cd..eb663bbc72c1d61d41fc25f9ee95917643c8ae48 100644
--- a/pykat/utilities/maps.py
+++ b/pykat/utilities/maps.py
@@ -1,24 +1,9 @@
from pykat.utilities.romhom import makeReducedBasis, makeEmpiricalInterpolant, makeWeights
-import numpy
+import numpy as np
import math
-
-class tiltmap(surfacemap):
-
- def __init__(self, name, size, step_size, tilt, center=(0,0) ):
- surfacemap.__init__(name, "phase", size, center, step_sizes)
- self.tilt = tilt
-
- @property
- def tilt(self):
- return self.__tilt
-
- @tilt.setter
- def tilt(self, value):
-
- self.__tilt = value
-class surfacemap:
+class surfacemap(object):
def __init__(self, name, maptype, size, center, step_size, scaling, data=None):
self.name = name
@@ -26,7 +11,12 @@ class surfacemap:
self.center = center
self.step_size = step_size
self.scaling = scaling
- self.data = data
+
+ if data == None:
+ self.data = np.zeros(size)
+ else:
+ self.data = data
+
self._rom_weights = None
def write_map(self, filename):
@@ -48,11 +38,11 @@ class surfacemap:
@property
def x(self):
- return self.step_size[0] * (numpy.array(range(1, self.data.shape[0]+1)) - self.center[0])
+ return self.step_size[0] * (np.array(range(1, self.data.shape[0]+1)) - self.center[0])
@property
def y(self):
- return self.step_size[1] * (numpy.array(range(1, self.data.shape[1]+1))- self.center[1])
+ return self.step_size[1] * (np.array(range(1, self.data.shape[1]+1))- self.center[1])
@property
def size(self):
@@ -60,7 +50,7 @@ class surfacemap:
@property
def offset(self):
- return numpy.array(self.step_size)*(self.center - numpy.array(self.size)/2)
+ return np.array(self.step_size)*(self.center - np.array(self.size)/2)
@property
def ROMWeights(self):
@@ -70,7 +60,7 @@ class surfacemap:
if "phase" in self.type:
k = math.pi * 2 / wavelength
- return numpy.exp(2j * k * self.scaling * self.data)
+ return np.exp(2j * k * self.scaling * self.data)
else:
raise BasePyKatException("Map type needs handling")
@@ -87,7 +77,7 @@ class surfacemap:
if len(xp) > 0: EI["xp"] = makeEmpiricalInterpolant(makeReducedBasis(xp, isModeMatched=isModeMatched))
if len(ym) > 0: EI["ym"] = makeEmpiricalInterpolant(makeReducedBasis(ym, isModeMatched=isModeMatched))
if len(yp) > 0: EI["yp"] = makeEmpiricalInterpolant(makeReducedBasis(yp, isModeMatched=isModeMatched))
-
+
#if len(x0) > 0: EI["x0"] = makeEmpiricalInterpolant(makeReducedBasis(x0, isModeMatched=isModeMatched))
#if len(y0) > 0: EI["y0"] = makeEmpiricalInterpolant(makeReducedBasis(y0, isModeMatched=isModeMatched))
@@ -122,8 +112,27 @@ class surfacemap:
return fig
+
+class tiltmap(surfacemap):
+
+ def __init__(self, name, size, step_size, tilt):
+ surfacemap.__init__(self, name, "phase", size, (np.array(size)+1)/2.0, step_size, 1)
+ self.tilt = tilt
+
+ @property
+ def tilt(self):
+ return self.__tilt
+
+ @tilt.setter
+ def tilt(self, value):
+ self.__tilt = value
+
+ xx, yy = np.meshgrid(self.x, self.y)
+
+ self.data = xx * self.tilt[0] + yy * self.tilt[1]
+
+
-
def read_map(filename):
with open(filename, 'r') as f:
@@ -137,7 +146,7 @@ def read_map(filename):
- data = numpy.loadtxt(filename, dtype=numpy.float64,ndmin=2,comments='%')
+ data = np.loadtxt(filename, dtype=np.float64,ndmin=2,comments='%')
return surfacemap(name,maptype,size,center,step,scaling,data)