diff --git a/pykat/maths/__init__.py b/pykat/math/__init__.py similarity index 100% rename from pykat/maths/__init__.py rename to pykat/math/__init__.py diff --git a/pykat/maths/hermite.py b/pykat/math/hermite.py similarity index 95% rename from pykat/maths/hermite.py rename to pykat/math/hermite.py index f300ad42f7746b24c88fdee764d38adba1a97a8e..9a295b85c47ef88a0024e683e725b27da647ae0c 100644 --- a/pykat/maths/hermite.py +++ b/pykat/math/hermite.py @@ -1,3 +1,8 @@ +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): diff --git a/pykat/math/jacobi.py b/pykat/math/jacobi.py new file mode 100644 index 0000000000000000000000000000000000000000..51dc17f095d38214b657c23bdd44394157daafec --- /dev/null +++ b/pykat/math/jacobi.py @@ -0,0 +1,31 @@ +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 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/maths/zernike.py b/pykat/math/zernike.py similarity index 100% rename from pykat/maths/zernike.py rename to pykat/math/zernike.py diff --git a/pykat/optics/gaussian_beams.py b/pykat/optics/gaussian_beams.py index 8be799f5e36849e172b7ba19edc95fb83c019f3f..eef4786cb2201ec6296dd2259fe5f5bc04eb10c3 100644 --- a/pykat/optics/gaussian_beams.py +++ b/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 diff --git a/pykat/optics/knm.py b/pykat/optics/knm.py index 9a412baa4ba4b6417e9a9a1f6b21bf55e3161ce6..c36aaf5b44979f1496990799a1b29885567a223f 100644 --- a/pykat/optics/knm.py +++ b/pykat/optics/knm.py @@ -1,5 +1,5 @@ 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() diff --git a/pykat/optics/maps.py b/pykat/optics/maps.py index d086ee10dc32613bbd330370367d8b23ec56584a..21c908b6c677e1b1d9db7c633a7e64d668b18353 100644 --- a/pykat/optics/maps.py +++ b/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 diff --git a/pykat/optics/pdtype.py b/pykat/optics/pdtype.py index 82f46a96fa6d2e04a50d687505685174127b9551..4bd75e3983f448ef7fcfc5b842305b8e5648e2c9 100644 --- a/pykat/optics/pdtype.py +++ b/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') diff --git a/pykat/optics/romhom.py b/pykat/optics/romhom.py index 238616e6748aa49add0028f8be58be0ae3dc6ef6..392523581af8502b4999b75e23ab326d86efb004 100644 --- a/pykat/optics/romhom.py +++ b/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)