ExamplesRefraction.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import copy\n",
    "np.set_printoptions(precision=2)\n",
    "from pathlib import Path\n",
    "\n",
    "import wave as w\n",
    "import constants as c\n",
    "import material as m\n",
    "import style"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example usage of simulated materials\n",
    "This notebook is ment as an reference on how to use `PlaneWave` and `Material` classes,\n",
    "and for visual test of the implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figures_dir = \"figures/\"\n",
    "Path(figures_dir).mkdir(parents=True, exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate an example spectrum\n",
    "Just an example how to create and plot wave."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = w.GaussianPlanarWave(mean=c.GREEN, std=(5*c.MICRO))\n",
    "s_replay = w.WhitePlanarWave()\n",
    "z = np.linspace(-6*c.MICRO, 6*c.MICRO, 1000)\n",
    "s.delay(0)\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))\n",
    "s.plot(ax=ax1, spectrum_axis=ax1.twinx(), wavelength=True)\n",
    "s.plot_amplitude(z, ax=ax2)\n",
    "ax2.legend()\n",
    "fig.tight_layout()\n",
    "plt.savefig(figures_dir + \"spectrum.pdf\", transparent=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare different materials\n",
    "\n",
    "#### Analitical vs constant vs matrix theory\n",
    "Create the same material (empty space) using three different material models in oder to inspect that they give the same resut."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air = m.EmptySpace(z=z, name=\"analytic\")\n",
    "common_n = 1\n",
    "fixed_dielectric = m.ConstantMaterial(z=z, n0=common_n, name=\"phase shift\")\n",
    "dielectric = m.SimpleDielectric(z=z, n0=1, name=\"transfer matrices\")\n",
    "\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 5))\n",
    "\n",
    "air.plot(ax1, imaginary_axis=ax1.twinx())\n",
    "fixed_dielectric.plot(ax2, imaginary_axis=ax2.twinx())\n",
    "dielectric.plot(ax3, imaginary_axis=ax3.twinx())\n",
    "fig.tight_layout() \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate interference pattern in those three materials: again we can see that the patterns match. \n",
    "The index of refraction after recording would be different in those materials because ConstantMaterial's index of refraction cannot be modified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fig, (ax1) = plt.subplots(1, 1, figsize=(8,4))\n",
    "\n",
    "air.record(s,s)\n",
    "air.plot_recent_energy(ax1, label=\"air\")\n",
    "fixed_dielectric.record(s,s)\n",
    "fixed_dielectric.plot_recent_energy(ax1, label=\"correlation\")\n",
    "dielectric.record(s,s)\n",
    "dielectric.plot_recent_energy(ax1, label=\"matrix theory\")\n",
    "ax1.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### An example recording\n",
    "\n",
    "Left: interference pattern \n",
    "created in the material. Middle: index of \n",
    "refraction after recording. Right: reflected pulse, when the \n",
    "incoming pulse was the same as the recorded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = np.linspace(-6*c.MICRO, 6*c.MICRO, 1000)\n",
    "dielectric = m.SimpleDielectric(z=z, n0=1.45, name=\"transfer matrices\")\n",
    "s0 = w.GaussianPlanarWave(mean=c.GREEN, std=(5*c.MICRO))\n",
    "dielectric.record(s0, s0)\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(9, 3))\n",
    "dielectric.plot_recent_energy(ax1)\n",
    "ax1.set_title(\"interference\")\n",
    "dielectric.plot(ax=ax2)\n",
    "dielectric.reflect(s_replay).plot(ax3, wavelength=True)\n",
    "plt.tight_layout()\n",
    "plt.savefig(figures_dir + \"recording.pdf\", transparent=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Layered material\n",
    "\n",
    "Left: interference \n",
    "pattern \n",
    "created in the material. Middle: index of refraction after recording. Right: \n",
    "reflected pulse. All inner implementations are the same as for dielectric.\n",
    "Additionally, on the first two plots shading is \n",
    "added to indicate layers of different index of \n",
    "refraction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_layered = np.linspace(-2*c.MICRO, 2*c.MICRO, 2000)\n",
    "layered = m.LayeredMaterial(z=z_layered, n_layers=20, layer_width=80*c.NANO, n0=1.5, n1=2.0)\n",
    "\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(9, 3))\n",
    "layered.shade_plot(ax1, alpha=0.2)\n",
    "layered.plot(ax1)\n",
    "layered.shade_plot(ax2, alpha=0.2)\n",
    "layered.reflect(s_replay).plot(ax3, label=\"reflected\", wavelength=True)\n",
    "layered.record(s, s)\n",
    "layered.plot_recent_energy(ax2)\n",
    "fig.tight_layout() \n",
    "plt.savefig(figures_dir + \"layered.pdf\",transparent=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Composite Material\n",
    "\n",
    "Left: \n",
    "interference pattern \n",
    "created in the material. Middle: index of refraction after recording. Right: \n",
    "reflected pulse. All inner implementations are the same as for dielectric. Again, the shading is \n",
    "added to indicate different components. The strange reflected wave shape comes mostly from interfaces with air which haven't been modeled before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_layered = np.linspace(-1*c.MICRO, 1*c.MICRO, 2000)\n",
    "z_air = np.linspace(-0.5*c.MICRO, 0, 500)\n",
    "z_clean = np.linspace(-1*c.MICRO, 0, 1000)\n",
    "air = m.EmptySpace(z=z_air)\n",
    "air2 = m.EmptySpace(z=np.copy(z_air))\n",
    "clean = m.SimpleDielectric(z=z_clean, n0=1.5)\n",
    "layered = m.LayeredMaterial(z=z_layered, n_layers=10, layer_width=100*c.NANO/1.5, n0=1.5, n1=2.0)\n",
    "composite = m.CompositeMaterial([air, layered, clean, air2])\n",
    "\n",
    "\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(9, 3))\n",
    "composite.shade_plot(ax1, alpha=0.2)\n",
    "composite.plot(ax1, color=\"C0\")\n",
    "composite.shade_plot(ax2, alpha=0.2)\n",
    "composite.reflect(s_replay).plot(ax3, label=\"reflected\", wavelength=True)\n",
    "composite.record(s, s)\n",
    "composite.plot_recent_energy(ax2, color=\"C0\")\n",
    "fig.tight_layout() \n",
    "plt.savefig(figures_dir + \"composite.pdf\",transparent=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random phase example\n",
    "\n",
    "Random phase. One pulse is kept the same, and the phase of the second \n",
    "one is choose uniformly at random form interval $(0, \\pi)$. On the top left \n",
    "plot we can see that the interference patter indeed get blurred, and on the top \n",
    "right we can see that there is indeed a Gaussian shape emerging at the basis of \n",
    "the interference pattern. Fortunately, the patter still has a lot of variance at \n",
    "the right frequency, and thus with each pulse more gets reflected, as depicted \n",
    "at the bottom plots. We can see that there gain is not the same with every \n",
    "pulse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "dielectric = m.SimpleDielectric(z=z, n0=1)\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=(9, 3))\n",
    "pulse = s0\n",
    "\n",
    "start = time.time()\n",
    "for pulse, color in enumerate(style.colors):\n",
    "    phase = np.random.uniform(0, np.pi)\n",
    "    s2 = s0 * np.exp(1j*phase)\n",
    "    dielectric.plot(ax=ax2, color=color, change_only=True)\n",
    "    dielectric.record(s0, s2)\n",
    "    dielectric.plot_recent_energy(ax1, label=\"pulse {}\".format(pulse + 1), c=color)\n",
    "    dielectric.reflect(s_replay).plot(ax3, color=color)\n",
    "\n",
    "ax1.set_title(\"interference\")\n",
    "ax3.set_title(\"reflected\")\n",
    "ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))\n",
    "\n",
    "print(\"{} recordings took {}s\".format(steps, time.time()-start))\n",
    "ax3.legend([f\"p. {nr + 1}\" for nr in range(steps)])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig(figures_dir + \"random_phase.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sweeping phase example\n",
    "\n",
    "Random phase. One pulse is kept the same, and the phase of the second \n",
    "one is changing by $\\pi/20$ with each pulse. On the top left \n",
    "plot we can see that the interference patter shifts to the right, and on the top \n",
    "right we can see the refractive index pattern is much sharper than in the previous example. The total reflection is however very similar. It seems \n",
    "to reach saturation faster, and be slightly sharper in wavelength."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dielectric = m.SimpleDielectric(z=z, n0=1)\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=(9, 3))\n",
    "\n",
    "for pulse, color in enumerate(style.colors):\n",
    "    phase = 0.5*np.pi*pulse/steps\n",
    "    s2 = s0 * np.exp(1j*phase)\n",
    "    dielectric.plot(ax=ax2, color=color, change_only=True)\n",
    "    dielectric.record(s0, s2)\n",
    "    dielectric.plot_recent_energy(ax1, label=\"pulse {}\".format(pulse + 1), c=color)\n",
    "    dielectric.reflect(s_replay).plot(ax3, color=color)\n",
    "    \n",
    "ax1.set_title(\"interference\")\n",
    "ax3.set_title(\"reflected\")\n",
    "ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))\n",
    "ax3.legend([f\"p. {nr + 1}\" for nr in range(steps)])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig(figures_dir + \"sweeping_phase.pdf\")      \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "dielectric = m.SimpleDielectric(z=z, n0=1)\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=(9, 3))\n",
    "pulse = s0\n",
    "interference = dielectric.deposited_energy\n",
    "\n",
    "start = time.time()\n",
    "for pulse, color in enumerate(style.colors):\n",
    "    phase = np.random.uniform(0, 2*np.pi)\n",
    "    s2 = s0 * np.exp(1j*phase)\n",
    "    dielectric.plot(ax=ax2, color=color, change_only=True)\n",
    "    dummy_dielectric =  m.SimpleDielectric(z=z, n0=1)\n",
    "    dummy_dielectric.record(s0, s2)\n",
    "    dielectric.deposited_energy += dielectric.energy_response(dummy_dielectric.recent_energy)\n",
    "    dummy_dielectric.plot_recent_energy(ax1, label=\"pulse {}\".format(pulse + 1), c=color)\n",
    "    dielectric.reflect(s_replay).plot(ax3, color=color)\n",
    "    \n",
    "print(\"{} recordings took {}s\".format(steps, time.time()-start))\n",
    "\n",
    "ax1.set_title(\"interference\")\n",
    "ax3.set_title(\"reflected\")\n",
    "ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))\n",
    "ax3.legend([f\"p. {nr + 1}\" for nr in range(steps)])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig(figures_dir + \"random_phase_worst_case.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Self organisation example\n",
    "**Warning: this example is not realistic.** It is a proof of concept that with small inintial modification (made by interference) much bigger modification can be build with *just one* pulse of the right frequency. The parameters are choosen for this simulation to give desired resutls and are not realistic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dielectric = m.SimpleDielectric(z=z, n0=0.5)\n",
    "dielectric.max_dn = 0.5\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=(9, 3))\n",
    "pulse = s0\n",
    "\n",
    "start = time.time()\n",
    "for pulse, color in enumerate(style.colors):\n",
    "    phase = np.random.uniform(0, np.pi)\n",
    "    if pulse > 0:\n",
    "        s2 = s0 * 0\n",
    "        s1 = 1.4 * s0\n",
    "    else:\n",
    "        s2 = 0.8 * s0\n",
    "        s1 = 0.8 * s0\n",
    "    dielectric.record(s1, s2)\n",
    "    dielectric.plot(ax=ax2, color=color, change_only=True)\n",
    "    if pulse > 0:\n",
    "        dielectric.plot_recent_energy(ax1, c=color)\n",
    "    dielectric.reflect(s_replay).plot(ax3, color=color)\n",
    "    \n",
    "ax1.set_title(\"interference\")\n",
    "ax3.set_title(\"reflected\")\n",
    "ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))\n",
    "ax3.legend([f\"p. {nr + 1}\" for nr in range(steps)])\n",
    "    \n",
    "plt.tight_layout()\n",
    "plt.savefig(figures_dir + \"self_organisation.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Chirps example\n",
    "\n",
    "#### Different chirps interfering with their copies. \n",
    "Left: chirping Gaussian spectrum leads to spreading of the pulse amplitude in \n",
    "time (and thus in space). Middle and right: the recorded patterns and reflected \n",
    "spectra look the same across the experiments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = np.linspace(-6*c.MICRO, 6*c.MICRO, 2000)\n",
    "probing_wave = w.GaussianPlanarWave(mean=c.GREEN, std=(5*c.MICRO))\n",
    "for chirp in np.linspace(0, 0.002*c.NANO, 3):\n",
    "    dielectric = m.SimpleDielectric(z=z, n0=1.45, name=\"transfer matrices\")\n",
    "    chirped_wave = w.ChirpedPlanarWave(mean=c.GREEN, std=(5*c.MICRO), skew=chirp)\n",
    "    chirped_wave.shift(chirp*2000)\n",
    "    dielectric.record(chirped_wave, chirped_wave)\n",
    "    fig, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4, figsize=(12, 3))\n",
    "    chirped_wave.plot_amplitude(z, ax=ax1)\n",
    "    chirped_wave.plot_amplitude(z, ax=ax0)\n",
    "    ax1.set_title(\"no-chirp\" if chirp == 0 else \"chirp\")\n",
    "    ax0.set_title(\"no-chirp\" if chirp == 0 else \"chirp\")\n",
    "    dielectric.plot(ax=ax2, change_only=True)\n",
    "    dielectric.reflect(s_replay).plot(ax3)\n",
    "    ax2.ticklabel_format(axis='y', style=\"sci\", scilimits=(0,0))\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(figures_dir + \"chirp_only_{:2.0e}.pdf\".format(chirp))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Chirps interfering with a non chirped pulse.\n",
    "\n",
    "Left: interference \n",
    "patterns look change if only one of the pulses is chirped. Right and middle:\n",
    "the recorded pattern and the reflected spectrum differ from plot to plot. \n",
    "It is hard to tell how significant this change is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = np.linspace(-6*c.MICRO, 6*c.MICRO, 2000)\n",
    "probing_wave = w.GaussianPlanarWave(mean=c.GREEN, std=(5*c.MICRO))\n",
    "for chirp in np.linspace(0, 0.002*c.NANO, 3):\n",
    "    dielectric = m.SimpleDielectric(z=z, n0=1.45, name=\"transfer matrices\")\n",
    "    chirped_wave = w.ChirpedPlanarWave(mean=c.GREEN, std=(5*c.MICRO), skew=chirp)\n",
    "    chirped_wave.shift(chirp*2000)\n",
    "    dielectric.record(chirped_wave, probing_wave)\n",
    "    fig, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4,  figsize=(12, 3))\n",
    "    chirped_wave.plot_amplitude(z, ax=ax1)\n",
    "    probing_wave.plot_amplitude(z, ax=ax0)\n",
    "    ax1.set_title(\"no-chirp\" if chirp == 0 else \"chirp\")\n",
    "    ax0.set_title(\"no-chirp\")\n",
    "    dielectric.plot(ax=ax2, change_only=True)\n",
    "    dielectric.reflect(s_replay).plot(ax3)\n",
    "    ax2.ticklabel_format(axis='y', useOffset=False, style=\"sci\", scilimits=(0,0))\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(figures_dir + \"chirp_mixed_{:2.0e}.pdf\".format(chirp))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gratings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "period = 230*c.NANO\n",
    "N = 20\n",
    "Z_max = N*period/2\n",
    "z_layered = np.linspace(-Z_max, Z_max, 2000)\n",
    "\n",
    "z_clean = np.linspace(-0.5*c.MICRO, 0, 3)\n",
    "clean = m.SimpleDielectric(z=z_clean, n0=1.5)\n",
    "clean2 = m.SimpleDielectric(z=z_clean, n0=1.5)\n",
    "\n",
    "layered = m.LayeredMaterial(z=z_layered, n_layers=N, layer_width=0.5*period, n0=1.5, n1=1.51)\n",
    "composite = m.CompositeMaterial([clean, layered, clean2])\n",
    "\n",
    "fig, axes = plt.subplots(2, 2, figsize=(10, 5))\n",
    "composite.plot(axes[0,0], color=\"C0\")\n",
    "composite.reflect(s_replay).plot(axes[0,1], label=\"reflected\", wavelength=True)\n",
    "axes[0,1].set_ylim(0, 0.02)\n",
    "\n",
    "layered = m.LayeredMaterial(z=z_layered, n_layers=N, layer_width=0.5*period, n0=1.505, n1=1.51)\n",
    "composite = m.CompositeMaterial([clean, layered, clean2])\n",
    "\n",
    "composite.plot(axes[1,0], color=\"C0\")\n",
    "composite.reflect(s_replay).plot(axes[1,1], label=\"reflected\", wavelength=True)\n",
    "axes[1,1].set_ylim(0, 0.02)\n",
    "\n",
    "fig.tight_layout() \n",
    "plt.savefig(figures_dir + \"grating.pdf\",transparent=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}