{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bilayer Sonophore model: computation of intermolecular pressure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from PySONIC.utils import logger, rmse, rsquared\n",
    "from PySONIC.neurons import getPointNeuron\n",
    "from PySONIC.core import BilayerSonophore, PmCompMethod\n",
    "from PySONIC.constants import *\n",
    "\n",
    "# Set logging level\n",
    "logger.setLevel(logging.INFO)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotPmavg(bls, Z, fs=15):\n",
    "    fig, ax = plt.subplots(figsize=(5, 3))\n",
    "    for skey in ['right', 'top']:\n",
    "        ax.spines[skey].set_visible(False)\n",
    "    ax.set_xlabel('Z (nm)', fontsize=fs)\n",
    "    ax.set_ylabel('Pressure (kPa)', fontsize=fs)\n",
    "    ax.set_xticks([0, bls.a * 1e9])\n",
    "    ax.set_xticklabels(['0', 'a'])\n",
    "    ax.set_yticks([-10, 0, 40])\n",
    "    ax.set_ylim([-10, 50])\n",
    "    for item in ax.get_xticklabels() + ax.get_yticklabels():\n",
    "        item.set_fontsize(fs)\n",
    "    ax.plot(Z * 1e9, bls.v_PMavg(Z, bls.v_curvrad(Z), bls.surface(Z)) * 1e-3, label='$P_m$')\n",
    "    ax.plot(Z * 1e9, bls.PMavgpred(Z) * 1e-3, label='$P_{m,approx}$')\n",
    "    ax.axhline(y=0, color='k')\n",
    "    ax.legend(fontsize=fs, frameon=False)\n",
    "    fig.tight_layout()\n",
    "\n",
    "def plotZprofiles(bls, f, A, Q, fs=15):\n",
    "    # Run simulation with integrated intermolecular pressure\n",
    "    data, _ = bls.simulate(f, A, Qm, Pm_comp_method=PmCompMethod.direct)\n",
    "    Z1 = data.loc[-NPC_FULL:, 'Z'].values * 1e9  # nm\n",
    "\n",
    "    # Run simulation with predicted intermolecular pressure\n",
    "    data, _ = bls.simulate(f, A, Qm, Pm_comp_method=PmCompMethod.predict)\n",
    "    Z2 = data.loc[-NPC_FULL:, 'Z'].values * 1e9  # nm\n",
    "    \n",
    "    # Plot figure    \n",
    "    t = np.linspace(0, 1 / f, NPC_FULL) * 1e6  # us\n",
    "    fig, ax = plt.subplots(figsize=(5, 3))\n",
    "    for skey in ['right', 'top']:\n",
    "        ax.spines[skey].set_visible(False)\n",
    "    ax.set_xlabel('time (us)', fontsize=fs)\n",
    "    ax.set_ylabel('Deflection (nm)', fontsize=fs)\n",
    "    ax.set_xticks([t[0], t[-1]])\n",
    "    for item in ax.get_xticklabels() + ax.get_yticklabels():\n",
    "        item.set_fontsize(fs)\n",
    "    \n",
    "    ax.plot(t, Z1, label='$P_m$')\n",
    "    ax.plot(t, Z2, label='$P_{m,approx}$')\n",
    "    ax.axhline(y=0, color='k')\n",
    "    ax.legend(fontsize=fs, frameon=False)\n",
    "    fig.tight_layout()\n",
    "    \n",
    "    return fig, Z1, Z2 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "pneuron = getPointNeuron('RS')\n",
    "bls = BilayerSonophore(32e-9, pneuron.Cm0, pneuron.Qm0)\n",
    "f = 500e3\n",
    "A = 100e3\n",
    "Qm = pneuron.Qm0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Profiles comparison over deflection range "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAADQCAYAAAA53LuNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xt8VOW56PHfs2ZygSRAgICgFETQioAgEQhqFVtarbXeTjliLdZe4LgtWi+7Vndra0+1lHqw3Z56Kq2t3VpblaqtooiyxcsuVBESgkiKW6iIYLgK5D4zz/ljzQyTMZcVsiYzk3m+n8/6rLXedXvI55OHN+961/uKqmKMMabnOekOwBhjcpUlYGOMSRNLwMYYkyaWgI0xJk0sARtjTJpYAjbGmDSxBGyMMWliCdgYY9LEErAxxqRJziTg8847T4GjWvbXNSs/7K+Vv//Xo76HLbbYklOLJzmTgPfs2XPU1zoihFXows/VGGM6lTMJuDvEgQgOaCTdoRhjehFLwB44IkQQS8DGGF9ZAvbAEVBLwMYYn1kC9sCtAVsThDHGX5aAPRDBmiCMMb6zBOyBtQEbY1LBErAHbgc0AZs9xBjjI0vAHlgbsDEmFSwBe2BtwMaYVLAE7IGIoAhiCdgY4yNLwB65TRDWBmyM8Y8lYI/c1Gs1YGOMfywBe6T2Es4Y4zNLwB5FrBuaMcZnloA9shqwMXDDDTcgIvGlrKyMOXPmUFtbm+7QslIw3QFkiwiCWBuwyXHV1dVUVFSwePFiVJU33niDm2++mVAoxOOPP57u8LKOJWCP7Es4Y9wEPGfOHKZPnw5ARUUFq1ev5rnnnktzZNnJmiA8iohj/YBNTqutraW2tpZTTjmlVXlJSQkNDQ1piiq7WQL2yMYDNrluw4YNAIwbNy5epqqsXr2aSZMmpSusrJaxTRAiMh14DfiMqq6Kln0WWAScBGwBblHVHvnbR3EQmxPO+OSOp99i0wcH0/LsccP78YMLT+n8xCSxBHzSSScRCoXYsWMHCxcuZOPGjTz55JN+h5kTMjIBi0gR8BAQSCgbB/wV+N/An4EvA0+JyGmq+laqY7IasMl11dXVAJSVlcXLRowYwUMPPcTFF1+crrCyWkYmYGAx8D4wJqHsemCNqt4Z3f++iJwZLZ+X6oBsMB7jp6OpgaZbdXU1M2fOZNGiRTiOQ2lpKccff3y6w8pqGZeAReTzwAXA+cCGhENnAY8lnb4KuLwn4lJxrBuayVmRSIRNmzZx4403Ul5e3u558+fPB2DLli1s3ryZu+66i9raWp544gl2797NsmXL+OQnP9lTYWe8jHoJJyKDgd8A3wD2Jx0+DtiRVPYBMKKD+80TkbUisnb37t3dis1GQzO5bMuWLTQ0NDBx4sQOz6usrKRPnz6sXLmSRx55hAULFjBhwgTWrFnDl770JZ544okeijg7ZFoN+H7gaVVdLiLHJR3rCzQmlTUBhe3dTFWXAEsAysvLu/UGzUZDM7ks1v7bUQKORCK88847rFy5Mv6l3PTp0zn//PPjxwcOHNgj8WaLjKkBi8hVwGTgpnZOaQAKksoKgLpUxnWEfQlncld1dTWFhYWMHTu23XNqamoYM2YMxcXFAFRVVTFjxoz48aqqqk5r0LnGcwIWkVEicrmIXC8i3xKRy0Sk3T//j8JXcZsZdonIYaAmWv6ciPwK2A4MS7pmOB9vlkgJFUGsBmxy1B133EFDQwOBQKDdcyorK5k8eXJ8f/369a32q6urLQEn6TABi0hQRL4qItXAu8CDwC3A7cCfgG0i8qaIzBWR7tamrwTGAZOiy+ei5d+IPu814Oyka2YCr3TzuZ4oDlg/YGPaVVVV1eqDjPXr18f3d+3aRVFRUbx2bFzttgGLSDluwj0E/A5YBmxRdd9EiYgApwLnADcAt4nIVar696MJRFVb1WRFJNbeu0NVa0XkXuBNEbkD+CNwBTANuOZontfl+KwbmjEdWrhwYav9ysrK+PYxxxzDli1bejqkjNfRS7h7gHmq+re2DqqqApXR5ecicjbwf4AzfY/SfV61iFyC+yXcLcBm4EJVfTsVz0tmY0EYY/zWbgJW1bO6ciNVfRkfk6+qvg9IUtky3Jp4j1PEPkU2xviqW+22IpIvImf4FUwmswHZjTF+85SARaRcRNaJSIuIhGMLbtewHnkJlm5q3dCMMT7zWgP+OW6ynQc04774WoT7IcT/TE1omUXFsW5oxhhfef0SbjLwKVV9U0TmATWqukREPsBNxktTFmGGsBqwMcZvXmvAAsQGU9gCTIhuP43bFa3XU3FwNJzuMIwxvYjXBLwR+Hx0exMQe/E2lIQxe3uzCAEcqwEbY3zktQnip8Cj0RdvfwR+ICJP4dZ+X0pVcJnE+gEbY/zmqQasqn8GpgOvq+o/cWvDzcCzwDdTF17msCYIY4zfOq0Bi0g/oFlV18bKVPUlcqTmG6ME7CWcMcZX7daARaRURJbhDox+SET+IiJl7Z3f20VwcKwJwuS4G264IT7Wr4hQVlbGnDlzqK2tTXdoWamjGvDdwGnAvwFh4Drcwc0v6YG4Mo6KYy/hTM6rrq6moqKCxYsXo6q88cYb3HzzzYRCIR5//PF0h5d1OkrA5wFfUdUXAUTkNeAlEclT1ZYeiS6DhCWAg7UBm9xWXV3NnDlzmD59OgAVFRWsXr2a5557Ls2RZaeOXsINARKne/87bpezoSmNKEOpBKwJwuS02tpaamtrOeWU1jM6l5SU0NDQkKaosltHCTgAR6p80XGAG4H8VAeVidwmCKsBm9y1YYM7Sfm4cePiZarK6tWrWw3EbrzLtEk5M5ZiNWDjo+e+C7uq0/PsYybA+Qs7Py9JLAGfdNJJhEIhduzYwcKFC9m4cSNPPvmk31HmhM4S8Okikjg9vACnicgxiSe1N2h7bxIRx7qhmZwWmxm5rOxIZ6gRI0bw0EMPcfHFF6crrKzWWQL+K0mDogOPJe0rufA5svWCMH46ihpoulVXVzNz5kwWLVqE4ziUlpZy/PHHpzssT8LhcIcTiqZLR23AxwOjo+uOltEpjjEj2Es4k8sikQibNm1ixowZlJeXc9ppp7WZfOfPn8/8+fM599xzGT58OA8++CCLFi1i+vTpnHDCCWzevLndZ6xYsYKKigomT57MKaecwquvvho/du2113LZZZdRUVHByJEjuf322zssj8WyYMECZs2axezZs9m+fTsXXHABJ598MpMnT+b111+Pn3vRRRfxwAMPAPDMM88wbdo0Dh8+3O2fW6dUtdMFGNLBsXle7pHuZcqUKdodzy/+hjb8YHC37mFMttq8ebMC+uijj3Z43tSpU/X666/XSCSiL730khYXF+uzzz6rqqq33HKL3nnnne1eu2fPHo1EIqqq+sILL+iFF14YPzZ9+nS98sorNRQK6d69e3XAgAH64Ycftlsei+Xqq6/WUCikoVBITz31VF2xYoWqqr788staXl4ev391dbWecMIJ+uqrr+qECRN0165dR/eDOsJTXvL6Eu4/ReQcVd0TKxCR0cADuPPALfH5/4WMo441QZjcFWv/nThxYrvnRCIR3nnnHVauXBn/Um769Omcf/758eMDBw5s9/qlS5fy+9//nrq6Og4dOsTpp58ev+7tt9/m6aefJhAIMHDgQI455hj279/fbvngwYOpqalh+fLlBAIBli1bxujRo5k1axYA48eP58MPP4w/e/z48UydOpXZs2fz8ssvM3Roz/S29Toc5fvAKhEZLK4bgGqglCNDU/ZqKkECloBNjqqurqawsJCxY8e2e05NTQ1jxoyhuLgYgKqqKmbMmBE/XlVV1W4CX7p0KU899RTLli2jqqqKSy+9lFNPdYca37x5M2VlZQwePBiAffv2sXfvXhobG9ssHz16NDU1NZx44omUlpbG40989oYNG5gwYUJ8f+vWrVRWVhIMBhk2bNjR/IiOitcEfCHwDu4APH8DfgzcAUxR1dc7urDXkICbgG1aIpOD7rjjDhoaGjp8kVVZWcnkyZPj++vXr2+1n5gE586d26rrWlVVFVOnTqW0tJR169Zx//33x69dv349tbW17N+/H1XlO9/5DgsWLGDjxo1tlufl5X0slmOPPZZNmzYB8NFHH3Hrrbfy7W9/G4A9e/Zw8cUXs2TJEi655BLuueceH35i3ngdjrIFuAz3y7jTgU+r6iLVHBqf0Yn+qOxFnDFtqqqqavVBxvr16+P7u3btoqioKF47XrduHccdd1z83Llz5/Loo48ydepUHnvsMQYPHhxPoOvWrWP+/Pl84QtfYNy4cZSVlXHbbbe1Wx6LJTEBX3755QQCAcaPH88ZZ5zBNddcw6xZs6ivr+eiiy7i9ttv58wzz+TWW2/lvvvuY+/evSn/eQGItlOjE5G22nUDwBXAe8DLsUJVnZeS6HxUXl6ua9eu7fzEdjx3302cX/sb+N5uCObkx4DG+OLAgQPMnj2bFStWeDr/3HPP5Wc/+xlTpkzxVJ4hkrvvtqmjl3DtNfasSTqeG3+TS/RPrxyq9BuTCgMGDPCcfAHeeuutVu21nZVnk3YTsKrO7MlAMp4TTcARS8DG9KTE3gpeyrNJRwOyf7UrN4r2jvh6tyPKVFYDNsb4rKOXcGeISKWIXCUiA9o7SUT6icg1uDMnz2jvvKwXewlnNWBjjE86aoL4poici9vl7H4RWYPbC2IPbuIuAyYDk4B1wLdV9YXUh5wmsSYI6wVhjPFJh1/Cqep/AjNEpAK3L/A03AHZFdgJvAAsyIm+wE70RxUJpTcOY0yv4elTZFVdDaxOcSyZTewlnDHGX16/hDPxDzEsARtj/GEJ2COxGrAxxmeWgL2KtgFHwpaAjTH+sATskUSbIEJhewlnjPFHlxKwiEwTka+KSImIjBOR3JnUM9oNTa0XhDHGJ54SsIj0F5EXcXtCPIDbB/inQLWIHNfhxb1E7P+asDVBGGN84rUG/FMgDxgB1EfLrgM+Au5OQVwZR6LjoEasCcIY4xOvTQgXAJep6g4Rd5Q1Vd0qIt8CnktVcBkl+hJOLQEbY3zitQY8ENjXRnkj0Me/cDKXRNuAw5aAjTE+8ZqAXwW+kbCv4naM/S458oVcLAGrtQEbY3zitQniZtxJOWcCBcC9wMm4k3J+JkWxZZRYAo5YLwhjjE+8zgm3EZgIrIguzcAjwMmqui514WWOeAK2GrAxxieeasAishK4TlW/n+J4MpbYl3DGGJ95bQOeBDSkMpBM5wSsCcIY4y+vbcCLcQdl/ymwlaRkrKof+B1YxrFuaMYYn3lNwN/Dffn2aVrPgizR/YDPcWUcsQRsjPGZ1wR8XkqjyAISjLUBt6Q5EmNMb+F1RoyXUx1Ixgvku+twc3rjMMb0Gl57Qazo6LiqftafcDKXE8wDQK0GbIzxidcmiB1tXDcWmADc42tEGUoctwZsTRDGGL94bYK4uq1yEfkBcKyvEWUoidaAsQRsjPFJd2fE+D0w249AMp0TbwO2BGyM8Ud3E/A0IDf6ZcVrwPYSzhjjj+68hOsHnAbc52tEGSpWA7aXcMYYvxztSzhwv4j7FfCwf+FkrlgvCCKWgI0x/ujWS7hcEggEiaigodxocTHGpJ7nNmAROVNEBke3LxeRZ0TkeyKSE1Pb5wUcWghaG7AxxjdeZ0X+F2AVMF5ETgUewh0H4lrgRymLLoPkBYQWAtYGbIzxjdfa6/XAPFVdBXwZ2KCqFwBXAl9JUWwZJS/gECJgbcDGGN94TcAjgRei258Dno1ubwGG+B1UJgpGa8DWD9gY4xevCXgHcIKInID7+fHz0fIzge2pCCzT5MfbgC0BG2P84bUb2hJgKdAEbFTV16LtwnfjjhXc6wUDDo1qTRDGGP947Yb2UxHZBJwA/CFavAeYr6oPpSq4TOK+hLMasDHGP15rwKjq07FtEckD3gVqUhFUJoq9hAvanHDGGJ947YY2UkReFJGpIlIIrAZeB7aJSHlKI8wQQcetAYs1QRhjfOL1JdzPgXxgF3AFMAZ3IJ5HgZ+lJrTMEnCEEAFLwMYY33htgpgJnKmq74nIF4BnVfUNEdkPVKYuvMwhIoQJItYEYYzxidcasAB1IhIAzuVIn+A+uD0jckJYApaAjTG+8VoDXgN8F7fnQ1/gGRE5FrgL+FuKYss4YQniWBOEMcYnXhPwAuAR3LbfG1V1t4j8O/BJ4POpCi7ThCQPRxvTHYYxppfw2g/4H0Byb4cfAterqvodVKayGrAxxk9dGY6yVERuEZHficgQ3LbgE1MXWuZpkQKCERuO0hjjD6/9gE8ENgNfwx0NrRi4DFgrIjNSF15maXHyCWrOvHM0xqSY1xrwPcBSVT2JI70evgw8DixMRWCZKCT55FkN2BjjE68JeDpwb2KBqkZwk+9kv4PKVCGngKBaAjbG+MNrAlbcPr/JhpBD/YBDTj752gy5897RGJNCXhPwX4Efi0hxdF9FZDTuJ8rLUhJZBgo7hThEwD7GMMb4wGsCvhEYCOwDinAH4tkCNAM3pya0zBN28t2Nlob0BmKM6RW8foiRr6oVIvIZYBJu4n1LVVemLrTMEw4UuhuhnGl1McakkNcEvFZELlXVF4EXUxlQJosEojXgkH0NZ4zpvq4MxpPz1T4NRt9DWgI2xvjAaw34t8ByEfkdsBVo1Qiqqo/4HVgmkmCBu2EJ2BjjA68J+PvR9W1tHFPcgXp6PcmLtgG3WAI2xnSf18F4PI8Z0ZtJXl93w2rAxhgfdJiARaQId9CdRmC1qh7ukagylJMXa4LI+eZwY4wP2k3AIjIReB4YGi36INoT4vUeiSwDBfLdGnC4uZ5AmmMxxmS/jpoWFgL/DczAnYCzBvhlTwSVqZwCNwG3NOb0HwLGGJ901ARRAcxU1UoAEfkm8A8RKVLVuh6JLsM4hf0ACNV/lOZIjDG9QUc14BLcaegBUNV3gRAwKNVBdUREAiLyExHZKSKHRWSpiAzt/MruCxSWABBuPNQTjzPG9HIdJWAHiCSVteC961qq/BC4CpgLfAo4DvhzTzw4WNCXkDqEGw72xOOMMb1cupNpl4hIPnA9cJ2qvhAtuxzYKiIzVDWlMzQX5gepoxBtshqwMab7OkvA14tIYntvEPgXEdmXeJKq3uV7ZG2bhNs0sirh2dtEZBtwFtBuAq6pqeGcc87p1sMPNLTQ58MDOH0eI++XVd26lzHG/YoLjW23HmdbNeGcds5rff2Rwo+dEy9PfLZ+/FpaD/etKIJQUti1uuqqVas8ndfRXd8Drkgq24U7F1wiBXoqAR8XXe9IKv8AGJF8sojMA+YBFBQUdPvhjkAEBycS7va9TObT6C+yqrsdiW+7v96x47Hf4/b2VY8kisT92ITiiftHzmv9HBL3oxdpQpzQ+hmt949cn1h+ZKVt3CNa/rF7kHik3WcnHkt4XKu4e4IkPF0SopF2zmnrvIAjnDpqSEriazcBq+qolDyxe/oCEVVNnhu+CShMPllVlwBLAMrLy9Xr/0rtqdp+gPCvP83IYUMZdM2z3bqXOSISUZpCEZpCYXfdkrDdUXlL0jkt7nZL2F1CYaU5aTsUjtAS1vg5LWElFI7QHFZCkQgtoQgtEfd4YgIRIgSJ4BAhSJhAfNtd5xHGkSPHjiyt94MSbnVdkCP78fPF3c4TdwmKki8RAo4SQAmIuwRFcXCXgERwiJVFcBKOBaPHEsvc+BUndiyhXIjeI3aNurFK7DyJ4Ki7n3iukFgeLUvYFxRR978Qif5MW5clbkfcRKjudvwajW4nXIPGzj0Sk59aCkrJu3WVr/eMyao2YNxBgBwRCapq4rQUBUDKu8b1yQ+wS/vgNOduG7CqmywPNrZwsCHEocYWDjWGoksLdc1hGppD1DeHqW8O09Acpr4lTH1TtKzlyPGG5jB1zSGaWkLkE6KAZgpooUBa3HVsie83R8+LlTXHzxkQiFDghCkUd50nYQrEXecTJk9C5BEmKGHyCJGnYYISIqghN2kGQgSdEIFAiIC6i5Ow9vuX+uh++EnrROIkLIHW+wg4SWUdLfFzpf17igNOe/cIJF0rrWOJ7RMrJ2Fb2jkveVs8XOMknNfZ/ZKvIb6dF2xrNjZ/ZFsC3h5dD0vYBhjOx5slfFdc4L6Ec1r2pvpRPaI5FGF/fTN7Dze767pm9h1uYl99C/vqmthf18L++uZ4cj3Y0EJLUx19wnWUSD39qaNEGiihnj7SRF+aKKKRPtJEmTRREmim2GmmWJopkib6ShN9aKJQGynURvK1iWCwmWAg+Q+aoyAOBPLByYNAsJPtInfbyXP3293OAyd4ZO0E3OQS23aC0USUuB+IbifuB91k1Wo/0PX7tZkEA0cShsk62ZaAq4BDwNnAwwAiMgoYBbyS6oeXFAY5oMXkNb+b6kd1W31ziPf21bPzQCM7P2pk18FGPoyud33UyP6PPiLYuI+BcpBBcohSDjFIDjJQDjJcDnFqoI7SQCP9pZ5i6inSw/SJ1BHIC0Nex89WBMkvgry+kN8X8oqi6/6QXxzdji2FECyEYAEECtx1bD8YO5afVBZdBxLKHfs43GSfrErAqtokIvcBd4vIHqAWuA94WVXXpPr5RflB9lFCYfMB9w1DmmseoXCEbXvr2LzrEFt317Ftbz3v7atj354P6VP3PsNkH8dEl1Gyj7OCHzHM2cfgyF4KaWyj1RxUAtB3ENJ3EBT2h8Jj3HVBPyjsl7Ddv3V5flE80UqwMO0/G2OyQVYl4Kjv4dbBHo6ulwPX9sSDHUeoCwzAIQyNB6BPaU88FnCbCzbtPMi6f+7nrQ8O8u7O3bC7huMj7zHS+ZBR8iHnBmsZyS6K9bDbKh6lThBKjkFKhkO/06FkGBSXQd9B0HcwFA2OrgchhQMseRrTQ7IuAUdfvt0UXXpcQ94A94Psur0pTcDNoQhv/nM/r2zZzVvvbsfZuZ5TIls42XmPmYHtjGQnTtD9UFHFQfuPwBk4GgaeDQNHQ+ko6Dcc+h2LFJW5bZDGmIySdQk43ZryS90EXL8XGOPrvQ82tvD8xl28Vv0OztZVTIlU80VnCzc77xMIRCAAof4jCQ6bDEOvhKGnwJBxyICRSDDf11iMMalnCbiLWgoHQT1Qv8eX+4Ujykuba3ll9WpKty3jLKnkUucdAk6EUH4RjJhKYOSXYcRUOHYKwcL+vjzXGJN+loC7qKWwzN04tLNb9znU2MJf/quKg2se4lNNq/iRsw0NCPWDJ+KMuwnGzCJ47BS3W5Qxpley3+4uihQPpYUgeQfeO6rrG1vCLH9+GcG1S/iSrqZAQhwYNJHw6XcRGH8pRf2G+RyxMSZTWQLuokElhXyggxl5FAn4zf9aQWjlnVwcqaRB+nJw3JWUzbyWAWUnpiBSY0ymswTcRWUlBfwzUsZx+7Z5nhfu4K6tbH14AVMOv8oB6ce2025h1OcW0KegJKWxGmMymyXgLiorLuA9HcKZe9/o/GOMcAs7l9/NgDfuYazC30b+L6ZcfhsD+tqLNGOMJeAuKysp4HkdhdO8EvZvdfvctmX/Nvb9x1yG7a9ilTONwf9jMTPGje/ZYI0xGc1653dRWUkB1ZHj3Z0P1rd9UvVSmn95BsF9/+Ce/rcy8aZnGG/J1xiTxBJwFx07oA+b9RM0BUtgy4utDzbXw1++BX/+OtXNw7hzxG+45lv/ysAi+0jCGPNx1gTRRQP65tO/qC9vFc/gtM3LoGG/+0ly7duw9Gto7dv8MnQRb5/0LX5xRTnBgP0fZ4xpmyXgozC6rIg/tlzEaU0r4E9fdsdd2PAozXn9mNfyXZpGnsODc6ZY8jXGdMgyxFEYO7SE5XsGE/7CL2DXRtj0V3aPnc3M+oXsHnoGS+ZOoSBo49MaYzqWkQlYRApEpEpErmzj2A0i8k8RqReRF0RkbE/HVzF6EIeaQmwY8kW4ZRubr97IrC2XEOw3hAevnkpJYScjlhtjDBmYgEWkBHgSmNjGsa8Dd+AORTkNd4645SLS/SmPu+DMMYPJDzj86fXtVH9wiCsfWEthMMDDX59GWUmPhmKMyWIZlYBF5DNAJTC0nVO+AyxW1aWqWg1cAQwBLuuhEAEoLcrnimmf4NG127nw/75GXkB4+BtTGTGwb0+GYYzJcpn2Eu7zwG+Bu4HGxAMiMgQ4EVgVK1PVwyKyFjgLeKTnwoR/u+BkPjGwLw0tYeZM/YR1NTPGdFlGJWBVvTG2LR//xPe46Dp59uMPgBEpDKtNeQGHr515fE8/1hjTi/RYAo7OXry1ncNNqtrGFJGtxP6+b0wqb6LN6SVBROYB86K7h0WkpvNIOzQY8GckdmNMtjia3/vlqnpeZyf1ZA14B3ByO8ciHq5viK6T33IVAHVtXaCqS4AlnqLzQETWqmq5X/czxmS+VP7e91gCVtUWYHM3brE9uh4GvJNQPhx4uxv3NcaYtMioXhAdUdVaYAtwdqxMRIqBcuCVdMVljDFHK6NewnmwGLhbRN4BNgJ3ATuBJ3ro+b41ZxhjskbKfu+zKgGr6q9EZABuIu4HvAacp6rNPfR8S8DG5JhU/t6Lqqbq3sYYYzqQNW3AxhjT21gC9kBEAiLyExHZKSKHRWSpiLT3ubQxxnhiCdibHwJXAXOBT+F+lffndAZkjMl+1gbcCRHJx/0K5jpVfTBaNgr3q74zVPVvaQvOGJPVrAbcuUlACa0HAdoGbMMdBMgY0wuJyKki8qyIHBCRZhGpEZG5fj7DEnDnMmoQIGNM6olIEbAC9/d8Gu745K8Av/bz/Y8l4M71BSLRT6kTtTsIkDEm6xXhfm9wnarWqOpm3A+/8nGHxfVFVn2IkSYNgCMiQVUNJZS3OwiQMSa7qWqtiPw/YK6ITAbG4jZHAvg24aMl4M4lDgK0PaF8OB9vljDG9AIiMgxYg/s7/jTwDG5zxFo/n2MJuHNVwCHcQYAehngviFHYIEDG9FaX4r58P0tVwwAi8rnosY/NFnG0LAF3QlWbROQ+3EGA9gC1wH3Ay6q6Jr3RGWNSZDfueDOXicjfgVOBf48e823mXUvA3nwPyMOtAecBy4Fr0xqRMSaVHgemAvcCxbhjkP8IuBU4HTcHdJt9iGGMMWli3dCMMSZNLAEbY0xQyPw4AAACpElEQVSaWAI2xpg0sQRsjDFpYgnYGGPSxBKwMcakiSVg0yuIyCgR0Q6WBzu5/jci8rUUxrdQRG5K1f1NdrJ+wKZXEJEAUNbGoWtxP6S5VFWfbOfac4GfA5NUNZKi+PoBm4CZqrolFc8w2ccSsOm1RORs4EXgXlW9sYPzVgO/VtXfpjiehcBwVfV1UG+TvawJwvRKIjIE+CPu6FW3dHBeBTAFeCqhbFV0EtaHROSgiLwvIveKSDB6/KsisllE5orIeyJSLyKPichwEfmDiNSJyHYRuSrpcY8Dl4vIcN//wSYrWQI2vY6IOMAfcAfPnt3GYPqJvgisUdV9SeU3AjVAOfAT3KaMyxOOjwauBD4PXAZcAlQDrwOn4Y4VcL+IlMYuUNU3gb3AeUf9jzO9iiVg0xt9D/g08BVV3d7JudOAt9oor1TVH6vqP1T1l8AGoCLheB5wrapuVNXngEpgo6r+QlVrcGdTKADGJN33LWB61/9JpjeyBGx6FRE5B7gdWBhNjJ0ZijvrdbJ/JO0fwK1RJ/rvhO26pP2G6Dp56MLd0WcaYwnY9B7Rdt9HgNeA73u8TGn796CprUckbIfb6DHhpQdFEAh7jM30cpaATa+Q0O7rAHNisxh4sJO2u6+lyuDoM42xAdlNr3ErbrvvHEBF5Jik481tvGgD96XZZ1MdHICICO7MCr/rieeZzGc1YNNbzMJtIvgTbg0zeXminev+CkwSkQE9EOMk3HnGlvXAs0wWsA8xTM6Lzvn1H9HeDql8zi+AAaqa3D/Y5ChLwCbnichngV8A47vQdtzVZwwA3gbOVtXkHhYmR1kThMl5qroCWAOkbDAe4Dbg55Z8TSKrARtjTJpYDdgYY9LEErAxxqSJJWBjjEkTS8DGGJMmloCNMSZN/j/pwJz+RlgMDAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x20efe1fd668>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Z = np.linspace(-0.4 * bls.Delta_, bls.a, 1000)\n",
    "fig = plotPmavg(bls, Z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Error quantification over a typical acoustic cycle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z-error: R2 = 0.9996, RMSE = 0.0454 nm (0.8212% dZ)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAADQCAYAAAA53LuNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl8VOX1+PHPmckkkz0hG4Swb7KJaETAFVusVC0u1Wqr1m5Srdraqq1drLa2tbS1v/bbWrXaaq1tXepSBRVcwA1FCoSw7zshIQvZl5k5vz/uBAMlZIAkdyY579frvpK5c5dDCIdnnnue5xFVxRhjTPfzuB2AMcb0VpaAjTHGJZaAjTHGJZaAjTHGJZaAjTHGJZaAjTHGJZaAjTHGJZaAjTHGJZaAjTHGJXFuBxCJ888/X1999VW3wzDGmEhJJAfFRAt43759bodgjDGdLiYSsDHG9ESWgI0xxiWWgI0xxiWWgI0xxiUxUQVhjk55bRPL126ketNiAmUbiK/dRXJLOZ5gI/EEaBEfjZ5kAv4+1KcMIiFvBNkjJjF26EAyk+PdDt+YXsMScA+gqqzfVUbxOy/h3zKPCY3/5ROesgPvN0kC1XFZBH2JBD0+vKH9xAe2klJXRXxdM+yFUJGwXgt4L3kijUPP56QzP83wvpku/qmM6fkkFlbEKCws1CVLlrgdRtRpCgR56513aFz0CNOa3iBd6mkUP7uzJhM/ZCo5J0wloe8YSOoDcpiyxFAIanZTu3stFWvfJbRtEf2qlpFAE+WayntJn8B76lf4xJmn4/d5u/8PaEzsiqgO2BJwDGoOhHh1/qtkLP4NZ+kSWohjZ99P0mfqF0kffS74/Mdx8Tr2F79C+Yf/YmDpm8QR5H2ZSMlJt3D+jM+QFG8fmoyJgCXgnujtJUWEXr2TcwLvUSsp7Bv/VQaddxOSktPp99KaEra//iCZxX8hLbSf92UiNWf8iE9OOxevJ6LfL2N6q+hLwCIyFlh5mLfOVNV32zvPEjBU1DTw1t9+yqdKH8UnIXaPncXgC29DEjO6/uZNteyc93vSl/6RxFA9Lyd+htFX/ZxRg/p3/b2NiU1RmYA/B/wBGH/IW+Wq2tLeeb09AX9YtBLPC7M4VVeyNXMq/b/wR3zZQ7s9Dq0rZ9sz32Xg1mfZq5m8O+5eLr70Knxeq2Y05hBRORfEOGC1qpYcsrWbfHszVeWl555g2HMzGK8b2H32rxl8y1xXki+AJGcx+LpHqP3CK3gTkrls5TeY+5uvsqdivyvxGBPr3EjAa7r5njGpORDimQd/woyiW2hOyEKvf4v8aV87fDVDN0sbMYXc2z5k+9ArmFn/b0p+fx4frrC/VmOOlhsJeJCIfCAiJSLyuohM6uYYol59Uwvzfn8DV+y9n51ZU+j3nXdIzB/rdlgHi09m8BcfpuS8PzGaLQz896d5bd5ct6MyJqZ0WwIWkURgKJAO3A58BtgNLBSR0Yc5/noRWSIiS8rKyg59u8dqaArw5u++xoXV/2LjgMsZ/I3/IAmpbofVrr5TP0/oy/OIi4vnzPeu4+l//oVYqKwxx+bWW29FRA5sOTk5XHXVVZSWlrodWkzq7odwaUCTqjaFX3uAYuBNVb25vfN6y0O45pYgb/7+a5xf8282D7uGoVf/X1R0OUQiUL2XvX+6kNz6TTw36Adcft238VipWo/zyU9+kvr6eu6//35UlY8++ojbbruNmTNn8swzz7gdXjSJvodwqlrdmnzDr0PAKmBAd8YRjUIhZe5D3+X8mn+zfvDnYyr5AsSl5ZH/zdfZnT6Bz23/Cc8+8nOCIWsJ9zTFxcVMmjSJyZMnM2XKFG655RYuu+wy5s+f73ZoMak7uyBOEZFqETm5zT4vcBJOEu7V5v7rj1y8789szP0UI6/9Y0wl31biT2fgzXPZmjGZz+76FU//5dfWHdGDlJaWUlpaytixBz+PSE1NpaGhwaWoYlt3toCLgK3AwyJyWnhQxl+BbOB33RhH1Hl/4Tymr7ubLUnjGfa1x8ETu3W14ktk8I3PszP9ZK7Y8TNe+Mef3A7JdJIVK1YAMGbMmAP7VJVFixZx0kknuRVWTOu2gf2qGhCRGcBs4CUgGXgPOEtVe20P/o7dexj45o3s92aQ//XnEF+i2yEdv/gkBnzjRbb/bgafXn8XL73cl4suvNTtqKLKPS+tYvXualfuPSY/jR9fdPRVNa0JeNSoUQQCAXbt2sV9993HypUref755zs7zF6hW2dWUdVdwBe6857RLBAIsuvxL3OKVFB5+X9ISMt1O6ROIwmpFNz4Ivt+dxanf3Qz7+YUcMZpVnEYy4qLiwHIyfl43pEBAwbwxBNPcPHFF7sVVkyzyXhcNP/vv2b6xp+yatztjP3sD90Op0s07d1A04PnUh5KJvSV+Qwb2Ouft8aswsJC0tLSmD17Nh6Ph8zMTIYMGeJ2WNEq+qogzMe2bNnAaRt+w8bECYy99Ptuh9NlEvJG0Hz5ExRIGSWPf4n9dU0dn2SiTigUYvXq1UydOpXCwkJOPvnkwybfWbNmMWvWLM4991zy8/N57LHHmD17NpMnT2bYsGGsXbvWheijlyVgF2goRNlTtxAvLfT5/EMx/dAtEtljzmHPpO9zevAj5j/yfauMiEEbNmygoaGBE0888YjHLV++nMTERN544w3+8Y9/cPPNNzN+/Hg++OADLr/8cp577rluijg29Ox/+VFq0bx/ManxfdaM+gZ9BvzPIMAeaeCMb7M57zwuqXiUeXOfdTscc5Ra+3+PlIBDoRAbN27k3nvvPTBSbvLkycyYMePA+3369OmWeGNFxAlYRBJEZLKIzBSRC8N1vQldGVxP1NjURL8PfsYuTz4TPttzux7+hwhDvvQopb4CJiy+g3Vbt7sdkTkKxcXF+P1+RowY0e4x69atY/jw4aSkpABQVFTE1KlTD7xfVFTUYQu6t+kwAYvIOSLyH6ASeB94HvgP8BGwT0T+LSJndW2YPcdH//4tQ9hJ7Vk/xuPrXf9/iT+NxCsfJVuq2P7kzTS2BN0OyUTonnvuoaGhAa+3/bUBly9fzsSJEw+8XrZs2UGvi4uLLQEfot0ELCIDRWQe8CTOAIpLcIYM+3FqeAcD1wA7gGdE5A0RGdy14ca26upKxq7/I2sSTmTU2Z9zOxxXZAw/jZ3jb2J6ywJee+Yht8MxnaioqOigARnLli078LqkpITk5OQDrWPjaLcMTURWAL8A/qUdPDURkTic+t7bVXVcZwfZU8rQ3nv8Lk7f8js2zXyRYRPPcTsc9wRb2PGrM0hu2Enpte9wwjB3Jpg3pgsddxnaqar6z46SLzij3FT1ceCUSKPrbRrraxm15XFW+k/u3ckXwOsj/fOPkCYN7Hj6dgLBkNsRGeOKdhNw21nLInUs5/QWK1/+A9lUoWfe5nYoUSFt4Hi2jfoy05te55U5VppkeqeIqiBEZIiI/FNEikVk/aFbVwcZ6zQUpN+av7Im7gTGTZnhdjhRY+hl97DPm8cJ/72bsqpat8MxpttFWob2N+A0YA7OQ7lDN3MEa999kf5aQs2EryA9fNDF0ZD4ZAKfuo8RsoPF/7zX7XCM6XaRTsZzMnCmqi7tymB6quDiP7OPdE6cfo3boUSdvpMuZcN7j3BmyWOs2ng9Y4fbAznTe0TaHNsAJHXWTcMDOgIick5nXTNaVe/ZzJiaRazMuxi/vwdMNdkF8i+fTbI0sfW5e2yYsulVIk3ANwF/EJFrReQMEZnadjuaG4pIMvAE0H5Fdw+y6fWHAcj/xNddjiR6JReMY3PBJUyve4nFS//rdjjGdJtIE/AJwGjgMeBt4N022ztHec/7gZ1HeU5sUiVn60sU+8YxcuSYjo/vxQZddi8B8dH02t3WCja9RqQJ+B7gUWAcMOSQLeJOOxH5NHABcMvRhRmbyjYspiC4k/KhM90OJerFZ+azbcS1nNX8Du++97bb4RjTLSJNwOnAbFVdrarbDt0iuYCIZAOPAF/FmVeixyt59280q5fhZ33e7VBiwsiZ36OeRFoW/IqQrahseoFIE/CzwPGuOfIQ8JKqvhrJwSJyvYgsEZElZWVlx3lrF4RC5O98hf/Gn8rAgv5uRxMTvClZ7Bp5Nee0vMv7H77vdjjGdLlIE/BW4Oci8o6I/FVEHm67dXSyiHwRmAh8J9LAVPVhVS1U1cK2a1DFioqNi8gKlVM99AK3Q4kpQy68gyaJJ7DgV9YXHIVuvfXWA3P9igg5OTlcddVVlJb22nV1j0ukdcBnAx+Gvx98DPe5DigASkQEPp6o4hUReVxVe1yJwN6PXiBNPQycZP2/RyMuLZd1gz/HmVueYNmKFZw8YYLbIZk2iouLmTJlCvfffz+qykcffcRtt91GIBDgmWeecTu8mBNRAlbVacd5n6uBtkWwfXGqJ74KzD/Oa0el1O1vsMJzAhOHDnQ7lJgz7KLb0P/7OyWv/x4mPOp2OKaN4uJirrrqKiZPngzAlClTWLRoEa+88orLkcWmiJelF5FEYAwQz8FTramqLjrSueHl6NteqzH87S5V7XGfXYKV2ylo2sSyvBs4WSKalc604c8axPqc6ZxROof12/cwcmA/t0MyQGlpKaWlpYwdO/ag/ampqTQ0NLgUVWyLdDKemcBuYDHwHgfXAb/bZdHFqJ2L/wNAyjjr/z1WeZ/6NmnSwPrXHnQ7FBO2YsUKAMaM+bimXVVZtGjRQROxm8hF2gK+G2cAxl1A1fHeVFV3EuGExbGoeeMC9mgfTpo4ye1QYlb68MlsSRzLiTv/SU39naQm+d0OqXO98j0oKXbn3n3Hw4z7jvq01gQ8atQoAoEAu3bt4r777mPlypU8//zznR1lrxBpFcRI4E5VLTrWOuBeQ5W88sWsip9AZkrvWvOt002+kYGyl8Wv28OdaNC6MnJOTg4+n4/BgwczZ84cnnjiCS6++HirVHunSFvAa4H+wOoujKVHCJSsIi20n+r+U9wOJeYNPv1yKhf8gIQVf0cvuhrpSf3px9ACdVtxcTHTpk1j9uzZeDweMjMzGTJkiNthRSQYDB5xQVG3RNoC/hnwp/DgiGnHMxlPT1da5BR1pI451+VIYp/EJbB3yCWc1vIRy1avczucXi0UCrF69WqmTp1KYWEhJ5988mGT76xZs5g1axbnnnsu+fn5PPbYY8yePZvJkyczbNgw1q5d2+495s2bx5QpU5g4cSJjx47lnXc+nmbmG9/4BpdddhlTpkxh0KBB3HXXXUfc3xrLzTffzPTp07niiivYsWMHF1xwAaNHj2bixIksXrz4wLEzZ87k0UedipuXX36Z0047jdrablgkQFU73IDQEbZgJNc4nu2UU07RWLHlDzN124+G6Z6qBrdD6RHqd69R/XGaznngdrdD6dXWrl2rgD711FNHPG7SpEn6zW9+U0OhkL711luakpKic+fOVVXV7373u/qzn/2s3XP37dunoVBIVVXnz5+vF1100YH3Jk+erFdffbUGAgEtLy/XjIwM3bt3b7v7W2P50pe+pIFAQAOBgE6YMEHnzZunqqoLFy7UwsLCA9cvLi7WYcOG6TvvvKPjx4/XkpKSY/tBfSyi3BZpF0RsfM5wmyqZFctZ5DuJ89N72EMjlyT2O4FNSRMYW/IiDU0/IzHB53ZIvVJr/++JJ57Y7jGhUIiNGzfyxhtvHBgpN3nyZGbMmHHg/T59+rR7/rPPPsvjjz9OXV0dNTU1nHrqqQfOW7NmDS+99BJer5c+ffrQt29fKisr292fnZ3NunXrePXVV/F6vcyZM4ehQ4cyffp0AMaNG8fevXsP3HvcuHFMmjSJK664goULF5KXl3fcP7NItNsFISKZrd/rYR686WEewolI+z/d3qBqO+nBSmqybfRWZwpNvJZBUsKSdyOaRsR0geLiYvx+PyNGjGj3mHXr1jF8+HBSUlIAKCoqYurUj3soi4qK2k3gzz77LC+88AJz5syhqKiISy+9lAnhUZBr164lJyeH7OxsACoqKigvL6exsfGw+4cOHcq6desYOXIkmZmZB+Jve+8VK1Ywfvz4A6+3bNnC8uXLiYuLo1+/7qs7P1If8Psi8m0R6XAlDBHJEJE7gV49g8r+jc54lLiBVn7WmYadcQVNxNO47Gm3Q+m17rnnHhoaGo74IGv58uVMnDjxwOtly5Yd9LptErz22msPKl0rKipi0qRJZGZmsnTpUh566KED5y5btozS0lIqKytRVe644w5uvvlmVq5cedj9Pp/vf2Lp378/q1c7NQT79+/nzjvv5Fvf+hYA+/bt4+KLL+bhhx/mkksu4be//W0n/MQic6QEfDowAdgrIk+LyHUicmp4heRh4WWFbhCRp3EmWB8DnNEdQUer/RsW0aDx9B91ituh9CiexDS2Zp3FxJoF7K2scTsc046ioqKDBmQsW7bswOuSkhKSk5MPtI6XLl1KQUHBgWOvvfZannrqKSZNmsTTTz9Ndnb2gQS6dOlSZs2axYUXXsiYMWPIycnh+9//frv7W2Npm4CvvPJKvF4v48aN4/TTT+eGG25g+vTp1NfXM3PmTO666y7OOOMM7rzzTh544AHKy8u7/OcFdPwQDmcljAeBXYQfuoW3ELADeAAYE2mn87FssfIQbtevz9DFPzpVaxpb3A6lx9m96BnVH6fpvBeecDsUc5wqKyt1+vTpER8/bdo0XbJkScT7o0REua3DMjRVXaOqX1fV/sAgnOXpJwEFqjpAVW9UVasPDjSTU7uGzQmjSUmIeIoNE6F+hRdRSzLxa23EVazLyMhg3rx5ER+/atWqg/prO9ofS44qU6jqDpxWrzlU2Rp82kJjrj2A6xJxCWzN/QSn7J1HWeV+cjLT3Y7IdJO21QqR7I8lkQ7EMB2o3roMgOSBEzs40hyrjMLPkiKNFL/7ktuhGNMpLAF3kppty2nQePoNsdWPu0r/iZ+iHj+6do7boRjTKbo1AYtIgYg8IyIVIlIlIv8SkfzujKGryN6VrNMBjOhnH427ivj87OgzlfG171NZ29jxCcZEuW5LwOLMpDIHyASm4Sxz1A+I/c+TqmRUr2OTZzA5NgNal/KPv4hcqaJ48Ztuh2LMcYt0QvZkEfm+iLwkIq+JyLy2W4T3ygPWAF9VZ1rLIuB+4OS2o+5iUvVukoLVVKWO6lkzdkWhAZMuJoCHplUvux2KMcct0iqIB4HPAq8C+47lRqpaAlzZ+lpECoBZwEeqWnks14wWWlKMAMHcsR0ea46PJ7kPm5ImMLR8IcGQ4vXYf3gmdkWagM8DrlHVZzvjpiLyAjATqATO6Yxruql2RxGpQMpAK0HrDk1Dz2Pcyl9SvHY148fYf3omdkXaB+wFijrxvnfhDOh4F3hdRPofekB47uElIrKkrKysE2/d+ep3rWGP9mFwgS0e2R0GnvppAHb/d67LkRhzfCJNwE8Ct0gndXCq6gpVXYzTJeEFvniYYx5W1UJVLczJyemM23YZT/kGNoX6MTIv1e1QeoW0gROo9GSSuONtt0Mx5rhEmoCTga8BO0Rk4bE8hBORPBG5su0+Va0HNuEsdxSbVEmt3cJObwHZVgHRPUTYkzWFsU3LqG5ocjsaY47Z0XRB/BOYD2zGmZin7RaJQcA/RaSwdYeIpAOjiOW15mr34g/VUZ082O1IehXfyE+SJTWsXfqe26EYc8wieginql/qhHstAd4BHhGR64EW4D6gDHi8E67vjn0bAGjObH+iatP5BhTOgPe+Tc3qeXC6rb9nYlPEAzHCcwE/JSKrRGS5iDwpIhHPPK6qIeBSYDnwMrAQqAbOVtVuWP2uawRKncUi4/NGuhxJ7+LPzGdH3CDS9y7u+GBjolRELWARORenBvi/OKPZvDgTtr8rItNVdWEk11HVfcB1xxZqdKrbvQafJtCnny2b190qswsZuWcuVbUNZKQkuh2OMUct0hbwz4EHVHWKqt6hqt9R1cnAH4B7uy686BcoXc9m7cfg7BS3Q+l1UkaeSZo0sGrZIrdDMeaYRJqAJ+CsfHGoh4BePf9ifNUmtmhfBmUlux1Kr1Mwwen7rd1g5WgmNkWagEuAgYfZPxCI2f7b4xZsIbmhhN2efmSnxLsdTa8TnzWIMk8uqXs/cjsUY45JpAn4KeBBEfmEiCSKSJKITAf+BHTK8OSYtH8nHoI0JA+wSXhcUpIxkRGNxTS3BN0OxZijFmkCvgenVnc+Tou3Bueh3GLgjq4JLQZUbgVAMwa5G0dvNmgKObKfjetWuB2JMUctogSsqg2q+hlgHM7w4YuBEar6+fBotl5JK7cBEJdtFRBu6TfubABK19iADBN72i1DE5F8Vd3d+n14dxXwXttjAFqP620aSzcRp15Sc60F7JbswSfSSDy6e6nboRhz1I5UB7xDRPqpaimwE9DDHCPh/d6uCC7aNZVtZq9m0y/TStBc441jl38kfapWuR2JMUftSAn4XKAi/P20bogl9lRtY7vm0j/DBgG4qSHnREZuf5ay/XXkpFs5oIkd7fYBq+pCVQ2EX56Ns3LFwrYbsAxnYvVeyV+7gx2aS74lYFclDCwkUZrZuta6IUxsaTcBi0i2iAwUkYHAj4ETWl+32X8ucEN3BRtVGvfjb6lityePzCSf29H0an1HTwGgerPVA5vYcqQuiBk4s5S19v0e7rdbgH93dlAxIVwBUZdUYDXALkvNP4FakvCVLHM7FGOOSrsJWFWfEJFNOK3kt3G6GiraHoJTDxy7c/kejyonAQfSDjdA0HQrj4ed/pHk1vTOX0UTu45YB6yq76vqu8AQ4AOgVFXfU9X3gAKgRFUjHoIUXhXjcRHZIyJV4SXuxx3Xn8At4UEY3iyrAY4G9VnjGBzcRk19g9uhGBOxSEfC9QXW4SxL1OqnwEoROSmSC4iIB3geGInTmp4K7AfeEJGsiCOOEsHyLezXJDKzct0OxQD+gnH4pYUt61e6HYoxEYs0Af8G+AdwZ5t9J+D0/94f4TUmAFOAL6vqYlVdDVwDpAAXRHiNqNG8bwvbrQIiauSNcFa6Kt9slRAmdkQ0ITtwEnBt2+4GVVURuR+nFC0S24ELcVrSrUI4D/IyI7xG1NDKbezUHKsBjhJ9Bo0jgAfda/3AJnZE2gKuAMYcZv8InAdxHVLVclWdE16aqNUtgB+IaGXlqKFKfO0udmqOtYCjhPgSKfH2J7VqXccHGxMlIk3AfwMeEpFrROSE8HY18CDw5LHcWEQ+A/wCuF9V1xzm/etFZImILCkrKzuWW3Sdun3EhRrZqTn0S/e7HY0Jq0gZQd+mzagebtS8MdEn0gR8N/Af4M/AKpzSs78AzwE/ONqbish1OP3HT9HOdJaq+rCqFqpqYU5OztHeomtVbQegJqEvfl+vnAYjKrVkjWYAeymvqOj4YGOiQKTTUQZU9QYgG5iE0yecoaq3qmrz0dxQRH4A/BWn9XztIV0SseFADfAAlwMxbfkLnIrG3RvsQZyJDREvSx82FhgPbAEGi0ikD/EAEJE7cBbxvEtVb9ZY/ay4fwcAngwbhBFN8oY5FZHVO+xBnIkNkS5Ln47T3TANZwTc28AvgeHhZel3RnCNE3FWV/4L8GcR6dvm7RpVrTva4N2ildup1mT6ZEVZ10gvl1Uwgha8BMvWux2KMRGJtAX8S5xkPQBoXQHjFpyBFL+O8BpX4swb/GVgzyHbrRFeIyoEKrayU7PJz7AHcNFEvD5KvPkk7t/sdijGRCTSLoQLgMtUdVfrxDOqukVEbgJeieQCqvp94PvHFGWUCVZutxrgKFWVNJjsmi1uh2FMRCJtAffh4Il4WjUCvSsLqRJXs5Ndmm01wFGoOWMoBbqH+sZGt0MxpkORJuB3gK+2ea0i4gW+Byzq9KiiWUMlcYF6G4QRpby5o4iXILu32IAME/0iTcC3AV8VkQ+BBOD/gPXAp4HvdlFs0SlcglYiuWQlx7scjDlUav/RAFTtsDXiTPSLtA54Jc5kOvNxhg0340zOM1pVe1fRZXgQRnNKfzwem4g92uQOcWqBm/ZaC9hEv4jreFV1F/DDLowlNlQ5NcCSaTXA0Sg1M5cK0vBVbHQ7FGM61G4CFpH5HH4p+v+hqud1WkTRrmo7tSSSkWk1wNGqzNefxLoOS9ONcd2RWsD2G3wYwcpt7Ahl079PktuhmHbUJhXQb3+R22EY06EjJeBE4CZV3SciZwGLVLWlm+KKWsEKqwGOdsH0QeRVvU59QwNJifb3ZKLXkR7CzcSp/wV4C8jo+nCinCqe6u3s0mxLwFHMlz0Uryh7tm9wOxRjjuhILeBiYIGIrMNZteJ5ETnszGeqem5XBBd16suJa6llu+YxLdMScLRK7jsMgOrdG2HUiS5HY0z7jpSAPwvchNPyPRvYCvTuJWcrnDkGtpFHv3RLwNGqT/+RADSWbXI5EmOOrN0ErKrbCU+WLiIjcPqDq7orsKgUTsDVSQOJjzvamTxNd+nTdxDNGgeVW90OxZgjinQgxjRVrRKR00TkOhFJFZExRzsfcMwr30QID5IxyO1IzBF4vF5KPHn4a7a7HYoxRxRRAhaRNBF5HWfeh0eBHJwpKotFpOBYbiwiD4nII8dyrmsqNlMiOeRmprkdielARUI+aY273A7DmCOK9HP0bMDH8c0HDIA4fgJcfzTnRQOt2MSmYB797QFc1KtLKiAvsAdidNEV0ztEmoAvAG4PD0cGnPmAcR7SfSLSm4nIUOBN4AYgtj4fqqL7NrEllMeQrGS3ozEdCKYNIoV6GmtsgU4Tvbp7PuApwGY+XlcudtSX42muZqv2ZUi2JeBo58ty+unLd1ktsIle3TofsKo+qapfUdWSo4gxOpQ6Cz2u1wJLwDEgOS9cC7zHStFM9Iq0iuE2nEEZ0/h4PuDRQCbwya4ITESuJ9xPPHBgFMw8ttdJwDviBpOTmuByMKYjWQXDAWgsi60PWqZ3idr5gFX1YVUtVNXCnJwomHmsdDU1njRSsvJpXRfPRK/c3L7UaCJaFVuPGkzvYvMBR6p0NRt0AMNyU92OxETAF+dlmyeXhNodbodiTLuO2AIWkTgRuV5E3hSRMhFpEpG9IvKKiFwtvaUpGAyge1dR1FLAmHyrAY4VFb5+pDYFTLU4AAAP+ElEQVTscTsMY9rVbgIWkRScWdD+BHiBp3Dqgf8FpACPA6+LiL8b4nRX2RqkpZ7loWGM6WcJOFbUJeWTFSixWmATtY7UBfFjYDBQqKrLDn1TRCYAc4BvAfd1SXTRYucSAJbpCH5oCThmBFILSK5qQBsqkaQ+HZ9gTDc7UhfEpcC3D5d8AVS1CLgTuOpYbqyq56jqVzs+MgrsXEKtN52G5AFWARFDPH0GA7DfStFMlDpSAu4PLOng/HeBnj0zjSpsfZulOorCwdaKiiX+7CEAVO+xBTpNdDpSAo4H6jo4vx7o2WUBFZuhajvzmsZyqiXgmJKe7wzGaCi1WmATnWxS245sfAOAt0MnMmmIJeBY0jevL9WaRKhym9uhGHNYHdUBf1NEjtQKTunMYKLS6hcp8Q2gzjuA0fYALqZkJcezlhziaqwW2ESnIyXg7cDnI7hGzx1qtH8Xuu09ng19lvNOysPr6R1lzz2FiLAvLo9h9bvdDsWYwzrSkkSDuzGO6LT8SQTlmZYp3DO2r9vRmGNQ4+9PZv0K52FqLxk3ZGKH9QG3p6URFv+Z5f5TaUkbzBnDs92OyByDppQCErUR6m1eYBN9LAG3Z/FDUFfKL6vP4wuTBxHntR9VTMpwZtIL2oM4E4UsqxxOTQm68Fcs9Z/GGv9JXDOlZ5c692QJVgtsopgl4EOFgvDc1wgGWvjO/iv41idGkOb3uR2VOUYpfYcCULd3s8uRGPO/LAG3pQpzb4Mtb3N34Dr6DhnHtVMGux2VOQ55OblUaTItFdYFYaJPxPMB93iBZif5Ln2cx2QmC5LO47krT8JjpWcxrV+Gn+2aQ9r+nlstaWKXJWCAqu3os19Bdi7modBM/p78RZ74ymRy03r+TJs9XZrfR4nkklu3q+ODjelm3doFISJeEfmFiOwRkVoReVZE8rozhoPUV8CCXxL6wySadq3gpuabmd/v6zx34xm28GYPsj+hH+lNe0CVZdsrKatpcjskY4Du7wO+G/gicC1wFlAA/LubY3BavK98j+D9Y2HBz5nXNI6L9H7Gf+pLPDVrik052cPUJxeQoI3sL9/DbX96hi2/ORfd9j6Vdc3MLd5DKGQTtht3dFsXhIjEA98EblHV+eF9VwJbRGSqqr7fpQFUbkXXz6Nu+fMk71lEEA8vBqfwD+9Mzjz7HJ49fQjpiVbt0BOF0gdAJaxbu4qve19iEitpfO0e7kq9j90r3uK0vOfoc+FPmL0xn+yUBL4ysgmyR1BeHyAzKd6eA5gu0519wCfhTF25oHWHqm4Vka3AmUDnJuDGatjxIQ3r3iC47jVSajYjwN5QP14OXcKK3M8wbdLJPD6xPykJ1hXek3lzRsJWKF37Aed5nSmu/bs/YC1r+anvabKq1lA/507+tOcuPuVZzFde/3/sPeFaphSdzxdOzuGnaS8SyBnNk41nMLpfms2K15sEA+DxdtkwdtFuWi9LRC7F6W6IV9WWNvvfA5ap6k3tnZuamqqnnHJKxPdq2LedxFrnqXcIoUaT2C+phPyZJCenkJHkw2cj23qNirpmUsv+SwgPCbSwW3LJ11J2ajYFso8W4vARYEVoKMM9u0jC6SNeGhpBX6kgX8oBWBUaRJ0kc1I/P/GhekhIB699auoRNAQtDWhzPYGmOrS5HgnUExdsItS/EK/v6LolFyxYEFHG7s6mXxIQapt8w5qA/yk3EJHrgesBEhKO7g/f6EmkWnIIxKfhSUwj2R/PgASfzcXSSyXFe6nRJPpIDSA0+PNoaqikQPYBsD5UwBjPNk7w7iRemynx5JEX2ssIXxnJwf3U+zLwttQxzLuXzaF++ErWAEoAL2UpJ9AvzUdTYwO+1Gw89ksWnYItEGwm2NJEsLkBDTRCSxMSbMIbasKrAQAEJyk2EU+DxtPsSSE9GOqy/2e7MwE3AB4RiVMN/2kdCRxm5Q1VfRh4GKCwsFAXLFgQ8Y1UFbF/CCYsGFJ+/qOb+JHv79SkjeDFqf8m8PJtXBc3j0C/Uxi+5Tvc4ft/XOBdTKV/KIVVP+d23x+50PsBLaRxVuP9jPds5uH43wKl7NSB3NVyHb/xPYiPEhIkiI8W3vaP5KTQKnaRS+jUrxFa9xo5595Exfr3GDjlMmrKdpA7/GQ0GCDOnwpe6/qKWLAFmmqguZZgYw31NVU01O6nsW4/TXX7aWmoJthQAw0VxNXvI76pnMTmCpIDFaQE9+MldNDlmtXLLs1mFzlU+vpSn5hPS8YQJOcE0gpOYGi/LIZmp5AY7+3SP1Z3/ga0zordr833APlApxZpWvI1bXk9wub+FzF/z2omn/9DTs/N5prgpzktdR+jp9/FiXPj+L9dlzI9v4X6U28n+Cy8mnUNF+ouVhV8nj3/zWLU8FGQsIrQ1neZVXsrq3QIlVN/wNBFd1Km6WyVAs5qfJN6TWAY+4h/7xbn5k/PoS/A8ntJAWo0kSQaKYkvIBBUqhPzIdhMTXwuWVpFtS+X1KR4quJySElOxhNqIejPxO9VGuLSiPMl4I9PIJCUS6KnhZS0TGolhcSkZHyhekLeROK9QnxSuvOJLxSEhFQEQIMQF/6w2bbrUUNOPyc4fZ3BQLjPU5yvIs7xGoJQAILNztfwvmAwSEugBQ00OS3MlhaCwQChQAuhQCOhlmYCLY0EmpsItTQQbG4k1FSDNtVCUy3SXIc01+IN1OFtqSUuUEdcoB5fsI74YD3+UAM+Pv7g7MV5mHS4tdDqNYF9pFEpGez29KE2bjiNiX0IJeXgSc3Bm1GAP2coGbkF5KUnMik9kfg497oju7MPOAEoA25U1b+H9w0GtgBTVPWD9s4tLCzUJUs6Wh/UmPZV1DVT1xRgQJ8kAD7cXM6EARn4fV5qmwKUVjcyNMdZ4GXNnmryMxJJT/Shqry5tpTCwX1I94Ug0MT8zQ00B0JcML4vG9/+J8/tyeJrnzqV9a/8gRWJkxmR0kjZ8rkkjTmP5KUPsyP7LAaVvklzUh55taup9PUlt3kHVaQwXHaBeMimii2hPAZJKUE8eAnhkc7/txlUwRu+but92n7fdl93CKpQRyK1+KnTROrw04CfBk8SLd4kAnHJBH3JaHwKxKcg/lS8/lTi/GnEJ6eRkJRGYko6SSkZJKelk5aaSkp8XDRUrkQUQLclYAARuQ+4LryVAg8Ajap6zpHOswRseopAMESc10MgGEJEqKhrJis5nuamBnbVhshI8LBn3z7S4pTyqv2E4lNpqi6jOiCkBaudxmhzA57GcmpaPASb6kinjqZAEPXE4wk20BIMQXO9cyxCXKABFQ8qQlywkZYQeDxeNBQkhMf5xBhsIShevAQJeeIBwUMIVSXOIyCgeAiJl6D4UI8X8XjxejyIJw6Px4N641FvPOKNQzxxiMcHvnjEG4/Hl0BcfCJeXwJxPj++5DR8iWkk+JPxx8fh93nw+7zEez3RkDw7Q9Q9hAP4IeAD/h7++irwjW6OwRjXtM4r3fq1ddCPPzGJYYnOMVlpzhzGAwe2njWkO0M03ahbE3D44dt3wpsxxvRqVgxrjDEusQRsjDEusQRsjDEusQRsjDEu6dYytGMlImXA0a4pkw3s64JwjDG9x7HmkX2qen5HB8VEAj4WIrJEVQvdjsMYE7u6Oo9YF4QxxrjEErAxxrikJyfgh90OwBgT87o0j/TYPmBjjIl2PbkFbIwxUc0SsDHGuKTHJWAR8YrIL0Rkj4jUisizIpLndlzGmNghInki8ng4j1SJyGsiMq6z79PjEjBwN/BF4FrgLKAAZzFQY4zpkIh4gOeBkcBMYCqwH3hDRLI69V496SGciMTjjFq5RVUfC+8bjLPqxumq+r5rwRljYoKITASWAmNUdU14XwJQAdygqn/rrHv1tBbwSThLRS1o3aGqW4GtwJmuRGSMiTXbgQuBdW32hXBWucjszBv1tGVZC8JfD13kczcwoJtjMcbEIFUtB+YcsvsWwA/M68x79bQWcBIQUtWWQ/Y34fzwjDHmqIjIZ4BfAPe3dkl0lp6WgBsAj4gc2rJPAOpciMcYE8NE5Dqch/hPAXd09vV7WgLeEf7a75D9+fxvt4QxxrRLRH4A/BV4ELhWVUOdfY+eloCLgBrg7NYd4SqIwcDbrkRkjIk5InIHcC9wl6rerF1ULtajytAAROQ+4LrwVgo8ADSq6jnuRWWMiRUiciJOGdrjwA8OebtGVTutO7OntYABfgg8CfwdeAtnJY3PuhqRMSaWXAl4gS8Dew7Zbu3MG/W4FrAxxsSKntgCNsaYmGAJ2BhjXGIJ2BhjXGIJ2BhjXGIJ2BhjXGIJ2BhjXGIJ2EQNERkjIhe0eb1VRH7oQhxnisjCTrrW+SLyUmdcy/Q8loBNNHkROLXN61OB33ZnACLiBx4Bbu+M66nqq4BfRL7QGdczPYslYBNNpO0LVS3rzGGfEboO2Kuqizvxmr8B7hURbyde0/QAloBNVBCRBcAw4McisjW870AXhIjcLSKvisgdIlIqIjUi8oCIDBSRuSJSLyLrReT8NtdMEJHfhBdWrBaRhSIyuYNQbgWebnON60QkcEisB+0Lv14jIk0isk1E7gmvK9bqDSAduOSYfjimx7IEbKLFpThLR/2Gg7sh2poGnAicgbNCwQ3Ahzhzf5yCs4TMY22O/xvOwqxXAIXAm8BbIjLycBcXkVE4CzHOjTTo8MQtD+FM2jIC+BZO98XVrceEFwiYB3wm0uua3qGnLUlkYpSqVohIEKhV1bJ2DhNgVrhbYr2IzAbmqeqTACLyADBXRHJwWpxXAONUdVX4/HtE5AzgO8Csw1z/NJxJ/bccRejDAAW2qep2YLuIfBLYechxK3FW6jbmAEvAJpbsOaRPuA7Y1OZ1Q/hrAjAx/P2HIgd1LSeEt8PJA8qPcu7XV3Fa4UtEZCPwGvBUOBm3VRa+vjEHWAI2seTQtf7AWa32cJrDX6fwcWJu1dTOOUpk3XIH/t2oagNwtogUAjOA84EbReR7qjr7kHOCEVzb9CLWB2yiSWfOjdra7ZCnqhtbN5yHbDPbOWcPkC0HN5mbAa+IJLXZN6L1GxGZLiI/UtUlqvpTVT0dZwmb6w65dnb4+sYcYAnYRJMaYKSI5B/vhcLJ9ingYRGZISLDRORnwNeB9la2XQzEA2Pa7PsA5z+Ge0RksIh8joOTazNO5cYtIjJERKbgPCz88JBrn3yYfaaXswRsosn9OB/jVxxSxnWsvopT0fBXnIdgM4BLVfWNwx2sqhuAtTgJtHXfZpxqi8vD732NNoM0VHUhzsoJ1wOrgReAhThVGgCEV+k+HWegiTEH2IoYxrQhIjcCX1LV9krhjuWalwC/BkapaqCj403vYS1gYw72KJAlIlM78Zq3APdY8jWHsgRsTBuq2oTTpfDLzrieiHwaaFDVv3XG9UzPYl0QxhjjEmsBG2OMSywBG2OMSywBG2OMSywBG2OMSywBG2OMS/4/SdU9LIJAGegAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x20efe343278>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, Z1, Z2 = plotZprofiles(bls, f, A, Qm)\n",
    "error_Z = rmse(Z1, Z2)\n",
    "r2_Z = rsquared(Z1, Z2)\n",
    "print('Z-error: R2 = {:.4f}, RMSE = {:.4f} nm ({:.4f}% dZ)'.format(\n",
    "    r2_Z, error_Z, error_Z / (Z1.max() - Z1.min()) * 1e2))"
   ]
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}