{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Frontiers: HFE Droplet in Water\n\nThis example corresponds to section 3.2 in\n`our publication <CitingOsaft>`.\nIn this example we compute the acoustic radiation force (ARF) on a HFE droplet\nsuspended in water subjected to a plane standing wave. We compare\nthe theories from King (1934), Yosioka & Kawasima (1955), and Gor'kov (1962).\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As always we start off by importing the nececassry Python modules. For this\nexample we are going to need the osaft library, and the third party packages\nNumPy and Matplotlib.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nfrom matplotlib import pyplot as plt\n\nimport osaft"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The next step is to define the properties for our example,\nthese include the material properties, the properties of the acoustic field\nand the radius of the particle. In the osaft library we are always\nassuming SI-units..\n\nThe wave type is set using the ``osaft.WaveType`` enum. Currently, there are\ntwo options:\n``osaft.WaveType.STANDING`` and ``osaft.WaveType.TRAVELLING`` for a plane\nstanding wave and a plane travelling wave, respectively.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# --------\n# Geometry\n# --------\n# Radius\nR_0 = 1e-6  # [m]\n\n# -----------------\n# Properties of HFE\n# -----------------\n# Speed of sound\nc_hfe = 659  # [m/s]\n# Density\nrho_hfe = 1_621  # [kg/m^3]\n\n# -------------------\n# Properties of Water\n# -------------------\n# Speed of sound\nc_w = 1_498  # [m/s]\n# Density\nrho_w = 997  # [kg/m^3]\n\n# --------------------------------\n# Properties of the Acoustic Field\n# --------------------------------\n# Frequency\nf = 5e6  # [Hz]\n# Pressure\np_0 = 1e5  # [Pa]\n# Wave type\nwave_type = osaft.WaveType.STANDING\n# Position of the particle in the field\nposition = np.pi / 4  # [rad]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Once all properties are defined we can initialize the solution classes.\nIn this example, we use the classes ``osaft.king1934.ARF()``,\n``osaft.yosioka1955.ARF()``, and ``osaft.gorkov1962.ARF()``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Theory of King\nking = osaft.king1934.ARF(\n    f=f,\n    R_0=R_0,\n    rho_s=rho_hfe,\n    rho_f=rho_w,\n    c_f=c_w,\n    p_0=p_0,\n    wave_type=wave_type,\n    position=position,\n)\n\n# Theory of Yosioka\nyosioka = osaft.yosioka1955.ARF(\n    f=f,\n    R_0=R_0,\n    rho_s=rho_hfe,\n    c_s=c_hfe,\n    rho_f=rho_w,\n    c_f=c_w,\n    p_0=p_0,\n    wave_type=wave_type,\n    position=position,\n)\n\n# Theory of Gor'kov\ngorkov = osaft.gorkov1962.ARF(\n    f=f,\n    R_0=R_0,\n    rho_s=rho_hfe,\n    c_s=c_hfe,\n    rho_f=rho_w,\n    c_f=c_w,\n    p_0=p_0,\n    wave_type=wave_type,\n    position=position,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First, want to evaluate the compressibility of the droplet and the\nfluid and the scattering coefficients $f_1$ and $f_2$.\nAll these quantities are properties of the solution classes and can easily be\nevaluated.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Compressibility\nprint(f\"{yosioka.scatterer.kappa_f = :.2e}\")\nprint(f\"{yosioka.fluid.kappa_f = :.2e}\")\n\n# Gor'kov scattering coefficients\nprint(f\"{gorkov.f_1 = :.2f}\")\nprint(f\"{gorkov.f_2 = :.2f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we want to compare the ARF from the different\ntheories in the small particle limit. To plot the ARF\nwe need to initialize an ``osaft.ARFPlot()`` instance.\nWith the method ``add_solutions()`` we can add our models to the instance.\nWith ``set_abscissa()`` we define the attribute that we want to\nplot the ARF against.  Here we select the radius `R_0`.\nFinally, ``osaft.ARFPlot.plot_solutions()`` will generate the plot.\n``plot_solutions()`` returns two Matplotlib objects, a ``Figure`` and\nan ``Axes`` object. These can be used to further manipulate the plot and to\nsave it.\nTo display the plot we need to call ``plt.show()``\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "arf_plot = osaft.ARFPlot()\n\n# Add solutions to be plotted\narf_plot.add_solutions(gorkov, yosioka, king)\n\n# Define independent plotting variable (in this case the radius)\narf_plot.set_abscissa(np.linspace(1e-6, 15e-6, 300), \"R_0\")\nfig, ax = arf_plot.plot_solutions()\n\n# Setting the axis labels and the title\nax.set_xlabel(r\"$R \\, \\mathrm{[m]}$\")\nax.set_title(r\"$\\mathrm{(a)} \\quad R \\ll  \\lambda$\")\n\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div class=\"alert alert-info\"><h4>Note</h4><p>It is only possible to plot the ARF  against\n      properties that wrap an underlying ``PassiveVariable``, i.e. an input\n      parameter of the model.\n      You can get a list of all input variables using the method\n      ``input_variables()``.</p></div>\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(f\"{gorkov.input_variables() = }\")\nprint(f\"{yosioka.input_variables() = }\")\nprint(f\"{king.input_variables() = }\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We redo our ARF plot, but this time the particle size is ranging from\n``1e-6`` to ``120e-6`` meters.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Plot\narf_plot.set_abscissa(np.linspace(1e-6, 120e-6, 300), \"R_0\")\nfig, ax = arf_plot.plot_solutions()\n\n# Additional Matplotlib commands for the publication\nax.set_xlabel(r\"$R \\, \\mathrm{[m]}$\")\nax.set_title(r\"$\\mathrm{(b)} \\quad R \\sim  \\lambda$\")\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Lastly, we want to plot the mode shapes of the particle. We can do this\nusing the ``osaft.ParticleWireframePlot()`` class. To plot a specific model\nwe have to pass this model when initializing ``ParticleWireframePlot``. We\ncan also pass a value for the ``scale_factor``. The displacements in the\nmode shape plot are exaggerated, the ``scale_factor`` is the ratio between\nthe maximal displacement and the particle radius in the exaggerated plot.\n\nIn our example we plot the mode shape for three different radii. Again,\nwe call ``plt.show()`` to display the plot.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Creating a figure with subplots\nfig, axes = plt.subplots(\n    1,\n    3,\n    figsize=(9, 2.5),\n    gridspec_kw={\n        \"width_ratios\": [1, 1, 1],\n    },\n)\n\n# List of radii\nradii = [1e-6, 30e-6, 90e-6]\n\n# We loop through the three radii and the three subplots\nfor radius, ax in zip(radii, axes):\n\n    # During each loop we change the radius in the model\n    yosioka.R_0 = radius\n\n    # We plot the wireframe plot in the respect\n    wireframe_plot = osaft.ParticleWireframePlot(yosioka, scale_factor=0.1)\n    wireframe_plot.plot(ax=ax)\n\n    # Making the plot prettier\n    um_r = int(yosioka.R_0 * 1e6)\n    ax.axis(False)\n    ax.set_aspect(1)\n    ax.set_title(rf\"$R_0 = {{{um_r}}}\\mathrm{{\\mu m}}$\")\n\nfig.tight_layout()\nplt.show()"
      ]
    }
  ],
  "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.9.15"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}