alchemyst/Skogestad-Python

{
  "cells": [
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline\n",
        "import numpy as np\n",
        "import scipy.signal as scs\n",
        "\n",
        "\n",
        "# Example 2.7 compares w_B and w_BT as indicators of performance\n",
        "# it illustrates that w_B might be a better indicator, but care\n",
        "# must be taken and the phase must also be considered.\n",
        "\n",
        "\n",
        "# Creates a general transfer function based on the coefficients\n",
        "# of the numerator and denominator\n",
        "def TF(num, den, s):\n",
        "    TF = np.polyval(num, s)/np.polyval(den, s)\n",
        "    return TF\n",
        "\n",
        "\n",
        "# Initialize freq range and other constants\n",
        "w = np.logspace(-2, 2, 100)\n",
        "z = 0.1\n",
        "tau = 1\n",
        "\n",
        "# Coeff found in sage of the T function\n",
        "Tnum = [-1, z]\n",
        "Tden = [tau, tau*z + 1, z]\n",
        "\n",
        "# Time domain output of a step response\n",
        "[Tim, yout, xout] = scs.lsim((Tnum, Tden), np.ones(np.size(w)), w)\n",
        "# Plot time domain response\n",
        "plt.subplot(2, 1, 1)\n",
        "plt.plot(Tim, yout)\n",
        "plt.ylabel('y(t)')\n",
        "plt.xlabel('Time [s]')\n",
        "plt.title('Time domain response')\n",
        "\n",
        "# Gain of both S and T transfer functions\n",
        "Tgain = np.abs([TF(Tnum, Tden, (1j*i)) for i in w])\n",
        "\n",
        "Sgain = np.abs([TF(Tnum, Tden, (1j*i)) - 1 for i in w])\n",
        "\n",
        "# Plot Gains\n",
        "plt.subplot(2, 1, 2)\n",
        "plt.loglog(w, Tgain)\n",
        "plt.loglog(w, Sgain)\n",
        "# TODO plot w180, wbt\n",
        "plt.legend([\"|T|\", \"|S|\"])\n",
        "plt.ylabel('Magnitude')\n",
        "plt.xlabel('Frequency [rad/s]')\n",
        "plt.title('Magnitude of S and T transfer function')\n",
        "plt.show()\n",
        "\n",
        "# TODO GM, PM, Ms, Mt, wb, wc\n"
      ],
      "outputs": [],
      "execution_count": null
    }
  ],
  "metadata": {
    "anaconda-cloud": {},
    "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.6.1"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 1
}