diff --git a/pykat/optics/gaussian_beams.py b/pykat/optics/gaussian_beams.py index 390cd522a38fa752102deb2d20b7bcb85cc9b1cb..8be799f5e36849e172b7ba19edc95fb83c019f3f 100644 --- a/pykat/optics/gaussian_beams.py +++ b/pykat/optics/gaussian_beams.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 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. @@ -316,9 +268,15 @@ 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() diff --git a/pykat/optics/pdtype.py b/pykat/optics/pdtype.py new file mode 100644 index 0000000000000000000000000000000000000000..82f46a96fa6d2e04a50d687505685174127b9551 --- /dev/null +++ b/pykat/optics/pdtype.py @@ -0,0 +1,72 @@ +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] +