diff --git a/docs/source/api/utils/index.rst b/docs/source/api/utils/index.rst index d65b1ae57..4194d2915 100644 --- a/docs/source/api/utils/index.rst +++ b/docs/source/api/utils/index.rst @@ -6,6 +6,7 @@ Utils :hidden: grid + mode_decomposition gerchberg_saxton refractive_index laser_utils diff --git a/docs/source/api/utils/mode_decomposition.rst b/docs/source/api/utils/mode_decomposition.rst new file mode 100644 index 000000000..8e61249e9 --- /dev/null +++ b/docs/source/api/utils/mode_decomposition.rst @@ -0,0 +1,6 @@ +Modal decomposition for Hermite- and Laguerre-Gaussians +======================================================= + +.. automodule:: lasy.utils.mode_decomposition + :members: + :undoc-members: diff --git a/docs/source/tutorials/denoised_laser.ipynb b/docs/source/tutorials/denoised_laser.ipynb index cc4d0a1b5..822b25837 100644 --- a/docs/source/tutorials/denoised_laser.ipynb +++ b/docs/source/tutorials/denoised_laser.ipynb @@ -2,15 +2,15 @@ "cells": [ { "cell_type": "markdown", - "id": "05402e6d", + "id": "f4e405ec", "metadata": {}, "source": [ - "# Laser from longitudinal and transverse profiles" + "# Denoised laser from longitudinal and transverse profiles" ] }, { "cell_type": "markdown", - "id": "1a41f95c", + "id": "d16542a9", "metadata": {}, "source": [ "In this tutorial, you will learn to create a Laser from an experimental measurement of the laser spectrum and a 2D image of the transverse intensity profile. We'll start by loading all packages required." @@ -19,28 +19,35 @@ { "cell_type": "code", "execution_count": null, - "id": "dcf32943-cf31-42c3-995c-2ed75c9761d2", + "id": "caec9fbc", "metadata": {}, "outputs": [], "source": [ + "from copy import deepcopy\n", + "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from mpl_toolkits.axes_grid1 import make_axes_locatable\n", "from PIL import Image\n", + "from scipy.constants import c, epsilon_0\n", "\n", "from lasy.laser import Laser\n", "from lasy.profiles.combined_profile import CombinedLongitudinalTransverseProfile\n", - "from lasy.profiles.longitudinal import LongitudinalProfileFromData\n", - "from lasy.profiles.transverse import TransverseProfileFromData\n", - "from lasy.profiles.transverse.hermite_gaussian_profile import (\n", - " HermiteGaussianTransverseProfile,\n", + "from lasy.profiles.longitudinal import (\n", + " ContinuousWaveProfile,\n", + " LongitudinalProfileFromData,\n", ")\n", - "from lasy.utils.mode_decomposition import hermite_gauss_decomposition" + "from lasy.profiles.transverse import TransverseProfileFromData\n", + "from lasy.utils.mode_decomposition import (\n", + " estimate_best_HG_waist,\n", + " hermite_gauss_composition,\n", + " hermite_gauss_decomposition,\n", + ")" ] }, { "cell_type": "markdown", - "id": "c1e972d0", + "id": "3e452d90", "metadata": {}, "source": [ "Next, let's define the physical parameters defining our laser pulse. " @@ -49,18 +56,18 @@ { "cell_type": "code", "execution_count": null, - "id": "e4a2fd86-a83f-490e-9e10-21d0e3ca7c4d", + "id": "122e1b52", "metadata": {}, "outputs": [], "source": [ - "polarization = (1, 0) # Linearly polarized in the x direction\n", + "pol = (1, 0) # Linearly polarized in the x direction\n", "energy_J = 1 # Pulse energy in Joules\n", "cal = 0.2e-6 # Camera pixel size in meters. Used for calibration" ] }, { "cell_type": "markdown", - "id": "bda26c20", + "id": "bc3c2088", "metadata": {}, "source": [ "In what follows, we do a few steps to create a LASY profile. Profile is the object used in LASY to describe the properties of the pulse. This profile is then used to create a Laser object, i.e., a regular mesh with the values of the vector potential of the pulse. This mesh can be 3D for a 3D (xyt) geometry, or can consist of a few 2D grids for cylindrical geometry with mode decomposition." @@ -68,7 +75,7 @@ }, { "cell_type": "markdown", - "id": "f88ab1e4-1f5a-4f07-8b7c-4c105b4aa9eb", + "id": "8509c82f", "metadata": {}, "source": [ "## Reconstruct longitudinal profile" @@ -76,7 +83,7 @@ }, { "cell_type": "markdown", - "id": "e6321294", + "id": "ac3edf8b", "metadata": {}, "source": [ "We start from an example dataset obtained with a [Frequency-resolved optical grating (FROG)](https://en.wikipedia.org/wiki/Frequency-resolved_optical_gating) measurement. This provides us with the spectrum and spectral phase." @@ -85,7 +92,7 @@ { "cell_type": "code", "execution_count": null, - "id": "7ee0575d", + "id": "73db82f3", "metadata": {}, "outputs": [], "source": [ @@ -101,7 +108,7 @@ }, { "cell_type": "markdown", - "id": "6b2b6ab4", + "id": "f9ce8dca", "metadata": {}, "source": [ "Now, we initialize a LASY LongitudinalProfile from experimentally measured spectrum. The central wavelength is calculated in this step, and user later on." @@ -110,7 +117,7 @@ { "cell_type": "code", "execution_count": null, - "id": "e70f703e-15a3-4751-a857-71b17651671a", + "id": "388f700b", "metadata": {}, "outputs": [], "source": [ @@ -131,7 +138,7 @@ }, { "cell_type": "markdown", - "id": "cb7c9541", + "id": "5428264b", "metadata": {}, "source": [ "Plot the longitudinal profile data." @@ -140,7 +147,7 @@ { "cell_type": "code", "execution_count": null, - "id": "e0085e86", + "id": "f17bf576", "metadata": {}, "outputs": [], "source": [ @@ -185,7 +192,7 @@ }, { "cell_type": "markdown", - "id": "1ba9ced4-cc19-48ba-9fec-ddcb6e77128b", + "id": "89d4e79f", "metadata": {}, "source": [ "## Reconstruct transverse profile" @@ -193,7 +200,7 @@ }, { "cell_type": "markdown", - "id": "70d133d0", + "id": "b385844f", "metadata": {}, "source": [ "For the following reconstruction, the data is provided in the .png image format, which can be read using pillow package." @@ -202,7 +209,7 @@ { "cell_type": "code", "execution_count": null, - "id": "bcfff530-c04b-4d28-977c-8143ffe87c62", + "id": "bdec6670", "metadata": {}, "outputs": [], "source": [ @@ -218,7 +225,7 @@ }, { "cell_type": "markdown", - "id": "a563b9c0", + "id": "5a77fa7e", "metadata": {}, "source": [ "Create a LASY TransverseProfile from our experimental data." @@ -227,23 +234,29 @@ { "cell_type": "code", "execution_count": null, - "id": "ac17f984", + "id": "0ed58977", "metadata": {}, "outputs": [], "source": [ + "dimensions = \"xyt\"\n", "nx, ny = intensity_data.shape\n", - "lo = (0, 0) # Lower bounds in x and y\n", - "hi = (ny * cal, nx * cal) # Upper bounds in x and y\n", + "lo = (0, 0, -200e-15) # Add the lower bounds in x and y to the temporal data\n", + "hi = (\n", + " ny * cal,\n", + " nx * cal,\n", + " 200e-5,\n", + ") # Add the upper bounds in x and y to the temporal data\n", + "num_points = (nx, ny, 1)\n", "\n", - "# Create the transverse profile. This also centers the data by default\n", + "# Create the transverse profile\n", "transverse_profile = TransverseProfileFromData(\n", - " intensity_data, [lo[0], lo[1]], [hi[0], hi[1]]\n", + " intensity_data, [lo[0], lo[1]], [hi[0], hi[1]], center_data=True\n", ")" ] }, { "cell_type": "markdown", - "id": "fe9c0b5f", + "id": "cd06c342", "metadata": {}, "source": [ "Plot the original transverse profile from the file. The laser occupies a small fraction of the image." @@ -252,13 +265,14 @@ { "cell_type": "code", "execution_count": null, - "id": "12b0d5eb-ad0a-4777-ae82-d612eba69fb0", + "id": "f21b6b31", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "cax = ax.imshow(\n", " intensity_data,\n", + " cmap=\"Reds\",\n", " aspect=\"auto\",\n", " extent=np.array([lo[0], hi[0], lo[1], hi[1]]) * 1e6,\n", ")\n", @@ -271,125 +285,126 @@ }, { "cell_type": "markdown", - "id": "dd87054b-ab2c-4591-be03-c6bfeae859a1", + "id": "e99cf479", "metadata": {}, "source": [ - "## Combine longitudinal and transverse profiles" + "Evaluate the image data on a sensible grid and calculate waist" ] }, { "cell_type": "code", "execution_count": null, - "id": "ab6ecbd2-0b4e-4251-8709-ce741dc4174e", + "id": "3fa83389", "metadata": {}, "outputs": [], "source": [ "laser_profile_raw = CombinedLongitudinalTransverseProfile(\n", - " wavelength=longitudinal_profile.lambda0,\n", - " pol=polarization,\n", - " laser_energy=energy_J,\n", - " long_profile=longitudinal_profile,\n", - " trans_profile=transverse_profile,\n", - ")" + " longitudinal_profile.lambda0,\n", + " pol,\n", + " ContinuousWaveProfile(longitudinal_profile.lambda0),\n", + " transverse_profile,\n", + " peak_fluence=energy_J,\n", + ")\n", + "\n", + "# Plot the original and denoised profiles\n", + "# Create a grid for plotting\n", + "x = np.linspace(-50e-6, 50e-6, 200)\n", + "X, Y = np.meshgrid(x, x)\n", + "\n", + "# Raw transverse profile\n", + "new_intensity_data = np.abs(transverse_profile.evaluate(X, Y)) ** 2\n", + "w0x, w0y = estimate_best_HG_waist(\n", + " x, x, transverse_profile.evaluate(X, Y), longitudinal_profile.lambda0\n", + ") # Make sure to pass the field not the intensity to estimate spot size\n", + "\n", + "fig, ax = plt.subplots()\n", + "cax = ax.imshow(\n", + " new_intensity_data,\n", + " cmap=\"Reds\",\n", + " aspect=\"auto\",\n", + " extent=np.array([x[0], x[-1], x[0], x[-1]]) * 1e6,\n", + ")\n", + "color_bar = fig.colorbar(cax)\n", + "color_bar.set_label(r\"Fluence (J/cm$^2$)\")\n", + "ax.set_xlabel(\"x ($ \\\\mu m $)\")\n", + "ax.set_ylabel(\"y ($ \\\\mu m $)\")\n", + "plt.show()" ] }, { "cell_type": "markdown", - "id": "3d55f0ea-a2ea-4157-895f-5d5681403851", + "id": "69dfce2b", "metadata": {}, "source": [ "## Clean/denoise the profile\n", "LASY functions can be used for denoising/cleaning. Here, the measured profile is decomposed into Hermite-Gauss modes. The cleaning is done by keeping only the first few modes. \n", - "Take a look at the following [example](https://github.com/LASY-org/lasy/blob/13f0e4515493deca36c1375be1d9e83c7e379d42/examples/example_modal_decomposition_data.py)." + "Take a look at the following [example](https://github.com/LASY-org/lasy/blob/13f0e4515493deca36c1375be1d9e83c7e379d42/examples/example_modal_decomposition_data.py). Note that this is done for a CW since only the spatial profile is important (and assumed to be independent)." ] }, { "cell_type": "code", "execution_count": null, - "id": "82a8a391-c4f1-473f-928e-cc18524e1cec", + "id": "b147af61", "metadata": {}, "outputs": [], "source": [ + "dimensions = \"xyt\"\n", + "nx, ny = new_intensity_data.shape\n", + "lo = (x[0], x[0], -200e-15) # Add the lower bounds in x and y to the temporal data\n", + "hi = (x[-1], x[-1], 200e-5) # Add the upper bounds in x and y to the temporal data\n", + "num_points = (nx, ny, 1)\n", + "\n", + "long_prof = ContinuousWaveProfile(longitudinal_profile.lambda0)\n", + "tran_prof = TransverseProfileFromData(\n", + " new_intensity_data, [lo[0], lo[1]], [hi[0], hi[1]]\n", + ")\n", + "laser_profile = CombinedLongitudinalTransverseProfile(\n", + " longitudinal_profile.lambda0, pol, long_prof, tran_prof, peak_fluence=energy_J\n", + ")\n", + "\n", + "laser_raw = Laser(dimensions, lo, hi, num_points, laser_profile)\n", + "\n", "# Maximum Hermite-Gauss mode index in x and y\n", "m_max = 10\n", "n_max = 10\n", "\n", "# Calculate the decomposition and waist of the laser pulse\n", - "modeCoeffs, waistX, waistY = hermite_gauss_decomposition(\n", - " transverse_profile,\n", - " longitudinal_profile.lambda0,\n", - " m_max=m_max,\n", - " n_max=n_max,\n", - " res=cal,\n", - " lo=[-2e-4, -2e-4],\n", - " hi=[2e-4, 2e-4],\n", - ")\n", + "modes = hermite_gauss_decomposition(laser_raw, w0x, w0y, m_max, n_max)\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", - " waistX, waistY, 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", - " laser_profile_cleaned = modeCoeffs[\n", - " mode_key\n", - " ] * CombinedLongitudinalTransverseProfile(\n", - " longitudinal_profile.lambda0,\n", - " polarization,\n", - " longitudinal_profile,\n", - " tmp_transverse_profile,\n", - " laser_energy=energy_J,\n", - " )\n", - " else: # All other modes\n", - " laser_profile_cleaned += modeCoeffs[\n", - " mode_key\n", - " ] * CombinedLongitudinalTransverseProfile(\n", - " longitudinal_profile.lambda0,\n", - " polarization,\n", - " longitudinal_profile,\n", - " tmp_transverse_profile,\n", - " laser_energy=energy_J,\n", - " )\n", - "\n", - "# Energy loss due to decomposition\n", - "energy_loss = 1 - energy_frac\n", - "print(f\"Energy loss: {energy_loss * 100:.2f}%\")" + "laser_cleaned = deepcopy(laser_raw) # Make a copy of the input grid\n", + "hermite_gauss_composition(laser_cleaned, w0x, w0y, modes, skipAsymmetricModes=False)" ] }, { "cell_type": "code", "execution_count": null, - "id": "fcb44ccf-ad2b-4a9b-bc19-9f615d057335", + "id": "a883e1b9", "metadata": {}, "outputs": [], "source": [ "# Plot the original and denoised profiles\n", - "# Create a grid for plotting\n", - "x = np.linspace(-5 * waistX, 5 * waistX, 500)\n", - "X, Y = np.meshgrid(x, x)\n", - "\n", - "# 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 m_max = %i, n_max = %i\" % (m_max, n_max)\n", ")\n", - "pltextent = np.array([np.min(x), np.max(x), np.min(x), np.max(x)]) * 1e6 # in microns\n", "cbar_labels = [\"Intensity (norm.)\", \"Intensity (norm.)\", \"|Intensity Error| (%)\"]\n", "title = [\"Original Profile\", \"Reconstructed Profile\", \"Error\"]\n", "scale = [\"linear\", \"log\"]\n", "lw = [1.0, 2.5]\n", - "\n", + "pltextent = (\n", + " laser_raw.grid.axes[0].min() * 1e6,\n", + " laser_raw.grid.axes[0].max() * 1e6,\n", + " laser_raw.grid.axes[1].min() * 1e6,\n", + " laser_raw.grid.axes[1].max() * 1e6,\n", + ")\n", "\n", "# Original profile\n", - "prof1 = np.abs(laser_profile_raw.evaluate(X, Y, 0)) ** 2\n", - "maxInten = np.max(prof1)\n", - "prof1 /= maxInten\n", + "intensity = epsilon_0 * c / 2 * np.abs(laser_raw.grid.get_temporal_field()) ** 2\n", + "prof1 = np.sum(intensity, axis=-1).T * laser_raw.grid.dx[-1]\n", "\n", "# Reconstructed profile\n", - "prof2 = np.abs(laser_profile_cleaned.evaluate(X, Y, 0)) ** 2\n", - "prof2 /= maxInten\n", + "intensity = epsilon_0 * c / 2 * np.abs(laser_cleaned.grid.get_temporal_field()) ** 2\n", + "prof2 = np.sum(intensity, axis=-1).T * laser_cleaned.grid.dx[-1]\n", "\n", "# Normalized error\n", "prof3 = (prof1 - prof2) / np.max(prof1)\n", @@ -403,7 +418,7 @@ "for i in range(3):\n", " divider = make_axes_locatable(ax[i])\n", " ax_cb = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n", - " pl = ax[i].imshow(100 * np.abs(prof[i]), cmap=\"magma\", extent=pltextent)\n", + " pl = ax[i].imshow(100 * np.abs(prof[i]), cmap=\"Reds\", extent=pltextent)\n", " cbar = fig.colorbar(pl, cax=ax_cb)\n", " cbar.set_label(cbar_labels[i])\n", " ax[i].set_xlim(pltextent[0], pltextent[1])\n", @@ -432,20 +447,30 @@ " ax2[i].set_xlabel(\"x ($ \\\\mu m $)\")\n", " ax2[i].set_ylabel(\"Intensity (norm.)\")\n", "\n", - "plt.show()" + "plt.show()\n", + "\n", + "# Update the transverse profile to be recentered and denoised\n", + "transverse_profile = TransverseProfileFromData(prof2, [lo[0], lo[1]], [hi[0], hi[1]])\n", + "laser_profile_cleaned = CombinedLongitudinalTransverseProfile(\n", + " longitudinal_profile.lambda0,\n", + " pol,\n", + " longitudinal_profile,\n", + " transverse_profile,\n", + " laser_energy=energy_J,\n", + ")" ] }, { "cell_type": "markdown", - "id": "103c9501-3e91-4f9b-a0bb-d3a5e2320d0c", + "id": "93c2dddd", "metadata": {}, "source": [ - "## Create a laser" + "## Create the laser" ] }, { "cell_type": "markdown", - "id": "db52cd0e", + "id": "63dfc88a", "metadata": {}, "source": [ "Now the hard part is done! From the cleaned profile, we create a LASY Laser object (3D cartesian or cylindrical geometry) and write to a file compliant with the openPMD standard." @@ -454,7 +479,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6b0f1ffd", + "id": "7c9b8f19", "metadata": {}, "outputs": [], "source": [ @@ -468,9 +493,10 @@ " 20,\n", ") # Number of points in each dimension. Use (300, 300, 200) for production.\n", "\n", + "\n", "# Constructing the object using 3D geometry might take a while to run depending on the hardware used.\n", "laser_xyt = Laser(dimensions, lo, hi, num_points, laser_profile_cleaned) # Laser\n", - "laser_xyt.normalize(energy_J * energy_frac) # Normalize the laser energy\n", + "laser_xyt.normalize(energy_J) # Normalize the laser energy\n", "laser_xyt.show()\n", "\n", "# Save the laser object to a file\n", @@ -480,29 +506,15 @@ { "cell_type": "code", "execution_count": null, - "id": "07be1792", + "id": "757ff37c", "metadata": {}, "outputs": [], - "source": [ - "# Then, cylindrical geometry. Here, we also propagate the pulse backwards by 0.5 mm.\n", - "dimensions = \"rt\" # Use cylindrical geometry\n", - "lo = (0, -150e-15)\n", - "hi = (40e-6, 150e-15)\n", - "num_points = (500, 400)\n", - "\n", - "laser = Laser(dimensions, lo, hi, num_points, laser_profile_cleaned)\n", - "laser.normalize(energy_J * energy_frac)\n", - "laser.propagate(-0.5e-3) # Propagate the laser pulse backwards\n", - "laser.show()\n", - "\n", - "# Save the laser object to a file\n", - "laser.write_to_file(\"Laser_rt_propagated\", \"h5\", save_as_vector_potential=True)" - ] + "source": [] } ], "metadata": { "kernelspec": { - "display_name": "venv", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -516,7 +528,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.2" + "version": "3.11.7" } }, "nbformat": 4, diff --git a/docs/source/tutorials/index.rst b/docs/source/tutorials/index.rst index 518188289..4f8191e1a 100644 --- a/docs/source/tutorials/index.rst +++ b/docs/source/tutorials/index.rst @@ -15,6 +15,7 @@ Additionally, a set of static (automatically tested) examples can be found below gaussian_laser.ipynb axiparabola.ipynb + modal_composition.ipynb denoised_laser.ipynb 1d_temporal_laser.ipynb 2d_spatial_laser.ipynb diff --git a/docs/source/tutorials/modal_composition.ipynb b/docs/source/tutorials/modal_composition.ipynb new file mode 100644 index 000000000..2822c6a56 --- /dev/null +++ b/docs/source/tutorials/modal_composition.ipynb @@ -0,0 +1,250 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "8340aefb-af9d-4d60-97e2-c4616504fefa", + "metadata": {}, + "source": [ + "# Constructing and deconstructing beams using modal representations" + ] + }, + { + "cell_type": "markdown", + "id": "47c36cea-b7ca-41fc-9429-276e77257b75", + "metadata": {}, + "source": [ + "The transverse intensity profile of any laser beam can be represented as a summation over infinitely many modes in a basis set. In this tutorial, we consider two different basis sets: the Hermite-Gaussians and the Laguerre-Gaussians. Other sets, such as the Zernike polynomials, are also possible but are not currently implemented in LASY. In this tutorial, beams are constructed using a complex dictionary of the modal coefficients. It is shown that a decomposition of these beams back into the basis sets returns the original mode coefficients, provided the assumed beam waist remains the same. First, we import the required modules:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "944abf27-8d51-4e93-a29b-8fd44830a78b", + "metadata": {}, + "outputs": [], + "source": [ + "from copy import deepcopy\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from scipy.constants import c, epsilon_0\n", + "\n", + "from lasy.laser import Laser\n", + "from lasy.profiles import CombinedLongitudinalTransverseProfile\n", + "from lasy.profiles.longitudinal import ContinuousWaveProfile\n", + "from lasy.profiles.transverse import (\n", + " GaussianTransverseProfile,\n", + ")\n", + "from lasy.utils.mode_decomposition import (\n", + " hermite_gauss_composition,\n", + " hermite_gauss_decomposition,\n", + " laguerre_gauss_composition,\n", + " laguerre_gauss_decomposition,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "c6d0576b-83c1-475d-a998-48d6d8dfc311", + "metadata": {}, + "source": [ + "Then we initialise the laser. We will just use a CW laser, since we only care about the spatial profile." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d3f803cc-6ff0-4739-93ed-2c6f255e6c03", + "metadata": {}, + "outputs": [], + "source": [ + "peak_fluence = 1.0e4 # J/m^2\n", + "spot_size = 10e-3\n", + "wavelength = 800e-9\n", + "omega0 = 2 * np.pi * c / wavelength\n", + "pol = (1, 0)\n", + "\n", + "dimensions = \"xyt\" # Use cylindrical geometry\n", + "lo = (-2.0 * spot_size, -2.0 * spot_size, None) # Lower bounds of the simulation box\n", + "hi = (2.0 * spot_size, 2.0 * spot_size, None) # Upper bounds of the simulation box\n", + "num_points = (256, 256, 1) # Number of points in each dimension\n", + "\n", + "long_prof = ContinuousWaveProfile(wavelength)\n", + "tran_prof = GaussianTransverseProfile(spot_size)\n", + "laser_profile = CombinedLongitudinalTransverseProfile(\n", + " wavelength, pol, long_prof, tran_prof, peak_fluence=peak_fluence\n", + ")\n", + "\n", + "laser = Laser(dimensions, lo, hi, num_points, laser_profile)\n", + "\n", + "# Plot the xt laser profile\n", + "laser.show(envelope_type=\"intensity\") # In this case this is the fluence (CW beam)\n", + "\n", + "# Plot the transverse laser profile\n", + "intensity = epsilon_0 * c / 2 * np.abs(laser.grid.get_temporal_field()) ** 2\n", + "fluence = np.sum(intensity, axis=-1).T * laser.grid.dx[-1]\n", + "extent = (\n", + " laser.grid.axes[0].min() * 1e3,\n", + " laser.grid.axes[0].max() * 1e3,\n", + " laser.grid.axes[1].min() * 1e3,\n", + " laser.grid.axes[1].max() * 1e3,\n", + ")\n", + "plt.figure()\n", + "plt.imshow(fluence, extent=extent, cmap=\"Reds\")\n", + "plt.xlabel(\"x (mm)\")\n", + "plt.ylabel(\"y (mm)\")\n", + "plt.colorbar(label=r\"Fluence (J/m$^2$)\")" + ] + }, + { + "cell_type": "markdown", + "id": "d5af5820-4420-468a-9147-d41b03c8dabe", + "metadata": {}, + "source": [ + "We now define some random coefficients. These can be used to construct arbitrary beams as superpositions of Hermite- and Laguerre-Gaussian modes (see Denoised Laser tutorial for an example of expressing an experimental laser profile as a summation of Hermite-Gaussian modes). " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e09b8234-f8c7-43ea-add0-380757bb5671", + "metadata": {}, + "outputs": [], + "source": [ + "# Modal coefficients used to construct the beam\n", + "modes = {}\n", + "modes[(0, 0)] = 1 - 0.2j\n", + "modes[(0, 1)] = -0.4 - 0.2j\n", + "modes[(1, 0)] = 0.5 + 1j\n", + "modes[(2, 2)] = -0.05 + 0.1j" + ] + }, + { + "cell_type": "markdown", + "id": "f0c7a309-3eb9-48fc-bb20-b60c0d579d89", + "metadata": {}, + "source": [ + "We can use these mode coefficients to generate laser pulses. Here is shown first for the Hermite-Gaussian modes..." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7f272dc2-679c-46d4-9963-89e81abdd601", + "metadata": {}, + "outputs": [], + "source": [ + "HG_laser = deepcopy(laser) # Make a copy of the input grid\n", + "hermite_gauss_composition(\n", + " HG_laser, spot_size, spot_size, modes, skipAsymmetricModes=False\n", + ")\n", + "\n", + "HG_laser.show(envelope_type=\"intensity\") # In this case this is the fluence (CW beam)\n", + "\n", + "intensity = epsilon_0 * c / 2 * np.abs(HG_laser.grid.get_temporal_field()) ** 2\n", + "fluence = np.sum(intensity, axis=-1).T * HG_laser.grid.dx[-1]\n", + "extent = (\n", + " HG_laser.grid.axes[0].min() * 1e3,\n", + " HG_laser.grid.axes[0].max() * 1e3,\n", + " HG_laser.grid.axes[1].min() * 1e3,\n", + " HG_laser.grid.axes[1].max() * 1e3,\n", + ")\n", + "plt.figure()\n", + "plt.imshow(fluence, extent=extent, cmap=\"Reds\")\n", + "plt.xlabel(\"x (mm)\")\n", + "plt.ylabel(\"y (mm)\")\n", + "plt.colorbar(label=r\"Fluence (J/m$^2$)\")\n", + "\n", + "calculatedModes = hermite_gauss_decomposition(HG_laser, spot_size, spot_size, 3, 3)\n", + "print(\"Mode coefficients for HG beam (up to third order)\")\n", + "for calcModeKey in calculatedModes.keys():\n", + " print(\n", + " \"%i,%i : %.2f , %.2f i\"\n", + " % (\n", + " calcModeKey[0],\n", + " calcModeKey[1],\n", + " np.real(calculatedModes[calcModeKey]),\n", + " np.imag(calculatedModes[calcModeKey]),\n", + " )\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "d6d3eff0-c99c-47fe-8ff8-d94da63fe199", + "metadata": {}, + "source": [ + "... where we have also confirmed that the recovered mode coefficients match the input mode coefficients. Similarly, we could use a Laguerre-Gaussian modal basis (which is particularly relevant for beams with orbital angular momentum):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e253de6d-c805-4ed2-9356-3865fc785013", + "metadata": {}, + "outputs": [], + "source": [ + "LG_laser = deepcopy(laser) # Make a copy of the input grid\n", + "laguerre_gauss_composition(LG_laser, spot_size, modes, skipAsymmetricModes=False)\n", + "\n", + "LG_laser.show(envelope_type=\"intensity\") # In this case this is the fluence (CW beam)\n", + "\n", + "intensity = epsilon_0 * c / 2 * np.abs(LG_laser.grid.get_temporal_field()) ** 2\n", + "fluence = np.sum(intensity, axis=-1).T * LG_laser.grid.dx[-1]\n", + "extent = (\n", + " LG_laser.grid.axes[0].min() * 1e3,\n", + " LG_laser.grid.axes[0].max() * 1e3,\n", + " LG_laser.grid.axes[1].min() * 1e3,\n", + " LG_laser.grid.axes[1].max() * 1e3,\n", + ")\n", + "plt.figure()\n", + "plt.imshow(fluence, extent=extent, cmap=\"Reds\")\n", + "plt.xlabel(\"x (mm)\")\n", + "plt.ylabel(\"y (mm)\")\n", + "plt.colorbar(label=r\"Fluence (J/m$^2$)\")\n", + "\n", + "calculatedModes = laguerre_gauss_decomposition(LG_laser, spot_size, 3, 3)\n", + "print(\"Mode coefficients for HG beam (up to third order)\")\n", + "for calcModeKey in calculatedModes.keys():\n", + " print(\n", + " \"%i,%i : %.2f , %.2f i\"\n", + " % (\n", + " calcModeKey[0],\n", + " calcModeKey[1],\n", + " np.real(calculatedModes[calcModeKey]),\n", + " np.imag(calculatedModes[calcModeKey]),\n", + " )\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b97d22fd-48ca-4331-b324-ee2f4ea805f2", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.7" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/example_gerchberg_saxton_algo.py b/examples/example_gerchberg_saxton_algo.py deleted file mode 100644 index a4e214305..000000000 --- a/examples/example_gerchberg_saxton_algo.py +++ /dev/null @@ -1,180 +0,0 @@ -import copy - -import matplotlib.pyplot as plt -import numpy as np -from mpl_toolkits.axes_grid1 import make_axes_locatable - -from lasy.laser import Laser -from lasy.profiles.gaussian_profile import GaussianProfile -from lasy.utils.phase_retrieval import gerchberg_saxton_algo -from lasy.utils.zernike import zernike - -# DEFINE PHYSICAL PARAMETERS & CREATE LASER PROFILE -wavelength = 800e-9 -pol = (1, 0) -laser_energy = 1e-3 -w0 = 25e-6 -tau = 30e-15 -t_peak = 0 -pulseProfile = GaussianProfile(wavelength, pol, laser_energy, w0, tau, t_peak) - -# DEFNIE COMPUTATION PARAMETERS AND CREATE LASER OBJECT -dim = "xyt" -lo = (-75e-6, -75e-6, -50e-15) -hi = (75e-6, 75e-6, 50e-15) -npoints = (100, 100, 100) -laser = Laser(dim=dim, lo=lo, hi=hi, npoints=npoints, profile=pulseProfile) - -# CALCULATE THE REQUIRED PHASE ABERRATION -x = np.linspace(lo[0], hi[0], npoints[0]) -y = np.linspace(lo[1], hi[1], npoints[1]) -X, Y = np.meshgrid(x, y) -pupilRadius = 2 * w0 -phase = -0.2 * zernike(X, Y, (0, 0, pupilRadius), 3) - -R = np.sqrt(X**2 + Y**2) -phaseMask = np.ones_like(phase) -phaseMask[R > pupilRadius] = 0 - -# NOW ADD THE PHASE TO EACH SLICE OF THE FOCUS -phase3D = np.repeat(phase[:, :, np.newaxis], npoints[2], axis=2) -laser.grid.set_temporal_field( - np.abs(laser.grid.get_temporal_field()) * np.exp(1j * phase3D) -) - -# PROPAGATE THE FIELD FIELD FOWARDS AND BACKWARDS BY 1 MM -propDist = 2e-3 -laserForward = copy.deepcopy(laser) -laserForward.propagate(propDist) -laserBackward = copy.deepcopy(laser) -laserBackward.propagate(-propDist) - -# PERFORM GERCHBERG-SAXTON ALGORITHM TO RETRIEVE PHASE -phaseBackward, phaseForward, amp_error = gerchberg_saxton_algo( - laserBackward, - laserForward, - 2 * propDist, - condition="amplitude_error", - max_iterations=50, - amplitude_error=1e-6, - debug=True, -) - -# GET THE FIELD AND PLOT IT -fig, ax = plt.subplots(2, 5, figsize=(15, 5), tight_layout=True) -tIndx = 50 -extent = (lo[0] * 1e6, hi[0] * 1e6, lo[1] * 1e6, hi[1] * 1e6) - - -def addColorbar(im, ax, label=None): - """Create a colorbar and add it to the plot.""" - divider = make_axes_locatable(ax) - cax = divider.append_axes("right", size="5%", pad=0.05) - plt.colorbar(im, cax=cax) - cax.set_ylabel(label) - - -field = laser.grid.get_temporal_field() -im0 = ax[0, 0].imshow(np.abs(field[:, :, tIndx]) ** 2, extent=extent, cmap="PuRd") -addColorbar(im0, ax[0, 0], "Intensity (norm.)") -ax[0, 0].set_title("Inten. z = 0.0 mm") -ax[0, 0].set_xlabel("x ($\mu m$)") -ax[0, 0].set_ylabel("y ($\mu m$)") - - -im1 = ax[0, 1].imshow( - np.angle(field[:, :, tIndx]) * phaseMask, extent=extent, cmap="coolwarm" -) -addColorbar(im1, ax[0, 1], "Phase (rad.)") -ax[0, 1].set_title("Phase z = 0.0 mm") -ax[0, 1].set_xlabel("x ($\mu m$)") -ax[0, 1].set_ylabel("y ($\mu m$)") - - -laser_calc = copy.deepcopy(laserBackward) -laser_calc.grid.set_temporal_field( - np.abs(laser_calc.grid.get_temporal_field()) * np.exp(1j * phaseBackward) -) -laser_calc.propagate(propDist) - -field_calc = laser_calc.grid.get_temporal_field() -im2 = ax[1, 0].imshow( - np.abs(np.abs(field[:, :, tIndx]) ** 2 - np.abs(field_calc[:, :, tIndx]) ** 2), - extent=extent, - cmap="PuRd", -) -ax[1, 0].set_title("Inten. Res. z = 0.0 mm") -ax[1, 0].set_xlabel("x ($\mu m$)") -ax[1, 0].set_ylabel("y ($\mu m$)") -addColorbar(im2, ax[1, 0], "Intensity (norm.)") - -phaseResidual = np.angle(field_calc[:, :, tIndx]) - np.angle(field[:, :, tIndx]) -phaseResidual -= phaseResidual[int(npoints[1] / 2), int(npoints[0] / 2)] -maxPhaseRes = np.max(np.abs(phaseResidual) * phaseMask) -im3 = ax[1, 1].imshow( - phaseResidual * phaseMask, - extent=extent, - cmap="coolwarm", - vmin=-maxPhaseRes, - vmax=maxPhaseRes, -) -ax[1, 1].set_title("Phase Res. z = 0.0 mm") -ax[1, 1].set_xlabel("x ($\mu m$)") -ax[1, 1].set_ylabel("y ($\mu m$)") -addColorbar(im3, ax[1, 1], "Phase (rad.)") - -field_bw = laserBackward.grid.get_temporal_field() -im4 = ax[0, 2].imshow(np.abs(field_bw[:, :, tIndx]) ** 2, extent=extent, cmap="PuRd") -addColorbar(im4, ax[0, 2], "Intensity (norm.)") -ax[0, 2].set_title("Inten. z = %.1f mm" % (-propDist * 1e3)) -ax[0, 2].set_xlabel("x ($\mu m$)") -ax[0, 2].set_ylabel("y ($\mu m$)") - -im5 = ax[0, 3].imshow(np.angle(field_bw[:, :, tIndx]), extent=extent, cmap="coolwarm") -addColorbar(im5, ax[0, 3], "Phase (rad.)") -ax[0, 3].set_title("Phase z = %.1f mm" % (-propDist * 1e3)) -ax[0, 3].set_xlabel("x ($\mu m$)") -ax[0, 3].set_ylabel("y ($\mu m$)") - -phaseResidual = np.angle(field_bw[:, :, tIndx]) - phaseBackward[:, :, tIndx] -phaseResidual -= phaseResidual[int(npoints[1] / 2), int(npoints[0] / 2)] -maxPhaseRes = np.max(np.abs(phaseResidual) * phaseMask) -im6 = ax[0, 4].imshow( - phaseResidual * phaseMask, - extent=extent, - cmap="coolwarm", - vmin=-maxPhaseRes, - vmax=maxPhaseRes, -) -addColorbar(im6, ax[0, 4], "Phase (rad.)") -ax[0, 4].set_title("Phase Res. z = %.1f mm" % (-propDist * 1e3)) -ax[0, 4].set_xlabel("x ($\mu m$)") -ax[0, 4].set_ylabel("y ($\mu m$)") - -field_fw = laserForward.grid.get_temporal_field() -im7 = ax[1, 2].imshow(np.abs(field_fw[:, :, tIndx]) ** 2, extent=extent, cmap="PuRd") -addColorbar(im7, ax[1, 2], "Intensity (norm.)") -ax[1, 2].set_title("Inten. z = %.1f mm" % (propDist * 1e3)) -ax[1, 2].set_xlabel("x ($\mu m$)") -ax[1, 2].set_ylabel("y ($\mu m$)") -im8 = ax[1, 3].imshow(np.angle(field_fw[:, :, tIndx]), extent=extent, cmap="coolwarm") -addColorbar(im8, ax[1, 3], "Phase (rad.)") -ax[1, 3].set_title("Phase z = %.1f mm" % (propDist * 1e3)) -ax[1, 3].set_xlabel("x ($\mu m$)") -ax[1, 3].set_ylabel("y ($\mu m$)") - -phaseResidual = np.angle(field_fw[:, :, tIndx]) - phaseForward[:, :, tIndx] -phaseResidual -= phaseResidual[int(npoints[1] / 2), int(npoints[0] / 2)] -maxPhaseRes = np.max(np.abs(phaseResidual) * phaseMask) -im9 = ax[1, 4].imshow( - phaseResidual * phaseMask, - extent=extent, - cmap="coolwarm", - vmin=-maxPhaseRes, - vmax=maxPhaseRes, -) -addColorbar(im9, ax[1, 4], "Phase (rad.)") -ax[1, 4].set_title("Phase Res. z = %.1f mm" % (propDist * 1e3)) -ax[1, 4].set_xlabel("x ($\mu m$)") -ax[1, 4].set_ylabel("y ($\mu m$)") -plt.show() diff --git a/lasy/utils/denoise.py b/lasy/utils/denoise.py deleted file mode 100644 index 5fc404c5f..000000000 --- a/lasy/utils/denoise.py +++ /dev/null @@ -1,68 +0,0 @@ -from lasy.profiles.transverse.hermite_gaussian_profile import ( - HermiteGaussianTransverseProfile, -) -from lasy.profiles.transverse.transverse_profile import TransverseProfile -from lasy.utils.mode_decomposition import hermite_gauss_decomposition - - -def hg_reconstruction( - transverse_profile, - wavelength, - resolution, - lo, - hi, - n_modes_x=10, - n_modes_y=10, -): - """ - Denoise the transverse profile by decomposing it into a set of Hermite-Gaussian modes. - - The profiles are weighted according to mode coefficients and then added. - - Parameters - ---------- - transverse_profile : TransverseProfile object - An instance of a class or sub-class of TransverseProfile. - Defines the transverse envelope of the laser. - - wavelength : float (in meter) - Central wavelength at which the Hermite-Gauss beams are to be defined. - - resolution : float - The number of grid points in x and y used during the decomposition calculation. - - lo, hi : array of floats - The lower and upper bounds of the spatial grid on which the - decomposition is be performed. - - n_modes_x, n_modes_y : ints (optional) - The maximum values of `n_x` and `n_y` up to which the expansion is performed. - - Returns - ------- - transverse_profile_cleaned : class instance - Denoised transverse profile after decomposition and recombination. - - w0x, w0y : floats - Beam waist for which the decomposition is calculated. - It is computed as the waist for which the weight of order 0 is maximum. - """ - assert isinstance(transverse_profile, TransverseProfile) - - # Calculate the decomposition and waist of the laser pulse - modeCoeffs, w0x, w0y = hermite_gauss_decomposition( - transverse_profile, wavelength, resolution, lo, hi, n_modes_x, n_modes_y - ) - - # Denosing the laser profile - for i, mode_key in enumerate(list(modeCoeffs)): - transverse_profile_temp = HermiteGaussianTransverseProfile( - w0x, w0y, mode_key[0], mode_key[1], wavelength - ) # Create a new profile for each mode - - if i == 0: # First mode (0,0) - transverse_profile_cleaned = modeCoeffs[mode_key] * transverse_profile_temp - else: # All other modes - transverse_profile_cleaned += modeCoeffs[mode_key] * transverse_profile_temp - - return transverse_profile_cleaned, w0x, w0y diff --git a/lasy/utils/exp_data_utils.py b/lasy/utils/exp_data_utils.py index ef861c406..d5ca326f3 100644 --- a/lasy/utils/exp_data_utils.py +++ b/lasy/utils/exp_data_utils.py @@ -10,11 +10,14 @@ def find_center_of_mass(img): img: 2Darray of floats The image on which to calculate the COM + cal: float (m) + Pixel calibration in meters + Returns ------- x0 , y0: floats The center of mass of the image along the horizontal - and the vertical. The units are in pixels. + and the vertical. The units are in pixels unless cal is specified """ rows, cols = np.shape(img) diff --git a/lasy/utils/mode_decomposition.py b/lasy/utils/mode_decomposition.py index aa9d23412..c0d896982 100644 --- a/lasy/utils/mode_decomposition.py +++ b/lasy/utils/mode_decomposition.py @@ -1,112 +1,263 @@ -import math - import numpy as np from lasy.profiles.transverse.hermite_gaussian_profile import ( HermiteGaussianTransverseProfile, ) -from lasy.profiles.transverse.transverse_profile import TransverseProfile +from lasy.profiles.transverse.laguerre_gaussian_profile import ( + LaguerreGaussianTransverseProfile, +) from lasy.utils.exp_data_utils import find_d4sigma +def get_hermite_mode(grid_in, w0x, w0y, i, j): + """Calculate the projection of a field onto a particular mode. + + Projects a laser field onto a Hermite-Gaussian mode to + obtain the complex mode coefficient. + + Parameters + ---------- + grid_in : Grid + Grid object at the input plane. + + w0x : float (in m) + Spot size in the x-direction + + w0y : float (in m) + Spot size in the y-direction + + i : integer + Order of the x-direction mode + + j : integer + Order of the y-direction mode + + Returns + ------- + coeff : complex float + The projected complex modal coefficient for the (i,j) Hermite-Gaussian mode + """ + X, Y, T = np.meshgrid( + grid_in.grid.axes[0], grid_in.grid.axes[1], grid_in.grid.axes[2], indexing="ij" + ) + + dx = np.mean(np.diff(grid_in.grid.axes[0])) + dy = np.mean(np.diff(grid_in.grid.axes[1])) + + hg = HermiteGaussianTransverseProfile( + w0x, w0y, i, j, grid_in.profile.lambda0 + ).evaluate(X, Y) + + coeff = np.sum(grid_in.grid.get_temporal_field() * np.conj(hg)) * dx * dy + + return coeff + + def hermite_gauss_decomposition( - laserProfile, - wavelength, - res, - lo, - hi, - m_max=12, - n_max=12, + grid_in, w0x, w0y, Mmax, Nmax, skipAsymmetricModes=False ): - """ - Decomposes a laser profile into a set of hermite-gaussian modes. + """Decompose a laser field onto a Hermite-Gaussian basis. - The function takes an instance of `TransverseProfile`. + Loops through the mode coefficients, calculating the mode coefficient + for each one via projection. Returns a dictionary of the complex mode + coefficients. Parameters ---------- - laserProfile : class instance - An instance of a class or sub-class of TransverseLaserProfile + grid_in : Grid + Grid object at the input plane. - wavelength : float (in meter) - Central wavelength at which the Hermite-Gauss beams are to be defined. + w0x : float (in m) + Spot size in the x-direction - res : float - The resolution of grid points in x and y that will be used - during the decomposition calculation + w0y : float (in m) + Spot size in the y-direction - lo, hi : array of floats - The lower and upper bounds of the spatial grid on which the - decomposition will be performed. + Mmax : integer + Maximum order of the x-direction mode - m_max, n_max : ints - The maximum values of `m` and `n` up to which the expansion - will be performed + Nmax : integer + Maximum order of the y-direction mode + + skipAsymmetricModes : Boolean + Allows the user to only consider symmetric modal coefficients Returns ------- - 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 (`m`,`n`) + cxy : dict + A dictionary of complex modal coefficients + """ + cxy = {} + for i in range(Mmax): + for j in range(Nmax): + if (i != j) and (skipAsymmetricModes): + cxy[(i, j)] = 0 + continue + cxy[(i, j)] = get_hermite_mode(grid_in, w0x, w0y, i, j) - w0x, w0y : floats - Beam waist for which the decomposition is calculated. - It is computed as the waist for which the weight of order 0 is maximum. + return cxy + + +def hermite_gauss_composition(grid_in, w0x, w0y, cxy, skipAsymmetricModes=False): + """Compose a laser field from the mode coefficients. + + Uses a dictionary of complex modal coefficients to generate a beam + from a Hermite-Gaussian basis. Updates the laser object. + + Parameters + ---------- + grid_in : Grid + Grid object at the input plane. + + w0x : float (in m) + Spot size in the x-direction + + w0y : float (in m) + Spot size in the y-direction + + cxy : dict + A dictionary of complex modal coefficients + + skipAsymmetricModes : Boolean + Allows the user to only consider symmetric modal coefficients """ - # Check if the provided laserProfile is a transverse profile. - assert isinstance(laserProfile, TransverseProfile), ( - "laserProfile must be an instance of TransverseProfile" + X, Y, T = np.meshgrid( + grid_in.grid.axes[0], grid_in.grid.axes[1], grid_in.grid.axes[2], indexing="ij" ) + field = np.zeros(X.shape, dtype="complex128") + + for nxy in list(cxy): + if (nxy[0] != nxy[1]) and (skipAsymmetricModes): + continue + field += cxy[nxy] * ( + HermiteGaussianTransverseProfile( + w0x, w0y, nxy[0], nxy[1], grid_in.profile.lambda0 + ).evaluate(X, Y) + ) + + grid_in.grid.set_temporal_field(field) - # Get the field, sensible spatial bounds for the profile - lo0 = lo[0] + laserProfile.x_offset - lo1 = lo[1] + laserProfile.x_offset - hi0 = hi[0] + laserProfile.y_offset - hi1 = hi[1] + laserProfile.y_offset - Nx = int((hi0 - lo0) // (2 * res) * 2) + 2 - Ny = int((hi1 - lo1) // (2 * res) * 2) + 2 +def get_laguerre_mode(grid_in, w0, i, j): + """Calculate the projection of a field onto a particular mode. - # Define spatial arrays - x = np.linspace( - (lo0 + hi0) / 2 - (Nx - 1) / 2 * res, - (lo0 + hi0) / 2 + (Nx - 1) / 2 * res, - Nx, + Projects a laser field onto a Laguerre-Gaussian mode to + obtain the complex mode coefficient. + + Parameters + ---------- + grid_in : Grid + Grid object at the input plane. + + w0 : float (in m) + Spot size + + i : integer + Order of the radial mode + + j : integer + Order of the azimuthal mode + + Returns + ------- + coeff : complex float + The projected complex modal coefficient for the (i,j) Laguerre-Gaussian mode + """ + X, Y, T = np.meshgrid( + grid_in.grid.axes[0], grid_in.grid.axes[1], grid_in.grid.axes[2], indexing="ij" ) - y = np.linspace( - (lo1 + hi1) / 2 - (Ny - 1) / 2 * res, - (lo1 + hi1) / 2 + (Ny - 1) / 2 * res, - Ny, + + dx = np.mean(np.diff(grid_in.grid.axes[0])) + dy = np.mean(np.diff(grid_in.grid.axes[1])) + + lg = LaguerreGaussianTransverseProfile(w0, i, j, grid_in.profile.lambda0).evaluate( + X, Y ) - X, Y = np.meshgrid(x, y) - dx = x[1] - x[0] - dy = y[1] - y[0] - # Get the field on this grid - field = laserProfile.evaluate(X, Y) + coeff = np.sum(grid_in.grid.get_temporal_field() * np.conj(lg)) * dx * dy - # Get estimate of w0 - w0x, w0y = estimate_best_HG_waist(x, y, field, wavelength) + return coeff - # Next we loop over the modes and calculate the relevant weights - weights = {} - for m in range(m_max): - for n in range(n_max): - HGMode = HermiteGaussianTransverseProfile(w0x, w0y, m, n, wavelength) - coef = np.real( - np.sum(field * HGMode.evaluate(X, Y)) * dx * dy - ) # modalDecomposition - if math.isnan(coef): - coef = 0 - weights[(m, n)] = coef - return weights, w0x, w0y +def laguerre_gauss_decomposition(grid_in, w0, Mmax, Nmax, skipAsymmetricModes=False): + """Decompose a laser field onto a Laguerre-Gaussian basis. + Loops through the mode coefficients, calculating the mode coefficient + for each one via projection. Returns a dictionary of the complex mode + coefficients. -def estimate_best_HG_waist(x, y, field, wavelength): + Parameters + ---------- + grid_in : Grid + Grid object at the input plane. + + w0 : float (in m) + Spot size + + Mmax : integer + Maximum order of the radial mode + + Nmax : integer + Maximum order of the azimuthal mode + + skipAsymmetricModes : Boolean + Allows the user to only consider symmetric modal coefficients + + Returns + ------- + cxy : dict + A dictionary of complex modal coefficients + """ + cxy = {} + for i in range(Mmax): + for j in range(Nmax): + if (i != j) and (skipAsymmetricModes): + cxy[(i, j)] = 0 + continue + cxy[(i, j)] = get_laguerre_mode(grid_in, w0, i, j) + + return cxy + + +def laguerre_gauss_composition(grid_in, w0, cxy, skipAsymmetricModes=False): + """Compose a laser field from the mode coefficients. + + Uses a dictionary of complex modal coefficients to generate a beam + from a Laguerre-Gaussian basis. Updates the laser object. + + Parameters + ---------- + grid_in : Grid + Grid object at the input plane. + + w0 : float (in m) + Spot size + + cxy : dict + A dictionary of complex modal coefficients + + skipAsymmetricModes : Boolean + Allows the user to only consider symmetric modal coefficients """ - Estimate the waist that maximises the weighting of the first mode. + X, Y, T = np.meshgrid( + grid_in.grid.axes[0], grid_in.grid.axes[1], grid_in.grid.axes[2], indexing="ij" + ) + field = np.zeros(X.shape, dtype="complex128") + + for nxy in list(cxy): + if (nxy[0] != nxy[1]) and (skipAsymmetricModes): + continue + field += cxy[nxy] * ( + LaguerreGaussianTransverseProfile( + w0, nxy[0], nxy[1], grid_in.profile.lambda0 + ).evaluate(X, Y) + ) + + grid_in.grid.set_temporal_field(field) + + +def estimate_best_HG_waist(x, y, field, wavelength): + """Estimate the waist that maximises the weighting of the first mode. Calculates a D4Sigma waist as a first estimate and then tests multiple gaussians with waists around this value to determine which has the best diff --git a/tests/test_denoise.py b/tests/test_denoise.py index 2f59af1d8..15f258186 100644 --- a/tests/test_denoise.py +++ b/tests/test_denoise.py @@ -6,35 +6,65 @@ value. """ +from copy import deepcopy + import numpy as np +from lasy.laser import Laser +from lasy.profiles.combined_profile import CombinedLongitudinalTransverseProfile +from lasy.profiles.longitudinal import ContinuousWaveProfile from lasy.profiles.transverse.super_gaussian_profile import ( SuperGaussianTransverseProfile, ) -from lasy.utils.denoise import hg_reconstruction +from lasy.utils.mode_decomposition import ( + hermite_gauss_composition, + hermite_gauss_decomposition, +) def test_denoise_hg_reconstruction(): # Parameters + dimensions = "xyt" + pol = (0, 1) + peak_fluence = 1.0 + lambda0 = 800e-9 w0 = 20e-6 shape_parameter = 3 - wavelength = 8e-7 - resolution = 0.2e-6 - lo = [-2e-4, -2e-4] - hi = [2e-4, 2e-4] + lo = [-2e-4, -2e-4, None] + hi = [2e-4, 2e-4, None] + num_points = [1000, 1000, 1] # Define the transverse profile + longitudinal_profile = ContinuousWaveProfile(lambda0) transverse_profile = SuperGaussianTransverseProfile( w0, shape_parameter ) # Super-Gaussian profile - transverse_profile_cleaned, w0x, w0y = hg_reconstruction( - transverse_profile, wavelength, resolution, lo, hi - ) # Denoised profile + laser_profile = CombinedLongitudinalTransverseProfile( + lambda0, + pol, + longitudinal_profile, + transverse_profile, + peak_fluence=peak_fluence, + ) + + laser_raw = Laser(dimensions, lo, hi, num_points, laser_profile) + + # Calculate the decomposition and waist of the laser pulse + modes = hermite_gauss_decomposition(laser_raw, w0, w0, 10, 10) + + laser_cleaned = deepcopy(laser_raw) # Make a copy of the input grid + hermite_gauss_composition(laser_cleaned, w0, w0, modes) # Calculate the error - x = np.linspace(-5 * w0x, 5 * w0x, 500) + x = np.linspace(-5 * w0, 5 * w0, 500) X, Y = np.meshgrid(x, x) - prof1 = transverse_profile.evaluate(X, Y) - prof2 = transverse_profile_cleaned.evaluate(X, Y) + + # Original profile + prof1 = laser_raw.grid.get_temporal_field() + + # Reconstructed profile + prof2 = laser_cleaned.grid.get_temporal_field() + error = np.sum(np.abs(prof2 - prof1) ** 2) / np.sum(np.abs(prof1) ** 2) + print(error) assert error < 0.02