.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        Click :ref:`here <sphx_glr_download_auto_examples_plot_RienstraCutOns.py>`     to download the full example code
    .. rst-class:: sphx-glr-example-title

    .. _sphx_glr_auto_examples_plot_RienstraCutOns.py:


How to compute wavenumbers and cut-ons
======================================

In this example, the wavenumbers and cut-on modes in a circular waveguide without flow are computed and plotted. The r
esults correspondto Rienstra "Fundamentals of Duct Acoustics" (Figures 7 and 8).

.. image:: ../../image/ripples.JPG
   :width: 800

1. Initialization
-----------------
First, we import the packages needed for this example.


.. code-block:: default

    import matplotlib.pyplot as plt
    import numpy
    import acdecom








In this example we compute the wavenumbers for arbitrary mode orders in a circular duct of radius 0.1 m
without background flow.


.. code-block:: default

    section = "circular"
    radius = 1  # [m]
    Mach_number = 0







We create an object for our test section. We use the standard conditions for ambient air inside the duct, which are
predefined.


.. code-block:: default

    td = acdecom.WaveGuide(cross_section = section, dimensions=(radius,), M=Mach_number)







From the *td* object, we can get the  :class:`.WaveGuide` specific properties.


.. code-block:: default

    print(td.get_domainvalues())




.. rst-class:: sphx-glr-script-out

 Out:

 .. code-block:: none

    {'density': 1.2002432611667804, 'dynamic_viscosity': 1.813e-05, 'specific_heat': 1002.8961443192019, 'heat_capacity': 1.401, 'thermal_conductivity': 0.02587, 'speed_of_sound': 343.35637592387417, 'Mach_number': 0, 'radius': 1, 'bulk-viscosity': 1.0878e-05}




2. Extract the Wavenumbers
--------------------
We define a non-dimensional frequency from which we compute the maximum frequency for our decomposition.


.. code-block:: default

    omega_nondimensional = 20
    omega = omega_nondimensional * td.speed_of_sound/radius
    f_decomb = omega/(2*numpy.pi)







We can conveniently compute the wavenumbers for an arbitrary set of modes by creating two arrays containing all
combinations of *m* and *n*. In this example we compute the first 24 *m* modes and the first 10 *n* modes.


.. code-block:: default

    m = numpy.arange(0,24)
    n = numpy.arange(0,10)
    nn, mm = numpy.meshgrid(n, m)







From the mode arrays we compute the wavenumbers in the left and in the right running direction. We can alter the
direction by setting the *sign* parameter. The  :math:`p_+` direction has the sign 1 and the :math:`p_-` direction
has the sign -1.


.. code-block:: default

    k_left = td.get_wavenumber(mm, nn, sign = -1, f=f_decomb)
    k_right = td.get_wavenumber(mm, nn, sign = 1, f=f_decomb)







3. Plot
----
Finally, we can plot the real and the imaginary part of the wavenumbers on the m - n grid.


.. code-block:: default

    text_params = {'ha': 'center', 'va': 'center', 'family': 'sans-serif',
                   'fontweight': 'bold'}
    fig, axs = plt.subplots(2,1,figsize = (10,10))
    ax=axs[0]
    limit = numpy.max(numpy.abs(numpy.imag(k_right)))
    kimag_left = ax.pcolor(-n, m, numpy.imag(k_left), cmap="rainbow",vmin=-limit,vmax=limit)
    kimag_right = ax.pcolor(n, m, numpy.imag(k_right), cmap="rainbow",vmin=-limit,vmax=limit)
    ax.set_xticks(numpy.append(-n,n))
    ax.set_yticks(m)
    ax.set_xticklabels(numpy.append(n,n))
    ax.set_xlabel("n")
    ax.set_ylabel("m")
    ax.set_title("$Im(k_{mn})$, $\omega$ = "+str(omega_nondimensional))
    fig.colorbar(kimag_left, ax=ax)
    ax.grid(b=True, which='major', color='#666666', linestyle='-')
    ax.text(numpy.max(n)*0.5,numpy.max(m)*0.5,"right running wave",bbox=dict(boxstyle="round",
                       ec="k",fc="w"),**text_params)
    ax.text(-numpy.max(n)*0.5,numpy.max(m)*0.5,"left running wave",bbox=dict(boxstyle="round",
                       ec="k",fc="w"),**text_params)
    ax=axs[1]
    limit = numpy.max(numpy.abs(numpy.real(k_right)))
    kimag_left = ax.pcolor(-n, m, numpy.real(k_left), cmap="rainbow",vmin=-limit,vmax=limit)
    kimag_right = ax.pcolor(n, m, numpy.real(k_right), cmap="rainbow",vmin=-limit,vmax=limit)
    ax.set_xticks(numpy.append(-n,n))
    ax.set_xticklabels(numpy.append(n,n))
    ax.set_yticks(m)
    ax.grid(b=True, which='major', color='#666666', linestyle='-')
    ax.text(numpy.max(n)*0.5,numpy.max(m)*0.5,"right running wave",bbox=dict(boxstyle="round",
                       ec="k",fc="w"),**text_params)
    ax.text(-numpy.max(n)*0.5,numpy.max(m)*0.5,"left running wave",bbox=dict(boxstyle="round",
                       ec="k",fc="w"),**text_params)
    ax.set_xlabel("n")
    ax.set_ylabel("m")
    ax.set_title("$Re(k_{mn})$, $\omega$ = "+str(omega_nondimensional))
    fig.colorbar(kimag_left, ax=ax)
    plt.show()



.. image:: /auto_examples/images/sphx_glr_plot_RienstraCutOns_001.png
    :alt: $Im(k_{mn})$, $\omega$ = 20, $Re(k_{mn})$, $\omega$ = 20
    :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 Out:

 .. code-block:: none

    /home/docs/checkouts/readthedocs.org/user_builds/acdecom/checkouts/latest/examples/plot_RienstraCutOns.py:64: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.
      kimag_left = ax.pcolor(-n, m, numpy.imag(k_left), cmap="rainbow",vmin=-limit,vmax=limit)
    /home/docs/checkouts/readthedocs.org/user_builds/acdecom/checkouts/latest/examples/plot_RienstraCutOns.py:65: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.
      kimag_right = ax.pcolor(n, m, numpy.imag(k_right), cmap="rainbow",vmin=-limit,vmax=limit)
    /home/docs/checkouts/readthedocs.org/user_builds/acdecom/checkouts/latest/examples/plot_RienstraCutOns.py:80: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.
      kimag_left = ax.pcolor(-n, m, numpy.real(k_left), cmap="rainbow",vmin=-limit,vmax=limit)
    /home/docs/checkouts/readthedocs.org/user_builds/acdecom/checkouts/latest/examples/plot_RienstraCutOns.py:81: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.
      kimag_right = ax.pcolor(n, m, numpy.real(k_right), cmap="rainbow",vmin=-limit,vmax=limit)




Additionally, we can compute the wavenumbers for a set of modes as a function of the frequency. Again, we create
two arrays that contain the combinations of modes. Here, we compute only the *m* modes and set *n* to 0.


.. code-block:: default

    m = numpy.arange(0, 17)
    n = numpy.zeros((m.shape[0],),dtype=int)







We also create an array for the dimensionless frequencies from 0.1 to 20 and allocate arrays for the left and right
running wavenumbers.


.. code-block:: default

    omegas = numpy.linspace(0.1,20,100)
    k_left = numpy.zeros((m.shape[0],omegas.shape[0]))
    k_right = numpy.zeros((m.shape[0],omegas.shape[0]))







We iterate through all the dimensionless frequencies and compute the wavenumbers for all modes.


.. code-block:: default

    for ii,omega_nondimensional in enumerate(omegas):
        omega = omega_nondimensional * td.speed_of_sound / radius
        f_decomb = omega / (2 * numpy.pi)
        k_left[:,ii] = numpy.real(td.get_wavenumber(m,n,f=f_decomb,sign=-1))
        k_right[:,ii] = numpy.real(td.get_wavenumber(m,n,f=f_decomb,sign=1))







Finally, we plot the results.


.. code-block:: default

    plt.figure(2)
    for mode in range(k_left.shape[0]):
        plt.plot(omegas,k_right[mode,:],"k")
        plt.plot(omegas, k_left[mode, :], "k")

    plt.text(0,0,"m = 0",bbox=dict(boxstyle="round",
                       ec="k",fc="w"),**text_params)
    plt.text(numpy.max(omegas)*0.95,0,"m = "+str(numpy.max(m)),
             bbox=dict(boxstyle="round", ec="k",fc="w"),**text_params)
    plt.text(numpy.max(omegas)*0.5,numpy.max(k_right)*0.2,"cut-on right running",
             bbox=dict(boxstyle="round", ec="k",fc="w"),**text_params)
    plt.text(numpy.max(omegas)*0.5,numpy.min(k_left)*0.2,"cut-on left running",
             bbox=dict(boxstyle="round", ec="k",fc="w"),**text_params)
    plt.xlabel("$\omega$")
    plt.ylabel("$Re(k_{m0})$")
    plt.grid()
    plt.xlim(0,numpy.max(omegas))
    plt.show()



.. image:: /auto_examples/images/sphx_glr_plot_RienstraCutOns_002.png
    :alt: plot RienstraCutOns
    :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** ( 0 minutes  1.550 seconds)


.. _sphx_glr_download_auto_examples_plot_RienstraCutOns.py:


.. only :: html

 .. container:: sphx-glr-footer
    :class: sphx-glr-footer-example



  .. container:: sphx-glr-download sphx-glr-download-python

     :download:`Download Python source code: plot_RienstraCutOns.py <plot_RienstraCutOns.py>`



  .. container:: sphx-glr-download sphx-glr-download-jupyter

     :download:`Download Jupyter notebook: plot_RienstraCutOns.ipynb <plot_RienstraCutOns.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_