adding in new integrator for knm using higher order newton-cotes methods and...

adding in new integrator for knm using higher order newton-cotes methods and analytic solution for square aperture

bin/test_square_knm.py

+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

pykat/optics/gaussian_beams.py

@@ -319,6 +319,14 @@ class HG_beam(object):
         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)))

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:
@@ -168,8 +173,32 @@ def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None
         U1 = Hg_in.Unm(x+delta[0], y+delta[1])
         U2 = Hg_out.Unm(x,y).conjugate()
 
+        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])
         stry = "u2[%i,%i]" % (mode_in[1], mode_out[1])
@@ -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
     

pykat/optics/romhom.py

@@ -177,6 +177,10 @@ def makeReducedBasis(x, isModeMatched=True, tolerance = 1e-12, sigma = 1, greedy
     
     greedyfile = os.path.join(pykat.__path__[0],'optics','greedypoints',greedypts)
     
+    # 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])