{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "393e0b32",
   "metadata": {},
   "source": [
    "# Nonlinear propagation using a split-step approach"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c87e770c",
   "metadata": {},
   "source": [
    "In this tutorial, we'll try to simulate the propagation of an intense laser pulse through a piece of fused silica glass, and observe the influence of the dispersion and self phase modulation that the pulse experiences in the material.\n",
    "\n",
    "This should give an overview of the functionality of the `AngularSpectrumDFFTPropagator`, the nonlinear phase applied by the `NonlinearKerrStep` and how do combine them in a simple, iterative split-step method.\n",
    "\n",
    "First, we'll import a few relevant functions and classes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35c997c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.constants import c\n",
    "\n",
    "from lasy.laser import Laser\n",
    "from lasy.profiles import GaussianProfile\n",
    "from lasy.propagators import AngularSpectrumPropagator, NonlinearKerrStep\n",
    "from lasy.utils.laser_utils import get_dispersion, get_spectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34232ba1",
   "metadata": {},
   "source": [
    "Next, we can define key parameters of the laser pulse and use them to create a Gaussian laser profile object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1bc884f",
   "metadata": {},
   "outputs": [],
   "source": [
    "wavelength = 800e-9  # wavelength in meters\n",
    "laser_energy = 0.5e-3  # laser pulse energy in Joule\n",
    "tau = 30e-15  # pulse duration in seconds\n",
    "w0 = 1e-3  # laser waist in meters\n",
    "polarisation = (0, 1)  # polarisation vector of the laser\n",
    "t_peak = 0  # laser peak position in time\n",
    "\n",
    "profile = GaussianProfile(\n",
    "    wavelength=wavelength,\n",
    "    laser_energy=laser_energy,\n",
    "    tau=tau,\n",
    "    w0=w0,\n",
    "    pol=polarisation,\n",
    "    t_peak=t_peak,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d6dacd3",
   "metadata": {},
   "source": [
    "To create the full laser object, we now define the parameters of the simulation grid, and combine the grid and profile to a laser object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0135763",
   "metadata": {},
   "outputs": [],
   "source": [
    "dim = \"xyt\"\n",
    "hi = (2.5e-3, 2.5e-3, 500e-15)\n",
    "lo = (-2.5e-3, -2.5e-3, -500e-15)\n",
    "npoints = (100, 100, 200)\n",
    "\n",
    "laser = Laser(dim=dim, hi=hi, lo=lo, npoints=npoints, profile=profile)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4bf0d5b",
   "metadata": {},
   "source": [
    "To check that the pulse looks as expected, we can show its profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44330e59",
   "metadata": {},
   "outputs": [],
   "source": [
    "laser.show(envelope_type=\"intensity\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4251a970",
   "metadata": {},
   "source": [
    "To see the influence of the self-phase modulation on the spectral properties of the laser, we can extract the initial spectrum and GDD curve of the laser pulse. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a964250",
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_spectrum, omega = get_spectrum(\n",
    "    grid=laser.grid, dim=laser.dim, omega0=laser.profile.omega0\n",
    ")\n",
    "\n",
    "initial_gdd, initial_gdd0 = get_dispersion(\n",
    "    grid=laser.grid, dim=laser.dim, omega0=laser.profile.omega0, order=2\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f25d99f4",
   "metadata": {},
   "source": [
    "Now let's define the properties of the fused silica material that the we want the pulse to propagate in. We will need the refractive index (ideally in form of the Sellmeier equation) and the intensity dependent refractive index, $n_2$. Both taken from [refractiveindex.info](https://www.refractiveindex.info).\n",
    "\n",
    "At this point we can also define the thickness of the fused silica window that we want to propagate through. In this example we'll use a 5mm thick window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b0dbdf1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# refractive index\n",
    "def n_fusedsilica(wavelength):\n",
    "    \"\"\"Sellmeier equation for fused silica.\"\"\"\n",
    "    x = wavelength * 1e6  #  convert wavelength to µm\n",
    "\n",
    "    return np.sqrt(\n",
    "        1\n",
    "        + 0.6961663 / (1 - (0.0684043 / x) ** 2)\n",
    "        + 0.4079426 / (1 - (0.1162414 / x) ** 2)\n",
    "        + 0.8974794 / (1 - (9.896161 / x) ** 2)\n",
    "    )\n",
    "\n",
    "\n",
    "# intensity dependent refractive index\n",
    "n2 = 2.7e-20\n",
    "\n",
    "# material thickness\n",
    "d_FS = 5e-3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b288ed6",
   "metadata": {},
   "source": [
    "A very common and proven approach for accurately simulating the nonlinear propagation of laser pulses is the so-called split-step method. In this method, the linear propagation, that entails effects such as dispersion and diffraction, and the non-linear propagation (containing the intensity dependent effects) into two different propagation steps. \n",
    "\n",
    "For our purpose, we can use the `AngularSpectrumDFFTPropagator` that describes the linear propagation in dispersive media, and the `NonlinearKerrStep` for the nonlinear propagation step. We will create them in the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51833efa",
   "metadata": {},
   "outputs": [],
   "source": [
    "linear_propagator = AngularSpectrumPropagator(\n",
    "    omega0=profile.omega0, n=n_fusedsilica, dim=laser.dim\n",
    ")\n",
    "nonlinear_phase = NonlinearKerrStep(n2=n2, k0=laser.profile.omega0 / c)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a138dcb5",
   "metadata": {},
   "source": [
    "To accurately simulate the interplay of the two propagation steps, the split-step approach devides total propagation distance into multiple sub-steps and alternately applies the two propagation steps. In general, the simulation result will become more accurate the smaller the propagation steps. In this specific case 10 steps with a size of $\\frac{d_{FS}}{N_{steps}}=\\frac{5mm}{10}=0.5mm$ seems like a good compromise between accurate results and a reasonably short simulation time.\n",
    "\n",
    "Lets do this iterative propagation next:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5cee4f82",
   "metadata": {},
   "outputs": [],
   "source": [
    "Nsteps = 10  # number of propagation steps\n",
    "dz = d_FS / Nsteps  # length of the individual propagation steps\n",
    "\n",
    "laser.add_propagator(linear_propagator)\n",
    "\n",
    "for i in range(Nsteps):\n",
    "    # linear propagation 1/2 step\n",
    "    laser.propagate(dz / 2)\n",
    "\n",
    "    # nonlinear propagation step\n",
    "    laser.grid = nonlinear_phase.apply(grid_in=laser.grid, distance=dz)\n",
    "\n",
    "    # linear propagation 1/2 step\n",
    "    laser.propagate(dz / 2)\n",
    "    # print progress\n",
    "    print(f\"Progress: {100 * (i + 1) / Nsteps:.1f}%\", end=\"\\r\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98788c5a",
   "metadata": {},
   "source": [
    "If we now again show the laser profile we can already see that the temporal shape of the pulse has changed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9fa237d",
   "metadata": {},
   "outputs": [],
   "source": [
    "laser.show(envelope_type=\"intensity\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70f71a92",
   "metadata": {},
   "source": [
    "To see the effects on the spectral properties, we can again extract the spectrum and GDD curve to compare to the initial ones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74bffd47",
   "metadata": {},
   "outputs": [],
   "source": [
    "output_spectrum, omega = get_spectrum(\n",
    "    grid=laser.grid, dim=laser.dim, omega0=laser.profile.omega0\n",
    ")\n",
    "\n",
    "output_gdd, output_gdd0 = get_dispersion(\n",
    "    grid=laser.grid, dim=laser.dim, omega0=laser.profile.omega0, order=2\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d97b198",
   "metadata": {},
   "source": [
    "Now lets plot them for comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d47001b",
   "metadata": {},
   "outputs": [],
   "source": [
    "wavelength_nm = 1e9 * 2 * np.pi * c / omega\n",
    "\n",
    "# spectrum\n",
    "plt.plot(wavelength_nm, initial_spectrum, label=\"Input spectrum\")\n",
    "plt.plot(wavelength_nm, output_spectrum, label=\"Output spectrum\")\n",
    "\n",
    "plt.ylabel(\"Spectral intensity\")\n",
    "plt.xlabel(\"Wavelength [nm]\")\n",
    "\n",
    "plt.ylim(\n",
    "    0,\n",
    ")\n",
    "plt.xlim(670, 950)\n",
    "\n",
    "plt.legend(loc=\"upper left\")\n",
    "plt.twinx()\n",
    "\n",
    "# GDD\n",
    "plt.plot(wavelength_nm, initial_gdd * 1e30, label=\"Input GDD\", linestyle=\"--\")\n",
    "plt.plot(wavelength_nm, output_gdd * 1e30, label=\"Output GDD\", linestyle=\"--\")\n",
    "\n",
    "plt.ylabel(\"GDD [fs$^2$]\")\n",
    "plt.xlabel(\"Wavelength [nm]\")\n",
    "\n",
    "plt.ylim(-5000, 5000)\n",
    "\n",
    "plt.legend(loc=\"upper right\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "533a174a",
   "metadata": {},
   "source": [
    "We can clearly see how the spectrum broadens and how the phase is affected by the intensity dependent refractive index and the resulting self-phase modulation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b5f65ef",
   "metadata": {},
   "source": [
    "Finally, we can let the laser propagate in vacuum for some more distance to show the influence of the self-focusing that the pulse experiences in the fused silica. The beam size decreases and the spatial profile is deformed from the ideal Gaussian that it is intially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ced6269",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define propagator for vacuum\n",
    "linear_propagator_vacuum = AngularSpectrumPropagator(\n",
    "    omega0=profile.omega0, n=1.0, dim=laser.dim\n",
    ")\n",
    "laser.add_propagator(linear_propagator_vacuum)\n",
    "\n",
    "# propagate in vacuum\n",
    "laser.propagate(0.2)\n",
    "\n",
    "# show laser after propagation\n",
    "laser.show(envelope_type=\"intensity\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "work",
   "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}