diff --git a/bin/test_square_knm.py b/bin/test_square_knm.py
new file mode 100644
index 0000000000000000000000000000000000000000..c1950c9807926489c0e05d185961065033c0baaf
--- /dev/null
+++ b/bin/test_square_knm.py
@@ -0,0 +1,20 @@
+import pykat
+import numpy as np
+from pykat.optics.knm import square_aperture_HG_knm, riemann_HG_knm
+import math
+
+R = 0.15
+q = pykat.beam_param(w0=6e-3, z=2000)
+N = 1201
+
+mode_i = (0,0)
+mode_o = (0,0)
+
+k = square_aperture_HG_knm(mode_i, mode_o, q, R)
+
+x = y = np.linspace(-R, R, N)
+Axy = np.ones((N,N))
+k_ = riemann_HG_knm(x, y, mode_i, mode_o, q, q, Axy=Axy, newtonCotesOrder=2)
+
+print "%15.15f + %15.15fi" % (k.real, k.imag)
+print "%15.15f + %15.15fi" % (k_.real, k_.imag), abs(k-k_)/1e-6, "ppm"
\ No newline at end of file
diff --git a/pykat/optics/gaussian_beams.py b/pykat/optics/gaussian_beams.py
index b7e524fee880b38e4ce03ebb0e124fadee9a5adb..6d8c0ed4f2879343ea251c0b929d80493b652f41 100644
--- a/pykat/optics/gaussian_beams.py
+++ b/pykat/optics/gaussian_beams.py
@@ -318,7 +318,15 @@ class HG_beam(object):
             self._qy = copy.deepcopy(value)
         else:
             self._qy = beam_param(q=complex(value))
-            
+    
+    @property
+    def constant_x(self):
+        return self.__xpre_const
+        
+    @property
+    def constant_y(self):
+        return self.__ypre_const
+        
     def _calc_constants(self):
         self.__xpre_const = math.pow(2.0/math.pi, 0.25)
         self.__xpre_const *= np.sqrt(1.0/(self._qx.w0 * 2**(self._n) * np.math.factorial(self._n)))
@@ -338,7 +346,7 @@ class HG_beam(object):
         
         self.__invqx = 1/ self._qx.q
         self.__invqy = 1/ self._qy.q
-        
+    
     def Un(self, x):
         return self.__xpre_const * self._hn(self.__sqrt2_wxz * x) * np.exp(-0.5j * self.__kx * x*x * self.__invqx)
     
diff --git a/pykat/optics/knm.py b/pykat/optics/knm.py
index 69435a818e210b90265c0394217e6da92db132f6..881efc4fad864c47603083e44b65e534929c8320 100644
--- a/pykat/optics/knm.py
+++ b/pykat/optics/knm.py
@@ -6,6 +6,10 @@ from pykat.external.progressbar import ProgressBar, ETA, Percentage, Bar
 from scipy.interpolate import interp2d
 from scipy.integrate import dblquad
 from pykat.optics.romhom import ROMWeights
+from math import factorial
+from pykat.maths.hermite import hermite
+from scipy.misc import comb
+from scipy.integrate import newton_cotes
 
 import time
 import pykat.optics.maps
@@ -144,15 +148,16 @@ def adaptive_knm(mode_in, mode_out, q1, q2, q1y=None, q2y=None, smap=None, delta
     
     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={}):
+def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None,
+                     Axy=None, cache=None, delta=(0,0), params={}, newtonCotesOrder=0):
 
-    if Axy == None:
+    if Axy is None:
         Axy == np.ones((len(x), len(y)))
     
-    if q1y == None:
+    if q1y is None:
         q1y = q1
         
-    if q2y == None:
+    if q2y is None:
         q2y = q2
         
     if len(mode_in) != 2 or len(mode_out) != 2:
@@ -164,11 +169,35 @@ def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None
     if cache == None:
         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])
-    
+        
         U1 = Hg_in.Unm(x+delta[0], y+delta[1])
         U2 = Hg_out.Unm(x,y).conjugate()
-        
-        return dx * dy * np.einsum('ij,ij', Axy, U1*U2)
+
+        if newtonCotesOrder > 0:
+                
+            W = newton_cotes(newtonCotesOrder, 1)[0]
+            
+            if newtonCotesOrder > 1:
+                if (len(x) - len(W)) % newtonCotesOrder != 0:
+                    raise ValueError("To use Newton-Cotes order {0} the number of data points in x must ensure: (N_x - ({0}+1)) mod {0} == 0".format(newtonCotesOrder) )
+
+                if (len(y) - len(W)) % newtonCotesOrder != 0:
+                    raise ValueError("To use Newton-Cotes order {0} the number of data points in y must ensure: (N_y - ({0}+1)) mod {0} == 0".format(newtonCotesOrder) )
+                
+            wx = np.zeros(x.shape, dtype=np.float64)    
+            wy = np.zeros(y.shape, dtype=np.float64)
+    
+            N = len(W)
+
+            for i in range(0, (len(wx)-1)/newtonCotesOrder): wx[(i*(N-1)):(i*(N-1)+N)] += W
+            for i in range(0, (len(wy)-1)/newtonCotesOrder): wy[(i*(N-1)):(i*(N-1)+N)] += W
+            
+            Wxy = np.outer(wx, wy)
+            
+        if newtonCotesOrder == 0:
+            return dx * dy * np.einsum('ij,ij', Axy, U1*U2)
+        else:
+            return dx * dy * np.einsum('ij,ij', Axy, U1*U2*Wxy)
     else:
         
         strx = "u1[%i,%i]" % (mode_in[0], mode_out[0])
@@ -438,6 +467,63 @@ 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 __sqaure_knm_int(n, _n, R):
+    # This uses the H_n(x) * H_m(x) product identity to reduce the overlap into
+    # a sum of factorial and an integral of a single Hermite with a gaussian function
+    # thus making it easier to solve
+    expR = math.exp(-(R**2))
+    S = 0
+    
+    for j in range(0, min(n, _n)+1):
+        _j1 = _n + n - 2*j - 1
+        
+        if _j1+1 == 0:
+            # for the zeroth order we just have the gaussian integral to solve
+            L = math.sqrt(math.pi) * math.erf(R)    
+        elif (_j1+1) % 2 == 1:
+            # if the Hermite is odd then the integral is always 0 as its anti-symmetric
+            L = 0
+        else:
+            L = 2 * hermite(_j1, 0) - expR * (hermite(_j1, R) - hermite(_j1, -R))
+        
+        I = 2**j * factorial(j) * comb(n, j) * comb(_n, j)
+        
+        S += I * L
+                
+    return S
+
+
+def square_aperture_HG_knm(mode_in, mode_out, q, R):
+    """
+    Computes the coupling coefficients for a square aperture.
+    """
+    # x modes
+    n = mode_in[0]
+    _n = mode_out[0]
+
+    # y modes
+    m = mode_in[1]
+    _m = mode_out[1]
+    
+    hg1 = HG_beam(q, n=n, m=m)
+    hg2 = HG_beam(q, n=_n, m=_m)
+        
+    kx = hg1.constant_x * hg2.constant_x.conjugate()
+    ky = hg1.constant_y * hg2.constant_y.conjugate()
+    
+    print hg1.constant_x, hg2.constant_x, hg1.constant_y, hg2.constant_y
+    
+    f = q.w / math.sqrt(2)
+    R = R / (q.w / math.sqrt(2))
+    
+    kx *= f
+    kx *= __sqaure_knm_int(n, _n, R)
+    
+    ky *= f
+    ky *= __sqaure_knm_int(m, _m, R)
+    
+    return kx * ky
+
 
 
 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={}):
@@ -547,6 +633,9 @@ def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riema
         return K.reshape(couplings.shape[:-1])
 
 
+
+
+
 def plot_knm_matrix(couplings, knm):
     import pylab as plt
     
diff --git a/pykat/optics/romhom.py b/pykat/optics/romhom.py
index b0b75f18f156300ae55aa661aa734bedf83000e7..1938d249f8e088c5de92cbd9e6933b9f259b418e 100644
--- a/pykat/optics/romhom.py
+++ b/pykat/optics/romhom.py
@@ -177,7 +177,11 @@ def makeReducedBasis(x, isModeMatched=True, tolerance = 1e-12, sigma = 1, greedy
     
     greedyfile = os.path.join(pykat.__path__[0],'optics','greedypoints',greedypts)
     
-    limits = np.loadtxt(greedyfile, usecols=(1,))[:5]
+    # Two potential different formats
+    try:
+        limits = np.loadtxt(greedyfile, usecols=(3,))[:5]
+    except IndexError:
+        limits = np.loadtxt(greedyfile, usecols=(1,))[:5]
     
     romlimits = ROMLimits(w0min=limits[0], w0max=limits[1], zmin=limits[2], zmax=limits[3], max_order=limits[4])