{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "\n",
    "import numpy as np\n",
    "import toolz\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.io import savemat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Frequency response\n",
    "\n",
    "Implicitly plot the function\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\gamma^2}{z^2} &= (\\omega^2 - \\alpha - \\frac{3}{4}\\beta z^2)^2 + (\\delta \\omega)^2 \\implies \\\\\n",
    "z &= \\frac{\\gamma}{\\sqrt{(\\omega^2 - \\alpha - \\frac{3}{4}\\beta z^2)^2 + (\\delta \\omega)^2}} \\implies\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def right(gamma, alpha, beta, delta, z, omega):\n",
    "    denom = (omega**2 - alpha - (3/4)*beta*z**2)**2 + (delta*omega)**2\n",
    "    return gamma/np.sqrt(denom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "gamma = 1\n",
    "alpha = 1\n",
    "beta_mul = 0.01\n",
    "delta = 0.1\n",
    "\n",
    "z, omega = np.meshgrid(\n",
    "    np.linspace(0, 15, 500),\n",
    "    np.linspace(0, 2, 500)\n",
    ")\n",
    "\n",
    "betas = (-0.3, 0.0, 10)\n",
    "for i, beta in enumerate(betas):\n",
    "    fig = plt.figure(i)\n",
    "    G = right(gamma, alpha, beta*beta_mul, delta, z, omega)\n",
    "    plt.contour(omega, z, (z-G), [0])\n",
    "    plt.xlabel(\"$\\omega$\")\n",
    "    plt.ylabel(\"$z/\\gamma$\")\n",
    "    fig.suptitle(\"Frequency response, beta=\" + str(beta))\n",
    "    fig.savefig('frequency-response-%d.png' % i)\n",
    "\n",
    "fig = plt.figure(len(betas))\n",
    "for beta in betas:\n",
    "    G = right(gamma, alpha, beta*beta_mul, delta, z, omega)\n",
    "    plt.contour(omega, z, (z-G), [0])\n",
    "    plt.xlabel(\"$\\omega$\")\n",
    "    plt.ylabel(\"$z/\\gamma$\")\n",
    "fig.suptitle(\"Frequency response, beta=\" + ', '.join(map(str, betas)))\n",
    "fig.savefig('frequency-response-%d.png' % len(betas))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma = 1\n",
    "alpha = 1\n",
    "beta_mul = 0.01\n",
    "delta = 0.1\n",
    "\n",
    "z, omega = np.meshgrid(\n",
    "    np.linspace(0, 15, 1000),\n",
    "    np.linspace(0, 2, 1000)\n",
    ")\n",
    "\n",
    "data = {'z': z, 'omega': omega}\n",
    "for beta in (-0.3, 0.0, 1, 4, 10):\n",
    "    data['G_beta%.1f' % beta] = right(gamma, alpha, beta*beta_mul, delta, z, omega)\n",
    "savemat('duffing-freq-response.mat', data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def right(gamma, alpha, beta, delta, z, f):\n",
    "    omega = 2*np.pi*f\n",
    "    denom = (omega**2 - alpha - (3/4)*beta*z**2)**2 + (delta*omega)**2\n",
    "    return gamma/np.sqrt(denom)\n",
    "\n",
    "f_resonance = 8e3\n",
    "omega_resonance = 2*np.pi*f_resonance\n",
    "\n",
    "gamma = 1e7\n",
    "alpha = omega_resonance**2\n",
    "beta_mul = 1e8\n",
    "delta = 200\n",
    "\n",
    "z, fs = np.meshgrid(\n",
    "    np.linspace(0, 1.2, 500),\n",
    "    np.linspace(7e3, 9e3, 500)\n",
    ")\n",
    "\n",
    "plt.figure()\n",
    "for beta in (-0.5, 0.0, 1, 4):\n",
    "    G = right(gamma, alpha, beta*beta_mul, delta, z, fs)\n",
    "    plt.contour(fs, z, (z-G), [0])\n",
    "    plt.xlabel(\"$\\omega$\")\n",
    "    plt.ylabel(\"$z$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import skew\n",
    "from scipy.interpolate import interp1d\n",
    "\n",
    "def get_pointmap(gamma, alpha, beta, delta, z, f):\n",
    "    fig = plt.figure()\n",
    "    G = right(gamma, alpha, beta, delta, z, f)\n",
    "    cs = plt.contour(f, z, (z-G), [0])\n",
    "    p = cs.collections[0].get_paths()[0]\n",
    "    v = p.vertices\n",
    "    plt.close(fig)\n",
    "    return v[:,0][::-1], v[:,1][::-1]\n",
    "\n",
    "def clip_pointmap(xs, ys):\n",
    "    # If softening spring, mirror the distribution across 0 and back\n",
    "    if skew(xs) > 0:\n",
    "        # Is order/direction different when softening?\n",
    "        raise NotImplementedError()\n",
    "        _xs, _ys = clip_pointmap(-1*xs, ys)\n",
    "        return -1*_xs, _ys\n",
    "    \n",
    "    # Find peak and continue s.t. x is monotonically increasing\n",
    "    indexes = []\n",
    "    N = len(xs)\n",
    "    max_x = xs[0] - 1e-9\n",
    "    for i in range(0, N):\n",
    "        if xs[i] > max_x:\n",
    "            indexes.append(i)\n",
    "            max_x = xs[i]\n",
    "    return xs[indexes], ys[indexes]\n",
    "\n",
    "def resample_pointmap(xs, ys, n):\n",
    "    lin = interp1d(xs, ys)\n",
    "    new_xs = np.linspace(xs[0], xs[-1], n)\n",
    "    return new_xs, lin(new_xs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "f_resonance = 8e3\n",
    "omega_resonance = 2*np.pi*f_resonance\n",
    "\n",
    "gamma = 1e7\n",
    "alpha = omega_resonance**2\n",
    "beta = 2e8\n",
    "delta = 200\n",
    "\n",
    "zs, fs = np.meshgrid(\n",
    "    np.linspace(0, 1.2, 500),\n",
    "    np.linspace(7e3, 9e3, 500)\n",
    ")\n",
    "\n",
    "plt.ioff()\n",
    "xs, ys = get_pointmap(gamma, alpha, beta, delta, zs, fs)\n",
    "plt.ion()\n",
    "\n",
    "print(\"Original\")\n",
    "print(len(xs))\n",
    "print(len(ys))\n",
    "print(skew(xs[::-1]))\n",
    "print(np.argmax(ys))\n",
    "plt.figure()\n",
    "plt.plot(xs, ys)\n",
    "\n",
    "xs, ys = clip_pointmap(xs, ys)\n",
    "print(\"Clipped\")\n",
    "print(len(xs))\n",
    "print(len(ys))\n",
    "plt.figure()\n",
    "plt.plot(xs, ys)\n",
    "\n",
    "xs, ys = resample_pointmap(xs, ys, 100)\n",
    "print(\"Resampled\")\n",
    "print(len(xs))\n",
    "print(len(ys))\n",
    "plt.figure()\n",
    "plt.plot(xs, ys)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rough method\n",
    "\n",
    "#### Simulating frequency response -- no coupling\n",
    "\n",
    " 1. Calculate linear stiffness $\\alpha$ from resonance frequencies.\n",
    " 2. Derive slope and half-width maximum by search on $\\beta$ and $\\delta$. Scale signals and apply MSE loss function.\n",
    " 3. Finally, scale amplitude by the inverse of scaling factor in 2).\n",
    " \n",
    "Step 3) outputs ballpark parameters for frequency scan simulations. See other notebook for coupling and further optimization."
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}