From 0529779df849daa096e9de3fcfdfc7362ae272cc Mon Sep 17 00:00:00 2001 From: Rob Shalloo Date: Thu, 14 Nov 2024 17:53:57 +0100 Subject: [PATCH] Updating Definition of Hermite-Gaussian Modes (#319) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Maxence Thévenet --- .../transverse/hermite_gaussian_profile.rst | 4 +- docs/source/tutorials/denoised_laser.ipynb | 63 +++++- examples/example_modal_decomposition_data.py | 110 ---------- .../transverse/hermite_gaussian_profile.py | 190 ++++++++++++++---- lasy/utils/mode_decomposition.py | 28 ++- tests/test_laser_profiles.py | 32 +-- 6 files changed, 242 insertions(+), 185 deletions(-) delete mode 100644 examples/example_modal_decomposition_data.py diff --git a/docs/source/api/profiles/transverse/hermite_gaussian_profile.rst b/docs/source/api/profiles/transverse/hermite_gaussian_profile.rst index 8f359b70..0992aaae 100644 --- a/docs/source/api/profiles/transverse/hermite_gaussian_profile.rst +++ b/docs/source/api/profiles/transverse/hermite_gaussian_profile.rst @@ -2,10 +2,8 @@ Hermite Gaussian Transverse Profile =================================== Used to define a Hermite-Gaussian transverse laser profile. -Hermite-Gaussian modes are a family of solutions to the paraxial wave equation written in cartesian coordinates. The modes are characterised by two transverse indices :math:`n_x` and :math:`n_y`. +Hermite-Gaussian modes are a family of solutions to the paraxial wave equation written in cartesian coordinates. The modes are characterised by two transverse indices :math:`m` and :math:`n`. -.. image:: https://user-images.githubusercontent.com/27694869/211018259-15d925bb-f123-42c6-b5ba-86488356ae70.png - :alt: Hermite-Gauss-Modes ------------ diff --git a/docs/source/tutorials/denoised_laser.ipynb b/docs/source/tutorials/denoised_laser.ipynb index c8c899bb..093d5992 100644 --- a/docs/source/tutorials/denoised_laser.ipynb +++ b/docs/source/tutorials/denoised_laser.ipynb @@ -309,19 +309,19 @@ "outputs": [], "source": [ "# Maximum Hermite-Gauss mode index in x and y\n", - "n_modes_x = 2\n", - "n_modes_y = 2\n", + "m_max = 10\n", + "n_max = 10\n", "\n", "# Calculate the decomposition and waist of the laser pulse\n", "modeCoeffs, waist = hermite_gauss_decomposition(\n", - " transverse_profile, n_x_max=n_modes_x, n_y_max=n_modes_y, res=cal\n", + " transverse_profile, longitudinal_profile.lambda0, m_max=m_max, n_max=n_max, res=cal\n", ")\n", "\n", "# Create a num profile, summing over the first few modes\n", "energy_frac = 0\n", "for i, mode_key in enumerate(list(modeCoeffs)):\n", " tmp_transverse_profile = HermiteGaussianTransverseProfile(\n", - " waist, mode_key[0], mode_key[1]\n", + " waist, waist, mode_key[0], mode_key[1], longitudinal_profile.lambda0\n", " )\n", " energy_frac += modeCoeffs[mode_key] ** 2 # Energy fraction of the mode\n", " if i == 0: # First mode (0,0)\n", @@ -365,13 +365,14 @@ "# Determine the figure parameters\n", "fig, ax = plt.subplots(1, 3, figsize=(12, 4), tight_layout=True)\n", "fig.suptitle(\n", - " \"Hermite-Gauss Reconstruction using n_x_max = %i, n_y_max = %i\"\n", - " % (n_modes_x, n_modes_y)\n", + " \"Hermite-Gauss Reconstruction using m_max = %i, n_max = %i\" % (m_max, n_max)\n", ")\n", "\n", "# Plot the original profile\n", "pltextent = np.array([np.min(x), np.max(x), np.min(x), np.max(x)]) * 1e6 # in microns\n", "prof1 = np.abs(laser_profile_raw.evaluate(X, Y, 0)) ** 2\n", + "maxInten = np.max(prof1)\n", + "prof1 /= maxInten\n", "divider0 = make_axes_locatable(ax[0])\n", "ax0_cb = divider0.append_axes(\"right\", size=\"5%\", pad=0.05)\n", "pl0 = ax[0].imshow(prof1, cmap=\"magma\", extent=pltextent, vmin=0, vmax=np.max(prof1))\n", @@ -383,6 +384,7 @@ "\n", "# Plot the reconstructed profile\n", "prof2 = np.abs(laser_profile_cleaned.evaluate(X, Y, 0)) ** 2\n", + "prof2 /= maxInten\n", "divider1 = make_axes_locatable(ax[1])\n", "ax1_cb = divider1.append_axes(\"right\", size=\"5%\", pad=0.05)\n", "pl1 = ax[1].imshow(prof2, cmap=\"magma\", extent=pltextent, vmin=0, vmax=np.max(prof1))\n", @@ -403,6 +405,51 @@ "ax[2].set_ylabel(\"y ($ \\\\mu m $)\")\n", "ax[2].set_title(\"Error\")\n", "\n", + "plt.show()\n", + "\n", + "fig2, ax2 = plt.subplots(2, 1, figsize=(8, 6), tight_layout=True)\n", + "fig2.suptitle(\n", + " \"Hermite-Gauss Reconstruction using m_max = %i, n_max = %i\" % (m_max, n_max)\n", + ")\n", + "ax2[0].plot(\n", + " x * 1e6,\n", + " prof2[int(len(x) / 2), :],\n", + " label=\"Reconstructed Profile\",\n", + " color=(1, 0.5, 0.5),\n", + " lw=2.5,\n", + ")\n", + "ax2[0].plot(\n", + " x * 1e6,\n", + " prof1[int(len(x) / 2), :],\n", + " label=\"Original Profile\",\n", + " color=(0.3, 0.3, 0.3),\n", + " lw=1.0,\n", + ")\n", + "ax2[0].legend()\n", + "ax2[0].set_xlim(pltextent[0], pltextent[1])\n", + "ax2[0].set_xlabel(\"x ($ \\\\mu m $)\")\n", + "ax2[0].set_ylabel(\"Intensity (norm.)\")\n", + "\n", + "ax2[1].plot(\n", + " x * 1e6,\n", + " prof2[int(len(x) / 2), :],\n", + " label=\"Reconstructed Profile\",\n", + " color=(1, 0.5, 0.5),\n", + " lw=2.5,\n", + ")\n", + "ax2[1].plot(\n", + " x * 1e6,\n", + " prof1[int(len(x) / 2), :],\n", + " label=\"Original Profile\",\n", + " color=(0.3, 0.3, 0.3),\n", + " lw=1.0,\n", + ")\n", + "ax2[1].legend()\n", + "ax2[1].set_xlim(pltextent[0], pltextent[1])\n", + "ax2[1].set_yscale(\"log\")\n", + "ax2[1].set_xlabel(\"x ($ \\\\mu m $)\")\n", + "ax2[1].set_ylabel(\"Intensity (norm.)\")\n", + "\n", "plt.show()" ] }, @@ -473,7 +520,7 @@ ], "metadata": { "kernelspec": { - "display_name": ".venv", + "display_name": "base", "language": "python", "name": "python3" }, @@ -487,7 +534,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.0" + "version": "3.12.7" } }, "nbformat": 4, diff --git a/examples/example_modal_decomposition_data.py b/examples/example_modal_decomposition_data.py deleted file mode 100644 index 2ddefade..00000000 --- a/examples/example_modal_decomposition_data.py +++ /dev/null @@ -1,110 +0,0 @@ -import matplotlib.pyplot as plt -import numpy as np -import skimage -from mpl_toolkits.axes_grid1 import make_axes_locatable - -from lasy.profiles.combined_profile import CombinedLongitudinalTransverseProfile -from lasy.profiles.longitudinal.gaussian_profile import GaussianLongitudinalProfile -from lasy.profiles.transverse.hermite_gaussian_profile import ( - HermiteGaussianTransverseProfile, -) -from lasy.profiles.transverse.transverse_profile_from_data import ( - TransverseProfileFromData, -) -from lasy.utils.mode_decomposition import hermite_gauss_decomposition - -# Define the transverse profile of the laser pulse -img_url = "https://user-images.githubusercontent.com/27694869/228038930-d6ab03b1-a726-4b41-a378-5f4a83dc3778.png" -intensityData = skimage.io.imread(img_url) -intensityData[intensityData < 2.1] = 0 -pixel_calib = 0.186e-6 -lo = ( - -intensityData.shape[0] / 2 * pixel_calib, - -intensityData.shape[1] / 2 * pixel_calib, -) -hi = ( - intensityData.shape[0] / 2 * pixel_calib, - intensityData.shape[1] / 2 * pixel_calib, -) -energy = 0.5 -pol = (1, 0) -transverse_profile = TransverseProfileFromData(intensityData, lo, hi) - -# Define longitudinal profile of the laser pulse -wavelength = 800e-9 -tau = 30e-15 -t_peak = 0.0 -longitudinal_profile = GaussianLongitudinalProfile(wavelength, tau, t_peak) - -# Combine into full laser profile -profile = CombinedLongitudinalTransverseProfile( - wavelength, pol, energy, longitudinal_profile, transverse_profile -) - -# Calculate the decomposition into hermite-gauss modes -n_x_max = 20 -n_y_max = 20 -modeCoeffs, waist = hermite_gauss_decomposition( - transverse_profile, n_x_max=n_x_max, n_y_max=n_y_max, res=0.5e-6 -) - -# Reconstruct the pulse using a series of hermite-gauss modes -for i, mode_key in enumerate(list(modeCoeffs)): - tmp_transverse_profile = HermiteGaussianTransverseProfile( - waist, mode_key[0], mode_key[1] - ) - if i == 0: - reconstructedProfile = modeCoeffs[ - mode_key - ] * CombinedLongitudinalTransverseProfile( - wavelength, pol, energy, longitudinal_profile, tmp_transverse_profile - ) - else: - reconstructedProfile += modeCoeffs[ - mode_key - ] * CombinedLongitudinalTransverseProfile( - wavelength, pol, energy, longitudinal_profile, tmp_transverse_profile - ) - -# Plotting the results -x = np.linspace(-5 * waist, 5 * waist, 500) -X, Y = np.meshgrid(x, x) - -fig, ax = plt.subplots(1, 3, figsize=(12, 4), tight_layout=True) - -pltextent = (np.min(x) * 1e6, np.max(x) * 1e6, np.min(x) * 1e6, np.max(x) * 1e6) -prof1 = np.abs(profile.evaluate(X, Y, 0)) ** 2 -divider0 = make_axes_locatable(ax[0]) -ax0_cb = divider0.append_axes("right", size="5%", pad=0.05) -pl0 = ax[0].imshow(prof1, cmap="magma", extent=pltextent, vmin=0, vmax=np.max(prof1)) -cbar0 = fig.colorbar(pl0, cax=ax0_cb) -cbar0.set_label("Intensity (norm.)") -ax[0].set_xlabel("x (micron)") -ax[0].set_ylabel("y (micron)") -ax[0].set_title("Original Profile") - -prof2 = np.abs(reconstructedProfile.evaluate(X, Y, 0)) ** 2 -divider1 = make_axes_locatable(ax[1]) -ax1_cb = divider1.append_axes("right", size="5%", pad=0.05) -pl1 = ax[1].imshow(prof2, cmap="magma", extent=pltextent, vmin=0, vmax=np.max(prof1)) -cbar1 = fig.colorbar(pl1, cax=ax1_cb) -cbar1.set_label("Intensity (norm.)") -ax[1].set_xlabel("x (micron)") -ax[1].set_ylabel("y (micron)") -ax[1].set_title("Reconstructed Profile") - - -prof3 = (prof1 - prof2) / np.max(prof1) -divider2 = make_axes_locatable(ax[2]) -ax2_cb = divider2.append_axes("right", size="5%", pad=0.05) -pl2 = ax[2].imshow(100 * np.abs(prof3), cmap="magma", extent=pltextent, vmin=0, vmax=2) -cbar2 = fig.colorbar(pl2, cax=ax2_cb) -cbar2.set_label("|Error| (%)") -ax[2].set_xlabel("x (micron)") -ax[2].set_ylabel("y (micron)") -ax[2].set_title("Error") - -fig.suptitle( - "Hermite-Gauss Reconstruction using n_x_max = %i, n_y_max = %i" % (n_x_max, n_y_max) -) -plt.show() diff --git a/lasy/profiles/transverse/hermite_gaussian_profile.py b/lasy/profiles/transverse/hermite_gaussian_profile.py index 8511361c..434caee3 100644 --- a/lasy/profiles/transverse/hermite_gaussian_profile.py +++ b/lasy/profiles/transverse/hermite_gaussian_profile.py @@ -10,31 +10,58 @@ class HermiteGaussianTransverseProfile(TransverseProfile): r""" A high-order Gaussian laser pulse expressed in the Hermite-Gaussian formalism. + Definition is according to Siegman "Lasers" pg. 646 eq. 60, as explicitly given + in https://doi.org/10.1364/JOSAB.489884 eq. 3, with the beam center upon the + optical axis, :math:`x_0,y_0 = (0,0)`. + More precisely, the transverse envelope (to be used in the :class:`.CombinedLongitudinalTransverseLaser` class) corresponds to: .. math:: - \mathcal{T}(x, y) = \, - \sqrt{\frac{2}{\pi}} \sqrt{\frac{1}{2^{n} n! w_0}}\, - \sqrt{\frac{1}{2^{n} n! w_0}}\, - H_{n_x}\left ( \frac{\sqrt{2} x}{w_0}\right )\, - H_{n_y}\left ( \frac{\sqrt{2} y}{w_0}\right )\, - \exp\left( -\frac{x^2+y^2}{w_0^2} \right) + \mathcal{H}_m (x) \, \mathcal{H}_n(y) \, \exp(i \left ( \Phi_x(z) + \Phi_y(z)\right ) ) + + with + + .. math:: + \mathcal{H}_p(q) = A_p h_p \left ( \frac{\sqrt{2}q}{w_q(z)} \right) \exp{\left( -\frac{q^2}{w_q^2(z)}\right)} \exp{ \left ( -i k_0 \frac{q^2}{2 R_q(z)} \right )} + + \Phi_q(z) = \left(p+\frac{1}{2}\right) \arctan\left({\frac{z}{Z_q}}\right) + + w_q(z) = w_{0,q} \sqrt{1 + \left( \frac{z}{Z_q}\right)^2} + + Z_q = \frac{\pi w_{0,q}^2}{\lambda_0} + + A_p = \frac{1}{\sqrt{w_q(z) 2^{p-1/2} p!\sqrt{\pi}}} + + R_q(z) = z + \frac{Z_q^2}{z} + + - where :math:`H_{n}` is the Hermite polynomial of order :math:`n`. + + + where :math:`h_{p}` is the Hermite polynomial of order :math:`p`, :math:`w_q(z)` is the + spot size of the laser along the :math:`q` axis (:math:`q` is :math:`x` or :math:`y`), :math:`Z_q` is the corresponding Rayleigh + length and, :math:`\lambda_0` and :math:`k_0` are the wavelength and central wavenumber of the + laser respectively. + + + The z-dependence shown in the above equations is required to correctly define the + electric field of the transverse profile relative to that of the pulse at the focus. + The absolute z position will be overwritten when creating a laser object. Parameters ---------- - w0 : float (in meter) - The waist of the laser pulse, i.e. :math:`w_0` in the above formula. - n_x : int (dimensionless) + w_0x : float (in meter) + The waist of the laser pulse in the x direction, + w_0y : float (in meter) + The waist of the laser pulse in the y direction, + m : int (dimensionless) The order of hermite polynomial in the x direction - n_y : int (dimensionless) + n : int (dimensionless) The order of hermite polynomial in the y direction wavelength : float (in meter), optional The main laser wavelength :math:`\lambda_0` of the laser. - (Only needed if ``z_foc`` is different than 0.) z_foc : float (in meter), optional Position of the focal plane. (The laser pulse is initialized at ``z=0``.) @@ -50,20 +77,91 @@ class HermiteGaussianTransverseProfile(TransverseProfile): Both methods are in principle equivalent, but note that the first method uses the paraxial approximation, while the second method does not make this approximation. + + Examples + -------- + >>> import matplotlib.pyplot as plt + >>> import numpy as np + >>> from lasy.profiles.transverse.hermite_gaussian_profile import ( + ... HermiteGaussianTransverseProfile, + ... ) + >>> # Create evaluation grid + >>> xy = np.linspace(-30e-6, 30e-6, 200) + >>> X, Y = np.meshgrid(xy, xy) + >>> # Create an array of plots + >>> fig, ax = plt.subplots(3, 6, figsize=(10, 5), tight_layout=True) + >>> extent = (1e6 * xy[0], 1e6 * xy[-1], 1e6 * xy[0], 1e6 * xy[-1]) + >>> for m in range(3): + >>> for n in range(3): + >>> transverse_profile = HermiteGaussianTransverseProfile( + ... w_0x = 10e-6, # m + ... w_0y = 15e-6, # m + ... m = m, # + ... n = n, # + ... wavelength = 0.8e-6, # m + ... ) + ... intensity = np.abs(transverse_profile.evaluate(X,Y))**2 + ... vmax_intensity = np.max(intensity) + >>> ax[m,n].imshow(intensity,extent=extent,cmap='bone_r',vmin=0,vmax=vmax_intensity) + >>> ax[m,n].set_title('Inten: m,n = %i,%i' %(m,n)) + >>> phase = np.angle(transverse_profile.evaluate(X,Y)) + ... vmax_phase = np.max(np.abs(phase)) + >>> ax[m,n+3].imshow(phase,extent=extent,cmap='seismic',vmin=-vmax_phase,vmax=vmax_phase) + >>> ax[m,n+3].set_title('Phase: m,n = %i,%i' %(m,n)) + >>> if m==2: + >>> ax[m,n].set_xlabel("x (µm)") + >>> ax[m,n+3].set_xlabel("x (µm)") + >>> else: + >>> ax[m,n].set_xticks([]) + >>> ax[m,n+3].set_xticks([]) + >>> if n==0: + >>> ax[m,n].set_ylabel("y (µm)") + >>> ax[m,n+3].set_yticks([]) + >>> else: + >>> ax[m,n].set_yticks([]) + >>> ax[m,n+3].set_yticks([]) + + + """ - def __init__(self, w0, n_x, n_y, wavelength=None, z_foc=0): + def __init__(self, w_0x, w_0y, m, n, wavelength, z_foc=0): super().__init__() - self.w0 = w0 - self.n_x = n_x - self.n_y = n_y - if z_foc == 0: - self.z_foc_over_zr = 0 - else: - assert ( - wavelength is not None - ), "You need to pass the wavelength, when `z_foc` is non-zero." - self.z_foc_over_zr = z_foc * wavelength / (np.pi * w0**2) + self.w_0x = w_0x + self.w_0y = w_0y + self.m = m + self.n = n + self.wavelength = wavelength + self.z_foc = z_foc + z_eval = -z_foc # this links our observation position to Siegmann's definition + + self.k0 = 2 * np.pi / wavelength + + # Calculate Rayleigh Lengths + Zx = np.pi * w_0x**2 / wavelength + Zy = np.pi * w_0y**2 / wavelength + + # Calculate Size at Location Z + wxZ = w_0x * np.sqrt(1 + (z_eval / Zx) ** 2) + wyZ = w_0y * np.sqrt(1 + (z_eval / Zy) ** 2) + + # Calculate Multiplicative Factors + Anx = 1 / np.sqrt(wxZ * 2 ** (m - 1 / 2) * factorial(m) * np.sqrt(np.pi)) + Any = 1 / np.sqrt(wyZ * 2 ** (n - 1 / 2) * factorial(n) * np.sqrt(np.pi)) + + # Calculate the Phase contributions from propagation + phiXz = (m + 1 / 2) * np.arctan2(z_eval, Zx) + phiYz = (n + 1 / 2) * np.arctan2(z_eval, Zy) + + self.z_eval = z_eval + self.Zx = Zx + self.Zy = Zy + self.wxZ = wxZ + self.wyZ = wyZ + self.Anx = Anx + self.Any = Any + self.phiXz = phiXz + self.phiYz = phiYz def _evaluate(self, x, y): """ @@ -81,22 +179,34 @@ def _evaluate(self, x, y): Contains the value of the envelope at the specified points This array has the same shape as the arrays x, y """ - # Term for wavefront curvature, waist and Gouy phase - diffract_factor = 1.0 - 1j * self.z_foc_over_zr - w = self.w0 * abs(diffract_factor) - psi = np.angle(diffract_factor) - - envelope = ( - np.sqrt(2 / np.pi) - * np.sqrt(1 / (2 ** (self.n_x) * factorial(self.n_x) * self.w0)) - * np.sqrt(1 / (2 ** (self.n_y) * factorial(self.n_y) * self.w0)) - * hermite(self.n_x)(np.sqrt(2) * x / w) - * hermite(self.n_y)(np.sqrt(2) * y / w) - * np.exp( - -(x**2 + y**2) / (self.w0**2 * diffract_factor) - - 1.0j * (self.n_x + self.n_y) * psi - ) - # Additional Gouy phase - * (1.0 / diffract_factor) + z_eval = self.z_eval + Zx = self.Zx + Zy = self.Zy + k0 = self.k0 + wxZ = self.wxZ + wyZ = self.wyZ + Anx = self.Anx + Any = self.Any + m = self.m + n = self.n + phiXz = self.phiXz + phiYz = self.phiYz + + # Calculate the HG in each plane + HGnx = ( + Anx + * hermite(m)(np.sqrt(2) * (x) / wxZ) + * np.exp(-((x) ** 2) / wxZ**2) + * np.exp(-1j * k0 * (x) ** 2 / 2 / (z_eval**2 + Zx**2) * z_eval) ) + HGny = ( + Any + * hermite(n)(np.sqrt(2) * (y) / wyZ) + * np.exp(-((y) ** 2) / wyZ**2) + * np.exp(-1j * k0 * (y) ** 2 / 2 / (z_eval**2 + Zy**2) * z_eval) + ) + + # Put it altogether + envelope = HGnx * HGny * np.exp(1j * (phiXz + phiYz)) + return envelope diff --git a/lasy/utils/mode_decomposition.py b/lasy/utils/mode_decomposition.py index 52ab6051..0c386db1 100644 --- a/lasy/utils/mode_decomposition.py +++ b/lasy/utils/mode_decomposition.py @@ -9,7 +9,7 @@ from lasy.utils.exp_data_utils import find_d4sigma -def hermite_gauss_decomposition(laserProfile, n_x_max=12, n_y_max=12, res=1e-6): +def hermite_gauss_decomposition(laserProfile, wavelength, m_max=12, n_max=12, res=1e-6): """ Decomposes a laser profile into a set of hermite-gaussian modes. @@ -20,8 +20,11 @@ def hermite_gauss_decomposition(laserProfile, n_x_max=12, n_y_max=12, res=1e-6): laserProfile : class instance An instance of a class or sub-class of TransverseLaserProfile - n_x_max, n_y_max : ints - The maximum values of `n_x` and `n_y` out to which the expansion + wavelength : float (in meter) + Central wavelength at which the Hermite-Gauss beams are to be defined. + + m_max, n_max : ints + The maximum values of `m` and `n` up to which the expansion will be performed res : float @@ -33,7 +36,7 @@ def hermite_gauss_decomposition(laserProfile, n_x_max=12, n_y_max=12, res=1e-6): weights : dict of floats A dictionary of floats corresponding to the weights of each mode in the decomposition. The keys of the dictionary are tuples - corresponding to (`n_x`,`n_y`) + corresponding to (`m`,`n`) waist : Beam waist for which the decomposition is calculated. It is computed as the waist for which the weight of order 0 is maximum. @@ -74,24 +77,24 @@ def hermite_gauss_decomposition(laserProfile, n_x_max=12, n_y_max=12, res=1e-6): field = laserProfile.evaluate(X, Y) # Get estimate of w0 - w0 = estimate_best_HG_waist(x, y, field) + w0 = estimate_best_HG_waist(x, y, field, wavelength) # Next we loop over the modes and calculate the relevant weights weights = {} - for i in range(n_x_max): - for j in range(n_y_max): - HGMode = HermiteGaussianTransverseProfile(w0, i, j) + for m in range(m_max): + for n in range(n_max): + HGMode = HermiteGaussianTransverseProfile(w0, w0, m, n, wavelength) coef = np.real( np.sum(field * HGMode.evaluate(X, Y)) * dx * dy ) # modalDecomposition if math.isnan(coef): coef = 0 - weights[(i, j)] = coef + weights[(m, n)] = coef return weights, w0 -def estimate_best_HG_waist(x, y, field): +def estimate_best_HG_waist(x, y, field, wavelength): """ Estimate the waist that maximises the weighting of the first mode. @@ -109,6 +112,9 @@ def estimate_best_HG_waist(x, y, field): field : 2D numpy array representing the field (not the laser intensity). the laser field profile in a 2D slice. + wavelength : float (in meter) + Central wavelength at which the Hermite-Gauss beams are to be defined. + Returns ------- w0 : scalar @@ -131,7 +137,7 @@ def estimate_best_HG_waist(x, y, field): for i, wTest in enumerate(waistTest): # create a gaussian - HGMode = HermiteGaussianTransverseProfile(wTest, 0, 0) + HGMode = HermiteGaussianTransverseProfile(wTest, wTest, 0, 0, wavelength) profile = HGMode.evaluate(X, Y) coeffTest[i] = np.real(np.sum(profile * field)) w0 = waistTest[np.argmax(coeffTest)] diff --git a/tests/test_laser_profiles.py b/tests/test_laser_profiles.py index 5bdf6f81..af1f059a 100644 --- a/tests/test_laser_profiles.py +++ b/tests/test_laser_profiles.py @@ -113,22 +113,28 @@ def test_transverse_profiles_rt(): def test_transverse_profiles_3d(): - npoints = 200 - w0 = 10.0e-6 + w0_x = 10.0e-6 + w0_y = 12.0e-6 # HermiteGaussianTransverseProfile print("HermiteGaussianTransverseProfile") - n_x = 2 - n_y = 2 - std_th = np.sqrt(5.0 / 4) * w0 - profile = HermiteGaussianTransverseProfile(w0, n_x, n_y) - x = np.linspace(-4 * w0, 4 * w0, npoints) - y = np.zeros_like(x) - field = profile.evaluate(x, y) - std = np.sqrt(np.average(x**2, weights=np.abs(field) ** 2)) - print("std_th = ", std_th) - print("std = ", std) - assert np.abs(std - std_th) / std_th < 0.01 + m = 2 + n = 2 + std_th_x = np.sqrt(5.0 / 4) * w0_x + std_th_y = np.sqrt(5.0 / 4) * w0_y + profile = HermiteGaussianTransverseProfile(w0_x, w0_y, m, n, wavelength=800e-9) + x = np.linspace(-4 * w0_x, 4 * w0_x, 200) + y = np.linspace(-4 * w0_y, 4 * w0_y, 150) + field_x = profile.evaluate(x, np.zeros_like(x)) + field_y = profile.evaluate(np.zeros_like(y), y) + std_x = np.sqrt(np.average(x**2, weights=np.abs(field_x) ** 2)) + std_y = np.sqrt(np.average(y**2, weights=np.abs(field_y) ** 2)) + print("std_th_x = ", std_th_x) + print("std_x = ", std_x) + print("std_th_y = ", std_th_y) + print("std_y = ", std_y) + assert np.abs(std_x - std_th_x) / std_th_x < 0.01 + assert np.abs(std_y - std_th_y) / std_th_y < 0.01 # TransverseProfileFromData print("TransverseProfileFromData")