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

pykat/maths/hermite.py → pykat/math/hermite.py

+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
 import numpy as np
 
 def hermite(n, x): 

pykat/math/jacobi.py

+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+import numpy as np
+import scipy.special
+
+def jacobi(n,a,b,x):
+  """
+  Jacobi Polynomial P_n^{a,b} (x)for real x.
+
+                  n    / n+a \ / n+b \ / x-1 \^(n-s) / x+1 \^s
+  P_n^{a,b}(x)= Sum    |     | |     | | --- |       | --- |
+                  s=0   \ s  / \ n-s / \  2  /       \  2  /
+  
+  n,a,b (int)
+  x (real)
+  P (real)
+
+  Implementation of Jacobi fucntion using binominal coefficients.
+  This can handle values of alpha, beta < -1 which the special.eval_jacobi
+  function does not.
+  Andreas Freise 15.05.2016
+  """
+
+  P=0.0
+  for s in np.arange(0,n+1):
+    P=P+scipy.special.binom(n+a,s) * scipy.special.binom(n+b,n-s) * (x-1.0)**(n-s) * (x+1.0)**s
+  P=P*0.5**n
+  return P

pykat/math/laguerre.py

+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+import numpy as np
+import scipy.special
+from math import factorial
+
+def laguerre(p,l,x):
+    """ function to evaluate associated Laguerre Polynomial L_p^l (x).
+    Usage: L = laguerre (p,l,x)
+
+                   p    1   / l+p \       
+    L_p^l(x)=   Sum    ---  |     |  (-x)^j
+                 j=0   j!   \ p-j /    
+    p,l (int)
+    x (real)
+    L (real)
+    Andreas Freise 15.05.2016
+    """
+
+    L=0.0
+    for j in np.arange(0,p+1):
+        L = L + scipy.special.binom(l+p,p-j) / factorial(j) * (-x)**j
+    
+    return L
+

pykat/optics/gaussian_beams.py

@@ -9,7 +9,9 @@ import math
 import copy
 import warnings
 import cmath
+from math import factorial
 from scipy.special import hermite
+from pykat.math.jacobi import jacobi
 from pykat.SIfloat import SIfloat
 
 class gauss_param(object):
@@ -269,12 +271,11 @@ 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:
+class HG_mode(object):
+    """ Hermite-Gauss mode 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=gb.HG_mode(qx,n=2,m=0)
     beam.plot()
     """    
     def __init__(self, qx, qy=None, n=0, m=0):
@@ -385,7 +386,7 @@ class HG_beam(object):
         return np.outer(_un, _um)
         
     def plot(self, ndx=100, ndy=100, xscale=4, yscale=4):
-        """ Make a simple plot the HG_beam """
+        """ Make a simple plot the HG_mode """
         import pykat.plotting 
         import matplotlib.pyplot as plt
         
@@ -404,4 +405,96 @@ class HG_beam(object):
         cbar = fig.colorbar(axes)
         plt.show()
         
+
+def HG2LG(n,m):
+    """A function for Matlab which returns the coefficients and mode indices of
+    the LG modes required to create a particular HG mode.
+    Usage: coefficients,ps,ls = HG2LG(n,m)
+    
+    n,m:          Indces of the HG mode to re-create.
+    coeffcients:  Complex coefficients for each order=n+m LG mode required to
+                  re-create HG_n,m.
+    ps,ls:        LG mode indices corresponding to coefficients.
+    """
+    # Mode order
+    N = n+m;
+    
+    # Create empty vectors for LG coefficients/ indices
+    coefficients = np.linspace(0,0,N+1,dtype=np.complex_)
+    ps = np.linspace(0,0,N+1)
+    ls = np.linspace(0,0,N+1)
+    
+    # Calculate coefficients
+    for j in np.arange(0,N+1):
+        
+        # Indices for coefficients
+        l = 2*j-N
+        p = int((N-np.abs(l))/2)
+        
+        ps[j] = p
+        ls[j] = l
+        
+        signl = np.sign(l)
+        if (l==0):
+            signl = 1.0
+
+        # Coefficient
+        c = (signl*1j)**m * np.sqrt(factorial(N-m)*factorial(m)/(2**N * factorial(np.abs(l)+p)*factorial(p)))
+        c = c * (-1.0)**p * (-2.0)**m * jacobi(m,np.abs(l)+p-m,p-m,0.0)
+
+        coefficients[j] = c
+        
+    return coefficients, ps, ls 
         
+
+    
+def LG2HG(p,l):
+    """ Function to compute the amplitude coefficients
+    of Hermite-Gauss modes whose sum yields a Laguerre Gauss mode
+    of order n,m.
+    Usage: coefficients, ns, ms = LG2HG(p,l)
+    p:    Radial LG index
+    l:    Azimuthal LG index
+    The LG mode is written as u_pl with 0<=|l|<=p.
+    The output is a series of coefficients for u_nm modes,
+    given as complex numbers and respective n,m indices
+    coefficients (complex array): field amplitude for mode u_nm
+    ns (int array): n-index of u_nm
+    ms (int array): m-index of u_nm
+
+    
+    The formula is adpated from M.W. Beijersbergen et al 'Astigmatic
+    laser mode converters and transfer of orbital angular momentum',
+    Opt. Comm. 96 123-132 (1993)
+    We adapted coefficients to be compatible with our
+    definition of an LG mode, it differs from
+    Beijersbergen by a (-1)^p factor and has exp(il\phi) rather
+    than exp(-il\phi).  Also adapted for allowing -l.
+    Andreas Freise, Charlotte Bond    25.03.2007"""
+
+    # Mode order
+    N=2*p+np.abs(l)
+
+    # Indices for coefficients
+    n = np.abs(l)+p
+    m = p
+
+    # create empty vectors
+    coefficients = np.linspace(0,0,N+1,dtype=np.complex_)
+    ns = np.linspace(0,0,N+1)
+    ms = np.linspace(0,0,N+1)
+
+    # l positive or negative
+    signl = np.sign(l)
+    if (l==0):
+        signl = 1.0
+    
+    # Beijersbergen coefficients
+    for j in np.arange(0,N+1):
+        ns[j]=N-j
+        ms[j]=j
+
+        c=(-signl*1j)**j * math.sqrt(factorial(N-j)*factorial(j)/(2**N * factorial(n)*factorial(m)))
+        coefficients[j] = c * (-1.0)**p * (-2)**j * jacobi(j,n-j,m-j,0.0)
+    
+    return coefficients, ns, ms

pykat/optics/knm.py

 from itertools import combinations_with_replacement as combinations
-from pykat.optics.gaussian_beams import beam_param, HG_beam
+from pykat.optics.gaussian_beams import beam_param, HG_mode
 from pykat.exceptions import BasePyKatException
 from pykat.optics.romhom import u_star_u
 from pykat.external.progressbar import ProgressBar, ETA, Percentage, Bar
@@ -7,10 +7,10 @@ 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 pykat.math.hermite import hermite
 from scipy.misc import comb
 from scipy.integrate import newton_cotes
-from pykat.maths import newton_weights
+from pykat.math import newton_weights
 
 import time
 import pykat.optics.maps
@@ -65,8 +65,8 @@ def adaptive_knm(mode_in, mode_out, q1, q2, q1y=None, q2y=None, smap=None, delta
     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])
+    Hg_in  = HG_mode(qx=q1, qy=q1y, n=mode_in[0], m=mode_in[1])
+    Hg_out = HG_mode(qx=q2, qy=q2y, n=mode_out[0], m=mode_out[1])
     
     Nfuncs = []
     Nfuncs.append(0)
@@ -168,8 +168,8 @@ def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None,
     dy = abs(y[1] - y[0])    
         
     if cache is 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])
+        Hg_in  = HG_mode(qx=q1, qy=q1y, n=mode_in[0], m=mode_in[1])
+        Hg_out = HG_mode(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()
@@ -505,8 +505,8 @@ def square_aperture_HG_knm(mode_in, mode_out, q, R):
     m = mode_in[1]
     _m = mode_out[1]
     
-    hg1 = HG_beam(q, n=n, m=m)
-    hg2 = HG_beam(q, n=_n, m=_m)
+    hg1 = HG_mode(q, n=n, m=m)
+    hg2 = HG_mode(q, n=_n, m=_m)
         
     kx = hg1.constant_x * hg2.constant_x.conjugate()
     ky = hg1.constant_y * hg2.constant_y.conjugate()

pykat/optics/maps.py

@@ -16,7 +16,7 @@ from __future__ import print_function
 from pykat.optics.romhom import makeWeightsNew
 from scipy.interpolate import interp2d, interp1d
 from scipy.optimize import minimize
-from pykat.maths.zernike import *        
+from pykat.math.zernike import *        
 from pykat.exceptions import BasePyKatException
 from copy import deepcopy
 

pykat/optics/pdtype.py

@@ -13,6 +13,9 @@ from math import factorial
 from scipy.integrate import quad
 from scipy.special import hermite
 from pykat.SIfloat import SIfloat
+from pykat.optics.gaussian_beams import HG2LG
+from pykat.math.jacobi import jacobi
+from pykat.math.laguerre import laguerre
 
 
 def finesse_split_photodiode(maxtem, xy='x'):
@@ -70,3 +73,80 @@ def HG_split_diode_coefficient_numerical(n1,n2):
     c_n1n2= quad(f, 0.0, np.Inf, epsabs=1e-10, epsrel=1e-10, limit=200)
     return A*c_n1n2[0]
     
+
+def LG_bullseye_coefficients_numerical(l1, p1, l2, p2, w, r):
+    """ Function to compute a beat coefficient for two LG modes on a
+    bullseye photo detector, used by finesse_bullseye_photodiode().
+    l1,p1 mode indices of first LG mode
+    l2,p2 mode indices of second LG mode
+    w: beam radius on diode [m]
+    r: radius of inner disk element [m]
+    w = beam radisu [m]
+    r = radius of inner element
+    """
+    if (l1 != l2):
+        return 0.0
+        
+    l = np.abs(l1)
+    # computer pre-factor of integral
+    A = math.sqrt(factorial(p1)*factorial(p2)/(factorial(l+p1)*factorial(l+p2)))
+
+    # definng integrand
+    f = lambda x: x**l * laguerre(p1,l,x) * laguerre(p2,l,x) * math.exp(-x)
+    # defining integration limit
+    int_limit = 2.0 * r**2 / w**2
+    # perform numerical integrations
+    val1, res= quad(f, 0.0, int_limit, epsabs=1e-10, epsrel=1e-10, limit=500)
+    val2, res= quad(f, int_limit, np.Inf, epsabs=1e-10, epsrel=1e-10, limit=500)
+    return A*(val2-val1)
+
+def HG_bullseye_coefficients_numerical(n1, m1, n2, m2, w, r):
+    """ Function to compute a beat coefficient for two HG modes on a
+    bullseye photo detector, used by finesse_bullseye_photodiode().
+    n1,m1 mode indices of first HG mode
+    n2,m2 mode indices of second HG mode
+    w: beam radius on diode [m]
+    r: radius of inner disk element [m]
+
+    """
+    mparity=np.mod(m1+m2,2)
+    nparity=np.mod(n1+n2,2)
+                
+    k = 0.0 + 0.0*1j
+    if (mparity==0 and nparity==0):            
+        a,p1,l1 = HG2LG(n1,m1)
+        b,p2,l2 = HG2LG(n2,m2)
+        for idx1, aa in enumerate(l1):
+            for idx2, bb in enumerate(l2):
+                if l1[idx1]==l2[idx2]:
+                    c = LG_bullseye_coefficients_numerical(l1[idx1],p1[idx1], l2[idx2], p2[idx2], w, r)
+                    k = k + a[idx1] * np.conj(b[idx2]) * c
+    return float(np.real(k))
+    
+def finesse_bullseye_photodiode(maxtem, w=1.0, r=0.5887050112577, name='bullseye'):
+    """Prints beat coefficients for HG modes on a bullseye photo detector
+    in the format for insertion into the kat.ini file of Finesse.
+
+    The default kat.ini numbers have been generated with:
+    finesse_bullseye_photodiode(6)
+
+    maxtem: highest mode order (int)
+    w: beam radius on diode [m]
+    r: radius of inner disk element [m]
+    name: name of entry in kat.ini (string)
+
+    The standard parameter of w=1 and r=0.5887050112577 have been chosen to achieve
+    a small a HG00xHG00 coefficient of <1e-13.
+    """
+    assert (isinstance(maxtem, int) and maxtem >=0), "maxtem must be integer >=0" 
+
+    print('PDTYPE {0}'.format(name))
+
+    for n1 in np.arange(0,maxtem+1):
+        for m1 in np.arange(0,maxtem+1):
+            for n2 in np.arange(n1,maxtem+1):
+                for m2 in np.arange(m1,maxtem+1):
+                    c = HG_bullseye_coefficients_numerical(n1,m1, n2, m2, w, r)
+                    if (np.abs(c)>1e-13 and not np.isnan(c)):
+                        print("{:2d} {:2d} {:2d} {:2d} {:.15g}".format(n1,m1,n2,m2, c))
+    print('END')

pykat/optics/romhom.py

@@ -13,11 +13,11 @@ import itertools
 from copy import copy
 from pykat.external.progressbar import ProgressBar, ETA, Percentage, Bar
 from itertools import combinations_with_replacement as combinations
-from pykat.optics.gaussian_beams import beam_param, HG_beam
+from pykat.optics.gaussian_beams import beam_param
 from scipy.linalg import inv
 from math import factorial
-from pykat.maths.hermite import *
-from pykat.maths import newton_weights
+from pykat.math.hermite import *
+from pykat.math import newton_weights
 from scipy.integrate import newton_cotes
 from multiprocessing import Process, Queue, Array, Value, Event
 from pykat.exceptions import BasePyKatException
@@ -658,4 +658,4 @@ def makeWeightsNew(smap, EIxFilename, EIyFilename=None, verbose=True, newtonCote
     if verbose:
         p.finish()
     
-    return ROMWeights(w_ij_Q1=w_ij_Q1, w_ij_Q2=w_ij_Q2, w_ij_Q3=w_ij_Q3, w_ij_Q4=w_ij_Q4, EIx=EIx, EIy=EIy, direction=direction)
\ No newline at end of file
+    return ROMWeights(w_ij_Q1=w_ij_Q1, w_ij_Q2=w_ij_Q2, w_ij_Q3=w_ij_Q3, w_ij_Q4=w_ij_Q4, EIx=EIx, EIy=EIy, direction=direction)