adding functions to compute bullseye coefficients

pykat/math/jacobi.py

@@ -7,9 +7,23 @@ import numpy as np
 import scipy.special
 
 def jacobi(n,a,b,x):
-  """ Implementation of Jacobi fucntion using binominal coefficients.
+  """
+  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 cannot."""
+  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

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

@@ -10,7 +10,7 @@ import copy
 import warnings
 import cmath
 from math import factorial
-from scipy.special import hermite, eval_jacobi
+from scipy.special import hermite
 from pykat.math.jacobi import jacobi
 from pykat.SIfloat import SIfloat
 
@@ -271,9 +271,8 @@ class gauss_param(object):
 class beam_param(gauss_param):
     pass
 
-# Should be renamed to HG_mode?    
 class HG_mode(object):
-    """ Hermite-Gauss beam profile. Example usage:
+    """ 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_mode(qx,n=2,m=0)
@@ -441,7 +440,7 @@ def HG2LG(n,m):
 
         # 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 * eval_jacobi(m,np.abs(l)+p-m,p-m,0.0)
+        c = c * (-1.0)**p * (-2.0)**m * jacobi(m,np.abs(l)+p-m,p-m,0.0)
 
         coefficients[j] = c
         

pykat/optics/pdtype.py

@@ -13,6 +13,10 @@ from math import factorial
 from scipy.integrate import quad
 from scipy.special import hermite
 from pykat.SIfloat import SIfloat
+#import scipy.special
+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 +74,60 @@ 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):
+    """
+    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):
+    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.5887050, name='bullseye'):
+    """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, 0.01, 0.03, name="bull1:3")
+    The default kat.ini numbers have been generated with:
+    ...
+    """
+    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)>10e-13 and not np.isnan(c)):
+                        print("{:2d} {:2d} {:2d} {:2d} {:.15g}".format(n1,m1,n2,m2, c))
+    print('END')