{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nThe empty circular duct with flow and higher-order modes\n========================================================\n\nIn this example we extract higher-order modes from measurement data in a circular duct with moderate mean flow.\nThe data is part of `this study <https://www.ingentaconnect.com/content/dav/aaua/2016/00000102/00000005/art00008>`_\\,\nwhich is referred here for further details.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "![](../../image/wave2.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": [
        "The duct has a circular cross section and is filled with air for which we use the standard properties.\nThe duct radius is 0.075 m. The distance from the microphones to the reference cross section is 0.55 m.\nThe highest frequency of interest is 3200 Hz.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "section = \"circular\"\nradius = 0.075  # m\ndistance = 0.55  # [m]\nf_max = 3200  # Hz"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We create objects for the upstream and the downstream section of the duct.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "td_us = acdecom.WaveGuide(dimensions=(radius,), cross_section=section, f_max=f_max, damping=\"dokumaci\",\n                          distance=distance, flip_flow=True)\ntd_ds = acdecom.WaveGuide(dimensions=(radius,), cross_section=section, f_max=f_max, damping=\"dokumaci\",\n                          distance=distance)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div class=\"alert alert-info\"><h4>Note</h4><p>The standard flow direction is in $P_+$ direction. On the upstream side, the Mach-number therefore must be\n  either set negative or the argument *flip_flow* must be set to *True*.</p></div>\n\n2. Sensor Positions\n-------------------\n\nThe microphone coordinates are saved in the\n`emptyUS.mic <https://github.com/ssackMWL/acdecom/blob/master/examples/data/emptyUS.mic>`_ and\n`emptyDS.mic <https://github.com/ssackMWL/acdecom/blob/master/examples/data/emptyDS.mic>`_ file.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "td_us.read_microphonefile(\"data/emptyUS.mic\", cylindrical_coordinates=True)\ntd_ds.read_microphonefile(\"data/emptyDS.mic\", cylindrical_coordinates=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div class=\"alert alert-info\"><h4>Note</h4><p>In this case, the microphone coordinates are defined in a cylindrical coordinate system with the circumferential\n  position in deg. Therefore, we set the argument *cylindrical_coordinates* to *True*. This will transform the\n  circumferential position from deg. to radians.</p></div>\n\n3. Decomposition\n----------------\n\nThe measurement must be pre-processed in a format that is understood by the *WaveGuide* object. Generally,\nthis must be a numpy.ndarray, wherein the columns contain the measurement data, such as the measured frequency, the\npressure values for that frequency, the bulk Mach-number, and the temperature. The rows can be different frequencies\nor different sound excitations (cases). In this example, the\nmeasurement was post-processed into the\n`higherOrderModes.txt <https://github.com/ssackMWL/acdecom/blob/master/examples/data/higherOrderModes.txt>`_ file and\ncan be loaded with the `numpy.loadtxt <https://numpy.org/doc/stable/reference/generated/numpy.loadtxt.html>`_\nfunction.\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>The pressure used for the decomposition must be pre-processed, for example to account for microphone.</p></div>\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pressure = numpy.loadtxt(\"data/higherOrderModes.txt\",dtype=complex, delimiter=\",\", skiprows=1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We examine the file header to understand how the data is stored in our input file.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "with open(\"data/higherOrderModes.txt\") as pressurefile:\n    print(pressurefile.readline().split(\",\"))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Mach-number, temperature, and frequency are stored in columns 0, 1, and 2. The upstream microphones 1-12 are in\ncolumns 3 - 14, the downstream microphones 13-24 are in columns 15-26, and the case number is in the last column.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Mach_number = 0\ntemperature = 1\nf = 2\nmics_us = range(3, 15)\nMics_ds = range(15, 27)\ncase = -1"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, we can decompose the sound-fields into the propagating modes. We decompose the sound-fields on the upstream\nand downstream side of the duct, using the two *WaveGuide* objects defined earlier.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "decomp_us, headers_us = td_us.decompose(pressure, f, mics_us, temperature_col=temperature, case_col=case,\n                                        Mach_col=Mach_number)\ndecomp_ds, headers_ds = td_ds.decompose(pressure, f, Mics_ds, temperature_col=temperature, case_col=case,\n                                        Mach_col=Mach_number)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. note ::\n  The decomposition may show warnings for ill-conditioned modal matrices. This typically happens for frequencies close\n  to the cut-on of a mode. However, it can also indicate that the microphone array is insufficient to separate the\n  modes. The condition number of the wave decomposition is stored in the data returned by *decompose* and\n  should be checked in case a warning is triggered.\n\n4. Further Post-processing\n--------------------------\nWe can print the *headers_ds* to see the names of the columns of the arrays that store the decomposed sound fields.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(headers_ds)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We use that information to extract the modal data for the different frequencies and cases.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plusmodes = [0,1,2,3,4,5]\nminusmodes = [6,7,8,9,10,11]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Furthermore, we can get the unique decomposed frequency points.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "frequs = numpy.abs(numpy.unique(decomp_us[:,headers_us.index(\"f\")]))\nnof = frequs.shape[0]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For each of the frequencies we can compute the scattering matrix by solving a linear system of equations\n$S = p_+ p_-^{-1}$\\, where $S$ is the scattering matrix and  $p_{\\pm}$ are matrices containing the\nacoustic modes palded in rows and the different test cases placed in columns.\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>Details for the computation of the Scattering Matrix and the procedure to measure the different test-cases can be\n  found in `this study <https://www.ingentaconnect.com/content/dav/aaua/2016/00000102/00000005/art00008>`_\\.</p></div>\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "S = numpy.zeros((12,12,nof), dtype=complex)\nfor fIndx, f in enumerate(frequs):\n    frequ_rows = numpy.where(decomp_us[:,headers_us.index(\"f\")] == f)\n    ppm_us = decomp_us[frequ_rows]\n    ppm_ds = decomp_ds[frequ_rows]\n    pp = numpy.concatenate((ppm_us[:,plusmodes].T, ppm_ds[:,plusmodes].T))\n    pm = numpy.concatenate((ppm_us[:,minusmodes].T, ppm_ds[:,minusmodes].T))\n    S[:,:,fIndx] = numpy.dot(pp,numpy.linalg.pinv(pm))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "5. Plot\n-------\nFinally, we can plot the transmission and reflection coefficients of the 6 propagating modes.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mode_names = td_us.mode_vector\nfig,axs=plt.subplots(6,1,figsize=(10,10))\naxs[0].set_title(\"Empty Circular Duct with Higher Order Modes\")\nfor mode in range (6):\n    axs[mode].plot(frequs, numpy.abs(S[mode,mode+6,:]),\n                   color=\"#67A3C1\", label = str(mode_names[mode]) + \"-Mode Transmission\")\n    axs[mode].plot(frequs, numpy.abs(S[mode,mode,:]), ls=\"--\",\n                   color=\"#D38D7B\", label = str(mode_names[mode]) + \"-Mode Reflection\")\n    axs[mode].set_xlim([0, 3200])\n    axs[mode].set_ylim([0, 1.1])\n    axs[mode].set_xticks([])\n    axs[mode].legend()\n\naxs[2].set_ylabel(\"Scattering Magnitude\")\naxs[5].set_xticks(range(0,3200,250))\nplt.xlabel(\"Frequency [Hz]\")\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
}