Merge branch 'master' of gitlab.aei.uni-hannover.de:finesse/pykat

# Please enter a commit message to explain why this merge is necessary,
# especially if it merges an updated upstream into a topic branch.
#
# Lines starting with '#' will be ignored, and an empty message aborts
# the commit.

pykat/optics/gaussian_beams.py

+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
 import pykat.exceptions as pkex
 import numpy as np
 import math
@@ -7,59 +12,6 @@ import cmath
 from scipy.special import hermite
 from pykat.SIfloat import SIfloat
 
-def fit_beam_size(z_data, w_data, w0guess, z0guess, lam=1064e-9, show_plot=True,
-                  plotpts=100, title='Beam scan fit', w2zero=True, filename=None):
-    """
-    Function to fit a beam size (w0 and z) to a 2D intensity array. Plotting can
-    be switched on or off.
-
-    Returns best fit of (w0, z)
-    
-    Contributed by Paul Fulda on 15/01/16
-    """
-    
-    import scipy.optimize
-    
-    def zw0z02w(z, w0, z0, lam):
-        z_R = pl.pi*w0**2/lam
-        w = w0*pl.sqrt(1+((z-z0)/z_R)**2)
-        return w
-    
-    popt,pcov = scipy.optimize.curve_fit(lambda z_data, w0, z0: zw0z02w(z_data, w0, z0, lam), z_data, w_data, p0=[w0guess,z0guess])
-    w0out=popt[0]
-    z0out=popt[1]
-    
-    z_fit = pl.linspace(min(z_data),max(z_data),plotpts)
-    w_fit = zw0z02w(z_fit, w0out, z0out, lam)
-
-
-    if show_plot or filename is not None:
-        import pylab as pl
-        
-        um=1e6
-        
-        pl.figure()
-        pl.plot(z_data,w_data*um,'bo', label = 'Data')
-        pl.plot(z_fit,w_fit*um,'b', label = 'Fit')
-        pl.tight_layout
-        pl.grid()
-        pl.xlabel('Position [m]')
-        pl.ylabel('Beam size [$\mu$m]')
-        pl.xlim((min(z_data),max(z_data)))
-        
-        if w2zero:
-            pl.ylim((0,max(w_data)*um))
-            
-        pl.legend(loc=0)
-        pl.title(title)
-        
-        if filename is not None:
-            pl.savefig(filename)
-            
-    return w0out, z0out      
-    
-    
-
 class gauss_param(object):
     """
     Use beam_param instead, will be the future name of this object.
@@ -317,8 +269,14 @@ class gauss_param(object):
 class beam_param(gauss_param):
     pass
 
+# Should be renamed to HG_mode?    
 class HG_beam(object):
-    
+    """ Hermite-Gauss beam profile. Example usage:
+    import pykat.optics.gaussian_beams as gb
+    qx=gb.beam_param(w0=1e-3,z=0)
+    beam=gb.HG_beam(qx,n=2,m=0)
+    beam.plot()
+    """    
     def __init__(self, qx, qy=None, n=0, m=0):
         self._qx = copy.deepcopy(qx)
         self._2pi_qrt = math.pow(2.0/math.pi, 0.25)
@@ -427,7 +385,9 @@ class HG_beam(object):
         return np.outer(_un, _um)
         
     def plot(self, ndx=100, ndy=100, xscale=4, yscale=4):
-        import pylab
+        """ Make a simple plot the HG_beam """
+        import pykat.plotting 
+        import matplotlib.pyplot as plt
         
         xrange = xscale * np.linspace(-self._qx.w, self._qx.w, ndx)
         yrange = yscale * np.linspace(-self._qy.w, self._qy.w, ndy)
@@ -437,11 +397,11 @@ class HG_beam(object):
 
         data = self.Unm(xrange,yrange)
 
-        fig = pylab.figure()
-        axes = pylab.imshow(np.abs(data), aspect=dx/dy, extent=[min(xrange),max(xrange),min(yrange),max(yrange)])
-        pylab.xlabel('x [m]')
-        pylab.ylabel('y [m]')
+        fig = pykat.plotting.figure()
+        axes = plt.imshow(np.abs(data.T), aspect=dx/dy, extent=[min(xrange),max(xrange),min(yrange),max(yrange)])
+        plt.xlabel('x [m]')
+        plt.ylabel('y [m]')
         cbar = fig.colorbar(axes)
-        pylab.show()
+        plt.show()
         
         

pykat/optics/pdtype.py

+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+import pykat.exceptions as pkex
+import numpy as np
+import math
+import copy
+import warnings
+import cmath
+from math import factorial
+from scipy.integrate import quad
+from scipy.special import hermite
+from pykat.SIfloat import SIfloat
+
+
+def finesse_split_photodiode(maxtem, xy='x'):
+    """Prints beat coefficients for HG modes on a split photo detector
+    in the format for insertion into the kat.ini file of Finesse.
+    Example use:
+    finesse_split_photodiode(5, xy='x')
+    The default kat.ini numbers have been generated with:
+    finesse_split_photodiode(40, xy='x')
+    finesse_split_photodiode(40, xy='y')
+    """
+    assert (xy=="x" or xy=="y"), "xy function parameter must be 'x' or 'y'" 
+    assert (isinstance(maxtem, int) and maxtem >=0), "maxtem must be integer >=0" 
+
+    if (xy=="x"):
+        print('PDTYPE x-split')
+    else:
+        print('PDTYPE y-split')
+
+    # running through even modes
+    for n1 in np.arange(0,maxtem+1,2):
+        # running through odd modes
+        for n2 in np.arange(1,maxtem+1,2) :
+            c_n1n2= HG_split_diode_coefficient_numerical(n1,n2)
+            if (xy=="x"):       
+                print("{:2d} x {:2d} x {:.15g}".format(n1,n2,c_n1n2))
+            else:
+                print("x {:2d} x {:2d} {:.15g}".format(n1,n2,c_n1n2))
+    print('END')
+
+                
+def HG_split_diode_coefficient_analytical(n1,n2):
+    """ Using an analytical equation to compute beat coefficients
+    betwween mode n1 and n2 on a split photo detector. Uses arbitrary
+    precision (from gmpy) because otherwise errors get very large
+    for n1,n2>30. This is for comparison with the numerical method
+    HG_split_diode_coefficient_numerical only. """
+    import gmpy
+    
+    temp=gmpy.mpq(0.0)
+    for l in np.arange(0,n1/2+1):
+        for k in np.arange(0,(n2-1)/2+1):
+            temp +=  gmpy.mpq(pow((-1.0/4.0),l+k)*gmpy.mpz(factorial((n2+n1-1)/2-l-k))/gmpy.mpz((factorial(l)*factorial(k)*factorial(n1-2*l)*factorial(n2-2*k))))
+
+    c_n1n2=gmpy.mpq(temp*gmpy.mpq(math.sqrt(2.0**(n1+n2)*gmpy.mpz(factorial(n1)*factorial(n2))/np.pi)))
+    return float(gmpy.mpf(c_n1n2))
+    
+    
+def HG_split_diode_coefficient_numerical(n1,n2):
+    """ Compute beam coefficient between mode n1 and n2 on a split
+    photo detector using numerical integration, tested up to n1,n2 = 40.
+    This is primarily used by finesse_split_photodiode()"""
+    A = 2.0 * np.sqrt(2.0/np.pi) * np.sqrt(1.0 / (2.0**(n1+n2) * factorial(n1) * factorial(n2)))
+    f = lambda x: hermite(n1)(np.sqrt(2.0)*x) * math.exp(-2.0*x*x) * hermite(n2)(np.sqrt(2.0)*x)    
+    c_n1n2= quad(f, 0.0, np.Inf, epsabs=1e-10, epsrel=1e-10, limit=200)
+    return A*c_n1n2[0]
+