diff --git a/pykat/math/jacobi.py b/pykat/math/jacobi.py index 2ff8fd8b7b656bd8d27c6887ffbe487b3b295b68..51dc17f095d38214b657c23bdd44394157daafec 100644 --- a/pykat/math/jacobi.py +++ b/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 diff --git a/pykat/math/laguerre.py b/pykat/math/laguerre.py new file mode 100644 index 0000000000000000000000000000000000000000..909197659362e51e840a870459838ee9967c9024 --- /dev/null +++ b/pykat/math/laguerre.py @@ -0,0 +1,28 @@ +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 + diff --git a/pykat/optics/gaussian_beams.py b/pykat/optics/gaussian_beams.py index 43d047008c7aa2c89e57a1aa26172d9e7e967d91..eef4786cb2201ec6296dd2259fe5f5bc04eb10c3 100644 --- a/pykat/optics/gaussian_beams.py +++ b/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 diff --git a/pykat/optics/pdtype.py b/pykat/optics/pdtype.py index 82f46a96fa6d2e04a50d687505685174127b9551..7a96dea2330923d8fe8aaac96fd050a44e42d317 100644 --- a/pykat/optics/pdtype.py +++ b/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')