{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# This file is part of nannos\n# License: GPLv3\n%matplotlib notebook"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Stress tensor\n\nOptical forces on bulk and thin layers of dielectric and metal.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\nfrom scipy.constants import c, e, h\n\nimport nannos as nn"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Results are compared to the reference\n:cite:p:`Antonoyiannakis1999`.\nFirst we define the material dielectric functions\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.close(\"all\")\nplt.ion()\n\n\nnwl = 81\nwl, n = np.loadtxt(\"GaP_Aspnes.csv\", skiprows=1, delimiter=\",\", max_rows=nwl - 1).T\nwl1, k = np.loadtxt(\n    \"GaP_Aspnes.csv\", skiprows=nwl + 2, delimiter=\",\", max_rows=nwl - 1\n).T\nassert np.all(wl1 == wl)\n\n\ndef epsilon_GaP(wli):\n    return (np.interp(wli, wl, n) + 1j * np.interp(wli, wl, k)) ** 2\n\n\ndef epsilon_Al(wli):\n    f = h * c / e / wli * 1e6\n    f_p = 15\n    gamma = 0.1\n    return 1 - f_p**2 / (f * (f + 1j * gamma))\n\n\nwls = np.linspace(wl[0], wl[-1], 500)\neps_Al = epsilon_Al(wls)\neps_GaP = epsilon_GaP(wls)\nfs = h * c / e / wls * 1e6\n\nplt.figure()\nplt.plot(fs, eps_Al.real, label=\"Re Ag\", c=\"#5ba865\")\nplt.plot(fs, eps_Al.imag, \"--\", label=\"Im Ag\", c=\"#5ba865\")\nplt.plot(fs, eps_GaP.real, label=\"Re GaP\", c=\"#778bdf\")\nplt.plot(fs, eps_GaP.imag, \"--\", label=\"Im GaP\", c=\"#778bdf\")\nplt.ylim(-50, 30)\nplt.xlabel(\"frequency (eV)\")\nplt.ylabel(\"relative permittivity\")\nplt.legend()\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Define the simulation\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "lattice = nn.Lattice(([1.0, 0], [0, 1.0]))\nsup = lattice.Layer(\"Superstrate\", epsilon=1)\nfreqs = np.linspace(1, 25, 500)\n\n\ndef simulation(mat, slab_flag=False):\n    ref = []\n    trans = []\n    pressure = []\n    for f in freqs:\n        w = h * c / e / f * 1e6\n        pw = nn.PlaneWave(wavelength=w, angles=(0, 0, 0))\n        eps_sub = epsilon_GaP(w) if mat == \"GaP\" else epsilon_Al(w)\n        if slab_flag:\n            sub = lattice.Layer(\"Substrate\", epsilon=1)\n            slab = lattice.Layer(\"Slab\", epsilon=eps_sub, thickness=0.4)\n            stack = [sup, slab, sub]\n        else:\n            sub = lattice.Layer(\"Substrate\", epsilon=eps_sub)\n            stack = [sup, sub]\n        sim = nn.Simulation(stack, pw, 1)\n        R, T = sim.diffraction_efficiencies()\n        Tx1, Ty1, Tz1 = sim.get_z_stress_tensor_integral(\"Superstrate\")\n        if slab_flag:\n            Tx3, Ty3, Tz3 = sim.get_z_stress_tensor_integral(\"Slab\")\n            Tz1 -= Tz3\n        # Tz -= Tz1\n        ref.append(R)\n        trans.append(T)\n        pressure.append(-Tz1)\n\n    return ref, trans, pressure"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Do the calculation\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ref_Al_slab, trans_Al_slab, pressure_Al_slab = simulation(\"Al\", slab_flag=True)\nref_GaP, trans_GaP, pressure_GaP = simulation(\"GaP\", slab_flag=False)\nref_Al, trans_Al, pressure_Al = simulation(\"Al\", slab_flag=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Figure 3\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure()\nplt.plot(freqs, ref_Al, label=\"$r$ (Al, bulk)\", lw=1, c=\"#bbdf77\")\nplt.plot(freqs, ref_GaP, label=\"$r$ (GaP, bulk)\", lw=1, c=\"#77addf\")\nplt.plot(freqs, ref_Al_slab, \"--\", label=\"$r$ (Al, 400nm)\", c=\"#635a5e\")\nplt.plot(freqs, trans_Al_slab, \"--\", label=\"$t$ (Al, 400nm)\", c=\"#e77d7d\")\nplt.xlabel(\"frequency (eV)\")\nplt.ylabel(\"Intensity\")\nplt.legend()\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Figure 4\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure()\nplt.plot(freqs, pressure_Al, label=\"Al, bulk\", lw=1, c=\"#bbdf77\")\nplt.plot(freqs, pressure_GaP, label=\"GaP, bulk\", lw=1, c=\"#77addf\")\nplt.plot(freqs, pressure_Al_slab, \"--\", label=\"Al, 400nm\", c=\"#635a5e\")\nplt.xlabel(\"frequency (eV)\")\nplt.ylabel(\"total pressure (SI units)\")\nplt.ylim(0, 2.1)\nplt.tight_layout()\nplt.legend()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Define the simulation\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def simulation_angle(eps_sup, eps_sub, angle):\n    lattice = nn.Lattice(([1.0, 0], [0, 1.0]))\n    sup = lattice.Layer(\"Superstrate\", epsilon=eps_sup)\n    slab = lattice.Layer(\"Slab\", epsilon=1, thickness=1)\n    sub = lattice.Layer(\"Substrate\", epsilon=eps_sub)\n    pressure = []\n    for theta in angle:\n        pw = nn.PlaneWave(wavelength=1 / 0.01, angles=(theta, 0, 90))\n        sim = nn.Simulation([sup, slab, sub], pw, 1)\n        sim.solve()\n        T1x, T1y, T1z = sim.get_z_stress_tensor_integral(\"Slab\")\n        pressure.append(-T1z)\n    return np.array(pressure)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Figure 6\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "angle = np.linspace(0, 90 * 0.99, 500)\n\nplt.figure()\nfor eps_sup, eps_sub in zip([8, 9, 10, 10], [9, 9, 9, 9 + 0.1j]):\n    pressure = (\n        simulation_angle(eps_sup, eps_sub, angle) * 3.5e9 / (c / eps_sup.real**0.5)\n    )\n    if np.imag(eps_sub) == 0:\n        label = f\"{eps_sup} | 1 | {eps_sub}\"\n    else:\n        label = f\"{eps_sup} | 1 | {eps_sub.real} + {eps_sub.imag}j\"\n    plt.plot(np.cos(angle * nn.pi / 180), pressure, label=label)\n\nplt.ylim(-60, 20)\nplt.legend()\nplt.xlabel(r\"$\\cos{\\theta}$\")\nplt.ylabel(\"pressure on III (SI units)\")\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import nannos.utils.jupyter\n%nannos_version_table"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "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.10.13"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}