{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nHow to compute wavenumbers and cut-ons\n======================================\n\nIn this example, the wavenumbers and cut-on modes in a circular waveguide without flow are computed and plotted. The r\nesults correspondto Rienstra \"Fundamentals of Duct Acoustics\" (Figures 7 and 8).\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "![](../../image/ripples.JPG)\n\n   :width: 800\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "1. Initialization\n-----------------\nFirst, we import the packages needed for this example.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy\nimport acdecom"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In this example we compute the wavenumbers for arbitrary mode orders in a circular duct of radius 0.1 m\nwithout background flow.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "section = \"circular\"\nradius = 1  # [m]\nMach_number = 0"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We create an object for our test section. We use the standard conditions for ambient air inside the duct, which are\npredefined.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "td = acdecom.WaveGuide(cross_section = section, dimensions=(radius,), M=Mach_number)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "From the *td* object, we can get the  :class:`.WaveGuide` specific properties.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(td.get_domainvalues())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "2. Extract the Wavenumbers\n--------------------\nWe define a non-dimensional frequency from which we compute the maximum frequency for our decomposition.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "omega_nondimensional = 20\nomega = omega_nondimensional * td.speed_of_sound/radius\nf_decomb = omega/(2*numpy.pi)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can conveniently compute the wavenumbers for an arbitrary set of modes by creating two arrays containing all\ncombinations of *m* and *n*. In this example we compute the first 24 *m* modes and the first 10 *n* modes.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m = numpy.arange(0,24)\nn = numpy.arange(0,10)\nnn, mm = numpy.meshgrid(n, m)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "From the mode arrays we compute the wavenumbers in the left and in the right running direction. We can alter the\ndirection by setting the *sign* parameter. The  $p_+$ direction has the sign 1 and the $p_-$ direction\nhas the sign -1.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "k_left = td.get_wavenumber(mm, nn, sign = -1, f=f_decomb)\nk_right = td.get_wavenumber(mm, nn, sign = 1, f=f_decomb)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "3. Plot\n----\nFinally, we can plot the real and the imaginary part of the wavenumbers on the m - n grid.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "text_params = {'ha': 'center', 'va': 'center', 'family': 'sans-serif',\n               'fontweight': 'bold'}\nfig, axs = plt.subplots(2,1,figsize = (10,10))\nax=axs[0]\nlimit = numpy.max(numpy.abs(numpy.imag(k_right)))\nkimag_left = ax.pcolor(-n, m, numpy.imag(k_left), cmap=\"rainbow\",vmin=-limit,vmax=limit)\nkimag_right = ax.pcolor(n, m, numpy.imag(k_right), cmap=\"rainbow\",vmin=-limit,vmax=limit)\nax.set_xticks(numpy.append(-n,n))\nax.set_yticks(m)\nax.set_xticklabels(numpy.append(n,n))\nax.set_xlabel(\"n\")\nax.set_ylabel(\"m\")\nax.set_title(\"$Im(k_{mn})$, $\\omega$ = \"+str(omega_nondimensional))\nfig.colorbar(kimag_left, ax=ax)\nax.grid(b=True, which='major', color='#666666', linestyle='-')\nax.text(numpy.max(n)*0.5,numpy.max(m)*0.5,\"right running wave\",bbox=dict(boxstyle=\"round\",\n                   ec=\"k\",fc=\"w\"),**text_params)\nax.text(-numpy.max(n)*0.5,numpy.max(m)*0.5,\"left running wave\",bbox=dict(boxstyle=\"round\",\n                   ec=\"k\",fc=\"w\"),**text_params)\nax=axs[1]\nlimit = numpy.max(numpy.abs(numpy.real(k_right)))\nkimag_left = ax.pcolor(-n, m, numpy.real(k_left), cmap=\"rainbow\",vmin=-limit,vmax=limit)\nkimag_right = ax.pcolor(n, m, numpy.real(k_right), cmap=\"rainbow\",vmin=-limit,vmax=limit)\nax.set_xticks(numpy.append(-n,n))\nax.set_xticklabels(numpy.append(n,n))\nax.set_yticks(m)\nax.grid(b=True, which='major', color='#666666', linestyle='-')\nax.text(numpy.max(n)*0.5,numpy.max(m)*0.5,\"right running wave\",bbox=dict(boxstyle=\"round\",\n                   ec=\"k\",fc=\"w\"),**text_params)\nax.text(-numpy.max(n)*0.5,numpy.max(m)*0.5,\"left running wave\",bbox=dict(boxstyle=\"round\",\n                   ec=\"k\",fc=\"w\"),**text_params)\nax.set_xlabel(\"n\")\nax.set_ylabel(\"m\")\nax.set_title(\"$Re(k_{mn})$, $\\omega$ = \"+str(omega_nondimensional))\nfig.colorbar(kimag_left, ax=ax)\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Additionally, we can compute the wavenumbers for a set of modes as a function of the frequency. Again, we create\ntwo arrays that contain the combinations of modes. Here, we compute only the *m* modes and set *n* to 0.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m = numpy.arange(0, 17)\nn = numpy.zeros((m.shape[0],),dtype=int)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We also create an array for the dimensionless frequencies from 0.1 to 20 and allocate arrays for the left and right\nrunning wavenumbers.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "omegas = numpy.linspace(0.1,20,100)\nk_left = numpy.zeros((m.shape[0],omegas.shape[0]))\nk_right = numpy.zeros((m.shape[0],omegas.shape[0]))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We iterate through all the dimensionless frequencies and compute the wavenumbers for all modes.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for ii,omega_nondimensional in enumerate(omegas):\n    omega = omega_nondimensional * td.speed_of_sound / radius\n    f_decomb = omega / (2 * numpy.pi)\n    k_left[:,ii] = numpy.real(td.get_wavenumber(m,n,f=f_decomb,sign=-1))\n    k_right[:,ii] = numpy.real(td.get_wavenumber(m,n,f=f_decomb,sign=1))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, we plot the results.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure(2)\nfor mode in range(k_left.shape[0]):\n    plt.plot(omegas,k_right[mode,:],\"k\")\n    plt.plot(omegas, k_left[mode, :], \"k\")\n\nplt.text(0,0,\"m = 0\",bbox=dict(boxstyle=\"round\",\n                   ec=\"k\",fc=\"w\"),**text_params)\nplt.text(numpy.max(omegas)*0.95,0,\"m = \"+str(numpy.max(m)),\n         bbox=dict(boxstyle=\"round\", ec=\"k\",fc=\"w\"),**text_params)\nplt.text(numpy.max(omegas)*0.5,numpy.max(k_right)*0.2,\"cut-on right running\",\n         bbox=dict(boxstyle=\"round\", ec=\"k\",fc=\"w\"),**text_params)\nplt.text(numpy.max(omegas)*0.5,numpy.min(k_left)*0.2,\"cut-on left running\",\n         bbox=dict(boxstyle=\"round\", ec=\"k\",fc=\"w\"),**text_params)\nplt.xlabel(\"$\\omega$\")\nplt.ylabel(\"$Re(k_{m0})$\")\nplt.grid()\nplt.xlim(0,numpy.max(omegas))\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.7.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}