make_ACTN2_Affimer_PERPL_plots.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Make ACTN2 Affimer PERPL analysis plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import modelling_general\n",
    "import plotting\n",
    "import zdisk_modelling\n",
    "import zdisk_plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set average estimated localisation precision for Affimer and PALM data.\n",
    "This is the mean after filtering for localisation precision < 5 nm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "af647_precision = 3.1 # Mean value after filtering for precision < 5 nm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set standard maximum distance over which to plot distances and fit models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitlength = 100."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Choose Affimer relative position data and combine files if necessary\n",
    "### Insert paths to the Affimer relative position data here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "actn_affimer_relpos_path = r\"..\\data-perpl\\ACTN2-Affimer_PERPL-relpos_200.0filter_len2440488.pkl\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select desired datasets to combine here:\n",
    "These contain all data attributes, including the relative positions in both directions for each pair of localisations. We will filter e.g. for cell-axial data later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['..\\\\data-perpl\\\\ACTN2-Affimer_PERPL-relpos_200.0filter_len2440488.pkl']\n"
     ]
    }
   ],
   "source": [
    "path_list = [actn_affimer_relpos_path]\n",
    "print(path_list)\n",
    "relpos = zdisk_modelling.read_relpos_from_pickles(path_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data attributes and number of data points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "axial          49.426791\n",
       "perpaxis2D     34.358846\n",
       "tocoverslip    15.602947\n",
       "z_nzcorr       16.072208\n",
       "transverse     37.932126\n",
       "Name: 0, dtype: float64"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "relpos.iloc[0, :] # This shows the first relative position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2440488"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(relpos) # This shows how many relative positions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get the axial (X) distances, without duplicates\n",
    "The YZ-distance limit for pairs of localisations to include can be set here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is the YZ-distance limit for X-distances to include:\n",
    "transverse_limit = 10.\n",
    "\n",
    "axial_distances = zdisk_modelling.getaxialseparations_no_smoothing(\n",
    "    relpos,\n",
    "    max_distance=relpos.axial.max(),\n",
    "    transverse_limit=transverse_limit\n",
    "    )\n",
    "axial_distances = zdisk_modelling.remove_duplicates(axial_distances)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get the 1-nm bin histogram data\n",
    "Up to distance = fitlength"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist_values, bin_edges = zdisk_plots.plot_distance_hist(\n",
    "    axial_distances,\n",
    "    fitlength\n",
    "    )\n",
    "bin_centres = (bin_edges[0:(len(bin_edges) - 1)]\n",
    "               + bin_edges[1:]\n",
    "               ) / 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get the KDE data\n",
    "Estimate every 1 nm, with kernel size based on localisation precision estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kde_x_values, kde = zdisk_plots.plot_distance_kde(\n",
    "    axial_distances,\n",
    "    af647_precision,\n",
    "    100.\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate the axial RPD with smoothing for Churchman 1D function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calculation_points = np.arange(fitlength + 1.)\n",
    "axial_rpd = plotting.estimate_rpd_churchman_1d(\n",
    "    input_distances=axial_distances,\n",
    "    calculation_points=calculation_points,\n",
    "    combined_precision=(np.sqrt(2) * af647_precision)\n",
    ")\n",
    "plt.plot(calculation_points, axial_rpd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Choose axial model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "axial_model_with_info = zdisk_modelling.set_up_model_4_peaks_fixed_ratio_with_fit_settings()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit model to histogram bin values, at bin centres"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(params_optimised,\n",
    " params_covar,\n",
    " params_1sd_error) = zdisk_modelling.fitmodel_to_hist(\n",
    "    bin_centres,\n",
    "    hist_values,\n",
    "    axial_model_with_info.model_rpd,\n",
    "    axial_model_with_info.initial_params,\n",
    "    axial_model_with_info.param_bounds,\n",
    "    )\n",
    "print('')\n",
    "print('Initial parameter guesses:')\n",
    "print(axial_model_with_info.initial_params)\n",
    "print('')\n",
    "print('Parameter bounds:')\n",
    "print(axial_model_with_info.param_bounds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot fitted model over histogram data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'axial_distances' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-32-7787b3f4a424>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      1\u001b[0m fig, axes = zdisk_plots.plot_distance_hist_and_fit(\n\u001b[1;32m----> 2\u001b[1;33m     \u001b[0maxial_distances\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      3\u001b[0m     \u001b[0mfitlength\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m     \u001b[0mparams_optimised\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m     \u001b[0mparams_covar\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'axial_distances' is not defined"
     ]
    }
   ],
   "source": [
    "fig, axes = zdisk_plots.plot_distance_hist_and_fit(\n",
    "    axial_distances,\n",
    "    fitlength,\n",
    "    params_optimised,\n",
    "    params_covar,\n",
    "    axial_model_with_info\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot fitted model over histogram data, with confidence intervals on the model\n",
    "### NOTE: IT TAKES A WHILE TO CALCULATE THE CONFIDENCE INTERVALS\n",
    "### Skip this if you don't need it right now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zdisk_plots.plot_distance_hist_and_fit(\n",
    "    axial_distances,\n",
    "    fitlength,\n",
    "    params_optimised,\n",
    "    params_covar,\n",
    "    axial_model_with_info,\n",
    "    plot_95ci=True\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Akaike weights for the models\n",
    "Typed in AICc values for the different models here, to obtain relative likelihood, summing to one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from modelstats import akaike_weights\n",
    "weights = akaike_weights([\n",
    "    364.35,\n",
    "    364.51,\n",
    "    366.38,\n",
    "    370.28,\n",
    "    374.67\n",
    "])\n",
    "print(weights)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot model components for best model (4 peaks with fixed peak ratio)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zdisk_plots.plot_model_components_4peaks_fixed_peak_ratio(\n",
    "    fitlength,\n",
    "    *params_optimised)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Transverse distances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get the transverse (YZ) distances, without duplicates\n",
    "The X-distance limit for pairs of localisations to include can be set here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2440488, 5)\n"
     ]
    }
   ],
   "source": [
    "# This is the YZ-distance limit for X-distances to include:\n",
    "axial_limit = 10.\n",
    "print(relpos.shape)\n",
    "\n",
    "trans_distances = zdisk_modelling.get_transverse_separations(\n",
    "    relpos,\n",
    "    max_distance=relpos.transverse.max(),\n",
    "    axial_limit=axial_limit\n",
    "    )\n",
    "trans_distances = zdisk_modelling.remove_duplicates(trans_distances)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Choose analysis lengthscale for transverse distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitlength = 50."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAPEklEQVR4nO3dUYhc133H8e8vqmuHJKV2vRKKpFYuKFA5NA4INeA+OHYau7GJTMFFgRQVDHpxqQMpsZSX0IBAeTF5aR5EEiJIHHdp4lq4tImsxLiBYHmVuHVk2VjEqr1IaNdOQ+IXtVL+fZhrOpZ2vbPamd2ds98PiHvvmXtn/gfhn47PPXMnVYUkqS3vWukCJEnDZ7hLUoMMd0lqkOEuSQ0y3CWpQb+10gUA3HjjjbV169aVLkOSxsqJEyder6qJuV5bFeG+detWpqamVroMSRorSf5rvteclpGkBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAatim+oSlIrtu77lznbzxy8e1nrcOQuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUoIHCPcmZJM8neS7JVNd2Q5KjSV7uttf3nb8/yekkLyW5c1TFS5LmtphH/n60ql7vO94HHKuqg0n2dccPJdkO7AZuBt4PPJnkA1V1aWhVS9KYWe5HAS9lWmYXcLjbPwzc29f+aFVdqKpXgNPAziV8jiRpkQYN9wK+n+REkr1d24aqOgfQbdd37ZuA1/qune7a3ibJ3iRTSaZmZ2evrnpJ0pwGnZa5tarOJlkPHE3y4jucmzna6oqGqkPAIYAdO3Zc8bok6eoNNHKvqrPddgZ4jN40y/kkGwG67Ux3+jSwpe/yzcDZYRUsSVrYgiP3JO8B3lVVv+72Pw58ETgC7AEOdtvHu0uOAI8keZjeDdVtwPER1C5JI7dafhN1sQaZltkAPJbkrfMfqap/S/IsMJnkfuBV4D6AqjqZZBJ4AbgIPOBKGUlaXguGe1X9HPjQHO1vAHfMc80B4MCSq5MkXZXFrHOXJHXmm65ZLQx3Sc0Z13nyYfLZMpLUIMNdkhrktIykNWO1z5MPkyN3SWqQ4S5JDTLcJalBhrskNcgbqpJWDdenD48jd0lqkOEuSQ0y3CWpQYa7JDXIG6qSVj1vtC6e4S5pbK2lxwksltMyktQgw12SGuS0jKRl53TK6Dlyl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQa5zl7RkPvtl9THcpTXMUG6X0zKS1KCBwz3JuiQ/TfJEd3xDkqNJXu621/eduz/J6SQvJblzFIVLkua3mJH7g8CpvuN9wLGq2gYc645Jsh3YDdwM3AV8Jcm64ZQrSRrEQOGeZDNwN/DVvuZdwOFu/zBwb1/7o1V1oapeAU4DO4dTriRpEIOO3L8MfA74TV/bhqo6B9Bt13ftm4DX+s6b7treJsneJFNJpmZnZxdduCRpfguGe5J7gJmqOjHge2aOtrqioepQVe2oqh0TExMDvrUkaRCDLIW8Ffhkkk8A1wG/k+SbwPkkG6vqXJKNwEx3/jSwpe/6zcDZYRYtSXpnC47cq2p/VW2uqq30bpT+oKo+DRwB9nSn7QEe7/aPALuTXJvkJmAbcHzolUuS5rWULzEdBCaT3A+8CtwHUFUnk0wCLwAXgQeq6tKSK5U0dvzFpZWzqHCvqqeAp7r9N4A75jnvAHBgibVJkq6S31CVpAYZ7pLUIMNdkhpkuEtSgwx3SWqQz3OXdAWf8z7+HLlLUoMMd0lqkNMykgbmN07HhyN3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CCXQkrLZKW+9enyxbXJkbskNchwl6QGGe6S1CDn3KVVyiczaikcuUtSgwx3SWqQ4S5JDTLcJalB3lCVGuGXldTPcJc6rk5RS5yWkaQGGe6S1CCnZaQx49y6BrFguCe5DngauLY7/5+q6gtJbgD+EdgKnAH+sqr+u7tmP3A/cAn426r63kiql1ahxYavYa1RGGTkfgG4vareTHIN8KMk/wr8BXCsqg4m2QfsAx5Ksh3YDdwMvB94MskHqurSiPogrQhDWavZgnPu1fNmd3hN96eAXcDhrv0wcG+3vwt4tKouVNUrwGlg51CrliS9o4FuqCZZl+Q5YAY4WlXPABuq6hxAt13fnb4JeK3v8umu7fL33JtkKsnU7OzsUvogSbrMQOFeVZeq6hZgM7AzyQff4fTM9RZzvOehqtpRVTsmJiYGq1aSNJBFLYWsql8CTwF3AeeTbATotjPdadPAlr7LNgNnl1ypJGlgg6yWmQD+t6p+meTdwMeALwFHgD3AwW77eHfJEeCRJA/Tu6G6DTg+gtqlq+KNUK0Fg6yW2QgcTrKO3kh/sqqeSPJjYDLJ/cCrwH0AVXUyySTwAnAReMCVMpK0vBYM96r6T+DDc7S/AdwxzzUHgANLrk6SdFV8/IAkNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkM9zlxbgl540jhy5S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAa5FFJjYb7liGcO3r3MlUjjwZG7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapBLITXWXCIpzc2RuyQ1yHCXpAYZ7pLUIMNdkhrkDVWtKv6knTQcjtwlqUGO3NUk/w9Aa50jd0lq0ILhnmRLkh8mOZXkZJIHu/YbkhxN8nK3vb7vmv1JTid5Kcmdo+yAJOlKg0zLXAQ+W1U/SfI+4ESSo8BfA8eq6mCSfcA+4KEk24HdwM3A+4Enk3ygqi6NpgtazfwGqbQyFgz3qjoHnOv2f53kFLAJ2AXc1p12GHgKeKhrf7SqLgCvJDkN7AR+POziNb6cE5dGa1Fz7km2Ah8GngE2dMH/1j8A67vTNgGv9V023bVd/l57k0wlmZqdnV185ZKkeQ0c7kneC3wH+ExV/eqdTp2jra5oqDpUVTuqasfExMSgZUiSBjBQuCe5hl6wf6uqvts1n0+ysXt9IzDTtU8DW/ou3wycHU65kqRBDLJaJsDXgFNV9XDfS0eAPd3+HuDxvvbdSa5NchOwDTg+vJIlSQsZZLXMrcBfAc8nea5r+zxwEJhMcj/wKnAfQFWdTDIJvEBvpc0DrpRpnzdIpdVlkNUyP2LueXSAO+a55gBwYAl1SZKWwG+oSlKDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIH+JSYvil5Wk8eDIXZIaZLhLUoMMd0lqkHPua5w/gye1yZG7JDXIcJekBjkts4qt5JSJSx6l8Wa4N8Y5dEngtIwkNcmR+xrhNIu0tjhyl6QGGe6S1CDDXZIaZLhLUoMMd0lqkKtlxpArXyQtxJG7JDXIcJekBhnuktQgw12SGuQN1VXAG6SShm3BkXuSryeZSfKzvrYbkhxN8nK3vb7vtf1JTid5KcmdoypckjS/QaZlvgHcdVnbPuBYVW0DjnXHJNkO7AZu7q75SpJ1Q6tWkjSQBcO9qp4GfnFZ8y7gcLd/GLi3r/3RqrpQVa8Ap4GdQ6pVkjSgq72huqGqzgF02/Vd+ybgtb7zpru2KyTZm2QqydTs7OxVliFJmsuwV8tkjraa68SqOlRVO6pqx8TExJDLkKS17WrD/XySjQDddqZrnwa29J23GTh79eVJkq7G1Yb7EWBPt78HeLyvfXeSa5PcBGwDji+tREnSYi24zj3Jt4HbgBuTTANfAA4Ck0nuB14F7gOoqpNJJoEXgIvAA1V1aUS1S5LmsWC4V9Wn5nnpjnnOPwAcWEpRkqSl8fEDktQgw12SGmS4S1KDDHdJapDhLkkN8pG/IzDfI3zPHLx7mSuRtFY5cpekBjlyX0b+KIek5eLIXZIaZLhLUoMMd0lqkOEuSQ3yhuoAvBEqadw4cpekBhnuktQgp2X6OP0iqRWO3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJatCaXOfuenZJrXPkLkkNMtwlqUFNT8s4/SJprXLkLkkNMtwlqUEjC/ckdyV5KcnpJPtG9TmSpCuNZM49yTrgH4A/A6aBZ5McqaoXRvF5zq1L0tuNauS+EzhdVT+vqv8BHgV2jeizJEmXGdVqmU3Aa33H08Cf9J+QZC+wtzt8M8lLS/i8G4HXl3D9uFlr/QX7vFasuT7nS0vq8x/M98Kowj1ztNXbDqoOAYeG8mHJVFXtGMZ7jYO11l+wz2uFfR6eUU3LTANb+o43A2dH9FmSpMuMKtyfBbYluSnJbwO7gSMj+ixJ0mVGMi1TVReT/A3wPWAd8PWqOjmKz+oMZXpnjKy1/oJ9Xivs85CkqhY+S5I0VvyGqiQ1yHCXpAaNdbivhUccJPl6kpkkP+truyHJ0SQvd9vrV7LGYUuyJckPk5xKcjLJg117s/1Ocl2S40n+o+vz33ftzfYZet9mT/LTJE90x63390yS55M8l2SqaxtJn8c23PsecfDnwHbgU0m2r2xVI/EN4K7L2vYBx6pqG3CsO27JReCzVfVHwEeAB7q/25b7fQG4vao+BNwC3JXkI7TdZ4AHgVN9x633F+CjVXVL39r2kfR5bMOdNfKIg6p6GvjFZc27gMPd/mHg3mUtasSq6lxV/aTb/zW9//g30XC/q+fN7vCa7k/RcJ+TbAbuBr7a19xsf9/BSPo8zuE+1yMONq1QLcttQ1Wdg14QAutXuJ6RSbIV+DDwDI33u5uieA6YAY5WVet9/jLwOeA3fW0t9xd6/2B/P8mJ7hEsMKI+j/MvMS34iAONtyTvBb4DfKaqfpXM9Vfejqq6BNyS5HeBx5J8cKVrGpUk9wAzVXUiyW0rXc8yurWqziZZDxxN8uKoPmicR+5r+REH55NsBOi2Mytcz9AluYZesH+rqr7bNTffb4Cq+iXwFL17La32+Vbgk0nO0JtSvT3JN2m3vwBU1dluOwM8Rm96eSR9HudwX8uPODgC7On29wCPr2AtQ5feEP1rwKmqerjvpWb7nWSiG7GT5N3Ax4AXabTPVbW/qjZX1VZ6/+3+oKo+TaP9BUjyniTve2sf+DjwM0bU57H+hmqST9Cbt3vrEQcHVrikoUvybeA2eo9CPQ98AfhnYBL4feBV4L6quvym69hK8qfAvwPP8//zsZ+nN+/eZL+T/DG9m2nr6A26Jqvqi0l+j0b7/JZuWubvquqelvub5A/pjdahNyX+SFUdGFWfxzrcJUlzG+dpGUnSPAx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1KD/A1CjXHUDWDxnAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "hist_1nm_bins = plt.hist(trans_distances, bins=np.arange(fitlength + 1.))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimate RPD using Churchman's function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x258bde32f08>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXhV1b3/8fcXCPMYhhiGEGYIElACqNQWHCqOqIjFaosVf9ReW7VOIA5ULYporbZXb4uWlmqtBhkdARGKA4igkJEEJAwhIYwJQ8h41v0jx/uLiBIgJ/sMn9fz8OScdfbJ+a6H5MNi7bXXNuccIiISXup5XYCIiNQ+hbuISBhSuIuIhCGFu4hIGFK4i4iEoQZeFwDQrl07Fx8f73UZIiIhZd26dXudc+2P91pQhHt8fDxr1671ugwRkZBiZtu+6zVNy4iIhCGFu4hIGFK4i4iEIYW7iEgYUriLiIQhhbuISBhSuIuIhCGFu4iIB3w+xxufb2dpRkFAvn9QXMQkIhJJNu8+xJR5aazZup8rEmO5OCGm1j9D4S4iUkdKyit5ccVX/M+KzTRt2IAZYxK5bnDngHyWwl1EpA6s+mofD85PZcveI1w9qCMPXZFAu+aNAvZ5CncRkQA6cKSMJ97NZM66XOKim/LPW4byw97H3eurVincRUQCwDnHog15PPZWBkVHy/nViB7ccUEvmjSsXyefr3AXEalluQeKeWhBGiuy9jCwcytevXUY/WJb1mkNCncRkVpS6XP849Ot/GFJFgCPXJHA+PPiqV/P6rwWhbuISC3IyDvIA/NS2JBbxMg+7Xn86jPp3KapZ/Uo3EVETkNJeSXPL9vEzJVbaN0kiufHDeKqgR0xq/vRenU1Cncz2wocAiqBCudckplFA28A8cBW4Hrn3AH/8Q8AE/zH3+GcW1zrlYuIeOyTzXuZMj+VbfuKGTu4Mw9e3o/WTRt6XRZwciP3kc65vdWeTwaWOeemm9lk//NJZpYAjAP6Ax2BD8yst3OustaqFhHx0IEjZUx7N5M31+US37Ypr906jPN6tvO6rG84nWmZ0cAI/+PZwApgkr/9dedcKZBjZpuBocCq0/gsERHPHbu88b9G9OCOC3vROKpuljeejJqGuwOWmJkD/uqcmwnEOOfyAZxz+WbWwX9sJ2B1tffm+tu+wcwmAhMB4uLiTrF8EZG6sX1fMQ8uSOWjTXsZ2KU1r147oM6XN56Mmob7cOdcnj/Al5rZxu859nhnEdy3Gqr+gZgJkJSU9K3XRUSCQXmlj1kf5/DHD7Kpb8ajV/XnpnO6erK88WTUKNydc3n+r7vNbD5V0ywFZhbrH7XHArv9h+cCXaq9vTOQV4s1i4jUiQ07Cpk8L5XM/INcnBDDo1f1p2PrJl6XVSMn3M/dzJqZWYuvHwM/BtKARcB4/2HjgYX+x4uAcWbWyMy6Ab2ANbVduIhIoBwureB3i9K55sVP2H+klL/cNJiXfp4UMsEONRu5xwDz/Ws2GwCvOefeN7PPgWQzmwBsB8YCOOfSzSwZyAAqgNu1UkZEQsXi9F1MXZhOwaESbhrWlftG9aFl4yivyzppJwx359wWYOBx2vcBF37He6YB0067OhGROpJfdJSpC9NZklFA3zNa8OJNZ3N2XBuvyzplukJVRCJapc/xyqqtPLMkmwqfj8mX9mXCD7oRVT+070KqcBeRiJW2s4gp81NJyS3iR73b8/urz6RLtHf7wdQmhbuIRJwjpRX8cWk2sz7JIbpZI/50w1lcmRjr+X4wtUnhLiIR5YOMAh5ZmEZeUQk3Dovj/lF9adUk9E6YnojCXUQiwq6iEn63KJ3303fRJ6YFc396FoO7RntdVsAo3EUkrFU/YVpe6WPSqL7cen7onzA9EYW7iISt9LwipsxPY8OOQs7v1Y5pVw8grm14nDA9EYW7iISd4rIKnvtgE3/7OIc2TYPnBhp1SeEuImFledZuHpqfxs7Co4wb0oXJl/YNmhto1CWFu4iEhT2HSnns7Qze2pBHj/bNSP7luQztFr4nTE9E4S4iIc05R/LaHUx7J5OSch+/vag3t43oTqMGwXcDjbqkcBeRkLVlz2EemJfKZzn7GRofzRPXDqBnh+ZelxUUFO4iEnLKKnz89T9f8eflm2ncoB7Trx3A9UldqBfkN9CoSwp3EQkp63cUMnluCht3HeLyAbFMvSqBDi0ae11W0FG4i0hIKC6r4A9Lsvn7Jzl0aNGYl36exMUJMV6XFbQU7iIS9FZm72HK/FRyDxzlpnOq9oMJxRto1CWFu4gErcLiMh5/O5O5X+TSvX0z5tx2LkPiI3d548lQuItIUHovNZ+HF6ZTWFzGr0f25NcX9KRxVGQvbzwZCncRCSq7D5UwdWE676Xt4sxOLZl9yxD6d2zldVkhR+EuIkHBOce8L3by2NsZHC2v5P5RfZh4fncahPnujYGicBcRz+UVHmXK/FRWZO1hcNc2PDUmURcjnSaFu4h4xjnHG59XbR1Q4XNMvTKBn58bT31djHTaFO4i4oncA8U8MC+Vjzbt5Zzu0cwYMzBi9lqvCwp3EalTPp/jtTXbefLdTAAev/pMbhwap60DapnCXUTqzI79xUyam8KnX+3jBz3b8eS1A+gSrdF6ICjcRSTgfD7Hv/yj9XpmPHHNAG4Y2iWi7oxU1xTuIhJQO/YXc/+bKazaso/ze7Vj+phEOrVu4nVZYU/hLiIBcexo/clrBzBuiEbrdaXG4W5m9YG1wE7n3BVmFg28AcQDW4HrnXMH/Mc+AEwAKoE7nHOLa7luEQli1efWNVr3xslc+nUnkFnt+WRgmXOuF7DM/xwzSwDGAf2BUcCL/n8YRCTM+XyOV1Zv45LnVpKSW8QT1wzgn7cMVbB7oEbhbmadgcuBl6s1jwZm+x/PBq6u1v66c67UOZcDbAaG1k65IhKsduwv5saXP+PhBWkM7tqGxb/9IT8dFqdpGI/UdFrmOeB+oEW1thjnXD6Acy7fzDr42zsBq6sdl+tv+wYzmwhMBIiLizvJskUkWGhuPTidMNzN7Apgt3NunZmNqMH3PN7fqPtWg3MzgZkASUlJ33pdRILf9n1Vc+taCRN8ajJyHw5cZWaXAY2Blmb2KlBgZrH+UXsssNt/fC7Qpdr7OwN5tVm0iHjL53P8c9VWnno/i/r1NFoPRiecc3fOPeCc6+yci6fqROmHzrmbgEXAeP9h44GF/seLgHFm1sjMugG9gDW1XrmIeCJn7xF+MnMVv3srg2Hdo1ny2x9yw1DNrQeb01nnPh1INrMJwHZgLIBzLt3MkoEMoAK43TlXedqVioinKn2OWR/n8MySLBo1qMczYwcy5uxOCvUgZc55P92dlJTk1q5d63UZIvIdsgsOcf+bKazfUchF/WKYds2ZxLRs7HVZEc/M1jnnko73mq5QFZHvVFpRyYvLv+LFFZtp3qgBz48bxFUDO2q0HgIU7iJyXOu2HWDy3BQ27T7M6EEdeeSKBNo2b+R1WVJDCncR+YYjpRU8vTiL2au2EtuyMbNuTuKCvjFelyUnSeEuIv/nw40FPLwgnZ2FR/n5uV25f1RfmjdSTIQi/a2JCAUHS3j0rXTeTd1Fzw7NmXPbuQyJj/a6LDkNCneRCFbpc7y6ehtPL86irNLHvT/uzcQf9qBhg5PZU1CCkcJdJEKl5xUxZV4qG3KLOL9XOx4ffSbx7Zp5XZbUEoW7SIQpOlrOH5dm88rqbbRpGqXljWFK4S4SIXw+x7wvdzL9vUz2HSnjp0PjuP+SvrRqGuV1aRIACneRCJCeV8QjC9NZt+0Ag7q05u83D2VA51ZelyUBpHAXCWNFxeU8uzSLV1Zvo3XThswYk8h1gztTr56mYMKdwl0kDFVU+nhtzXaeXZrNwaPl3HROV+65uI+mYCKIwl0kzHy0aQ+Pv51BdsFhzukezSNX9CehY0uvy5I6pnAXCRM5e48w7Z0MPsjcTVx0U/5y02Au6R+jVTARSuEuEuIOHCnjzx9u5pXVW2lYvx6TRvXlF8PjaRxV3+vSxEMKd5EQVVJeyexPt/LfyzdzpLSC65O6cPePe9OhhfZZF4W7SMjx+RxvpeQx4/0sdhYeZWSf9ky+tB99zmjhdWkSRBTuIiHk06/28uS7G0ndWUT/ji2ZcV0iw3u287osCUIKd5EQkLaziBmLs1iZvYeOrRrz7PUDuXpQJ61Xl++kcBcJYtv3FfPMkiwWbcijddMoHrysHz87t6tOlsoJKdxFgtDew6X8edkmXluznfr1jNtH9mDiD3vQqokuQpKaUbiLBJFDJeW8tHILL3+cQ2mFj58M6cKdF/YipqVWwMjJUbiLBIGS8kpeXb2NF5Zv5kBxOZcPiOXuH/emR/vmXpcmIUrhLuKhikof877cyXNLs8krKuH8Xu2475I+JHZu7XVpEuIU7iIecM6xNKOApxdnsWn3YQZ2bsXTYwdqWaPUGoW7SB1bu3U/T763kXXbDtC9XTNevPFsLj3zDO0BI7VK4S5SRzYVHOKp97P4ILOADi0a8cQ1Axib1Jmo+roZtdQ+hbtIgBUcLOHZJdnMWbeDZg0bcN8lffjF8HiaNtSvnwTOCX+6zKwxsBJo5D/+TefcVDOLBt4A4oGtwPXOuQP+9zwATAAqgTucc4sDUr1IEDtSWsFfV27hpZVbqPD5+MXwbvx6ZE/aNGvodWkSAWoydCgFLnDOHTazKOBjM3sPuBZY5pybbmaTgcnAJDNLAMYB/YGOwAdm1ts5VxmgPogElUqfI3ntDp5dms2eQ6VcnhjLpEv6Ete2qdelSQQ5Ybg75xxw2P80yv/HAaOBEf722cAKYJK//XXnXCmQY2abgaHAqtosXCQY/Sd7D0+8k0lWwSEGd23DX382mLPj2nhdlkSgGk36mVl9YB3QE3jBOfeZmcU45/IBnHP5ZtbBf3gnYHW1t+f62479nhOBiQBxcXGn3gORILBlz2F+/04mH27cTde2TfmfG89mlFbAiIdqFO7+KZVBZtYamG9mZ37P4cf7aXbH+Z4zgZkASUlJ33pdJBQcKinnzx9u5u+f5NCoQX2mXNaXm8/rRsMGWgEj3jqp0/XOuUIzWwGMAgrMLNY/ao8FdvsPywW6VHtbZyCvNooVCRY+n+PNdbnMWLyRfUfKGDu4M/dd0pf2LRp5XZoIULPVMu2Bcn+wNwEuAp4CFgHjgen+rwv9b1kEvGZmz1J1QrUXsCYAtYt4YsOOQh5emEZKbhGDu7Zh1s1DtF2ABJ2ajNxjgdn+efd6QLJz7m0zWwUkm9kEYDswFsA5l25myUAGUAHcrpUyEg6KisuZsXgjr63ZTvvmjXjuJ4MYPaij5tUlKFnVYhhvJSUlubVr13pdhshx+XyOuV/kMv29jRwoLuPm87rx24t70aKx9lYXb5nZOudc0vFe0yVyIt9j466DPLwgjc+3HuDsuNb8c8JQ+nds5XVZIiekcBc5jpLySp5ftomZK7fQsnEDZoxJ5LrBnXXPUgkZCneRY3y6eS9T5qeydV8xYwd3Zspl/bRlgIQchbuIX2FxGdPeyWTOuly6tm3Ka7cO4zztry4hSuEuEc85x9sp+Tz6VjoHisu57Uc9uOuiXjSOqu91aSKnTOEuEW33oRIemp/GkowCEju3YvYtOmEq4UHhLhHJOcf8L3fy6FsZHC2vZMplfblleDca6MYZEiYU7hJxdhWV8OD8VJZt3M3grm2YcV0iPdo397oskVqlcJeI4ZxjzrpcHn87g/JKHw9fkcDN58VTX8sbJQwp3CUiFBwsYfLcFJZn7WFot2hmjEkkvl0zr8sSCRiFu4Q15xxvpeTz8II0SisqecQ/WtfFSBLuFO4StvYfKePhBWm8k5rPWXGt+cPYgXTX3LpECIW7hKWlGQU8MC+FoqPl3D+qDxPP766VMBJRFO4SVg6VlPPoWxm8uS6XhNiWvHrrMPqe0dLrskTqnMJdwsbqLfu4J3kD+UVH+fXIntxxYS/d7k4ilsJdQl5JeSV/WJLFyx/n0DW6KW/+6jzOjmvjdVkinlK4S0hL21nE3cnryS44zE3nxDHlsn40bagfaxH9FkhIqvQ5/vKfr3jug2zaNG3IP34xhBF9OnhdlkjQULhLyNm+r5i7k9ezdtsBLhtwBtOuHqD91kWOoXCXkOGcI3ntDh57K4N69Uw3qBb5Hgp3CQl7D5cyeW4qH2QWcF6PtjwzdiAdWzfxuiyRoKVwl6C3NKOAyXNTOFRawcNXJPALbR8gckIKdwlah0rKefztDJLXVl2Q9O9xg+gd08LrskRCgsJdgtKanP3cnbyevMKj3D6yB3de2FsXJImcBIW7BJXSikqeXZLNzI+2EBfdlDm3ncvgrtFelyUSchTuEjQy8g5yd/J6Nu46xI3Dqi5IatZIP6Iip0K/OeK5ikoff125hec+yKZ104b8/eYhjOyrC5JETofCXTy1Zc9h7pmzgS+3F3J5Yiy/H32mLkgSqQUnPENlZl3MbLmZZZpZupnd6W+PNrOlZrbJ/7VNtfc8YGabzSzLzC4JZAckNPl8jn98ksNlf/qILXuO8KcbzuKFn56tYBepJTUZuVcA9zjnvjCzFsA6M1sK3Awsc85NN7PJwGRgkpklAOOA/kBH4AMz6+2cqwxMFyTU5B4o5v43U/j0q32M6NOep8YkEtOysddliYSVE4a7cy4fyPc/PmRmmUAnYDQwwn/YbGAFMMnf/rpzrhTIMbPNwFBgVW0XL6HFOccbn+/g9+9k4pzjyWsHMG5IF20fIBIAJzXnbmbxwFnAZ0CMP/hxzuWb2ddnwDoBq6u9LdffJhEsv+gok+amsjJ7D+d2b8uM6xLpEt3U67JEwlaNw93MmgNzgbuccwe/Z7R1vBfccb7fRGAiQFxcXE3LkBDjnGPuFzt59K10Kiodj43uz03Dumr7AJEAq1G4m1kUVcH+L+fcPH9zgZnF+kftscBuf3su0KXa2zsDecd+T+fcTGAmQFJS0rfCX0Lf7oMlPDAvlWUbdzM0PpqnxybStW0zr8sSiQgnDHerGqL/Dch0zj1b7aVFwHhguv/rwmrtr5nZs1SdUO0FrKnNoiW4fT1af+ytdMoqfTxyRQI3a7MvkTpVk5H7cOBnQKqZrfe3TaEq1JPNbAKwHRgL4JxLN7NkIIOqlTa3a6VM5MgvOsqUeaksz9rDkPg2PDUmke7tm3tdlkjEqclqmY85/jw6wIXf8Z5pwLTTqEtCzNc30vj925mU+3xMvTKB8edqtC7iFV2hKqdtZ+FRJs9N4aNNexnWLZoZ12luXcRrCnc5ZT6f419rtjP93Uwc8Pjo/tyolTAiQUHhLqckZ+8RJs1NYU3Ofn7Qsx1PXjtA69ZFgojCXU5Kpc8x6+McnlmSRcMG9XhqzACuT9JVpiLBRuEuNZZdcIj73kxhw45CLuoXw7RrztSeMCJBSuEuJ1RW4eN/VnzFfy/fRIvGUfzphrO4MjFWo3WRIKZwl++1YUchk+amsHHXIa4c2JHfXZlA2+aNvC5LRE5A4S7HVVJeyR+XZvPSR1to36IRL/08iYsTYrwuS0RqSOEu3/LZln1MmpvC1n3F3DC0C5Mv7UerJlFelyUiJ0HhLv/ncGkFT723kVdWb6NLdBNeu3UY5/Vs53VZInIKFO4CwH+y9zBlXip5RUe5ZXg37r2kN00b6sdDJFTptzfCFRWX8/t3MpizLpce7Zvx5m3nMbhrmxO/UUSCmsI9gi1J38WDC9LYf6SM20f24DcX9KJxVH2vyxKRWqBwj0AHjpQxdVE6izbk0S+2JX+/eQhndmrldVkiUosU7hHm/bRdPLQglcLicn57UW/+a2QPourX87osEallCvcIsd8/Wn9rQx79O7bkn7cMI6FjS6/LEpEAUbhHgPfT8nloQRpFR8u5++Le/GqERusi4U7hHsaKist5ZFEaC9dXjdZfmTCMfrEarYtEAoV7mFqetZvJc1PYd7hMc+siEUjhHmYOl1Yw7Z0M/r1mB71jmvO38VoJIxKJFO5hZPWWfdw7ZwM7C4/yyx9157cX9da6dZEIpXAPA6UVlTyzOIuXP84hLropc355Lknx0V6XJSIeUriHuI27DnLX6+vZuOsQNw6L48HL+2lPGBFRuIcqn8/xt49zeHpxFi2bRPH3m4cwsm8Hr8sSkSChcA9BOwuPck/yelZv2c+PE2J48toBujuSiHyDwj3ELNqQx4PzU/H5HDPGJDI2qbPuZSoi36JwDxGHSsqZujCdeV/u5Oy41jz3k7OIa9vU67JEJEgp3EPAF9sPcNfr68k9UMxdF/Xi1yN70kAXJInI91C4B7FKn+OF5Zt5ftkmYls1Zs5t5zK4q5Y4isiJnXD4Z2azzGy3maVVa4s2s6Vmtsn/tU211x4ws81mlmVmlwSq8HCXe6CYcTNX8ezSbK5MjOXdO89XsItIjdXk//b/AEYd0zYZWOac6wUs8z/HzBKAcUB//3teNDNdInmS3k7J49LnPyIz/xDP/WQQz407i5aNo7wuS0RCyAmnZZxzK80s/pjm0cAI/+PZwApgkr/9dedcKZBjZpuBocCq2ik3vBWXVfDoogzeWLuDs+Ja86dxZ9ElWidNReTkneqce4xzLh/AOZdvZl9fPdMJWF3tuFx/27eY2URgIkBcXNwplhE+0nYWccfrX5Kz9wi3j+zBXRf11i6OInLKavuE6vEWXLvjHeicmwnMBEhKSjruMZHAOcesT7by1HsbadMsin9NGMZ5Pdt5XZaIhLhTDfcCM4v1j9pjgd3+9lygS7XjOgN5p1NgONt3uJR752xgedYeLurXgRnXDSS6WUOvyxKRMHCq/+9fBIz3Px4PLKzWPs7MGplZN6AXsOb0SgxPn2zey6jnP+KTr/bx6FX9eennSQp2Eak1Jxy5m9m/qTp52s7McoGpwHQg2cwmANuBsQDOuXQzSwYygArgdudcZYBqD0nllT6eXZrNX/7zFT3aN+eftwzVre9EpNbVZLXMDd/x0oXfcfw0YNrpFBWuduwv5jf//pL1Owq5YWgXHr4iQdvzikhAKFnqyKINeTw4LxUMXvjp2VyeGOt1SSISxhTuAXaktIKpi9J5c10uZ8e15nmtXReROqBwD6C0nUX85t9fsnXfEe64oCd3XNhLG36JSJ1QuAeAz+eY9UkOT72/kbbNGvHaredwbo+2XpclIhFE4V7L9vrXrq/I2sPFCTHMGJNIGy1xFJE6pnCvRSuydnPvnA0cLKng8dH9uemcrrpLkoh4QuFeC0orKnnqvSxmfZJDn5gWvHrrMPqeobXrIuIdhftp2lRwiDteX09m/kFuPi+eyZf2pXGUdjkWEW8p3E+Rc45/fbadx9/OoHmjBsy6OYkL+sZ4XZaICKBwPyX7j5QxaW4KSzMK+GHv9jwzNpEOLRp7XZaIyP9RuJ+k/2Tv4d45GygqLuehy/txy/Bu1Kunk6YiElwU7jVUUl7J9Pc28o9Pt9InpoU2/BKRoKZwr4HM/IPc+fqXZBcc5hfD45k0SidNRSS4Kdy/x9dXms54P4tWTaOYfctQftS7vddliYickML9O+wsPMq9yRtYtWUfFyfE8NSYRN1MQ0RChsL9GM45Fq7P4+GFafh8jhljEhmb1FlXmopISFG4V1NYXMZDC9J4OyWfwV3b8MfrBxHXVtvzikjoUbj7fbSpaonjvsNl3HdJH277UQ/qa4mjiISoiA/3o2WVTH8vk9mrttGzQ3P+Nn4IZ3Zq5XVZIiKnJaLD/cvtB7gneQNb9h7hluHduH9UHy1xFJGwEJHhXlbh488fbuKF5ZuJbdWE1/7fMM7r0c7rskREak3EhXt2wSHuTl5P2s6DXDe4M49cmUDLxlFelyUiUqsiJtwrfY5ZH+fw9JIsWjRqwF9/NphL+p/hdVkiIgEREeG+bd8R7puTwpqt+7k4IYYnrhlA+xaNvC5LRCRgwjrcv95z/Yl3M6lfz/jD2IFce3YnXZAkImEvbMN9V1EJ989NYWX2Hs7v1Y6nxiTSsXUTr8sSEakTYRfuzjkWrN/J1IXplFc63ahaRCJSWIX73sOlPDg/lcXpBSR1bcMzYwcS366Z12WJiNS5gIW7mY0CngfqAy8756YH6rMA3kvN58EFaRwurWDKZX2Z8IPu2j5ARCJWQMLdzOoDLwAXA7nA52a2yDmXUdufVVhcxtRF6Sxcn8eATq149vqB9IppUdsfIyISUgI1ch8KbHbObQEws9eB0UCthntKbiG3zl7L/iNl3H1xb341ogdR9evV5keIiISkQIV7J2BHtee5wLDqB5jZRGAiQFxc3Cl9SFx0U/qc0YJJo/pqsy8RkWoCNcw93mS3+8YT52Y655Kcc0nt25/aretaN23IKxOGKdhFRI4RqHDPBbpUe94ZyAvQZ4mIyDECFe6fA73MrJuZNQTGAYsC9FkiInKMgMy5O+cqzOzXwGKqlkLOcs6lB+KzRETk2wK2zt059y7wbqC+v4iIfDetGxQRCUMKdxGRMKRwFxEJQwp3EZEwZM65Ex8V6CLM9gDbTuNbtAP21lI5oSDS+gvqc6RQn09OV+fcca8CDYpwP11mttY5l+R1HXUl0voL6nOkUJ9rj6ZlRETCkMJdRCQMhUu4z/S6gDoWaf0F9TlSqM+1JCzm3EVE5JvCZeQuIiLVKNxFRMJQSIe7mY0ysywz22xmk72uJxDMbJaZ7TaztGpt0Wa21Mw2+b+28bLG2mZmXcxsuZllmlm6md3pbw/bfptZYzNbY2Yb/H1+1N8etn2Gqvstm9mXZva2/3m493ermaWa2XozW+tvC0ifQzbcq92E+1IgAbjBzBK8rSog/gGMOqZtMrDMOdcLWOZ/Hk4qgHucc/2Ac4Db/X+34dzvUuAC59xAYBAwyszOIbz7DHAnkFntebj3F2Ckc25QtbXtAelzyIY71W7C7ZwrA76+CXdYcc6tBPYf0zwamO1/PBu4uk6LCjDnXL5z7gv/40NU/fJ3Ioz77aoc9j+N8v9xhHGfzawzcDnwcrXmsO3v9whIn0M53I93E+5OHtVS12Kcc/lQFYRAB4/rCRgziwfOAj4jzPvtn6JYD+wGljrnwr3PzwH3A75qbeHcX6j6B1sNzpMAAAGZSURBVHuJma0zs4n+toD0OWA366gDJ7wJt4Q2M2sOzAXucs4dNDveX3n4cM5VAoPMrDUw38zO9LqmQDGzK4Ddzrl1ZjbC63rq0HDnXJ6ZdQCWmtnGQH1QKI/cI/km3AVmFgvg/7rb43pqnZlFURXs/3LOzfM3h32/AZxzhcAKqs61hGufhwNXmdlWqqZULzCzVwnf/gLgnMvzf90NzKdqejkgfQ7lcI/km3AvAsb7H48HFnpYS62zqiH634BM59yz1V4K236bWXv/iB0zawJcBGwkTPvsnHvAOdfZORdP1e/uh865mwjT/gKYWTMza/H1Y+DHQBoB6nNIX6FqZpdRNW/39U24p3lcUq0zs38DI6jaFrQAmAosAJKBOGA7MNY5d+xJ15BlZj8APgJS+f/zsVOomncPy36bWSJVJ9PqUzXoSnbOPWZmbQnTPn/NPy1zr3PuinDur5l1p2q0DlVT4q8556YFqs8hHe4iInJ8oTwtIyIi30HhLiIShhTuIiJhSOEuIhKGFO4iImFI4S4iEoYU7iIiYeh/Aeo3cobd4d5UAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fitlength = 50.\n",
    "calculation_points = np.arange(fitlength + 1.)\n",
    "combined_precision= np.sqrt(2) * af647_precision\n",
    "transverse_rpd = plotting.estimate_rpd_churchman_2d(\n",
    "    input_distances=trans_distances[trans_distances < (fitlength + 3 * combined_precision)],\n",
    "    calculation_points=calculation_points,\n",
    "    combined_precision=combined_precision\n",
    ")\n",
    "plt.plot(calculation_points, transverse_rpd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normalise for increasing search circle with increasing distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x258bdea9e48>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD6CAYAAACvZ4z8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU9b3/8ddnJntCNpJANgg7si9Bwa2KVRGXuuFWvbbVeu/t3moXf7ftbX+tbW8Xa/dbf2qxLljqVqso7mIVwbBp2NdAAllYskL27++PDC2lCQkhkzMzeT8fjzySOTmTeX9F3vlyzpnvMeccIiISuXxeBxARkeBS0YuIRDgVvYhIhFPRi4hEOBW9iEiEU9GLiES4bovezB4ys0ozKz5m2wIzW29m7WZWeILnpprZk2a2ycw2mtmcvgouIiI9E9WDfRYCvwb+eMy2YuBq4PfdPPcXwEvOuWvNLAZI6EmojIwMV1BQ0JNdRUQEWLVq1X7nXGZn3+u26J1zy8ys4LhtGwHMrMvnmVkycC7wicBzmoHmngQuKCigqKioJ7uKiAhgZiVdfS+Yx+hHAlXAH8xsjZk9YGaJXe1sZneYWZGZFVVVVQUxlojIwBLMoo8CZgC/c85NBxqAb3S1s3PufudcoXOuMDOz0399iIhILwSz6EuBUufcisDjJ+kofhER6UdBK3rnXDmwx8zGBTZdAGwI1uuJiEjnenJ55SJgOTDOzErN7DYzu8rMSoE5wAtmtjSwb46ZLTnm6Z8HHjOzD4BpwA/6fggiInIiPbnq5sYuvvVMJ/vuBeYf83gt0OV19iIiEnx6Z6yISIRT0YeJ6sPNPPJeCe9u3097u24WIyI915N3xoqHquqaeOBvO3h0eQkNzW0ADB+cwHWF+Vw7M48hyXEeJxSRUKeiD1Fl1Ue4/63tPPH+Hlra2rl0Sg63nz2C7VX1/On9Pfxk6WbufWUL54/L5PpZwzh/XCZRfv0DTUT+lYo+xNQ3tfL95zfw1OpSnIOrZ+Tyn+eNZkRGx5uKp+ancvWMPHbub2Bx0R6eXFXKqxuLOGv0YB76xCxio/wej0BEQo2F4s3BCwsL3UBc68Y5x+cWreGl4nI+fsYw7jh3JHlpJ14HrqWtnSdW7uZbf1nPldNyuPe6afh8Xa9BJCKRycxWOec6vcpRM/oQ8uh7JbzwwT6+Nm8cnzlvdI+eE+33ccucAmobW/nJ0s1kp8bz9Xnjg5xURMKJij5EFJfV8L3nN3L+uEz+49xRJ/38z5w3ir3VR/jdm9vJSYnjljkFfR9SRMKSij4E1Da28JnHVjM4KYaf9fLQi5nx3SsmUlHbyLefW09WchwXTxwahLQiEm50mYbHnHN8/ckP2Ft9hF/fNJ30xJhe/6wov49f3jidKXmpfGHRGlaVHOrDpCISrlT0Hnv43V28WFzO1+aNY+bw9FP+eQkxUTx4ayFDU+K4/eH32VFV3wcpRSScqeg9tG5PNfcs2chHT8vi0+eM7LOfm5EUy8OfPB0z4/Y/FtHU2tZnP1tEwo+K3iM1h1v47OOryRoUx08XTD3hbRl7oyAjkZ9fP40dVQ3c/9aOPv3ZIhJeVPQe+enLmymvaeRXN00nNaH3x+VP5CNjM7l0Sja/fmMbuw8cDspriEjoU9F7YG/1Ef70/h4WFOYzY1haUF/r25dNINrv49vPFROKb44TkeBT0Xvgd29up905Pnv+yV8vf7KGJMfx5QvH8ubmKpauLw/664lI6FHR97NjZ/PdLW/QV26dM5wJ2cl8968baGhq7ZfXFJHQoaLvZ0dn8585L/iz+aOi/D6+f9Uk9tU0ct+rW/rtdUUkNKjo+9G+mqOz+Tzy0/tnNn/UjGFp3Hh6Pg+9s4uN+2r79bVFxFsq+n70j9l8zxYs62tfnzeelPhovvlsse5SJTKAqOj7yb6aIzyx0pvZ/FGpCTHcfcl4VpUc4slVpZ5kEJH+p6LvJ17P5o+6ZkYeswrS+OGLG6k53OJpFhHpHyr6fhAKs/mjfD7ju1dMovpIC799c5unWUSkf6jo+0GozOaPmpCTzNXT8/jDu7soqz7idRwRCTIVfZCF0mz+WF+5aCwAP3t5s8dJRCTYVPRB9vu3doTUbP6o3NR4PnlWAc+sKWPDXl1uKRLJVPRBVNfYwp+L9nDF1JyQms0f9ZmPjCY5LpofvbTJ6ygiEkQq+iB6Zk0ZDc1t/NuZBV5H6VRKQjSfnzuaZVuq+NvW/V7HEZEgUdEHiXOOR5aXMDk3hal5KV7H6dItc4aTmxrPD1/cqDdRiUQoFX2QrNh5kK2V9dwyZ3if31SkL8VG+fnqxeNYv7eW59bt9TqOiASBij5IHlleQkp8NJdPyfE6SreumJrDxJxkfrJ0M40tuu2gSKRR0QdBRW0jS9eXs2BmHvExfq/jdMvnM+6+5DTKqo/w6HslXscRkT6mog+CRSt309ruuHn2cK+j9NjZYzI4Z0wGv3p9m5ZGEIkwKvo+1tLWzqKVuzl3bCYFGYlexzkp37hkPLWNLfzy9a1eRxGRPtRt0ZvZQ2ZWaWbFx2xbYGbrzazdzAq7eb7fzNaY2fN9ETjUvbKhgoraJv4tjGbzR03MSeG6mfk8/O4udlTVex1HRPpIT2b0C4F5x20rBq4GlvXg+V8ENp5crPD1yPISclPjOX98ltdReuWui8cRF+3nnhcGzB+ZSMTrtuidc8uAg8dt2+ic63aRFDPLAy4FHuh1wjCyrbKO5TsO8PHZw/D7QveSyhPJHBTL5+aO5rVNlSzbUuV1HBHpA8E+Rn8f8DWgvbsdzewOMysys6KqqvAsmEeWlxDj93FdYb7XUU7JJ88qYFh6At97fgOtbd3+0YlIiAta0ZvZZUClc25VT/Z3zt3vnCt0zhVmZmYGK1bQNDS18tTqMuZPHkpGUqzXcU5JbJSf/zP/NLZW1vP4yt1exxGRUxTMGf1ZwBVmtgt4AphrZo8G8fU89cyaMuqbWrllToHXUfrExROHMGfkYO59ZYsutxQJc0Ereufc3c65POdcAXAD8Lpz7uZgvZ6XnHM8+l4JE7KTmTEs1es4fcLM+PblE6g90sJ9r23xOo6InIKeXF65CFgOjDOzUjO7zcyuMrNSYA7wgpktDeybY2ZLghs59KzeXc2m8jpunh3a69qcrNOyk7l+1jAeWV7CtkpdbikSrnpy1c2Nzrls51x0YIb+oHPumcDXsc65Ic65iwP77nXOze/kZ7zpnLssGAMIBYtW7iYxxs8V00J/XZuTdedFY4mP9nPPCxu8jiIivaR3xp6imiMtPP/BXq6YlktSbJTXcfpcRlIsX7hgDG9sruKNTZVexxGRXlDRn6Jn15TR2NLOx88Y5nWUoLn1zAJGZyVx99MfUn242es4InKSVPSnwDnHopW7mZybwqTc0L25yKmKifJx3/XT2F/fxLf+st7rOCJyklT0p+DoSdgbT4/c2fxRk3JT+PKFY/nrur38ZW2Z13FE5CSo6E9BJJ+E7cx/fGQUM4en8c1niymrPuJ1HBHpIRV9L0X6SdjO+H3Gz6+bRnu7467F63SPWZEwoaLvpYFwErYzwwYn8N+XT2T5jgM89M5Or+OISA+o6HvBOcfjKyL/JGxXFhTmcdGEIfz4pc1sKq/1Oo6IdENF3wurd1ezuaKOmwbYbP4oM+OHV08mOT6KLz2xlqZW3VBcJJSp6Hvh6EnYy6cOjJOwnRmcFMuPr53CpvI67nlhI87peL1IqFLRn6SjJ2E/Nn3gnITtytzxQ7j97BH8cXkJP3pxk8peJEQN7KbqhaMnYW8aANfO98R/XXoaja1t/H7ZDqDjBuORtLCbSCRQ0Z+EgX4StjNmxvc+NglAZS8SolT0J+H9XYfYXFHHD6+e7HWUkKKyFwltKvqTsPDdnaTER3PltFyvo4Qclb1I6FLR91BZ9RGWrq/g9nNGEB/j9zpOSDq+7Jta2/n6vPH67yXiMRV9Dz2yvATnHLfMHu51lJB2tOyjfD4WvruLF4v3ceeF47hmZh5+n2b3Il7Q5ZU9cKS5jUUrd3PxxKHkpSV4HSfkmRnfuWIif/6POeSkxvO1pz7gkl8s441NlboEU8QDmtH3wLNry6g50sInzizwOkpYmVWQztP/eSYvFpfz45c28cmF7zNn5GC+Om8cU/NSezXDb2hqZXtVfcdHZQMVtY00NLdS19hKQ1Mr9U2tNDS14fcZU/JSmDEsjRnD05iQnUxMlOY1MjBZKM6wCgsLXVFRkdcxgI5LKi++bxlRPh8vfOFsnVzspebWdh5fUcIvX9/GwYZmYqN8jMpMYuyQJMYMGcSYrCRGZibS2NLOocPNHGz4x8eBhmZ2HzjM9qp69tU0/v1n+n1GZlIsSXFRJMX+4yMxNoojLa2s3V3N3sD+sVE+JuemcPqIdG47ewSDk2K9+k8hEhRmtso5V9jp91T0J/butv3c9MAKfnztFK4rzPc6TtirbWzhpeJytpTXsbWynq0VdX8v486YQVpCDPlp8YzKTGJUVhKjMpMYnZXIsPTEbmfp5TWNrN59iNUlh1i9+xAflNYwKC6Kb102gaum5+oXt0QMFf0puP3hItbsPsQ735hLXLSuHgmGusYWtlXWs+tAA/HRftITY0lPjCE9MYaU+Og+PYm7paKObzz1Aat3V3POmAx+cNVk8tN13kXCn4q+l3YfOMxHfvoGnzt/NHdeNM7rONJH2tsdj60o4X9e2kxreztfuXAsnzprBFF+HcOX8HWiotf/2Sfw8PJd+M24WZdURhSfz7hlTgGvfOVczh6dwQ+WbOLK375DyYEGr6OJBIWKvgsNTa0sfn8P8ydnMyQ5zus4EgTZKfH8v38r5Dc3zaD00BFufWglB+qbvI4l0udU9F14anUpdU2tfOKsAq+jSBCZGZdOyebBW2exr6aR2/9YRGOLbqQikUVF34n2dsfCd3YxNT+VGcPSvI4j/WDm8DR+ccM01u6p5otPrKFNNz6XCKKi78QbmyvZsb+BT+oNUgPKvEnZfOvSCSxdX8E9L2z0Oo5In9E7Y4/T3u6495Ut5KfHM39yttdxpJ996uwRlB46wkPv7CQ3LZ7bzh7hdSSRU6aiP85L68tZv7eWny2YqrfMD1D/delp7K0+wvdf2EBuahzzJukXvoQ3Ndkx2gKz+dFZSVw5XWvOD1R+n3HfDdOYlp/KF59Yy6qSQ15HEjklKvpjPLumjG2V9XzlwrFaUneAi4v28+Cts8hKjuXOxWtpatWVOBK+VPQBza3t3PfaFibmJDNv4lCv40gISE+M4Z4rJ7PrwGEeeHun13FEeq3bojezh8ys0syKj9m2wMzWm1m7mXW+toJZvpm9YWYbA/t+sS+D97XFRXvYc/AId100Dp9m8xJw7thMLpk0lF+9vpXSQ4e9jiPSKz2Z0S8E5h23rRi4Glh2gue1Anc6504DZgOfNbMJvQkZbI0tbfzq9a3MHJ7GeeMyvY4jIeabl03AML7/vC65lPDUbdE755YBB4/bttE5t7mb5+1zzq0OfF0HbARC8gzno++VUFHbxF0XjdOytfIvclPj+dzc0by0vpy3tlR5HUfkpPXLMXozKwCmAytOsM8dZlZkZkVVVf33l6m+qZXfvrmds0dnMGfU4H57XQkvt58zghEZiXznufU6MSthJ+hFb2ZJwFPAl5xztV3t55y73zlX6JwrzMzsv8Mnf/jbTg42NHPXxVqGWLoWG+XnO1dMZOf+Bp2YlbAT1KI3s2g6Sv4x59zTwXyt3qg53ML9b+/go6cNYVp+qtdxJMR9ZGwm8yZ2nJgtqz7idRyRHgta0VvHwe4HgY3OuXuD9Tqn4uevbqG+qZU7LxrrdRQJE9+6vON6gu/9dYPHSUR6rieXVy4ClgPjzKzUzG4zs6vMrBSYA7xgZksD++aY2ZLAU88CbgHmmtnawMf8II3jpC3bUsXCd3dxy+zhnJad7HUcCRO5qfF8fu4YnZiVsDIgbyV4sKGZefctIzk+muc/f7buBSsnpam1jXn3vY3fZ7z8pXP1vgsJCbqV4DGcc9z99AccOtzML26YppKXkxYb5efLF45lW2U9L2+o8DqOSLcGXNEvLtrD0vUV3HXROCbmpHgdR8LU/ElDGZaewO/e2k4o/qtY5FgDquh37m/gu3/dwJyRg/n0OSO9jiNhLMrv445zR7JuTzXLdxzwOo7ICQ2Yom9pa+dLf1pLlM/42XVTdVxVTtm1M/PISIrlf9/a4XUUkRMaMEX/q9e2sm5PNT+4ejI5qfFex5EIEBft51NnF7BsSxXFZTVexxHp0oAo+lUlB/n1G9u4ekYul03J8TqORJCbZw9nUGwU//vWdq+jiHQp4ot+5c6D3PZwEblp8Xz3iolex5EIkxwXzU2zh7Hkw33s2t/gdRyRTkV00T+zppSbH1hBekIMj3zqDAbFRXsdSSLQbWeNIMrn4/63daxeQlNEFr1zjp+/soUv/2kdM4an8vRnzqQgI9HrWBKhspLjuGZmHk+uKqWyrtHrOCL/IuKKvqm1jS//aS2/eG0r187M44+fOoPUhBivY0mE+/dzR9La1s4f3tnldRSRfxFRRX+woZmbH1jBs2v38tWLx/GTa6cQExVRQ5QQVZCRyCWTs3l0eQm1jS1exxH5JxHTgocamrnqt++wrrSGX904nc+eP1p3i5J+9Z8fGUVdUyuPvbfb6ygi/yRiij41IZp5k4ay6NOzuXyqLqGU/jcpN4VzxmTw4N926i5UElIipujNjLsvOY2Zw9O8jiID2O3njGR/fRNL12uxMwkdEVP0IqHgnNEZ5KfH89h7JV5HEfk7Fb1IH/L5jBtPH8aKnQfZVlnvdRwRQEUv0ucWzMwn2m8sWqmTshIaVPQifSxzUCwXTRzKk6tKaWzRSVnxnopeJAg+fvowao60sOTDfV5HEVHRiwTDnFGDGZmRyGMrdPhGvKeiFwkCs46TsqtKDrGpvNbrODLAqehFguSamXnERPl4XLN68ZiKXiRI0hNjmD9pKM+sLuNwc6vXcWQAU9GLBNFNZwynrqmV59fppKx4R0UvEkSzCtIYk5XEYyv0TlnxjopeJIjMjJvOGMa60hrdQFw8o6IXCbKrp+cRF+3TpZbiGRW9SJClJERz2ZQcnltbRn2TTspK/1PRi/SDm84YRkNzG39dt9frKDIAqehF+sH0/FRGZyXx56I9XkeRAUhFL9IPzIwFM/NYvbtayxdLv1PRi/STq2bk4vcZT64q9TqKDDAqepF+kjUojvPHZfLU6lJa29q9jiMDiIpepB9dOzOfqromlm2t8jqKDCDdFr2ZPWRmlWZWfMy2BWa23szazazwBM+dZ2abzWybmX2jr0KLhKu547MYnBjDn4t0+Eb6T09m9AuBecdtKwauBpZ19SQz8wO/AS4BJgA3mtmE3sUUiQwxUT6unJ7LqxsrONjQ7HUcGSC6LXrn3DLg4HHbNjrnNnfz1NOBbc65Hc65ZuAJ4GO9TioSIRYU5tHS5nh2TZnXUWSACOYx+lzg2IuGSwPbOmVmd5hZkZkVVVXp+KVErvFDk5mSl8Lioj0457yOIwNAMIveOtnW5f/Vzrn7nXOFzrnCzMzMIMYS8d6CmXlsKq9j/V7dfUqCL5hFXwrkH/M4D9D7v0WAK6bmEhPl0ztlpV8Es+jfB8aY2QgziwFuAJ4L4uuJhI2UhGgunjiUZ9fupbGlzes4EuF6cnnlImA5MM7MSs3sNjO7ysxKgTnAC2a2NLBvjpktAXDOtQKfA5YCG4HFzrn1wRqISLhZMDOPmiMtvLqxwusoEuGiutvBOXdjF996ppN99wLzj3m8BFjS63QiEeys0RnkpMTx56JSLpuS43UciWB6Z6yIR/w+45qZeby9tYrymkav40gEU9GLeOjamXm0O3hqtd4pK8Gjohfx0PDBicwemc7ioj20t+uaegkOFb2Ix66flU/JgcO8t/OA11EkQqnoRTx2yaRsBsVFsfh9XVMvwaGiF/FYXLSfK6fl8mJxOTWHW7yOIxFIRS8SAq6flU9Tazt/WaeFzqTvqehFQsCk3BQm5iTzJx2+kSBQ0YuEiBtm5bN+by3FZTVeR5EIo6IXCRFXTMslNsqnWb30ORW9SIhIiY/mkklDeXZtmRY6kz6lohcJIdfPGkZdYysvFu/zOopEEBW9SAiZPTKd4YMTdPhG+pSKXiSEmBnXFebz3o6D7Nrf4HUciRAqepEQc+3MPHwGi3X3KekjKnqREDMkOY6547N4clUprW3tXseRCKCiFwlB1xXmU1nXxJubq7yOIhFARS8Sgs4fn0XmoFgWrdztdRSJACp6kRAU7fdx4+nDeH1zJSUHdFJWTo2KXiRE3XzGMPxm/HF5iddRJMyp6EVCVFZyHJdOyWbx+3uob2r1Oo6EMRW9SAj7xJkF1DW18rTuKSunQEUvEsKmD0tjan4qC9/dpXvKSq+p6EVC3CfPLGBHVQNvb9vvdRQJUyp6kRA3f3I2mYNiWfjOTq+jSJhS0YuEuJgoHzefMZw3NlexU+vfSC+o6EXCwE1nDCPabzz87i6vo0gYUtGLhIHMQbFcPiWHPxftoa6xxes4EmZU9CJh4tYzC2hobuPJVbrUUk6Oil4kTEzNT2XGsFQe1qWWcpJU9CJh5BNnjWDXgcO8tUWrWkrPqehFwsglk4YyJDmWh3SppZwEFb1IGIn2+/jEmSN4e+t+VpUc8jqOhIlui97MHjKzSjMrPmZbupm9YmZbA5/Tunjul81svZkVm9kiM4vry/AiA9GtZw4nIymGnyzdhHM6Vi/d68mMfiEw77ht3wBec86NAV4LPP4nZpYLfAEodM5NAvzADaeUVkRIiInis+eP5r0dB3ln2wGv40gY6LbonXPLgIPHbf4Y8HDg64eBK7t4ehQQb2ZRQAKwt5c5ReQYN50xjJyUOH7y8mbN6qVbvT1GP8Q5tw8g8Dnr+B2cc2XAT4HdwD6gxjn3clc/0MzuMLMiMyuqqtIVBSInEhvl54sfHcO6PdW8urHS6zgS4oJ2MjZw3P5jwAggB0g0s5u72t85d79zrtA5V5iZmRmsWCIR45oZeYzISORnL2/WdfVyQr0t+gozywYIfO5sSvFRYKdzrso51wI8DZzZy9cTkeNE+X186aNj2FRex18/0FFR6Vpvi/454NbA17cCf+lkn93AbDNLMDMDLgA29vL1RKQTl0/JYfzQQdz36lZa29q9jiMhqieXVy4ClgPjzKzUzG4DfgRcaGZbgQsDjzGzHDNbAuCcWwE8CawGPgy81v1BGYXIAOXzGXdeNI6d+xt4SrcblC5YKJ6xLywsdEVFRV7HEAkLzjmu/O27VNU28sZXzyM2yu91JPGAma1yzhV29j29M1YkzJkZX7t4HHtrGnl8xW6v40gIUtGLRICzRmcwZ+RgfvPGNmqOaL16+WcqepEIcff88Rw63MJ3nlvvdRQJMSp6kQgxJS+Vz50/mmfWlPG8LreUY6joRSLI5+aOZmp+Kv/1TDHlNY1ex5EQoaIXiSDRfh/3XT+N5tZ27vrzOr1jVgAVvUjEGZGRyLcum8Dftu1n4bu7vI4jIUBFLxKBbjw9nwvGZ/GjlzaxpaLO6zjiMRW9SAQyM350zRQGxUbxpSfW0tyq5REGMhW9SITKHBTL/1wzhQ37arn3lS1exxEPqehFIthHJwzhxtPz+f2y7bz44T6v44hHVPQiEe6bl05gen4qn318NYvf3+N1HPGAil4kwiXGRvHo7Wdw1ugMvvbUB9y/bLvXkaSfqehFBoCEmCgevHUWl07J5gdLNvHjlzbpXrMDSJTXAUSkf8RE+fjlDdNJiY/mt29u59DhFr5/5ST8PvM6mgSZil5kAPH7jHuunERaQjS/eWM7tUdauPf6qSGxhn17u2N/fRN7axqpb2zF5wO/GVF+w2dGlM9HfIyf/PT4kMgbTlT0IgOMmfHVi8eTGh/DPUs28kFZNf9+7iiunZlHXHRwC9Q5R3ltIx+W1lC8t5bdBxrYW9PI3uojVNQ20tLW/eEkn0FuWjwjMpIYMTiBERmJjMpKYlp+KoPiooOaP1zpDlMiA9ibmyv5+atbWbenmsxBsXz6nBHcdMZwkmL7Zg5YXtPIB6XVfFhWw4dlNRSX1bC/vhnoKOyc1HhyUuLJSY0jOzU+8DiOQXHRtLU72p2jtd3R3u5oa3fUNbWwc/9hdu5vYNf+Bnbub6C+qRXo+NfKlLwUzhw1mDNHZTBzeFrQf3GFkhPdYUpFLzLAOedYvv0Av3lzG+9sO0BKfDS3nlnA1dNzyU2LJ9rf/TUbzjkONDTzYVkNH+yp4cOyaj4oraGyrgnoKOExWUlMyk1hcm4Kk3KTOS07mYSYU/uF4pyjqr6JzeV1vLfjAMu3H2BdaQ1t7Y4Yv48Zw1O5YPwQ5k0aSn56wim9VqhT0YtIj6zdU81v39jGyxsqgI5Z99DkOPLSE8hLiycvLYGkWD9VdU2U1zZRUdtIZW0jFbVNHGlpA8AMRmUmMSU3hSl5KUzOS2VCdjLxMf0zu65vauX9nQd5d/t+3t66n03lHWv9TMxJZt7EoVwyeSijswb1S5b+pKIXkZOyrbKO1burKT10hNKDhzs+HzrMvtpGnIO4aB9Dk+PISo5jSHIcQwbFMjQljkm5KUzMSQ6pY+W7Dxxm6fpyXizex+rd1QCMykzkkknZzJ+czWnZgzAL/yuPVPQi0ieaW9tpbG1jUGxUWJZjRW0jL68v58Xict7bcYB217Gs8/zJQ5k/OZsJ2clhOS5Q0YuI/Iv99U0sXV/Okg/3sXx7R+kXDE5g3qRsLjgti+n5qUT14PxEqFDRi4icwIH6Jpaur2DJh/t4b8cBWtsdqQnRnDc2k/PHZ3He2CxSEkLncFRnVPQiIj1Uc6SFt7dW8fqmSt7cXMXBhmb8PmN6fiqzRqRTODyNmcPTSE2I8TrqP1HRi4j0Qlu7Y11pNa9vrOTtbftZX1ZDa+A+vKOzkigcnsaM4WlMyE5mdFaSp9ftq+hFRPrAkeY21pVWs6rkEEW7DrKq5BC1jR1v2PIZDB+cyNghSYwbMoixQwcxOiuJgsGJ/fIL4ERFryUQRER6KD7Gz+yRg5k9chzHw4EAAAVaSURBVDDQsT7Pjv31bC6vZ3NFHVsr6thcUccrGyoITPzxGeSnJzAqM4lRmYmMykxiQk4yE7KT++1kr4peRKSXfD5jdNYgRmcN4lKy/769saWN7VX1bK9qYHtl/d+/fmfbfpoC9+9NjPEzY3gaswrSmVWQzrT81KC9qUxFLyLSx+Ki/UzMSWFiTso/bW9vd5RVH2HNnmqKdh1k5c6D/PzVLTgH0X5jen4aT9wxG18fLx2tohcR6Sc+n5GfnkB+egJXTM0BoOZwC6t2H2TlzkNUH27u85IHFb2IiKdSEqKZO34Ic8cPCdprhM/bvkREpFdU9CIiEa7bojezh8ys0syKj9mWbmavmNnWwOe0Lp6bamZPmtkmM9toZnP6MryIiHSvJzP6hcC847Z9A3jNOTcGeC3wuDO/AF5yzo0HpgIbe5lTRER6qduid84tAw4et/ljwMOBrx8Grjz+eWaWDJwLPBj4Oc3OuepTSisiIiett8fohzjn9gEEPmd1ss9IoAr4g5mtMbMHzCyxqx9oZneYWZGZFVVVVfUyloiIHC+YJ2OjgBnA75xz04EGuj7Eg3PufudcoXOuMDMzM4ixREQGlt4WfYWZZQMEPld2sk8pUOqcWxF4/CQdxS8iIv2ot2+Yeg64FfhR4PNfjt/BOVduZnvMbJxzbjNwAbChJz981apV+82spJvdMoD9Jxc7ImjcA4vGPbCcyriHd/WNbpcpNrNFwHmBABXAfwPPAouBYcBuYIFz7qCZ5QAPOOfmB547DXgAiAF2AJ90zh3q5SCOz1XU1ZKckUzjHlg07oElWOPudkbvnLuxi29d0Mm+e4H5xzxeCwy4PywRkVCid8aKiES4cC76+70O4BGNe2DRuAeWoIw7JG8lKCIifSecZ/QiItIDKnoRkQgXdkVvZvPMbLOZbTOzLt9pG+5OZdXQcGZm+Wb2RmC10/Vm9sXA9ogeu5nFmdlKM1sXGPd3A9sjetxHmZk/sFTK84HHA2Xcu8zsQzNba2ZFgW19PvawKnoz8wO/AS4BJgA3mtkEb1MFzUJ6v2poOGsF7nTOnQbMBj4b+DOO9LE3AXOdc1OBacA8M5tN5I/7qC/yz6vbDpRxA5zvnJt2zPXzfT72sCp64HRgm3Nuh3OuGXiCjpU0I05vVw0Nd865fc651YGv6+j4y59LhI/ddagPPIwOfDgifNwAZpYHXErHmyuPivhxn0Cfjz3cij4X2HPM49LAtoGiJ6uGRgwzKwCmAysYAGMPHL5YS8faUa8E1omK+HED9wFfA9qP2TYQxg0dv8xfNrNVZnZHYFufjz3cbg7e2e3RdX1oBDKzJOAp4EvOuVqzzv7oI4tzrg2YZmapwDNmNsnrTMFmZpcBlc65VWZ2ntd5PHCWc26vmWUBr5jZpmC8SLjN6EuB/GMe5wF7PcrihZ6sGhr2zCyajpJ/zDn3dGDzgBg7QOAGPW/ScY4m0sd9FnCFme2i41DsXDN7lMgfN/D3ZWNwzlUCz9BxeLrPxx5uRf8+MMbMRphZDHADHStpDhRHVw2FLlYNDXfWMXV/ENjonLv3mG9F9NjNLDMwk8fM4oGPApuI8HE75+52zuU55wro+Pv8unPuZiJ83ABmlmhmg45+DVwEFBOEsYfdO2PNbD4dx/T8wEPOuXs8jhQUJ7NqqFcZg8HMzgbeBj7kH8ds/w8dx+kjduxmNoWOE29+OiZgi51z/9fMBhPB4z5W4NDNXc65ywbCuM1sJB2zeOg4jP64c+6eYIw97IpeREROTrgduhERkZOkohcRiXAqehGRCKeiFxGJcCp6EZEIp6IXEYlwKnoRkQj3/wGhwh20C1TnjQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "normalised_transverse_rpd = transverse_rpd[calculation_points > 0.] / calculation_points[calculation_points > 0.]\n",
    "norm_rpd_calculation_points = calculation_points[calculation_points > 0.]\n",
    "plt.plot(norm_rpd_calculation_points, normalised_transverse_rpd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1 nm-bin  histogram result for comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(hist_1nm_bins[1][0:-1] + 0.5, hist_1nm_bins[0]/(hist_1nm_bins[1][0:-1] + 0.5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optional save/load to save time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save('normalised_transverse_rpd_smoothed_Churchman-4p4', normalised_transverse_rpd)\n",
    "# normalised_transverse_rpd = np.load('normalised_transverse_rpd_smoothed_Churchman-4p4.npy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up an RPD model and fit\n",
    "I've tried a few smoothing kernel widths here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trans_model_with_info = zdisk_modelling.set_up_model_2d_onepeak_plus_replocs_flat_bg_with_fit_settings()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(params_optimised,\n",
    " params_covar,\n",
    " params_1sd_error) = zdisk_modelling.fitmodel_to_hist(\n",
    "    norm_rpd_calculation_points[0:31],\n",
    "    normalised_transverse_rpd[0:31],\n",
    "    trans_model_with_info.model_rpd,\n",
    "    trans_model_with_info.initial_params,\n",
    "    trans_model_with_info.param_bounds\n",
    "    )\n",
    "print('')\n",
    "print('Initial parameter guesses:')\n",
    "print(trans_model_with_info.initial_params)\n",
    "print('')\n",
    "print('Parameter bounds:')\n",
    "print(trans_model_with_info.param_bounds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(norm_rpd_calculation_points,\n",
    "    normalised_transverse_rpd)\n",
    "zdisk_plots.plot_fitted_model(\n",
    "    norm_rpd_calculation_points,\n",
    "    fitlength,\n",
    "    params_optimised,\n",
    "    params_covar,\n",
    "    trans_model_with_info,\n",
    "    plot_95ci=False\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitlength = 50.\n",
    "calculation_points = np.arange(fitlength + 1.)\n",
    "combined_precision=8.\n",
    "transverse_rpd_s8 = plotting.estimate_rpd_churchman_2d(\n",
    "    input_distances=trans_distances[trans_distances < (fitlength + 3 * combined_precision)],\n",
    "    calculation_points=calculation_points,\n",
    "    combined_precision=combined_precision\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXhV1d33//c3J/NASEjCFCBMMggyBRBERdQWUUvVG8WqRUVxQEXbatunta3et/3Zu1q16tOKioADSrHgWIU6FAemMGmYERnCEMIQ5szr+SNHfxFDEkJO9sk5n9d15dpnr+x9znfpxYfNOmuvbc45REQkdEV4XYCIiASWgl5EJMQp6EVEQpyCXkQkxCnoRURCXKTXBVQnLS3NZWVleV2GiEiTsXTp0j3OufTqfheUQZ+VlUVOTo7XZYiINBlmtuVEv9PQjYhIiFPQi4iEOAW9iEiIU9CLiIQ4Bb2ISIhT0IuIhLhag97MppjZbjPLrdI2xsxWmVmFmWXXcr7PzJab2dsNUbCIiJycusyjnwo8BUyv0pYLXA48U4fzJwFrgGYnW5xXnHMcKy3ncHEZR4rLOVJcVvlTUrlfXFZBSVkFpeWV2xL/FsAXYfgijEj/1hdhxEb5aB4XRXJcFMnx/m1cFIkxkZiZx70VkVBXa9A75+abWdZxbWuAWkPKzDKBi4GHgJ/Vt8i6yj9YxLGScorKyiu3pRUUlZZTVFrOkZLKwD78zU9R5fZQURmHikort8X+bVEZ5RWBX6c/OjKCdilxdGiRQPvUeNqnxtOhRTwdWiTQMS0BX4T+EhCRUxfoO2MfB+4Dkmo70MwmABMA2rdvX68PO+d/P6LYf2VdkyifkRgTSWJsJIkxUSTFRtKmeSxJsUkkxUaS5G9PjPGREBNJQkwkif5tQrSP2CgfUb4IoiMrf6J8RrSvchSsvMJR7hzlFY6yCkd5eeW/Dg4cK6XwaCkHjpVy4FgJhUdL2XO4mK37jrJ13zEWbdrLkZLyb2tMiPbROzOZPu2a0zezOX3bN6dVs1j9C0BETlrAgt7MLgF2O+eWmtnw2o53zk0GJgNkZ2fX63L6f37cC1+EERdVGcYxURHfvo6P9n0b7jGRvvq8fZ1E+ux7/1FTgDbN42o8zznH3iMlbN13lE0FR/gir5CV2wqZ8unXlJZX/ufISIrhrC5pnNc9g3O7ppMcHxWYTohISAnkFf1ZwI/MbBQQCzQzs5ecc9cG6gPHZLcL1FsHnJmRlhhDWmIM/dun8F8DMgEoLitn9Y6DrNxWyLKthXy8bjezl2/HF2EM6JDCiO4ZjOieQdeMRF3ti0i1rC7PjPWP0b/tnOt1XPvHwC+cczWuQOa/ov+Fc+6SuhSVnZ3ttKhZ9corHCu2FfLh2nw+XFvAmp0HAeicnsCY7HZc3q8tGc1iPa5SRBqbmS11zlU7C7LWoDezGcBwIA3IB34P7AOeBNKBQmCFc+6HZtYGeM45N+q49xiOgj4gdh44xgdrdjNn+XZytuzHF2Gce1o6YwZkcn6PlkRH6lYJkXBwSkHvBQV9/WwqOMyspXm8viyP/IPFpCZEc0X/ttxwVsdavyMQkaZNQR9myisc8zcU8I+cbby/Kh8DftyvLbee24kuGbVOgBKRJqimoA/KB4/IqfFFGOd1y+C8bhls23eU5z/9mleXbGXW0jwu6NGS24Z3ZkCHFK/LFJFGoiv6MLHvSAnTPt/MtAWbKTxayqCOqfxyZDcGdEj1ujQRaQAaupFvHS0p47Ul2/jbx1+x+1AxF/VqxS9HdicrLcHr0kTkFCjo5XuOlpTx7PyveWb+V5SUVXDtmR246/yupCZEe12aiNSDgl5OaPehIh7/9wZeW7KN+Cgft53XmfHDOgb07mERaXg1Bb0mWYe5jKRY/nhZb96/+2wGd0rlf99bx0VPfMKiTXu9Lk1EGoiCXgDokpHEc+MGMu3GQZSUVXDV5IX8+p9fcOBoqdelicgpUtDLd5x7Wjpz7zmHW87pxMycPM7/y394+4sdBOMQn4jUjYJevic+OpJfj+rBGxPPonVyLHe8spybpuWw88Axr0sTkXpQ0MsJ9WqbzOzbh3L/JT1ZsGkvIx//hHe/3Ol1WSJykhT0UqNIXwTjh3Xk3bvOJistgdtfXsa9/1jJ4eIyr0sTkTpS0EudZKUlMOvWIdw5oguvL8vj4r9+wopthV6XJSJ1oKCXOovyRfDzH3Tj1QlDKCt3XPG3z3nqww2N8nxdEak/Bb2ctEEdU3l30tmM6t2aR+au55rnFrLncLHXZYnICSjopV6S46L469i+/Pm/zmD51kIuffJTDeWIBCkFvdSbmTEmux2v3zYUX4Rx5d8X8NqSrV6XJSLHqTXozWyKme02s9wqbWPMbJWZVZhZ9YvomMWa2WIzW+k/9oGGLFyCR6+2ybx1xzAGd0rll69/ya//+SXFZeVelyUifnW5op8KjDyuLRe4HJhfw3nFwAjnXB+gLzDSzM6sT5ES/FISopl6wyBuH96ZGYu3cuUzC3WDlUiQqDXonXPzqXwYeNW2Nc65dbWc55xzh/27Uf4fTc8IYb4I476R3fn7tf3ZmH+IS5/8jJUatxfxXEDH6M3MZ2YrgN3APOfcohqOnWBmOWaWU1BQEMiyJMBG9mrNnIlnERsVwVWTFzB31S6vSxIJawENeudcuXOuL5AJDDKzXjUcO9k5l+2cy05PTw9kWdIIurZMYvbtZ9GtVTNueWkpUz792uuSRMJWo8y6cc4VAh/z/bF+CWHpSTG8evOZ/KBnSx58ezV/eHOVbq4S8UDAgt7M0s2suf91HHABsDZQnyfBKS7ax/+9ZgDjh3Vk6uebueXFpRwt0To5Io2pLtMrZwALgG5mlmdm483sMjPLA4YA75jZ+/5j25jZu/5TWwMfmdkXwBIqx+jfDkw3JJj5Ioz7L+nJg6NP58O1+YydrDtpRRqTnhkrjeqDNflMfGUZbZrH8fJNg2mdHOd1SSIhQc+MlaBxfo+WTL9xMAUHi/mvvy1g854jXpckEvIU9NLoBnVMZcaEMzlWWs6YZxawdtdBr0sSCWkKevFEr7bJzLzlTCIMrnpmoRZEEwkgBb14pktGErNuHUpyXBTXPLuQz7/a43VJIiFJQS+eapcazz9uHUKb5nFc/8ISPlq32+uSREKOgl4817JZLDNvGULXjERueXEp/1mvJTBEGpKCXoJCSkI0L40fTOf0RCZMz+HTDRrGEWkoCnoJGikJ0bx802A6piVw0/QlGrMXaSAKegkqqf6wb58az/ipOSzctNfrkkSaPAW9BJ0WiTG8fNOZtE2J48apS1iyeV/tJ4nICSnoJSilJ8Xwys2DaZUcy/VTFrN0y36vSxJpshT0ErQykmKZcfOZpCfFcMMLi1m365DXJYk0SQp6CWotm8Xy4vjBxEX7uO75RWzbd9TrkkSaHAW9BL12qfFMv3EwxWUVXPf8IgoOaYljkZOhoJcmoVurJKZcP5D8g8Vc/8JiDhaVel2SSJOhoJcmY0CHFP52bX/W7TrEzdNyKCot97okkSahLk+YmmJmu80st0rbGDNbZWYVZlbtQvdm1s7MPjKzNf5jJzVk4RKehnfL4NEr+7B48z7unLGcsvIKr0sSCXp1uaKfyvcf6p0LXA7Mr+G8MuDnzrkewJnARDPrWZ8iRaoa3bctf7j0dOatzuc3s3MJxqekiQSTyNoOcM7NN7Os49rWAJhZTeftBHb6Xx8yszVAW2B1/csVqTRuaBZ7Dhfz5IcbaZcaxx0junpdkkjQqjXoG4L/L4p+wKLG+DwJDz+78DTy9h/jkbnraZcaz+i+bb0uSSQoBfzLWDNLBF4H7nbOnfCZcWY2wcxyzCynoEDL1ErtzIyHr+jNoI6p3PuPL7RUgsgJBDTozSyKypB/2Tn3z5qOdc5Nds5lO+ey09PTA1mWhJCYSB+TrxtAZkocN0/P4Ws9bFzkewIW9FY5gP88sMY595dAfY5I8/hoXrhhIBFm3PDCYvYdKfG6JJGgUpfplTOABUA3M8szs/FmdpmZ5QFDgHfM7H3/sW3M7F3/qWcB1wEjzGyF/2dUgPohYa5DiwSe/ekAdhwoYsJ0zbEXqcqCcWpadna2y8nJ8boMaYLe+WInE19ZxqV92vDXsX1rnBkmEkrMbKlzrtr7mnRnrISUi89ozX0ju/HWyh08/dFGr8sRCQqNMr1SpDHddm5n1u86xCNz19OtVTMu7NnS65JEPKUregk5ldMuz+CMzGTufnU56/O1jr2ENwW9hKTYKB/PXDeAuOhIbp6eQ+FRzcSR8KWgl5DVOjmOZ64bwM7CIu54RQugSfhS0EtIG9Ahhf+5rBefbtzDH99d63U5Ip7Ql7ES8q7MbseanQeZ8tnX9GidxJjsdl6XJNKodEUvYeE3o3owrEsav5mdy4pthV6XI9KoFPQSFiJ9ETz1k35kNIvh9peWapkECSsKegkbzeOj+fu1A9hzpIS7ZiynvCL47goXCQQFvYSVXm2T+e/Rp/Ppxj08Nm+91+WINAoFvYSdqwa256rsdjz10Ub+vTrf63JEAk5BL2HpgdGn06ttM+6ZuYIte7WGvYQ2Bb2EpdgoH3+7ZgARZtz60jKOlWhZYwldCnoJW+1S43l8bF/W7jrIb+fkEoxLdos0BAW9hLXzumVw14iuvL4sjxmLt3ldjkhAKOgl7E06vytnd03jgbdWsWbnCZ9fL9Jk1eVRglPMbLeZ5VZpG2Nmq8yswsyqfaLJic4VCTYREcZjV/UlOS6KiS8v43BxmdcliTSoulzRTwVGHteWC1wOzK/HuSJBJy0xhifG9mPz3iP8dvaXGq+XkFJr0Dvn5gP7jmtb45xbV59zRYLVkM4tuPuC05izYgczczReL6EjaMbozWyCmeWYWU5BQYHX5UiYmnheF4Z1SeN3b6xi7S6N10toCJqgd85Nds5lO+ey09PTvS5HwpTPP17fzD9ef0Tj9RICgiboRYJFelIMT1zVl017jnC/5tdLCFDQi1RjaJc0Jp3flX8u384/luZ5XY7IKanL9MoZwAKgm5nlmdl4M7vMzPKAIcA7Zva+/9g2ZvZuTecGphsiDe/OEV0Z0qkFv39jFV8VHPa6HJF6s2D8Z2l2drbLycnxugwR8g8WMfLx+bROjmP2xKHERPq8LkmkWma21DlX7X1NGroRqUHLZrE8MqYPq3ce5OF/6eHi0jQp6EVqcX6Pllw/NIsXPtvMB2u0fr00PQp6kTr41UXd6dG6GffO+oL8g0VelyNyUhT0InUQG+Xjyav7cayknHteW6HnzUqToqAXqaMuGYk88KPT+fyrvTwz/yuvyxGpMwW9yEkYk53JJWe05tG561m2db/X5YjUiYJe5CSYGQ9d1pvWybFMenU5h4pKvS5JpFYKepGTlBwXxRNj+7GjsIjfvbHK63JEaqWgF6mHAR1SuGtEV2Yv386c5du9LkekRgp6kXqaeF5nsjuk8Ns5uWzbd9TrckROSEEvUk+Rvggeu6ovBkx6dTll5RVelyRSLQW9yClolxrPQ5f3ZtnWQp78cKPX5YhUS0Evcop+1KcNl/dvy5MfbiBns56cKcFHQS/SAB4c3YvMlHgmvbqCA8c05VKCi4JepAEkxkTyxNi+7DpYpKdSSdBR0Is0kH7tU7jngq68uXIHszXlUoKIgl6kAd02vAuDslL53Rur2LpXUy4lONTlUYJTzGy3meVWaRtjZqvMrMLMqn2iif+4kWa2zsw2mtmvGqpokWDlizAeG9sXM7j7NU25lOBQlyv6qcDI49pygcuB+Sc6ycx8wNPARUBP4Goz61m/MkWajrbN4/jjZZVTLv+qKZcSBGoNeufcfGDfcW1rnHPrajl1ELDRObfJOVcCvAqMrnelIk3Ipf4pl09pyqUEgUCO0bcFtlXZz/O3VcvMJphZjpnlFBQUBLAskcbxwI9O/3bK5UGtcikeCmTQWzVtJ5xz5pyb7JzLds5lp6enB7AskcaRFBvF41WmXIp4JZBBnwe0q7KfCewI4OeJBJ3+7VOYdH5X3lixQ6tcimcCGfRLgK5m1tHMooGxwJsB/DyRoDTxvC4MzNIql+KdukyvnAEsALqZWZ6ZjTezy8wsDxgCvGNm7/uPbWNm7wI458qAO4D3gTXATOecntIgYccXYfzlysopl3e9upxSTbmURmbBeKt2dna2y8nJ8boMkQb11sod3DljOXec14Vf/LCb1+VIiDGzpc65au9r0p2xIo3k0j5tGDMgk6c/3siCr/Z6XY6EEQW9SCP6w49Op2OLBO55bQX7j5R4XY6ECQW9SCNKiInkr1f3Y++RYn75+hda5VIahYJepJH1apvMfT/sztzV+by8aKvX5UgYUNCLeGD8sI6c3TWN/357NRvyD3ldjoQ4Bb2IByIijEev7ENiTCR3zlhOUWm51yVJCFPQi3gkIymWR8b0Ye2uQ/zx3TVelyMhTEEv4qHzumdw07COTF+whfdyd3pdjoQoBb2Ix+4b2Z0+mcncO+sLLZEgAaGgF/FYdGQET/2kPwB3zFhOSZmWSJCGpaAXCQLtUuP50xVnsHJbIY/Mre2ZPiInR0EvEiRG9W7NtWe2Z/L8TXy4Nt/rciSEKOhFgshvL+5J91ZJ/HzmSnYeOOZ1ORIiFPQiQSQ2ysfT1/SnuKyCSTNWUKYljaUBKOhFgkzn9EQeuqwXizfv4/F/b/C6HAkBCnqRIHRZv0yuzM7kqY82arxeTpmCXiRIPTi6Fz1aN+Oe11Zqfr2ckro8SnCKme02s9wqbalmNs/MNvi3KSc4d5KZ5ZrZKjO7uyELFwl1sVE+/n5tfyqc47aXl2o9HKm3ulzRTwVGHtf2K+AD51xX4AP//neYWS/gZmAQ0Ae4xMy6nlK1ImGmQ4sE/nJlX3K3H+SBt/TIZamfWoPeOTcf2Hdc82hgmv/1NODH1ZzaA1jonDvqf1D4f4DLTqFWkbB0Yc+W3D68MzMWb2Nmzjavy5EmqL5j9C2dczsB/NuMao7JBc4xsxZmFg+MAtqd6A3NbIKZ5ZhZTkFBQT3LEglNP7vwNIZ2bsH9c3JZteOA1+VIExOwL2Odc2uAPwHzgPeAlUBZDcdPds5lO+ey09PTA1WWSJMU6Yvgr1f3IyU+mtteWsaBY6VelyRNSH2DPt/MWgP4t7urO8g597xzrr9z7hwqh380KVikntISY3j6mv7sKDzGz15bQUWFnjcrdVPfoH8TGOd/PQ54o7qDzCzDv20PXA7MqOfniQgwoEMK91/Skw/W7uYv89Z7XY40EXWZXjkDWAB0M7M8MxsPPAxcaGYbgAv9+5hZGzN7t8rpr5vZauAtYKJzbn+D90AkzPx0SAfGDmzHUx9t5K2VO7wuR5qAyNoOcM5dfYJfnV/NsTuo/NL1m/2z61+aiFTHzHhwdC++KjjMvbNWktUigd6ZyV6XJUFMd8aKNEHRkRH87doBpMZHM+HFHHYfKvK6JAliCnqRJiotMYZnx2VTeLSUW19cSnGZ7pyV6inoRZqw09sk88iYPizbWshvZufinGbiyPcp6EWauIvPaM1dI7owa2keUz7b7HU5EoQU9CIh4O4LTuMHPVvy0Dur+fdqLWss36WgFwkBERHGY1f15fQ2ydw5YzkrthV6XZIEEQW9SIhIiIlkyvUDSUuKZvzUJWzZe8TrkiRIKOhFQkh6UgzTbhhEhXOMm7KYvYeLvS5JgoCCXiTEdEpP5Llx2ew8UMRN03M4VqJpl+FOQS8SggZ0SOWJsf1Ysa2Qu15dTrkWQAtrCnqREDWyVyv+cOnpzFudzx/eXKU59mGs1rVuRKTpGjc0ix2Fx3hm/iZSE6K558LTvC5JPKCgFwlxvxzZnX1HSnjigw3ERvm4bXhnr0uSRqagFwlxERHGw1ecQXFZBX96by0xkRHcOKyj12VJI1LQi4QBX4Tx6JV9KC4r58G3VxMb5eMng9t7XZY0En0ZKxImonwRPHl1f87rls5v5nzJ60vzvC5JGkldnjA1xcx2m1lulbZUM5tnZhv825QTnHuPma0ys1wzm2FmsQ1ZvIicnG/WsR/auQX3zlqpJ1SFibpc0U8FRh7X9ivgA+dcV+AD//53mFlb4C4g2znXC/ABY0+pWhE5ZbFRPp79aTbZHVK5+7UV/OvLnV6XJAFWa9A75+YD+45rHg1M87+eBvz4BKdHAnFmFgnEA7p8EAkC8dGRTLlhIH0yk5n4yjJm5mzzuiQJoPqO0bd0zu0E8G8zjj/AObcdeATYCuwEDjjn5ta3UBFpWIkxkbx002DO6pLGfbO+4LlPNnldkgRIwL6M9Y/bjwY6Am2ABDO7tobjJ5hZjpnlFBQUBKosEakiPjqS58ZlM6p3K/7nnTU8Oned7qANQfUN+nwzaw3g3+6u5pgLgK+dcwXOuVLgn8DQE72hc26ycy7bOZednp5ez7JE5GTFRPp48ur+XJXdjic/3Mjv31xFhdbGCSn1Dfo3gXH+1+OAN6o5ZitwppnFm5kB5wNr6vl5IhJAvgjj4St6M+GcTkxfsIWfzVxBaXmF12VJA6nL9MoZwAKgm5nlmdl44GHgQjPbAFzo38fM2pjZuwDOuUXALGAZ8KX/syYHpBcicsrMjF9f1J17f9iNOSt2cNO0HA4WlXpdljQAC8bxuOzsbJeTk+N1GSJh69XFW/ntnFyy0hJ47qfZZKUleF2S1MLMljrnsqv7ne6MFZHvGTuoPS+OH8yew8WMfvozPt+4x+uS5BQo6EWkWkM6t+DNicNo2SyG66Ys5sUFm70uSepJQS8iJ9S+RTyv3zaU4aelc/8bq/jN7C/1JW0TpKAXkRolxUYx+afZ3HpuZ15etJVrn1vErgNFXpclJ0FBLyK18kUYv7qoO49d1Ycv8g5w0RPz+ffqfK/LkjpS0ItInV3WL5O37xpG6+Q4bpqew+/eyKWotNzrsqQWCnoROSmd0xOZPXEo44d1ZPqCLYx+6jPW5x/yuiypgYJeRE5aTKSP+y/pydQbBrL3SDGXPvkpLy7conVygpSCXkTqbXi3DP416RwGd2rB/XNyufrZhXxVcNjrsuQ4CnoROSXpSTFMvX4g/9/lvVm94yAXPf4Jj81bT3GZxu6DhYJeRE5ZRIRx9aD2fPDz4VzUuxVPfLCBix7/hAVf7fW6NEFBLyINKD0phifG9mPajYMorajg6mcX8vOZK9l9SPPuvaSgF5EGd+5p6cy9+1xuH96ZN1ZsZ/ifP+bRues4pNUwPaGgF5GAiIv2cd/I7sy95xzO65bBkx9u5Nw/f8zzn36t8ftGpqAXkYDqlJ7I09f05807zqJ7qyT+++3VjHjkP7y+NI9yPcmqUWg9ehFpNM45Ptmwhz+9t5ZVOw7SMS2Bm87uyBX9M4mN8nldXpNW03r0CnoRaXQVFY5/5e7i7//5ii+3H6BFQjTjhmZx3ZkdSEmI9rq8JumUgt7MpgCXALudc738banAa0AWsBm40jm3/7jzuvmP+UYn4HfOucdrK1hBLxIenHMs2LSXyfM38fG6AuKifFyZncmNwzrSoYWeanUyTjXozwEOA9OrBP3/Avuccw+b2a+AFOfcL2t4Dx+wHRjsnNtSW8EKepHws27XIZ79ZBNvrNhOabljSKcWXDWwHSN7tdKwTh2c8tCNmWUBb1cJ+nXAcOfcTjNrDXzsnOtWw/k/AH7vnDurLgUr6EXCV/7BIv6Rs43Xcraxbd8xmsVGMrpvW64a2I5ebZO9Li9oBSLoC51zzav8fr9zLqWG86cAy5xzT9VwzARgAkD79u0HbNlS64W/iISwigrHwk17eS1nG//K3UVJWQU9WjfjkjNaM6p3azrqgeXf4WnQm1k0sAM43TlXpycV6IpeRKo6cLSUOSu2M2fFdpZvLQSge6skLu7dmlFntKZzeqLHFXqvpqCPrOd75ptZ6ypDN7trOPYiKq/m9TgaEamX5Pgoxg3NYtzQLLYXHuO93F28++VOHp23nkfnradbyyRG9Mhg+Gnp9O+QQpRPtwhVVd8r+j8De6t8GZvqnLvvBOe+CrzvnHuhrkXpil5E6mLngcrQfy93F0u37KeswpEUG8nZXdMY3q0y+DOaxXpdZqM41Vk3M4DhQBqQD/wemAPMBNoDW4Exzrl9ZtYGeM45N8p/bjywDejknDtQ14IV9CJysg4WlfL5xj18tLaAj9fvJv9gMQCntUxkSKcWDOncgsEdW4TsPH3dMCUiYcU5x9pdh/jP+gI+/2ovOZv3cbSkcn2d7q2Svg39gVkptEiM8bjahqGgF5GwVlpewRd5hSz4ai8LNu0lZ/N+issqAOiUnsDADqlkZ6UwMCuVDi3iMTOPKz55CnoRkSqKy8rJ3X6AJZv3k7N5H0s27+fAscollNMSYxjQoTnZHVIZkJVCrzbJREcG/5e7gZh1IyLSZMVE+hjQIZUBHVLh3M5UVDg2Fhxm8df7WLZlPzlb9vP+qnz/sRH0yWzOgKwUBmVVhn+z2CiPe3BydEUvIlKN3QeLWLplP0v9wZ+7/QBlFQ4z6NGqGYM6pjKoYyoDs1JJT/J+nF9DNyIip+hYSTnLt+5n8eZ9lVf+W/dTVFo5zt+tZRLDuqYxrGsagzumEh/d+IMlCnoRkQZWWl5B7vYDLNy0j8827mHx5n2UlFUQ7YtgQIcUhnVN44IeLTmtZWKjfLmroBcRCbCi0nIWf72PTzfuYf76AtbuOgRAp7QERvZqxUW9WtOrbbOAhb6CXkSkke0+WMTc1fm8l7uLBZv2Ul7hyEyJY+TprRjdty29Mxt2JU4FvYiIh/YfKWHe6nzeW7WLTzfsoaS8gjMyk7lmcHsu7dOmQcb0FfQiIkHiwLFS5izfzsuLtrA+/zBJsZFc0T+Tnwxuz2ktk+r9vgp6EZEg45wjZ8t+Xlq4hX99uYuS8goGd0xl+vhBxESe/BO1dMOUiEiQMTMGZlXOw//dJcXMWprH13uO1Cvka6OgFxHxWIvEGG45t3PA3j/4F3AQEZFToqAXEQlxCnoRkRCnoBcRCXG1Br2ZTTGz3WaWW6Ut1czmmdkG/zblBOc2N7NZZrbWzNaY2ZCGLF5ERGpXlyv6qcDI49p+BXzgnM2DXosAAAQ/SURBVOsKfODfr84TwHvOue5AH2BNPesUEZF6qjXonXPzgX3HNY8GpvlfTwN+fPx5ZtYMOAd43v8+Jc65wlOqVkRETlp9x+hbOud2Avi3GdUc0wkoAF4ws+Vm9pyZJZzoDc1sgpnlmFlOQUFBPcsSEZHj1WkJBDPLAt52zvXy7xc655pX+f1+51zKcedkAwuBs5xzi8zsCeCgc+7+OnxeAbCllsPSgD21Fh961O/won6Hl1PpdwfnXHp1v6jvnbH5ZtbaObfTzFoDu6s5Jg/Ic84t8u/P4sRj+d9xomKrMrOcE63rEMrU7/CifoeXQPW7vkM3bwLj/K/HAW8cf4Bzbhewzcy6+ZvOB1bX8/NERKSe6jK9cgawAOhmZnlmNh54GLjQzDYAF/r3MbM2ZvZuldPvBF42sy+AvsAfG7oDIiJSs1qHbpxzV5/gV+dXc+wOYFSV/RVAoP75NTlA7xvs1O/won6Hl4D0OyjXoxcRkYajJRBEREKcgl5EJMQ1uaA3s5Fmts7MNppZnaZrNkWnssZQU2Zm7czsI//aSKvMbJK/PaT7bmaxZrbYzFb6+/2Avz2k+/0NM/P5b6x8278fLv3ebGZfmtkKM8vxtzV435tU0JuZD3gauAjoCVxtZj29rSpgplL/NYaasjLg5865HsCZwET//+NQ73sxMMI514fKGWojzexMQr/f35jEd9fCCpd+A5znnOtbZf58g/e9SQU9MAjY6Jzb5JwrAV6lct2dkFPfNYaaOufcTufcMv/rQ1T+4W9LiPfdVTrs343y/zhCvN8AZpYJXAw8V6U55Ptdgwbve1ML+rbAtir7ef62cFGXNYZChn/pjX7AIsKg7/7hixVU3mk+z39Xecj3G3gcuA+oqNIWDv2Gyr/M55rZUjOb4G9r8L43tYeDWzVtmh8agswsEXgduNs5d9Csuv/1ocU5Vw70NbPmwGwz6+V1TYFmZpcAu51zS81suNf1eOAs59wOM8sA5pnZ2kB8SFO7os8D2lXZzwR2eFSLF/L9awtRwxpDTZ6ZRVEZ8i875/7pbw6LvgP4l/P+mMrvaEK932cBPzKzzVQOxY4ws5cI/X4D395kinNuNzCbyuHpBu97Uwv6JUBXM+toZtHAWCrX3QkXta4x1NRZ5aX788Aa59xfqvwqpPtuZun+K3nMLA64AFhLiPfbOfdr51ymcy6Lyj/PHzrnriXE+w1gZglmlvTNa+AHQC4B6HuTuzPWzEZROabnA6Y45x7yuKSA8K8xNJzKZUvzgd8Dc4CZQHtgKzDGOXf8F7ZNmpkNAz4BvuT/H7P9P1SO04ds383sDCq/ePNReQE20zn3oJm1IIT7XZV/6OYXzrlLwqHfZtaJyqt4qBxGf8U591Ag+t7kgl5ERE5OUxu6ERGRk6SgFxEJcQp6EZEQp6AXEQlxCnoRkRCnoBcRCXEKehGREPf/ADlhQDo/YZYeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "normalised_transverse_rpd_s8 = transverse_rpd_s8[calculation_points > 0.] / calculation_points[calculation_points > 0.]\n",
    "norm_rpd_calculation_points = calculation_points[calculation_points > 0.]\n",
    "plt.plot(norm_rpd_calculation_points, normalised_transverse_rpd_s8)\n",
    "np.save('..\\data-perpl\\ACTN2-Affimer\\\\normalised_transverse_rpd_smoothed_Churchman-8', normalised_transverse_rpd_s8)\n",
    "# normalised_transverse_rpd = np.load('normalised_transverse_rpd_smoothed_Churchman-4p4.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(params_optimised,\n",
    " params_covar,\n",
    " params_1sd_error) = zdisk_modelling.fitmodel_to_hist(\n",
    "    norm_rpd_calculation_points[0:31],\n",
    "    normalised_transverse_rpd_s8[0:31],\n",
    "    trans_model_with_info.model_rpd,\n",
    "    trans_model_with_info.initial_params,\n",
    "    trans_model_with_info.param_bounds\n",
    "    )\n",
    "print('')\n",
    "print('Initial parameter guesses:')\n",
    "print(trans_model_with_info.initial_params)\n",
    "print('')\n",
    "print('Parameter bounds:')\n",
    "print(trans_model_with_info.param_bounds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(norm_rpd_calculation_points,\n",
    "    normalised_transverse_rpd_s8)\n",
    "zdisk_plots.plot_fitted_model(\n",
    "    norm_rpd_calculation_points,\n",
    "    fitlength,\n",
    "    params_optimised,\n",
    "    params_covar,\n",
    "    trans_model_with_info,\n",
    "    plot_95ci=False\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitlength = 50.\n",
    "calculation_points = np.arange(fitlength + 1.)\n",
    "combined_precision=5.\n",
    "transverse_rpd_s5 = plotting.estimate_rpd_churchman_2d(\n",
    "    input_distances=trans_distances[trans_distances < (fitlength + 3 * combined_precision)],\n",
    "    calculation_points=calculation_points,\n",
    "    combined_precision=combined_precision\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXhV5b328e9vZ55IGMKUQMJkGCJjmARRUCtS59ZWKqjVljprB23t6enk8bQ97zlqtR4tCmIda52qqFhEFEEFEkAIJAwyhilhTAhkft4/EjmUAgkhO2sP9+e6uLL32mtn34/KvZdreJY55xARkdDl8zqAiIj4l4peRCTEqehFREKcil5EJMSp6EVEQlyk1wFOpEOHDi4zM9PrGCIiQSMvL2+Pcy71RK8FZNFnZmaSm5vrdQwRkaBhZltO9pp23YiIhDgVvYhIiFPRi4iEOBW9iEiIU9GLiIQ4Fb2ISIhT0YuIhLiAPI9e/llVTR0riw6wfOsBMjskMLZ3B+KiI7yOJSJBQkUfgL4q9s837uXzjfvI3bKPiuq6o6/HRPo4p1d7JvTtyPi+HUlvG+9hWhEJdCr6APNK7jZ+89ZqyqtqAejbOYlrh3dndK/2DOmewvrdh5hXUMy8wt3M//tq+Ptq+nZO4rtjMvlWTjfMzOMRiEigsUC8w1ROTo4LtykQqmrqeGD2Gp77fAuje7bnhnMyGdmjHW0Tok+4vnOOjXvK+bCgmNkrd/BF0UGuGpLGg1dlEx+t72+RcGNmec65nBO9pkYIAMVlFdz2/DJyt+znB+N6cu/FWURGnPo4uZnRKzWRXqmJ3Dy2B4/P38BDH6xj9Y6DPDFlGL1SE1spvYgEOp1147FlW/dz2WMLyd9xkEcnD+H+Sf0aLfnj+XzGnRf04bmbRrLnUBWXP7aQ2St3+CmxiAQbFb2HXl6ylWv//DnRkT5ev3UMlw/qeka/b2yfDrxz11iyOidxx4vL+fVbq6mqqWv8jSIS0lT0Hnl8/gZ+9voqRvZsx9t3jKV/1zYt8nu7JMfx8rTRfHdMJrM+3cyUGYs50nBgV0TCk4reA/PXFvPf/1jL5YO6Muu7I0iJP/EB1+aKjvTxq8sG8PC3B7F08z7ueHEZNbXashcJVyr6VrZlbzl3v7Scvp3b8IdvDCTC57/TIa8aks5vr8hmXmExP39jFYF4hpWI+J/OumlFR6pqueX5ZQD8ecqwVrm6deqoDEpKK3j0ww2kJsVw78V9/f6ZIhJYVPStxDnH/a+vpHBXKTNvHE739q13NesPLzqLkkOVPD7/S1ITY7hxTI9W+2wR8Z6KvpU8++lm3lyxgx9fdBbjszq26mebGQ9ckc2eQ1X8ZvYa2ifGcNkZnuEjIsFD++hbwZJN+/iPdwq4sF8nbh/f25MMkRE+Hps8hJyMtvzolRUs2rDHkxwi0vpU9H62u7SC215YRrd28Tz07UH4/HjwtTGxURE8ff1wenZI5AfP5bGhuMyzLCLSelT0fuSc40evrOBwVQ1PThlGm9goryORHB/FrJuGExvlY9pzeZRVVHsdSUT8TEXvR28s386iDXv5+aR+ZHVO8jrOUV2S43hs8lC27D3Mfa+u1GmXIiFORe8n+8ur+I93ChjaPYXvjOjudZx/MbpXe342sS/v5e9i+oKNXscRET9S0fvJ794roPRINf959dme7pc/le+d24Ovn92FP8wp5FMdnBUJWSp6P/h8415eyS3ie+f2pG/nlpnDxh/MjD98cyA9UxO546Xl7DhwxOtIIuIHKvoWVllTy7+9sYr0tnHcfUEfr+M0KjEmkienDKOqpo5bX1hGZY0mQBMJNSr6FvbnjzfyZUk5D1yZHTQ38O7dMZH/vmYgX2w7wG/eXuN1HBFpYY0WvZnNNLNiM8s/Ztk1ZrbazOrM7IS3rjpm3QgzW25ms1sicCDbtKecP83fwKUDu7T61a9namJ2F245rxcvLt7KK7nbvI4jIi2oKVv0s4CJxy3LB64GFjTh/XcDBacXK/g45/jFm6uIifDxy0v7ex2nWX7ytbMY07s9v3gzn/ztB72OIyItpNGid84tAPYdt6zAObe2sfeaWTrwdeDpZicMEm+uqD9n/r5L+tKxTazXcZolMsLHo9cOoX1CNLc8n8f+8iqvI4lIC/D3PvpHgPuARu96YWbTzCzXzHJLSkr8HKtllVZU8x+zCxjcLYXrAvCc+dPRPjGGJ6YMo7i0krv/uoLaOl1MJRLs/Fb0ZnYpUOycy2vK+s656c65HOdcTmpqqr9i+cXjH25g3+EqHrgiO2DPmT8dg7ul8KvL+7NgXQl/nLfe6zgicob8uUU/BrjczDYDLwMTzOx5P36eJ7bsLeeZRZv5xtB0zk5P9jpOi/nOiO5cMyydR+etZ17Bbq/jiMgZ8FvRO+fud86lO+cygWuBD51zU/z1eV753buFREYY916c5XWUFmVmPHBlNgO6tuGHf13Blr3lXkcSkWZqyumVLwGfAVlmVmRmN5vZVWZWBIwG3jGz9xvW7Wpm7/o3cuD4fONe5qzexa3n9aJTkB6APZXYqAienDIMM+MHz+VxpEoXU4kEIwvEmQtzcnJcbm6u1zFOqbbOcfmfFrK/vIoPf3I+sVHBcXFUc8xfW8xNs5YyKbsLj00eEhLHIURCjZnlOedOeF2TroxtpteWFbF6Ryk/vaRvSJc8wPisjvz8kn68s2onv59T6HUcETlNumdsMxyqrOH/vb+WId1TuDxM7r36vXN7ULT/MNMXbCQtJY4bzsn0OpKINJGKvhme+GgDJWWVTJ9av/86HJgZv7xsANsPVPCbt1fTJTmWrw3o7HUsEWkC7bo5TUX7D/PUJ5u4cnBXhnRv63WcVhXhMx6bPISz01O46+XlrNh2wOtIItIEKvrT9Pv3CvEZ3Dexr9dRPBEXHcGMG3JITYrhe88uZevew15HEpFGqOhPQ+7mfcxeuZNp43rRNSXO6zie6ZAYw6zvjqCmznHjrCWaE0ckwKnom6iuzvHrt1fTuU0st5zX0+s4nuuVmshT1+dQtP8INz6zhH0qe5GApaJvolfzisjfXsr9k/oSH61j2ADDM9vxxHVDKdhVxjVPfqpbEYoEKBV9E5RVVPNf7xcyLKNt2JxO2VQX9OvEX24aQXFpJd984lM2FB/yOpKIHEdF3wR/+nADe8ur+NVl/cPmdMrTMapne16aNoqq2jq+9efPWFmks3FEAomKvhGb9pQzc9EmrhmWzsD0FK/jBKzstGT+dss5xEVFMHn653z65R6vI4lIAxV9Ix58Zw0xkRH8JMRmp/SHHh0SeO3Wc0hrG8eNM5fy3qqdXkcSEVT0p/TxuhI+KCjmzgm96ZgUerNT+kPn5Fhe+cFostPacOsLy7j/9VUcqqzxOpZIWFPRn0R1bR0PzF5DZvt4bhyT6XWcoJISH82L3x/FD8b15OWlW7n44QUs2qBdOSJeUdGfxPOfb2FD8SF+8fX+xESG9uyU/hAbFcH9k/rx6i2jiY70cd3Ti/n3N/Mp19a9SKtT0Z/AvvIqHp67jnP7dOCCfh29jhPUhmW04927zuWmMT14fvEWJv5xAZ9v3Ot1LJGwoqI/gQdmr+FwVS2/vFSnU7aEuOgIfnlZf17+/igM49rpnzN1xmI+WV9CIN74RiTUqOiPM7+wmDeWb+e28b3p0ynJ6zghZWTP9sy551zum5hF4a4yps5YwqRHF/LG8iKqa+u8jicSsnQrwWOUVVTztYcXkBQbyew7zyU6Ut+D/lJZU8vfl+9g+icb2VB8iK7Jsdw0tgdXDkmjQ2KM1/FEgs6pbiWooj/Gz99YxctLtvL6bWMY3E0XR7WGujrHR+uK+fPHG1m8aR8+g5zMdkwc0JmJ2Z3DepZQkdOhom+CT7/cw3eeWsy0cT35+aR+rfrZUq9wVynvrdrFnPxdrN1dBsCg9GQuzu7M1/p3oldqoo6ZiJyEir4Rh6tqmPjIJ/gM3rt7HHHROp3SaxtLDvH+6t3MWb2LLxruZNW9XTwX9OvIhf06MTyznXatiRxDRd+IB2avYcbCTfx12ihG9mzfap8rTbPjwBHmFRbzYcFuFn25l6qaOpJiIhl3VioTszvztQGddK2DhL1TFX3YT6y+bOt+Zi7axNRRGSr5ANU1JY6pozKYOiqDw1U1LFy/hw8Li5lXWMw7q3bSPiGab+akM3l4dzI7JHgdVyTghPUWfWVNLV9/dCFHqmp5/4fjSIwJ+++9oFJX51iwvoQXF29lXmExtXWOsb078J2R3bmofyeiIrRrR8KHtuhP4g/vrWVD8SGevWmESj4I+XzG+VkdOT+rI7tLK/jr0m28vGQrt72wjC7Jsfzqsv5cPKCzDuBK2AvbTZ6Xl2xl5qJN3HhOJuedlep1HDlDndrEctcFffjkpxOYeWMOKfHR3PL8Mr7/l1yK9h/2Op6Ip8Ky6D/9cg+/eDOfcWel8ouv61TKUBLhMyb07cTbd4zh3yb1Y9GGvVz00AKeWrCRGl19K2Eq7Ip+055ybn1+GT06JPCn7wwhUvtxQ1JkhI/vj+vJ3B+NY0zv9jz4bgGX/WkRy7fu9zqaSKsLq5Y7eLiam2ctJcJnzLhhOG1io7yOJH6W3jaep67P4ckpw9hfXsU3nviUl5Zs9TqWSKsKm6Kvrq3jthfz2Lb/ME9OGUb39vFeR5JWYmZMzO7M3B+N47yzUrn/9VU8+fGXXscSaTVhUfTOOX791moWbdjLf151NiN6tPM6knggKTaKP0/N4bJBXfn9e4X8YU6hpkmWsBDy5xQ653jy4428sHgrPzivJ9fkdPM6kngoOtLHI98eTFJsJE989CWlR6p54IpsfD6dgimhK6SL/lBlDfe/voq3v9jBpLM789OL+3odSQJAhM948MpskuOieOKjLymrqOF/vjVIF1hJyArZoi/YWcrtLyxj895y7r04i1vP66WtNjnKzPjpxL4kxUbyX3PWcqiyhv+9biixUZozR0JPyG3COOd4eclWrnx8EYcqa3jx+6O4fXxvlbyc0G3n9+bBq7KZv7aY+19fpX32EpJCaou+vLKGX7yZzxvLtzO2dwce/vZgUpN0tyI5tetGZrDvUBX/M3cdQzPaMnVUhteRRFpUo1v0ZjbTzIrNLP+YZdeY2WozqzOzE89/bNbNzOabWUHDune3ZPDjHThcxRWPL+LNFdv54YVn8exNI1Ty0mS3j+/N+KxUHnh7DSsa5r8XCRVN2XUzC5h43LJ84GpgwSneVwP82DnXDxgF3G5m/ZsTsimS46IY27sDL9w8krsv7EOEdtXIafD57Oj/Ad7+wjL2l1d5HUmkxTRa9M65BcC+45YVOOfWNvK+nc65ZQ2Py4ACIO0Msp6SmfHrywdwTu8O/voICXEp8dE8MWUoJWWV3PPXFdTVaX+9hIZWORhrZpnAEGDxKdaZZma5ZpZbUlLSGrFE/sXA9BR+dXl/Pl5XwmMfbvA6jkiL8HvRm1ki8Bpwj3Ou9GTrOeemO+dynHM5qamaNli8850R3bl6SBqPzFvHx+u00SHBz69Fb2ZR1Jf8C8651/35WSItxcx48KqzyeqUxD0vL2f7gSNeRxI5I34requ/rc8MoMA595C/PkfEH+KiI3hiyjCqax33/u0LnV8vQa0pp1e+BHwGZJlZkZndbGZXmVkRMBp4x8zeb1i3q5m92/DWMcBUYIKZrWj4M8lP4xBpcT06JPDTiVl8+uVeZq/c6XUckWYL65uDizSmts5xxeMLKSmrZN6Pz9e9hSVgnerm4CE3BYJIS4rwGQ9ckc3u0koenbfe6zgizaKiF2nEkO5tuXZ4N2Yu3MS63WVexxE5bSp6kSa4b2JfEmIi+eXf83VgVoKOil6kCdolRHPfxCw+37iPt77Y4XUckdOiohdpomuHd2dgejIPvlNAWUW113FEmkxFL9JEXx2YLTlUySMf6MCsBA8VvchpGNQthckjujPr080U7jrpjB4iAUVFL3Ka7v1aFm1iI/n1W6t1YFaCgope5DS1TYjmngvP4vON+/j0y71exxFplIpepBm+PbwbXZJjeWjuOm3VS8BT0Ys0Q2xUBLeP703elv18sn6P13FETklFL9JM38rpRlpKnLbqJeCp6EWaKTrSxx0TerNi2wE+WqsblEjgUtGLnIFvDkunW7s4Hv5AW/USuFT0ImcgKsLHneP7sLLoIPMKir2OI3JCKnqRM3TV0DQy2sdrq14Clope5AxFRfi4c0IfVu8o5R9rdnsdR+RfqOhFWsCVg7vSo0MCD89dR12dtuolsKjoRVpAZISPuy/oQ+GuMuas3uV1HJF/oqIXaSGXDepKr9QEHvlAW/USWFT0Ii0kwmfcdUEf1u0+xIeFOgNHAoeKXqQFTTq7C2kpcUxfsNHrKCJHqehFWlBUhI+bxvZgyeZ9LN+63+s4IoCKXqTFXTu8G21iI3nqE23VS2BQ0Yu0sISYSKaMymBO/i627C33Oo6Iil7EH248J5NIn4+nP9nkdRQRFb2IP3RsE8uVQ7ryt7xt7Cuv8jqOhDkVvYiffP/cnlRU1/HcZ1u8jiJhTkUv4id9OiUxoW9H/vLZZiqqa72OI2FMRS/iR9PG9WRveRWvLSvyOoqEMRW9iB+N7NGOQenJPP3JJmo1LYJ4REUv4kdmxrRxvdi0p5y5msJYPKKiF/Gziwd0olu7OF1AJZ5R0Yv4WWSEj++N7Unelv3kbdG0CNL6VPQireCanHTaxEYyc5EuoJLWp6IXaQXx0ZFMHtmdOfm72H7giNdxJMw0WvRmNtPMis0s/5hl15jZajOrM7OcU7x3opmtNbMNZvazlgotEoyuH50JwF8+2+xlDAlDTdminwVMPG5ZPnA1sOBkbzKzCOBx4BKgPzDZzPo3L6ZI8EtLiWPigM68tHgrh6tqvI4jYaTRonfOLQD2HbeswDm3tpG3jgA2OOc2OueqgJeBK5qdVCQE3DQ2k9KKGl5btt3rKBJG/LmPPg3YdszzooZlJ2Rm08ws18xyS0pK/BhLxDtDu7dlUHoyzyzapPvKSqvxZ9HbCZad9L9s59x051yOcy4nNTXVj7FEvGNm3DS2BxtLyvl4vTZopHX4s+iLgG7HPE8Hdvjx80SCwiXZXejUJoaZC3WqpbQOfxb9UqCPmfUws2jgWuAtP36eSFCIjvRx/ehMPlm/h3W7y7yOI2GgKadXvgR8BmSZWZGZ3WxmV5lZETAaeMfM3m9Yt6uZvQvgnKsB7gDeBwqAV5xzq/01EJFgMnlEd2IifTyzaLPXUSQMRDa2gnNu8kleeuME6+4AJh3z/F3g3WanEwlR7RKiuXpoGq8vK+K+i7NomxDtdSQJYboyVsQjN57Tg8qaOl5cstXrKBLiVPQiHsnqnMTY3h34y2ebqa6t8zqOhDAVvYiHbhqbye7SSt5dtdPrKBLCVPQiHjr/rI707JDAjIWbcE4XUIl/qOhFPOTz1V9AtbLoIEs27Wv8DSLNoKIX8dg3hqaTEh/F07qASvxERS/isbjoCKaMzOCDgt1s2lPudRwJQSp6kQBw/TkZRPl8PKM7UIkfqOhFAkDHpFguH9yVv+UWceBwlddxJMSo6EUCxM1je3CkupYXFusCKmlZKnqRANGvSxvO7dOBZz/dTFWNLqCSlqOiFwkgN4/tQXFZJbNXakZvaTkqepEAct5ZqfTpmMjTn+gCKmk5KnqRAGJm3Dy2B2t2lvLZxr1ex5EQoaIXCTBXDkmjfUI0T3+iUy2lZajoRQJMbFQEU0dn8GFhMRuKD3kdR0KAil4kAE0ZlUF0pI8ZCzd6HUVCgIpeJAB1SIzhmmHpvJa3nZ0Hj3gdR4Kcil4kQN1yXi9qnWP6Am3Vy5lR0YsEqG7t4rlqSBovLdlKSVml13EkiKnoRQLYbef3orKmjhmawljOgIpeJID1TE3k0oFdee6zzZrsTJpNRS8S4G4f34vyqlqeWbTZ6ygSpFT0IgGub+c2XDygE88s2kRZRbXXcSQIqehFgsAd4/tQWlHDXz7b4nUUCUIqepEgcHZ6MudnpTJj4SYOV9V4HUeCjIpeJEjcOaE3+8qreFE3JpHTpKIXCRLDMtoxumd7pi/YSEV1rddxJIio6EWCyJ0TelNcVsnf8oq8jiJBREUvEkRG92rP0O4pPDF/g7bqpclU9CJBxMy49+K+7DhYoatlpclU9CJBZnSv9lzUvxP/O3+D5sCRJlHRiwSh+y/pS2VNHQ/NXed1FAkCKnqRINQzNZGpozP469KtrN1V5nUcCXAqepEgdfcFfUiMieTBdwu8jiIBTkUvEqRS4qO564I+LFhXwkdri72OIwGs0aI3s5lmVmxm+ccsa2dmc81sfcPPtid57w/NbLWZ5ZvZS2YW25LhRcLd9aMzyWwfz4PvFFBTW+d1HAlQTdminwVMPG7Zz4B5zrk+wLyG5//EzNKAu4Ac51w2EAFce0ZpReSfREf6+Nkl/VhffIiXl27zOo4EqEaL3jm3ANh33OIrgGcbHj8LXHmSt0cCcWYWCcQDO5qZU0RO4uIBnRjRox0Pz12naYzlhJq7j76Tc24nQMPPjsev4JzbDvw3sBXYCRx0zv2juUFF5MTMjH//en/2llfxvx996XUcCUB+OxjbsN/+CqAH0BVIMLMpp1h/mpnlmlluSUmJv2KJhKSz05O5emgaMxZuYkPxIa/jSIBpbtHvNrMuAA0/T3TI/0Jgk3OuxDlXDbwOnHOyX+icm+6cy3HO5aSmpjYzlkj4+tklfUmIjuCevy6nqkYHZuX/NLfo3wJuaHh8A/D3E6yzFRhlZvFmZsAFgE74FfGTjkmx/O7qgeRvL+WP83TFrPyfppxe+RLwGZBlZkVmdjPwe+AiM1sPXNTwHDPrambvAjjnFgOvAsuAVQ2fNd0voxARACZmd+aaYek88dGX5G4+/hwKCVfmnPM6w7/Iyclxubm5XscQCUqHKmu45I8LAHj3rnNJio3yOJG0BjPLc87lnOg1XRkrEmISYyJ5+FuD2b7/CL99e43XcSQAqOhFQlBOZjtuPb8Xf8srYk7+Tq/jiMdU9CIh6u4LzuLstGTuf30VxaUVXscRD6noRUJUdKSPh789mCPVtdz76koC8XictA4VvUgI690xkZ9P6sfH60r4n3/olMtwFel1ABHxr6mjMlizo5Q/zd9AYmwkt5zXy+tI0spU9CIhzsx48KqzOVRZw+/fKyQhJpKpozK8jiWtSEUvEgYifFa/v76qln9/M5+E6AiuHprudSxpJdpHLxImoiJ8PH7dUM7p1Z57X13JnPxdXkeSVqKiFwkjsVERPHV9DgPTk7nrpeUsWKeZYsOBil4kzCTERDLrxhH06pjItOdy+VhlH/JU9CJhKDk+iuduHkFGuwRumLmEB2avoaK61utY4icqepEw1SExhjdvH8PUURnMWLiJKx9fxNpdZV7HEj9Q0YuEsbjoCB64MptnbhzOnkOVXPanhcxcuIm6Ol1FG0pU9CLC+L4dmXPPOM7t3YHfzl7DDc8sYefBI17HkhaiohcRoH5XztM35PDgVdks3byPc/8wn9tfXMZnX+7VPDlBThdMichRZsZ1IzMY1yeVWZ9u5tW8It5ZuZNeqQlcNzKDbwxLJzlONzIJNrrDlIicVEV1LW9/sYPnF2/li20HiI3ycVH/zgxKT6Z/1zYM6JJMcryKPxCc6g5TKnoRaZL87Qd5YfEW5hUUU1xWeXR5ets4BnRtQ1bnNrSNjyIxJpKk2CiSYiNJio0kISaSSJ8dXd845rGBz2f4DCLMMKt/HBMVQXxUBL5j3iendqqi164bEWmS7LRkfnf1QABKyipZs7OU1TsOsnpHKWt2lPKPNbtpye1GM0iMrv+ySIyt//JIjouia0osaSnxpLWNIy0ljvS2caQmxuhL4RRU9CJy2lKTYjgvKZXzzko9uqy6to5DFTWUVdRQVll99PGhyhrqGr4Bjv0icIBzDueg1jnqnKPO1S+rqK6tf39lzT/9nl0HK8jbsp+DR6r/KU90pI8+HRPp27kN/bok0a9LG/p2TqJ9Ykxr/OMIeCp6EWkRURE+2iZE0zYh2u+fdaiyhu37j7D9wGG27z/C1n2HKdxVxoL1Jby2rOjoeqlJMQxKT2FoRgpDu7dlYHoy8dHhV3vhN2IRCXqJMZFkdU4iq3PSv7y251AlhTvLKNxVv0tpxbYDfFCwG6ifrrlflySGdm/LsIy2DM9sR9eUuNaO3+p0MFZEQt7+8iqWb9vPsi0HWLZ1P19sO0B5Vf3cPmkpcYzo0Y7hme0YntmW3h0TMQu+/f06GCsiYa1tQjQT+nZiQt9OANTWOQp2lrJ08z6Wbt7HJ+v38Mby7fXrxkfVb/FntiUnox0D05OJjYrwMv4Z0xa9iIQ95xyb9x5m6ab64s/bup+NJeUAREUYA7omMyyjLYO6pTA4PYVu7eICbqtf59GLiJymfeVV5G3ZT+6WfeRt3s/K7QepqqkD6rf6B6anMKhbCoPSkzk7LZmObWI9zauiFxE5Q1U1dazbXcYXRQf4YtsBVhYdZN3uMr6a6DM1KYbsrm3ITktmQNdkstPakJbSelv+2kcvInKGoiN9ZKclk52WzHUjMwA4XFXD6h2l5G8/SP72+gvIFqzfQ21D+7dPiGZwt5T6P93r/w+gTWzrTxmhohcRaab46MiGs3XaHV1WUV1L4a4yVhUdYMW2g6zYtp95hcVHX++VmsConu25sH8nRvds3yoHerXrRkTEzw4eqWZl0QFWbD3A8m0H+HzjXg5X1RIfHcG4Pqlc2L8T47NSz+hKXu26ERHxUHJcFOf2SeXcPvVTRlRU1/L5xr18ULCbD9YUM2f1LsxgeGY7XvzeSCIjWvZWISp6EZFWFhsVwflZHTk/qyMPXOFYvaOUuWt2s7u0osVLHlT0IiKeMrOjB3n9RbcSFBEJcSp6EZEQp6IXEQlxjRa9mc00s2Izyz9mWTszm2tm6xt+tj3Je1PM7FUzKzSzAjMb3ZLhRUSkcU3Zop8FTDxu2c+Aec65PsC8hucn8kdgjnOuLzAIKGhmThERaaZGi945twDYd9ziK4BnGx4/C1x5/PvMrA0wDpjR8HuqnHMHziitiKJipVYAAAQKSURBVIictubuo+/knNsJ0PCz4wnW6QmUAM+Y2XIze9rMEk72C81smpnlmlluSUlJM2OJiMjx/HkwNhIYCjzhnBsClHPyXTw456Y753KcczmpqaknW01ERE5Tcy+Y2m1mXZxzO82sC1B8gnWKgCLn3OKG569yiqI/Vl5e3h4z29LIah2APU1OHDo07vCicYeXMxl3xsleaG7RvwXcAPy+4effj1/BObfLzLaZWZZzbi1wAbCmKb/cOdfoJr2Z5Z5sAp9QpnGHF407vPhr3E05vfIl4DMgy8yKzOxm6gv+IjNbD1zU8Bwz62pm7x7z9juBF8xsJTAY+M+WHoCIiJxao1v0zrnJJ3npghOsuwOYdMzzFUDYfSuLiASSYL4ydrrXATyicYcXjTu8+GXcAXnjERERaTnBvEUvIiJNoKIXEQlxQVf0ZjbRzNaa2QYza9J5+cHoTCaTC2Zm1s3M5jdMgrfazO5uWB7SYzezWDNbYmZfNIz7Nw3LQ3rcXzGziIYr6Gc3PA+XcW82s1VmtsLMchuWtfjYg6rozSwCeBy4BOgPTDaz/t6m8ptZNH8yuWBWA/zYOdcPGAXc3vDvONTHXglMcM4Nov5U5IlmNorQH/dX7uafJz0Ml3EDjHfODT7m/PkWH3tQFT0wAtjgnNvonKsCXqZ+grWQ09zJ5IKdc26nc25Zw+My6v/ypxHiY3f1DjU8jWr44wjxcQOYWTrwdeDpYxaH/LhPocXHHmxFnwZsO+Z5UcOycNGUyeRChpllAkOAxYTB2Bt2X6ygfkqRuQ3Th4T8uIFHgPuAumOWhcO4of7L/B9mlmdm0xqWtfjYg+3m4HaCZTo/NASZWSLwGnCPc67U7ET/6kOLc64WGGxmKcAbZpbtdSZ/M7NLgWLnXJ6Zne91Hg+Mcc7tMLOOwFwzK/THhwTbFn0R0O2Y5+nADo+yeGF3wyRynGIyuaBnZlHUl/wLzrnXGxaHxdgBGu7b8BH1x2hCfdxjgMvNbDP1u2InmNnzhP64gaOzCeCcKwbeoH73dIuPPdiKfinQx8x6mFk0cC31E6yFi68mk4OTTCYX7Kx+030GUOCce+iYl0J67GaW2rAlj5nFARcChYT4uJ1z9zvn0p1zmdT/ff7QOTeFEB83gJklmFnSV4+BrwH5+GHsQXdlrJlNon6fXgQw0zn3oMeR/KJhMrnzqZ+2dDfwK+BN4BWgO7AVuMY5d/wB26BmZmOBT4BV/N8+259Tv58+ZMduZgOpP/AWQf0G2CvOud+aWXtCeNzHath18xPn3KXhMG4z60n9VjzU70Z/0Tn3oD/GHnRFLyIipyfYdt2IiMhpUtGLiIQ4Fb2ISIhT0YuIhDgVvYhIiFPRi4iEOBW9iEiI+//k2hgPCrz6YwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "normalised_transverse_rpd_s5 = transverse_rpd_s5[calculation_points > 0.] / calculation_points[calculation_points > 0.]\n",
    "norm_rpd_calculation_points = calculation_points[calculation_points > 0.]\n",
    "plt.plot(norm_rpd_calculation_points, normalised_transverse_rpd_s5)\n",
    "np.save('normalised_transverse_rpd_smoothed_Churchman-5', normalised_transverse_rpd_s5)\n",
    "# normalised_transverse_rpd = np.load('normalised_transverse_rpd_smoothed_Churchman-4p4.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitlength = 50.\n",
    "calculation_points = np.arange(fitlength + 1.)\n",
    "combined_precision=6.\n",
    "transverse_rpd_s6 = plotting.estimate_rpd_churchman_2d(\n",
    "    input_distances=trans_distances[trans_distances < (fitlength + 3 * combined_precision)],\n",
    "    calculation_points=calculation_points,\n",
    "    combined_precision=combined_precision\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXRV1d3G8e8vIRBCEgJkIJCEMQSQmSgyqFQcEKk4z0pbfXFC0da26msH27etta1aq+KIaFWoRXFAHHBEEdEwCIEEAsgQhiQMgYQAmfb7R6420kBiyOXce/N81sq6ufuek/Pbi8WTk33O2ducc4iISOgK87oAERHxLwW9iEiIU9CLiIQ4Bb2ISIhT0IuIhLgWXhdQl/j4eNe1a1evyxARCRqLFy/e4ZxLqOuzgAz6rl27kpWV5XUZIiJBw8w2Hu4zDd2IiIQ4Bb2ISIhT0IuIhDgFvYhIiFPQi4iEOAW9iEiIU9CLiIS4gLyPvjk7UFHFV5uLKT1YSVl5FfsrqjhQUUVZeRUHK6rpGh/FkLR2pLRrjZl5Xa6IBAEFfQAoKjnIh7mFzMsp4JO8Ig5UVNe7T3x0K4akxTGkSzuGpLVjYGpbWrUIPwbVikiwUdB7pLDkALMW5/PeqgKWbi7GOejUNpKLM1M5pVcCCTGtaB0RTuuW4d++RoSHsXp7CUs37WbJpmKWbtrNu6sKgJrg//HIrlw5rAttoyI87p2IBBILxBWmMjMzXahOgVByoIIn5q/nqU++Zn9FFQNS2jKmdxKn9U2kb3Ls9x6O2Vl6kC837ObFLzYxf00RUS3DueT4VK4Z1Y2UdlF+6oWIBBozW+ycy6zzMwX9sVFeWc0Lizbyjw/WsmtfOeMHJPOzMzLoFt+myY6Rs20vT85fz+tfbcUB4/onc8MpPejbKbbJjiEigUlB76Hqascby7fy13dXs3nXfkb06MAdZ/VmQEqc3465tXg/zyz4mhlfbKasvJJrRnXjttN7EdVSI3UioUpB75H83WXc9q9lfLlhN32SY7njrN6cnB5/zO6W2VNWwb1v5zLji02ktm/Nn84bwKj0+GNybBE5thT0HpizfCt3vrIC5+DX4/ty4dAUwsK8uR3y8/U7ufOVFXy9Yx8XDk3h7rP7EBfV0pNaRMQ/jhT0emCqiZWVV/LLWcuZ/OJSuidEM/eWk7j4+FTPQh7gxO4deGvKSdw4ugezl27htPs/Zs7yrZ7VIyLHloK+CWVv2cP4hz7lpcWbuekHPZh1/XDSOgTGnS+REeH8Ymxv3pg8iuS2rZn84lLueHk5ByqqvC5NRPys3qA3s2lmVmhm2bXaLjKzlWZWbWZ1/qlQa9twM1tqZnOaouBA5Jxj2qdfc96jC9hXXskL1w7j52f2JiI88H6P9u0Uy+wbR3Dj6B7M/HIzF0z9jE07y7wuS0T8qCFJNB0Ye0hbNnA+ML8B+08Bcr5fWcHjYGUVt/97Ob+bs4pTeiXy9pSTGdEjsC94tggP4xdje/P0xEw27yrj7H98wjzfg1ciEnrqDXrn3Hxg1yFtOc651fXta2YpwNnAU42uMIAVlRzk8icX8fKSfG49LZ0nrx5KuzbBc5FzTJ8k3rzlJLp0iOJ/nsviT2/lUFlV//QLIhJc/D228CDwC6De9DCzSWaWZWZZRUVFfi7r6K3aupdzH1nAyq17eOTyIdx6Wq+gnGQstX0Us64fwRXD0nj84/Vc/tQiikoOel2WiDQhvwW9mY0HCp1zixuyvXPuCedcpnMuMyEhwV9lNYm3s7dzwdTPqKp2zLp+BGcPSPa6pKMSGRHOH87rzwOXDGR5fjETHv6UlVv3eF2WiDQRf57RjwTOMbMNwEzgVDN73o/H8zvnHA9/kMf1zy8mo2MMr08eSb/Obb0uq8mcNziFWdePwAEXTl3I29nbvS5JRJqA34LeOXency7FOdcVuBT4wDl3pb+O528HK6v46Utf8dd313DuoE7MnHQiibGRXpfV5Pp1bstrN40ko2MM1z+/mIc/yCMQH6oTkYZryO2VM4CFQIaZ5ZvZNWZ2npnlA8OBN83sHd+2ncxsrn9LPvZ2lh7kiicXMXvpFm4/oxcPXDKIyIjQnfs9MTaSmZNO5NxBnfjru2uYMnOZ7rcXCWKaAqEeawtL+Mn0LAr2HuBvFw9k/IBOXpd0zDjnmPrxOv7yzmoGpMTx5FVDQ/KvGJFQoCkQGumTvCLOe/QzysqrmDnpxGYV8gBmxo2je/L4lUPJKyhhwiMLyN6ii7QiwUZBfxjPf76RHz3zJZ3jWvPqTSMYnNbO65I8c8ZxHZl1/QgMuOixhbyzUhdpRYKJgv4QByuruPvVFdz9ajan9Epg1g0jtFITNVMnvDp5JL18F2mnfrROF2lFgoSCvpYtxfu5+PHPef7zTVx3SneevDqT6FZarOMbiTGR/Ms3hPXnt3O5/d/LOVipi7QigU4p5jN/TRFTZi6lssrx+FVDOfO4jl6XFJAiI8J56NJB9Ehow4Pv5bF5VxmPXTWU9kE09YNIc9Psz+irqx0PvZ/HxGe+IDEmktdvHqWQr4eZcetpvXjossEsyy/mnIc/JXf7Xq/LEpHDaNZBX1xWzrXPZXH/vDWcO6gzs28a0aSLdYe6cwZ24qXrhlNeWc35j36mJ2lFAlSzDHrnHK8syWfM3z7mk7wifj/hOO6/eKAWz26EQalxvHHzKNKTai7S/v29PKqrdZFWJJA0u2RbX1TK3a9m89m6nQxKjeOP5/Wnb6dYr8sKakmxNRdp73plBQ+8t4bc7Xv5m35xigSMZvM/8UBFFVM/WsfUj9bRKiKM/zu3H5efkObpWq6hJDIinL9dPJC+nWL549wcNkwt48mrh+rWVJEAEPJBX1Xt+CC3kD/OzeHrHfs4Z2An7h7fh8QYPcrf1MyMa0/qTnpSDJNfXMI5Dy/gH5cNZmTPwF5xSyTUhexcN9v3HOBfX27mX19uYuueA3TpEMXvJ/Tj5F6BPdd9qFhfVMr1zy9mbWEpt5+ZwQ2n9AjKhVlEgsWR5roJqTP6qmrH/DVFvLBoEx/kFlDt4KT0eH41vi9j+iTRskWzvPbsie4J0cy+cSR3vLKC+95ezbJNxfz14oHERkZ4XZpIsxMyQV96sJKxD84nf/d+4qNbct0pPbjs+DTSOmiM2CttWrXgoUsHMSQtjj+8mcOEhxfw2JVDyegY43VpIs1KyAR9dKsWnD0gmYEpcZyms/eAYWb8eGQ3+nVuy40vLOHcRxZw7wX9mTCos9eliTQbDVl4ZJqZFZpZdq22i8xspZlVm1mdY0JmFmlmX5jZV75t72nKwuty51l9GNc/WSEfgI7v2p43bx5F/85tmTJzGb96NVvz5IgcIw1JxOnA2EPasoHzgflH2O8gcKpzbiAwCBhrZic2pkgJDYmxkbzwP8O47uTu/PPzjVww9TM27SzzuiyRkFdv0Dvn5gO7DmnLcc6trmc/55wr9b2N8H0F3i0+ckxFhIdx57g+PHV1Jpt37efsf3yiqRNE/MyvYxxmFm5my4BCYJ5zbpE/jyfB47S+Scy5eRTdE6K5/vnF3PPGSsorq70uSyQk+TXonXNVzrlBQApwgpn1O9y2ZjbJzLLMLKuoqMifZUmASG0fxb+vG86PR3blmQUbuOjxheTv1lCOSFM7JlctnXPFwEf891h/7W2ecM5lOucyExL0UFNz0bJFGL/54XFMvWII6wtLOfuhT/kgt8DrskRCit+C3swSzCzO931r4DQg11/Hk+B2Vv9k5twyipR2rfnJ9Cz+/HYulVUayhFpCg25vXIGsBDIMLN8M7vGzM4zs3xgOPCmmb3j27aTmc317ZoMfGhmy4EvqRmjn+Ofbkgo6NKhDS/fMILLTkhj6kfruPypRRTuPeB1WSJBL2TnupHgNntpPne9kk2bVuE8dOlgRmhiNJEjOtJcN3qySALSeYNTeH3ySOKiWnLl04t45MO1BOJJiUgwUNBLwEpPiuG1m0Zy9oBO/OWd1Ux+cSll5ZVelyUSdBT0EtC+mRjtzrN681b2Ns5/9DM279ItmCLfh4JeAp6Zcd0pPXjmxyewtXg/5zz8KZ+t2+F1WSJBQ0EvQeOUXgm8NnkUHaJbcdXTX/DMgq81bi/SAAp6CSrd4tsw+8YR/CAjkXveWMUdL6/Q/fYi9VDQS9CJiYzgiauGMvkHPflX1mau++di9pdrymORw1HQS1AKCzNuPzOD35/bjw9WF3LV04soLiv3uiyRgKSgl6B21YldeOTyISzP38PFjy9k2579XpckEnAU9BL0xvVPZvpPjmdb8QEuePQz1haWeF2SSEBR0EtIGNEjnpnXnUh5lePCxxaydNNur0sSCRgKegkZx3Vqyys3jKBt6wiueGoRizcq7EVAQS8hJq1DFP++fjhJsZH86JkvyN6yx+uSRDynoJeQkxgTyfPXDiM2MoKrp31BXoHG7KV5U9BLSOoc15oXrh1GeJhxxVOL2Lhzn9cliXhGQS8hq2t8G164dhgVVdVc/uQithbr1ktpnhqywtQ0Mys0s+xabReZ2UozqzazOie6N7NUM/vQzHJ8205pysJFGqJXUgzP/WQYe/dXcOVTiygqOeh1SSLHXEPO6Kfz34t6ZwPnA/OPsF8l8DPnXB/gROAmM+vbmCJFjkb/lLY88+Pj2bbnAFc9vYg9+yu8LknkmKo36J1z84Fdh7TlOOdW17PfNufcEt/3JUAO0PkoahVptMyu7Xny6kzWFpYyZeZSqqo166U0H8dkjN7MugKDgUVH2GaSmWWZWVZRUdGxKEuamVHp8dwz4Tg+Wl3EX9454nmKSEjxe9CbWTTwMnCrc27v4bZzzj3hnMt0zmUmJCT4uyxppq4Y1oUrhqXx2MfreG3ZFq/LETkm/Br0ZhZBTci/4Jx7xZ/HEmmo3/zwOE7o1p5fzFrO8vxir8sR8Tu/Bb2ZGfA0kOOcu99fxxH5vlq2CGPqFUOIj27FpOcWU1hywOuSRPyqIbdXzgAWAhlmlm9m15jZeWaWDwwH3jSzd3zbdjKzub5dRwJXAaea2TLf1zg/9UPke+kQ3Yonrh7Knv0V3PD8Eg5WauESCV0WiGtuZmZmuqysLK/LkGbgzeXbuOnFJVySmcq9F/Sn5g9RkeBjZoudc3U+16QnY6VZO3tAMjefWrMk4YwvNntdjohfKOil2bvttF6clB7P7+asZF1RqdfliDQ5Bb00e2Fhxl8vGkjriHBunbmM8spqr0sSaVIKehEgKTaSP50/gBVb9vDge2u8LkekSSnoRXzG9uvIJZmpTP14HYvW7/S6HJEmo6AXqeXXP+xLl/ZR/PSlrzT5mYQMBb1ILW1ateCBSwaxfe8Bfv1adv07iAQBBb3IIQantWPKmHReW7ZV8+FISFDQi9ThxtE9GNqlHXe/mk3+7jKvyxE5Kgp6kTq0CA/jwUsG4Rzc8fIKAvEJcpGGUtCLHEZq+yh+eVZvPl27g1c1hCNBTEEvcgRXnJDG4LQ4fj8nh937yr0uR6RRFPQiRxAWZvzp/P7s3V/Bn97K8bockUZR0IvUo3fHWK49qTsvZeXzuR6kkiCkoBdpgClj0klt35q7Zq/Q3PUSdBqy8Mg0Mys0s+xabReZ2UozqzazOuc/Pty+IsGodctw/u/c/qwv2sfUj9Z5XY7I99KQM/rpwNhD2rKB84H5jdhXJCid0iuBcwZ24tEP12k6Ywkq9Qa9c24+sOuQthzn3OrG7CsSzH41vi+REWH872zdWy/BI2DG6M1skpllmVlWUVGR1+WI1CkhphV3juvD5+t3MWtxvtfliDRIwAS9c+4J51ymcy4zISHB63JEDuuSzFSO79qOP72VqxkuJSgETNCLBIuwMOO35xzH7rJyHno/z+tyROqloBdphOM6teXS41N59rMNujArAa8ht1fOABYCGWaWb2bXmNl5ZpYPDAfeNLN3fNt2MrO5R9rXP90QOfZ+dkYGrSPC+cObemJWAluL+jZwzl12mI9m17HtVmBcA/YVCXrx0a24eUxP/jg3l4/XFHFKL11bksCkoRuRo/CjEd3o2iGK389ZRUVVtdfliNRJQS9yFFq2CON/z+7L2sJSXvh8o9fliNRJQS9ylE7rk8ionvE88F6epjKWgKSgFzlKZsbd4/tQcqCCB99b43U5Iv9FQS/SBHp3jOXyYWk8v2gTawpKvC5H5DsU9CJN5KenZ9CmZTi/n7NK8+BIQFHQizSR9m1aMuW0XnySt4P5eTu8LkfkWwp6kSZ01YldSGsfxb1v5VJVrbN6CQwKepEm1LJFGLefmUHOtr28unSL1+WIAAp6kSY3vn8y/Tu35f55azhQoWUHxXsKepEmFhZm3DmuN1uK9/Pcwg1elyOioBfxhxE94hmdkcDDH6yluEwPUYm3FPQifvLLsb0pOVjJo1pMXDymoBfxkz7JsVwwJIXpn20gf3eZ1+VIM6agF/Gjn57eCwPuf1dTI4h3FPQiftQprjU/HtmN2cu2sHLrHq/LkWaqIStMTTOzQjPLrtV2kZmtNLNqM8s8wr5jzWy1ma01szuaqmiRYHLD6B60bR3BvW/lel2KNFMNOaOfDow9pC0bOB+Yf7idzCwceAQ4C+gLXGZmfRtXpkjwats6gsk/6MkneTv4VFMjiAfqDXrn3Hxg1yFtOc651fXsegKw1jm33jlXDswEJjS6UpEgdtXwLnSOa82f386lWlMjyDHmzzH6zsDmWu/zfW11MrNJZpZlZllFRUV+LEvk2GvVIpzbTu/Fii17mJu9zetypJnxZ9BbHW2HPZVxzj3hnMt0zmUmJGiRZQk95w3uTEZSDH99Z7XWl5Vjyp9Bnw+k1nqfAmz14/FEAlp4mPHzMzPYsLOMf325uf4dRJqIP4P+SyDdzLqZWUvgUuB1Px5PJOCN6ZNIZpd2/P39PMrKK70uR5qJhtxeOQNYCGSYWb6ZXWNm55lZPjAceNPM3vFt28nM5gI45yqBycA7QA7wknNupb86IhIMzIw7zupNUclBnlmwwetypJmwQFzyLDMz02VlZXldhojfXPvslyxav4v5v/gB7dq09LocCQFmttg5V+dzTXoyVsQDPz+zN6XllTz60VqvS5FmQEEv4oGMjjGcPziFZxduZEvxfq/LkRCnoBfxyG2np4ODB+dpwjPxLwW9iEdS2kVx1fAuvLwknzUFJV6XIyFMQS/ioZt+0JM2LVtw39v1zSgi0ngKehEPtW/TkutH9+C9nAIWrd/pdTkSohT0Ih77ychudIyN5I9v5RKItztL8FPQi3isdctwfnpGL77aXMzcFdu9LkdCkIJeJABcMCSFjKQY7nsnl/JKTXgmTUtBLxIAwsOMO8b1ZuPOMl5ctNHrciTEKOhFAsToXgmM6NGBhz5Yy94DFV6XIyFEQS8SIMyMO8/qw6595Tz+8Tqvy5EQoqAXCSD9U9oyYVAnnvrka7bt0dQI0jQU9CIB5vYzMnAOHtDUCNJEFPQiASa1fRRXD+/CrMX55G7f63U5EgIasvDINDMrNLPsWm3tzWyemeX5XtsdZt8pZpZtZivN7NamLFwklE0+tSfRrVrwhzdz9BCVHLWGnNFPB8Ye0nYH8L5zLh143/f+O8ysH/A/wAnAQGC8maUfVbUizURcVEtuGZPOJ3k7+CC30OtyJMjVG/TOufnArkOaJwDP+r5/Fji3jl37AJ8758p8ywp+DJx3FLWKNCsTR3SlR0Ibfj9nFQcrq7wuR4JYY8fok5xz2wB8r4l1bJMNnGxmHcwsChgHpB7uB5rZJDPLMrOsoqKiRpYlEjoiwsP41fi+bNhZxnStLytHwW8XY51zOcCfgXnA28BXwGGXvXfOPeGcy3TOZSYkJPirLJGgMjojkTG9E/nHB2spLDngdTkSpBob9AVmlgzge61zENE597Rzbohz7mRqhn/yGnk8kWbr7vF9OVhZpTnrpdEaG/SvAxN9308EXqtrIzNL9L2mAecDMxp5PJFmq1t8G34yqhuzFufz1eZir8uRINSQ2ytnAAuBDDPLN7NrgHuB080sDzjd9x4z62Rmc2vt/rKZrQLeAG5yzu1u8h6INAM3n5pOQkwrfvvGSqqrdbulfD8t6tvAOXfZYT4aU8e2W6m56PrN+5MaX5qIfCO6VQt+cWYGP5+1nFeXbeH8ISlelyRBRE/GigSJC4akMDA1jnvfymXfwcPe1yDyXxT0IkEiLMz47Q/7UlhykEc+XOt1ORJEFPQiQWRwWjsuGJLCk5+sJ6+gxOtyJEgo6EWCzF3jehPdqgW/fHk5VbowKw2goBcJMh2iW/HrH/ZlyaZinv9cyw5K/RT0IkHo3EGdOblXAve9ncuWYi1QIkemoBcJQmbGH87thwPunr1CUxnLESnoRYJUavsobj8jgw9XF/H6V1u9LkcCmIJeJIhNHNGVgalx3PPGKnbtK/e6HAlQCnqRIBYeZvz5gv7s3V/B/81Z5XU5EqAU9CJBrnfHWG4Y3YNXlm7h4zVay0H+m4JeJARMPrUnPRLacNcrKyg5UOF1ORJgFPQiIaBVi3Duu3AA2/bs539nZ+suHPkOBb1IiBjapT0/Pb0Xr3+1lZeyNntdjgQQBb1ICLlhdE9G9uzAb15fyRrNhSM+DVl4ZJqZFZpZdq229mY2z8zyfK/tDrPvbWa20syyzWyGmUU2ZfEi8l3hYcYDlwwiulULbnphCfvLq7wuSQJAQ87opwNjD2m7A3jfOZcOvO97/x1m1hm4Bch0zvUDwoFLj6paEalXYkwk9188iLzCUu55Y6XX5UgAqDfonXPzqVnYu7YJwLO+758Fzj3M7i2A1mbWAogC9PieyDFwcq8Ebhjdg5lfbtZTs9LoMfok59w2AN9r4qEbOOe2AH8FNgHbgD3OuXcbW6iIfD8/Pb0XQ9LiuOuVFWzYsc/rcsRDfrsY6xu3nwB0AzoBbczsyiNsP8nMsswsq6hID32IHK2I8DAeumww4WHGzTOWcrBS4/XNVWODvsDMkgF8r4V1bHMa8LVzrsg5VwG8Aow43A90zj3hnMt0zmUmJCQ0siwRqS2lXRR/uXAAK7bs4ZezllOthUqapcYG/evARN/3E4HX6thmE3CimUWZmQFjgJxGHk9EGumM4zry8zMzeHXZVu57Z7XX5YgHGnJ75QxgIZBhZvlmdg1wL3C6meUBp/veY2adzGwugHNuETALWAKs8B3rCb/0QkSO6MbRPbhiWBqPfbyO5xZu8LocOcYsEB+VzszMdFlZWV6XIRJSqqod1/1zMe/nFvDYlUM587iOXpckTcjMFjvnMuv6TE/GijQT4WHGPy4bzMCUOG6ZsZTFGw+9a1pClYJepBlp3TKcpydmktw2kmuezWJdUanXJckxoKAXaWY6RLfi2Z+cQIswY+K0Lyjce8DrksTPFPQizVCXDm14euLx7NpXzsWPL2TTzjKvSxI/UtCLNFMDU+P45zXDKN5fwflTPyN7yx6vSxI/UdCLNGNDu7Rj1vXDaRluXPrE53y2dofXJYkfKOhFmrmeiTG8cuNIOse1ZuIzX/CGJkELOQp6EaFj20heum44g1PbccvMpTyz4GuvS5ImpKAXEQDaRkXw3DUncEbfJO55YxW/e2OVJkILEQp6EflWZEQ4j14xlB+N6Mq0BV8z4eEF5Gzb63VZcpQU9CLyHeFhxm/POY5pP8pkR2k5Ex5ewGMfr6NKM18GLQW9iNTp1N5JvHvbyZzaO5F738rlsic+Z/Mu3W8fjBT0InJY7du0ZOqVQ/jbRQPJ2baXsQ/O58VFm3R2H2QU9CJyRGbGBUNTeOvWk+if0pa7Zq/g9Ac+5o2vtmohkyChoBeRBklpF8WL157IY1cOoYVvecJxD33Cuyu3E4jTnct/KOhFpMHCwoyx/ZJ5a8rJ/P3SQRysrGbSPxcz4ZEFfJhbqDP8ANWQFaammVmhmWXXamtvZvPMLM/32q6O/TLMbFmtr71mdmtTd0BEjr3wMGPCoM7Mu+1k/nLhAHbtK+fH07/kpPs+5P53V7Nx5z6vS5Ra6l1hysxOBkqB55xz/Xxt9wG7nHP3mtkdQDvn3C+P8DPCgS3AMOfcxvqK0gpTIsGlvLKat7K3MWtxPp+u3YFzcHzXdlw4NIVx/ZOJiYzwusSQd6QVphq0lKCZdQXm1Ar61cBo59w2M0sGPnLOZRxh/zOA3zjnRjakYAW9SPDatmc/s5duYdbifNYX7SMyIozh3TtwUnoCJ6XH0zMxGjPzusyQ44+gL3bOxdX6fLdz7r+Gb2p9Pg1Y4px7+AjbTAImAaSlpQ3duLHeE38RCWDOOZZtLua1ZVuZv6aI9TtqhnM6xkZyUno8o9LjGdatAx3bRnpcaWjwNOjNrCWwFTjOOVfQkIJ1Ri8SevJ3l/Fp3g4+ydvBp2t3sGd/BQCJMa0YmBrHwJS2DEiJY2BKHG2jNNTzfR0p6Fs08mcWmFlyraGbwiNsexY1Z/MNCnkRCU0p7aK49IQ0Lj0hjapqR/aWPSzZtJvl+Xv4Kr+Yeav+ExGp7VuTkRRLRsdoMjrG0rtjDN3i2xARrhsFG6OxQf86MBG41/f62hG2vQyY0cjjiEgICg+zmrP41G8HBtizv4LsLXtYtrmYnG17Wb29hA9XF377FG5EuNE9PpqeidH0SKx5TU+Mplt8GyIjwr3qSlBoyF03M4DRQDxQAPwGeBV4CUgDNgEXOed2mVkn4Cnn3DjfvlHAZqC7c67B65Rp6EZEAA5WVrG+aB9rCkrI3V5CXkEJawtL2bSrjG9u2Q8zSG0fRc8E3y+AWq/NaQjoqMfojzUFvYgcyYGKKr7esY+1haU1X0WlrCssZf2OfZRXVn+7XXx0K9K/OftPqgn/nknRJES3Crk7f/wxRi8i4pnIiHD6JMfSJzn2O+1V1Y783WX/+QXg+yXw6rItlByo/Ha72MgWpCfF0CspmvTEGNKToumVFENiTOj9AgAFvYiEkPAwo0uHNnTp0IYxfZK+bXfOUVhykLyCUtYWlpBXWEpeYSlvZ29nRtnmb7f7zy+AGDJ84d+rYwzx0a286E6TUdCLSMgzM5JiI0mKjWRUevy37c45dpSWk1dYQl5BKWsKal7fyt7GjC8qvt2ufdCeAG8AAAX/SURBVJuW9EqKprfvDqDeybH0SoomqmVwRGhwVCki4gdmRkJMKxJiWjGix3d/ARSVHGRNQSmrC2ouAuduL+GlrM2UlVf59oUu7aPI6BhDv05t6de55ishJvDO/hX0IiKHMDMSYyNJPOQvgOpqR/7u/eRs30vuthJWF+xl1da9vLPyP88AdIyN9IV+LEO7tGNol3aen/kr6EVEGigszEjrEEVahyjOPK7jt+0lBypYuXUv2Vv2kL1lDyu27OH93AKcgxZhRr/ObRnWvT3DurUns2t7Yo/xJG+6vVJExA9KD1ayeONuvvh6J4vW7+Kr/GIqqhxmMDg1jrMHdGJc/44kt23dJMfTffQiIh47UFHF0k3FfL5+J/NWFbBq214AMru04+wByZzVL/moJnhT0IuIBJj1RaXMXbGNOcu3kbu9BDM4oWt7nr92WKPm9NEDUyIiAaZ7QjSTT01n8qnprC2sCf2txfv9MnGbgl5ExGM9E6O5ZUy6336+5vwUEQlxCnoRkRCnoBcRCXEKehGREFdv0JvZNDMrNLPsWm3tzWyemeX5Xg+3Xmycmc0ys1wzyzGz4U1ZvIiI1K8hZ/TTgbGHtN0BvO+cSwfe972vy9+Bt51zvYGBQE4j6xQRkUaqN+idc/OBXYc0TwCe9X3/LHDuofuZWSxwMvC07+eUO+eKj6paERH53ho7Rp/knNsG4HtNrGOb7kAR8IyZLTWzp8yszeF+oJlNMrMsM8sqKipqZFkiInKoBk2BYGZdgTnOuX6+98XOubhan+92zrU7ZJ9M4HNgpHNukZn9HdjrnPtVA45XBGysZ7N4YEe9xYce9bt5Ub+bl6PpdxfnXEJdHzT2ydgCM0t2zm0zs2SgsI5t8oF859wi3/tZHH4s/zsOV2xtZpZ1uHkdQpn63byo382Lv/rd2KGb14GJvu8nAq8duoFzbjuw2cwyfE1jgFWNPJ6IiDRSQ26vnAEsBDLMLN/MrgHuBU43szzgdN97zKyTmc2ttfvNwAtmthwYBPyxqTsgIiJHVu/QjXPussN8NKaObbcC42q9Xwb468+vJ/z0cwOd+t28qN/Ni1/6HZDz0YuISNPRFAgiIiFOQS8iEuKCLujNbKyZrTaztWbWoNs1g9HRzDEUzMws1cw+9M2NtNLMpvjaQ7rvZhZpZl+Y2Ve+ft/jaw/pfn/DzMJ9D1bO8b1vLv3eYGYrzGyZmWX52pq870EV9GYWDjwCnAX0BS4zs77eVuU302n8HEPBrBL4mXOuD3AicJPv3zjU+34QONU5N5CaO9TGmtmJhH6/vzGF786F1Vz6DfAD59ygWvfPN3nfgyrogROAtc659c65cmAmNfPuhJzGzjEU7Jxz25xzS3zfl1Dzn78zId53V6PU9zbC9+UI8X4DmFkKcDbwVK3mkO/3ETR534Mt6DsDm2u9z/e1NRcNmWMoZPim3hgMLKIZ9N03fLGMmifN5/meKg/5fgMPAr8Aqmu1NYd+Q80v83fNbLGZTfK1NXnfg21xcKujTfeHhiAziwZeBm51zu01q+ufPrQ456qAQWYWB8w2s35e1+RvZjYeKHTOLTaz0V7X44GRzrmtZpYIzDOzXH8cJNjO6POB1FrvU4CtHtXihQLf3EIcYY6hoGdmEdSE/AvOuVd8zc2i7wC+6bw/ouYaTaj3eyRwjpltoGYo9lQze57Q7zfw7UOmOOcKgdnUDE83ed+DLei/BNLNrJuZtQQupWbeneai3jmGgp3VnLo/DeQ45+6v9VFI993MEnxn8phZa+A0IJcQ77dz7k7nXIpzris1/58/cM5dSYj3G8DM2phZzDffA2cA2fih70H3ZKyZjaNmTC8cmOac+4PHJfmFb46h0dRMW1oA/AZ4FXgJSAM2ARc55w69YBvUzGwU8Amwgv+M2d5FzTh9yPbdzAZQc+EtnJoTsJecc78zsw6EcL9r8w3d3O6cG98c+m1m3ak5i4eaYfQXnXN/8Effgy7oRUTk+wm2oRsREfmeFPQiIiFOQS8iEuIU9CIiIU5BLyIS4hT0IiIhTkEvIhLi/h8hFUQmCwebbQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "normalised_transverse_rpd_s6 = transverse_rpd_s6[calculation_points > 0.] / calculation_points[calculation_points > 0.]\n",
    "norm_rpd_calculation_points = calculation_points[calculation_points > 0.]\n",
    "plt.plot(norm_rpd_calculation_points, normalised_transverse_rpd_s6)\n",
    "np.save('normalised_transverse_rpd_smoothed_Churchman-6', normalised_transverse_rpd_s6)\n",
    "# normalised_transverse_rpd = np.load('normalised_transverse_rpd_smoothed_Churchman-4p4.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(params_optimised,\n",
    " params_covar,\n",
    " params_1sd_error) = zdisk_modelling.fitmodel_to_hist(\n",
    "    norm_rpd_calculation_points[0:31],\n",
    "    normalised_transverse_rpd_s6[0:31],\n",
    "    trans_model_with_info.model_rpd,\n",
    "    trans_model_with_info.initial_params,\n",
    "    trans_model_with_info.param_bounds\n",
    "    )\n",
    "print('')\n",
    "print('Initial parameter guesses:')\n",
    "print(trans_model_with_info.initial_params)\n",
    "print('')\n",
    "print('Parameter bounds:')\n",
    "print(trans_model_with_info.param_bounds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(norm_rpd_calculation_points,\n",
    "    normalised_transverse_rpd_s8)\n",
    "zdisk_plots.plot_fitted_model(\n",
    "    norm_rpd_calculation_points,\n",
    "    fitlength,\n",
    "    params_optimised,\n",
    "    params_covar,\n",
    "    trans_model_with_info,\n",
    "    plot_95ci=False\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up another model and fit\n",
    "Tried a few smoothing widths again\n",
    "### This model fits much better"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "trans_model_with_info = zdisk_modelling.set_up_model_2d_twopeaks_flat_bg_with_fit_settings()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "normalised_transverse_rpd_s6 = np.load('./normalised_transverse_rpd_smoothed_Churchman-6.npy')\n",
    "calculation_points = np.arange(fitlength + 1.)\n",
    "norm_rpd_calculation_points = calculation_points[calculation_points > 0.]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.10972642e+01 8.35992331e-03]\n",
      " [2.38092488e+01 2.03679789e-02]\n",
      " [6.59552871e+00 7.85793535e-03]\n",
      " [1.20984124e+02 4.88657166e-01]\n",
      " [1.84130391e+02 8.07139153e-01]\n",
      " [1.06337374e+01 1.97937077e-03]]\n",
      "SSR = 3.1998682683694416e-06\n",
      "AIC = -484.67803170839477\n",
      "AICcorr = -479.80846649100346\n",
      "\n",
      "Initial parameter guesses:\n",
      "[20, 40.0, 5, 1, 1, 20]\n",
      "\n",
      "Parameter bounds:\n",
      "([0, 0, 0, 0, 0, 0], [50, 50.0, 20, 1000, 1000, 100])\n"
     ]
    }
   ],
   "source": [
    "(params_optimised,\n",
    " params_covar,\n",
    " params_1sd_error) = zdisk_modelling.fitmodel_to_hist(\n",
    "    norm_rpd_calculation_points[0:31],\n",
    "    normalised_transverse_rpd_s6[0:31],\n",
    "    trans_model_with_info.model_rpd,\n",
    "    trans_model_with_info.initial_params,\n",
    "    trans_model_with_info.param_bounds\n",
    "    )\n",
    "print('')\n",
    "print('Initial parameter guesses:')\n",
    "print(trans_model_with_info.initial_params)\n",
    "print('')\n",
    "print('Parameter bounds:')\n",
    "print(trans_model_with_info.param_bounds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x at zero bg is -0.05775112568529942.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Code\\perpl-python3\\models_2d_distances_normalised.py:101: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  / x_values # Normalisation\n",
      "C:\\Code\\perpl-python3\\models_2d_distances_normalised.py:104: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  / x_values # Normalisation\n",
      "C:\\Users\\fbsacu\\Miniconda3\\envs\\PERPL\\lib\\site-packages\\numdifftools\\extrapolation.py:489: RuntimeWarning: invalid value encountered in less_equal\n",
      "  converged = err <= tol\n",
      "C:\\Users\\fbsacu\\Miniconda3\\envs\\PERPL\\lib\\site-packages\\numdifftools\\limits.py:173: RuntimeWarning: invalid value encountered in less\n",
      "  outliers = (((abs(der) < (a_median / trim_fact)) +\n",
      "C:\\Users\\fbsacu\\Miniconda3\\envs\\PERPL\\lib\\site-packages\\numdifftools\\limits.py:174: RuntimeWarning: invalid value encountered in greater\n",
      "  (abs(der) > (a_median * trim_fact))) * (a_median > 1e-8) +\n",
      "C:\\Users\\fbsacu\\Miniconda3\\envs\\PERPL\\lib\\site-packages\\numdifftools\\limits.py:175: RuntimeWarning: invalid value encountered in less\n",
      "  ((der < p25 - 1.5 * iqr) + (p75 + 1.5 * iqr < der)))\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZxO9f//8cdrNvturFkiWwnVZCuJyh5RKn1EUiJkK7SQpUVZEhGK+ERIiOyU0kIZ+yCRKFLIkm3M9v79YT7fn4+PZRpzzbmua573221uM3PmOs7z3TFPp/d1FnPOISIiwSvE6wAiIuJbKnoRkSCnohcRCXIqehGRIKeiFxEJcmFeB7iY/Pnzu5IlS3odQ0QkYKxbt+6wcy7yYj/zy6IvWbIk0dHRXscQEQkYZrb3Uj/T1I2ISJBT0YuIBDkVvYhIkFPRi4gEORW9iEiQU9GLiAQ5Fb2ISJDzy/PoxTun4xL48++zxMYncjYhidj4xOSPJM4mJHI2PonY+ATOnoklS+ZMVLgmDxUK5yBrhP4qifgr/XZmcLHxiWzYe5TNX6/n5PIV5N+0lvwnDnOxpxSEAlmTPxJCwohIjOeAGavzF+dI6fKEVrmZ/LfeRPmSkdxQJCf5s2dK38GIyEWp6DOYhMQkNu87xqZvN3F86Qpyr19NkSO/E56nCAmVq3Hq2RfIUak8mSPCyBweQqaw0PM+h/7XsmOn49n+619k+34DtjaaTAs/ocB7b3IqMYG5uQtzoHhZkipVJvK2ajStWZZiebN6PXyRDMn88QlTUVFRTrdASDun4xL4bPkGDs75jBzRqyl8eB+/5y7IgRtuJetddbmx1k3cWiofOTOHX/W2jp+OZ+v+o+xZu4UTP0QTsmUzhXdu5kCuguxt+AC3tmpC/YqFyRwemgYjE5H/MLN1zrmoi/5MRR+8Dp88y6fTlhL+7mhyHv+LLTffSZY6d3J9rSiql85HvnSaWvn92BmWf7ychMkfcM3eH/n+htvI1PpR7r37Jq4vkjNdMogEOxV9BvPLoZMsHjOdAtPe51R4Zn55qD33PnYvt5TI42mupCTH6m372fzuFAotmcvpsEzE1GnGDY8+QNNbiqXJ/1GIZFQq+gxiw88H+XbYeEotmMmuwtdx8vGneOj+WpSOzO51tP9x9FQcSxeu5vh773Pt1rVsKHUTsQ8/StuHa1Myfzav44kEHBV9EEtKcny1fjfbhrxNme+WsbZCDbJ36kir+pUpkCOz1/GuyDnHlr1/sea9GUTOnMKBnAU40a03T9xfPd2mlkSCgYo+SK3evIfNvfpRYvs6vq7emOuebk/L264je6bAPJnq4IlYZg+bQolJo9lSqgq5nutFm3o36hx9kRS4XNHrytgAdPDvWN7pOYyTDeqzO38JTn+6hAEz36Bd3fIBW/IABXJkptPApyi79juyli9HqQ6tGPJAT2Z+9zMJiUlexxMJWCr6AJKY5Jg15xu+vrUObtUqto+ZxsB/D6RFVHHCQ4NnV15XKCddx/cj78qVFAhLIk/LxrzU/hVWbP0Df/w/UBF/FzztEORidh9kQtMOZOneidUtO9Jw0Uc80zwqqM9HjypXmM6zR8OsTym5dzunG9en73Pj2PDrUa+jiQQUFb2fOxEbz4TBE9l/R232hWQlYc5Chg5uw3UF/O9MGl8wM+rVLE/7FdOIfXscN3wxl233NOHV8Uv5Ozbe63giAUFvxvop5xxLvtjE0Z69OJ1kHOzVj84P1SRXlox9rvnJswlMHTeXEsMG8P1Ndan1+gvcdUNhr2OJeE5vxgaY/UdOMbZ1b0LatCK69n3cuugTXnisdoYveYDsmcLo2K0lRb76kryJsZxs3JABw+fy18mzXkcT8VtXLHozm2RmB80s5rxlLc1sq5klmdlF/wU577WhZrbBzBakReBgt2jZOtZWq8Pf+//gj2mfMnRkFyoXy+11LL9TuVQkHeeN56/+Q7hl9EAm3NeJz6L36s1akYtIyRH9ZKDBBctigBbAqhSs3w3Y/s9iZTzHz8TzTo9hhLX7F6tadOBfs8fT9s5yhIaY19H8VkRYCI8/3oCy331Fpjw5CW3emP6DpvLH8Vivo4n4lSsWvXNuFXDkgmXbnXM7rrSumV0DNAbeT3XCDGDN5r3Mr9WUsJUr2DV+Bm++9jjF8+mWvilVrkguus18iyPD36XK1NF81LgtM1f9pKN7kWS+nqMfCfQGrni1i5l1MLNoM4s+dOiQj2P5h7iEJN4bMYO/G9Zn0w01qf7ZDJ5uUoWwIDonPr2EhhitH6zFLd99QeJ1Zcndqhkv9R7Pb0dOex1NxHM+axQzawIcdM6tS8nrnXMTnHNRzrmoyMhIX8XyGzv3H2VikyfIPXYE3w8ay6D3nqeK5uKvWsnI7PSaNIjj731IxYUf8Wmj1ny6ZpeO7iVD8+Wh421AUzPbA8wA6prZVB9uLyA45/j446/YcXsdDoRkIfen8+jXvq7u55KGQkKMBxtFcfvXizld7noyPdCMQW98zPHTOu9eMiafFb1z7nnn3DXOuZLAw8AXzrnWvtpeIPjz7zOMat+f7L2e5uv2fekyczj3VNQ54L5SLF82nv1gMH8MHUOV8UMYe38Xvtt50OtYIukuJadXTgdWA+XMbJ+ZtTez5ma2D6gBLDSzpcmvLWJmi3wbOTB99cNOvritEW7bVo5MncOQFx8KiNsIB7rQEKNdqzu59ssVZLIk/qjfkLenfMHZhESvo4mkG10Z62Ox8Yl8MGI6pd9+haVNn6DjwKcoUzCH17EypNNxCYwfPYcKbw1gTZ37eWRYb8pqX0iQ0JWxHtn5+zHea/YU+d8bxcYh43l1dA+VvIeyRoTRo9eDhM9fROGdW1hzV1OmLt6oN2ol6KnofcA5x6z5a9hW6y6OJoWSf/58erepHdR3mgwkd91ckhbLZ/JbveYUbP8Qg54by8G/dZGVBC8VfRo7ciqO4d2Hk6vTY6x55Bk6znqbOtcX8jqWXCAyRyZeGN6Vw5NmUGbRDGY2fJTPN/7qdSwRn1DRp6HVMb8x964HyLNyKfsnzeLVgW30hqsfMzNaNbiJqisXcrx4ac7e24jho+dxJk5v1EpwUdGngfjEJMaNm8+xBg3YXboS1RZ/Qrv6NxKi+9QEhOsK5uS5D19lZ7+hlB3Wj5EP92Lb/uNexxJJMyr6q7R131FGPdCdoq/3Y+MLQ3nxgwFULJrL61jyD2UKC6Vbh4bkWrSEzEcPE1O3Af+e/wNJSXqjVgKfij6VziYkMuGDZeyqVZe4U6eImL+Q559upCtcA9wdNxSh7YJJ/PjA4xTq2IbXe4zk4Am9USuBTefRp8KmvUdY1q0/FdZ/xaZnBtD56XvJnTXC61iShpxzzPg8hrM9uxOXJRvXjRtF3ZtKeh1L5JJ0Hn0aiY1P5J2JS9hz510knYkl26KlvPjs/Sr5IGRmtLr7Rm5f8Sn7y1Uh7t5GjHxrNrHxeqNWAo+O6FNo3S+HWdH9Za7f+DXbegykU8cm5MysR/tlBGcTEhn775VcO7g3P1esRv0xr1KxeF6vY4n8Fx3RX4UzcYmMHL+YfXXuxuLiyL1kKX26N1fJZyCZwkLp8fjd5F20iLj4BH6ufTcffLSShMQrPmZBxC+o6C/BOceqH//k7fufofRrffi510Ce/vRdalXQ3SYzqjsqFKLTnDFsfKoPBft0YUTbF9l7+KTXsUSuSEV/ERv2HmHQMyM4eU9dwhLjiFy+nJ5dm5E9k86oyehyZ43g5b4PkTD7M8L37GbN7fX5ZMl63S9H/JqK/jw//XmCgf0msffOuym49huOjJ1M1/njqF62oNfRxM80rVqKhxZMZnPTNuRs/whDuw3TaZjit/RmLPDbkdNMmbyEa8cPJyFzFuKfe4GHW96hI3i5oqQkx4crYjjbpzeZLInCY9+hXvUyXseSDEhvxl7CwROxvDlxBYvrNqfcuDc53qUXTb9bxBOt66rkJUVCQoy29W6kzpKP+TGqLrS4l5Evv8+JWD22UPxHhjyiP346nsmL1pE0YhjX7d/JnvbdaNntEQrl0g3IJPXiEpIYP/t78r78HKfyF+SGsW9xW6XiXseSDEJH9Jx7utCCzb/T77UZjL+nNeW7tSOsanUqRn9L1/6Pq+TlqkWEhdD1oRqUXzKfPSUqcKphPUa/MllH9+K5oD6ij41P5Msdh/hm+fdkmj+HqJ3RHC5QlNimLbjt8Qe5Xhe9iI/Exicy9uPvKPDKC8TmyU+5sW9Rq0pJr2NJELvcEX3QFX1cQhLf7DrElys3wdzZVP3xO07mzMvJxvdxQ9sHqVq+CKG6fbCkk3V7jzLrpbep+/kMfnzqWR7r01YX24lPZIiiPx2XwBvTvuPs7E+ovvUbEiMycaxeE8q2e4QalUsQFpphZqnEz8TGJzJ29hryD3qBhBy5KD12JLVvKeV1LAkyGWKOPkt4KOErllCuUE5yT/+Ie7d+S/vRL1Lr5mtV8uKpzOGh9Hz4Nm5c8ik/VYji7L2NGN1vHMfPaO5e0kfQHNHDudsWmGlaRvxXbHwi7879gVyDXiQkSxaufXcUtaNKex1LgkCGOKIHVPLi9zKHh9LjwRrcsngO2268jbNNGzPquZH8dfKs19EkiF2x6M1skpkdNLOY85a1NLOtZpZkZhef/DfLbGY/mNmm5NcOTMvgIoGscrHcDJ7Qhx1jphL++TKW1mzAgs836Z454hMpOaKfDDS4YFkM0AJYdZn1zgJ1nXOVgSpAAzOrnpqQIsEoU1goXZtHUXfRTDbUbUFEm4cY1b4/+46c8jqaBJkrFr1zbhVw5IJl251zO66wnnPO/eceruHJHzpcEblAuUI5GDKiC398+ClxP+1kbbU6zJi1ikQ9mFzSiE/n6M0s1Mw2AgeB5c657325PZFAFRpitKlbnlafTeL71s+QvVdnxjzQlR37j3kdTYKAT4veOZfonKsCXANUNbOKl3qtmXUws2gziz506JAvY4n4rWvyZOX1/v8icc4Cjp2OZ1fNO5g4bj5nE/SsWkm9dDnrxjl3DPiS/53rP/81E5xzUc65qMjIyPSIJeKXzIxmUSXoMns0a3u9Rq6hg5hUvy3RO373OpoEKJ8VvZlFmlnu5K+zAHcDP/pqeyLBJm+2CAY804T8CxfyW54i/FX3Lsa+/iEnzyZ4HU0CTEpOr5wOrAbKmdk+M2tvZs3NbB9QA1hoZkuTX1vEzBYlr1oYWGlmm4G1nJujX+CbYYgErzvLF+SFaa+x8bVx5JgygRl3PsBX63Z7HUsCSFBdGSsS7NbtPcqs/qOou3w629p2pc1LT5I3W4TXscQPZJgrY0WC3S0l8jBw4kv8NGoK2RbNY+FtjVn4VYwutJLLUtGLBJhMYaF0eaAatZZ+QkzNhoS3eoARXYbw+7HTXkcTP6WiFwlQ5Qvl5LV3erD/g1lErF3D1zXq8fGCtSTpQiu5gIpeJICFhhjt6t9Is8XT2dSsDbmeasPwDoP4/aiO7uX/U9GLBIHi+bLy6utPcvSjeWTZvIGVtRoxb+UWzd0LoKIXCRpmxsO1y3Lv4mlsvas5Wf7VkmG939EtkEVFLxJsSuTLxuC3nmHfhI8osHgOs+9qyYp1v3gdSzykohcJQqEhxuNNbqbainn8Vv5mEps2YvjrUzkRq8cXZkQqepEgVr5QTl6a2I+f3hhH8YmjmNSkPWt2/OF1LElnKnqRIJcpLJSurWtTasViTuTIw+F69Rg1biGx8bojZkahohfJIG4pmY+eM4fzY+9XKPVaX4Y91o9f/9JpmBmBil4kA8kaEcaznZuSZf5C8v+0haX1W/LFpl+9jiU+pqIXyYDuqlKcRktmcqDCTZy+tzFjPvxCjy4MYip6kQyqeL6s9J40kG3d+lOm79MM7vk2h3XOfVBS0YtkYJnDQ+ndqyUnpn5C+SWz+OC+jkTvPux1LEljKnoR4f46FblxxQJCwkL5ucG9TFm8UbdPCCIqehEB4IZr8vDE3Hf5+b7WFG7/MINen67HFgYJFb2I/J9cWcLpO+Rp9o98j1vGD2FomxfZ+ecJr2PJVVLRi8h/CQkx2j1Yi3yLl1B8dwyLm7fjh581bx/IVPQiclE1ri9CvSUfk5g7D9vue5DF6/d6HUlSSUUvIpdULF82Hps1hj+r1OBUyweYujzG60iSCip6EbmsPNkieGbiIH5u1po87Vvx9vRvdEZOgFHRi8gVZYkIpdfQruzs3JcKvZ5g8KjPiEtI8jqWpJCKXkRSJCw0hG69H2H/66Oo/kZv+g/8UKdfBggVvYikmJnRrm09Tk2cyh2ThzKox9scPBHrdSy5gisWvZlNMrODZhZz3rKWZrbVzJLMLOoS6xUzs5Vmtj35td3SMriIeKd5wyiyzZ1HzWXTeevJQew+dNLrSHIZKTminww0uGBZDNACWHWZ9RKAXs65CkB1oLOZXZ+akCLif2pHlabkkoVU2rqaaa17suHXo15Hkku4YtE751YBRy5Ytt05t+MK6x1wzq1P/voEsB0oehVZRcTPVLmuIDU+n0+Rk4dZ0bozm3475nUkuYh0maM3s5LATcD36bE9EUk/JQvkoOHCqeRJjGVBux5sP/C315HkAj4vejPLDswGujvnLvk3wMw6mFm0mUUfOnTI17FEJA0VyZ2V+vP+TbFjfzCzfV92HdScvT/xadGbWTjnSn6ac27O5V7rnJvgnItyzkVFRkb6MpaI+ECx/Nm5fcEMyu7fwQdP9WfvX6e8jiTJfFb0ZmbARGC7c26Er7YjIv6jVKFc3LJgFpV3rGXs06/y+7EzXkcSUnZ65XRgNVDOzPaZWXsza25m+4AawEIzW5r82iJmtih51duAR4G6ZrYx+aORj8YhIn6iXLG8VPhsNtU2fM6IbkN1nr0fMH+8Z0VUVJSLjo72OoaIXIV1P+5jf+NmfFX3IV4c2Z282SK8jhTUzGydc+6i1zXpylgR8Ylbyl9D5JzZ1Fk+jcF9x3H8TLzXkTIsFb2I+EyNyiXJMesT6s2bwMCXJuneOB5R0YuIT9W+tQxhU2dQf+bbvDx4KmfiEr2OlOGo6EXE5+65oyKJH0yl/gdDePmNWSQk6hbH6UlFLyLpotE9N3NizCTufncAw6Z943WcDEVFLyLp5v77arK758tU6P8MH339k9dxMgwVvYikqyd7Pszuei040/lpvtul252kBxW9iKSr0BCj/ajniS9QiC86Pc+ew7pVgq+p6EUk3eXMHE7D6e9Sdv8OxvQeyd+xOsfel1T0IuKJEpE5KD5jKnVXTOeVN3Umji+p6EXEM9UrleDMO+9x9/hBjPjoW6/jBC0VvYh4qkXTGuzq3p9y/boy85udXscJSip6EfFch54P8/M9zTn5dCfW/HzY6zhBR0UvIp4LCw2h/agXiI8syPJOz/PrX6e9jhRUVPQi4hdyZQmn/vRxlNu3ndG93+aEzsRJMyp6EfEb1xbIQbHpU6m7bBqvDp1FUpL/PS8jEKnoRcSv1KhckhOjJ1D33UG8t3Cj13GCgopeRPxOy2Y12Pqvp8nSpzvRv/zldZyAp6IXEb9jZjwxuBOnIwsxv+crHDsd53WkgKaiFxG/lCNzODUnv0PVDV8w7K05+OPzrQOFil5E/FalayM5/uZoao0dyNTPt3kdJ2Cp6EXErz3Sshabm7cjrld3YvYf9zpOQFLRi4hfMzPav9GdpGzZmdlriB4wngoqehHxe3mzRVBp4hiqf7eAEe/M13z9P6SiF5GAUK1CEQ4Mfouokf2Y/e0ur+MElCsWvZlNMrODZhZz3rKWZrbVzJLMLOqfrCsiklrtHr2LLQ0e4nC37uw6eMLrOAEjJUf0k4EGFyyLAVoAq1KxrohIqoSGGI8N702EOSY/N5zY+ESvIwWEKxa9c24VcOSCZdudcztSs66IyNUomDMzpSe+S80vZjHy/aVexwkIfjNHb2YdzCzazKIPHdKT4UXk0mpXLsGeF9+g4ut9WBS9x+s4fs9vit45N8E5F+Wci4qMjPQ6joj4uSefbMSW2k3Z3aUn+47q/vWX4zdFLyLyT4SHhtD67RfJffIoEwZN1C2NL0NFLyIBq1jerESMGMntM8cw/fMtXsfxWyk5vXI6sBooZ2b7zKy9mTU3s31ADWChmS1Nfm0RM1t0uXV9MwwRyaha3lOJtS07EPvcc+w5fMrrOH7J/PEKs6ioKBcdHe11DBEJEAeOn2FRneb8Vq0u/cY8R2iIeR0p3ZnZOufcRa9r0tSNiAS8wrmykPett6g5fxJTF633Oo7fUdGLSFC4747y/NC6GyF9erHr4Emv4/gVFb2IBAUzo8NL7TiZIw8f9RlOQmKS15H8hopeRIJGgRyZKTZyGNWXTWfy3DVex/EbKnoRCSqNq5Xihyf7kuOFXvx4QA8qARW9iAQZM6PTs604UKgEc3q+RrymcFT0IhJ88mXPRIWRb3DT1wuYNP0rr+N4TkUvIkGp/s0lWN+lHwX6P0vMb0e9juMpFb2IBK0uz9zPrlKVWPTMAM4mZNx716voRSRo5coaTtTIV6mwfiUTJy/zOo5nVPQiEtTqVCxCTM9BFBvclw2/HPY6jidU9CIS9Lo81ZiY66vxeY8BxCVkvLNwVPQiEvRyZA4n6o2XqbD+K6bO/NLrOOlORS8iGcI9lYqyruML5Hq5D7v+POF1nHSloheRDKNT1xbsLXIdc3sPyVBPpFLRi0iGEZkjE8VeH0SlFXP4ZMk6r+OkGxW9iGQoD9S8jlWPdCOhdy/+OB7rdZx0oaIXkQzFzOjQ91GOZ8rB1P7v4I9P2UtrKnoRyXBK5MtGxOBXqDTrfZat/snrOD6noheRDKlt/UqsbNKO/d16cvx0vNdxfEpFLyIZUlhoCP96uSMhp08zZei/vY7jUyp6EcmwKhbNxZF+r1BuwjBWb/3N6zg+o6IXkQytU4vqrKrVnE1dehMbH5x3uFTRi0iGliUilMaDe5Bn324+HPep13F8QkUvIhlezTKR/NRzAEXffJmte4PvDpdXLHozm2RmB80s5rxlLc1sq5klmVnUZdZtYGY7zGyXmfVNq9AiImmta9u7+L5ybVZ2eZGEIHvObEqO6CcDDS5YFgO0AFZdaiUzCwXGAA2B64FWZnZ96mKKiPhW7qwRVB38PMW3rGHmzJVex0lTVyx659wq4MgFy7Y753ZcYdWqwC7n3G7nXBwwA2iW6qQiIj7WqEpRVj/5PNn69eG3v055HSfN+HKOvihw/vlK+5KXXZSZdTCzaDOLPnTokA9jiYhcnJnRpXMzfipSlrnPDQma2yP4sujtIssu+V/NOTfBORflnIuKjIz0YSwRkUsrmjsLkQP6U27pJyz5MubKKwQAXxb9PqDYed9fA/zuw+2JiKSJNnXKsaxFR4706BEUt0fwZdGvBcqY2bVmFgE8DMz34fZERNJEaIjRru9jnE4ypg6Z5HWcq5aS0yunA6uBcma2z8zam1lzM9sH1AAWmtnS5NcWMbNFAM65BKALsBTYDnzsnNvqq4GIiKSlikVzcfz5gZSZ+BZrt+3zOs5VMX98syEqKspFR0d7HUNEMrjTcQkMfLQ/pQ/toe2yD8kUFup1pEsys3XOuYte16QrY0VELiFrRBgNBnQjz6+7mD5xoddxUk1FLyJyGXUqFGRjl5eIfO0lfv7juNdxUkVFLyJyBd3a12Nd2aos7j4gIM+tV9GLiFxBgRyZKTPgeUp9u4x5C3/wOs4/pqIXEUmBh2qUYsUjPYl/rieHTsR6HecfUdGLiKRASIjRqeeDHMien4/7v+N1nH9ERS8ikkJlCuaAF/tTZub7fLv+Z6/jpJiKXkTkH+jQuApL6z/Kzs49ORMXGI8eVNGLiPwDmcNDuf+lpwg9eoSPRs/0Ok6KqOhFRP6hmqXzs7tnf6556xVidh/0Os4VqehFRFKhe+s7+SqqPt907uv3jx5U0YuIpEKurOHcPvBZCv24kRlTl3kd57JU9CIiqdSwUhG+7/gCuQY8z56DJ7yOc0kqehGRVDIzunVswpZrKzO/+yC/vT2Cil5E5CoUypWZawe9RJlVi5i7wD9vj6CiFxG5Sg/VLM0XrXuQ0Lsnfxw/43Wc/6GiFxG5SiEhxtPPtuKPHJHMfGmU13H+h4peRCQNXJs/G5kHDqLCnA9Y/u12r+P8FxW9iEgaaVevIiuaPM7+br04fjre6zj/R0UvIpJGwkNDaNP/KYg9w5ShU7yO839U9CIiaahi0VwcfXEw5ScMY3XMb17HAVT0IiJprtMDNVh1+31seKaPX9zhUkUvIpLGMoeH0ujVXhTYu5PJE+Z7HUdFLyLiCzWvi+SnXgMp/mZ/Yn494mkWFb2IiI90bncPG2+oybKu/YhL8O4Ol1csejObZGYHzSzmvGV5zWy5me1M/pznEut2M7MYM9tqZt3TMriIiL/LlSWcam++TJmN3zBpylLPcqTkiH4y0OCCZX2Bz51zZYDPk7//L2ZWEXgSqApUBpqYWZmrSisiEmDuvrEIW3oMpMjg51n/y2FPMlyx6J1zq4ALJ5iaAf85SXQKcN9FVq0ArHHOnXbOJQBfAc2vIquISEDq2rEx2ypUZcUz/T05Cye1c/QFnXMHAJI/F7jIa2KAO8wsn5llBRoBxS71B5pZBzOLNrPoQ4cOpTKWiIj/yZE5nDveHkyFzd8y7v3F6b59n70Z65zbDrwBLAeWAJuAhMu8foJzLso5FxUZGemrWCIinqhZtiC7+rxC6Tdf5Nsdf6brtlNb9H+aWWGA5M8XfTquc26ic+5m59wdnJv+2ZnK7YmIBLyOTzRiR8UarOrRn79j0+9eOKkt+vlA2+Sv2wLzLvYiMyuQ/Lk40AKYnsrtiYgEvCwRodw9+hUqbl3D6PGL0m27KTm9cjqwGihnZvvMrD0wBLjHzHYC9yR/j5kVMbPz0882s23AZ0Bn59zRNB+BiEgAuena/Pze73VufKs/y7f8ni7bNH98xmFUVJSLjo72OoaIiE/EJSQx9v7OnDibxNOfjCJf9kxX/Wea2TrnXNTFfqYrY0VE0llEWAgNxr5GlR+/Z8S7C3z+UHEVvYiIB8oXzcOxQW9Sbcwg5q3z7e2MVfQiIh55pJpEG6oAAARUSURBVPU9/HZzLTa9MJg/jsf6bDsqehERj4SGGI3ffY1bdq5l6Oh5PpvCUdGLiHioZMGcxL42nDvfe41p3/3ik22o6EVEPHb/w3U5UPUOdrw6nNNxl7yBQKqFpfmfKCIi/4iZ0WT865w4fIysEWlfyyp6ERE/UCRvdsib3Sd/tqZuRESCnIpeRCTIqehFRIKcil5EJMip6EVEgpyKXkQkyKnoRUSCnIpeRCTI+eWDR8zsELD3gsX5gcMexElrwTIO0Fj8VbCMJVjGAekzlhLOuciL/cAvi/5izCz6Uk9PCSTBMg7QWPxVsIwlWMYB3o9FUzciIkFORS8iEuQCqegneB0gjQTLOEBj8VfBMpZgGQd4PJaAmaMXEZHUCaQjehERSQUVvYhIkPP7ojezBma2w8x2mVlfr/NcDTPbY2ZbzGyjmUV7neefMLNJZnbQzGLOW5bXzJab2c7kz3m8zJhSlxjLADPbn7xvNppZIy8zpoSZFTOzlWa23cy2mlm35OUBt18uM5aA2i9mltnMfjCzTcnjGJi83NN94tdz9GYWCvwE3APsA9YCrZxz2zwNlkpmtgeIcs4F3EUgZnYHcBL4t3OuYvKyN4Ejzrkhyf8I53HO9fEyZ0pcYiwDgJPOuWFeZvsnzKwwUNg5t97McgDrgPuAxwiw/XKZsTxIAO0XMzMgm3PupJmFA98A3YAWeLhP/P2Iviqwyzm32zkXB8wAmnmcKUNyzq0CjlywuBkwJfnrKZz7xfR7lxhLwHHOHXDOrU/++gSwHShKAO6Xy4wloLhzTiZ/G5784fB4n/h70RcFfjvv+30E4M4/jwOWmdk6M+vgdZg0UNA5dwDO/aICBTzOc7W6mNnm5Kkdv5/uOJ+ZlQRuAr4nwPfLBWOBANsvZhZqZhuBg8By55zn+8Tfi94ussx/55qu7Dbn3M1AQ6Bz8hSC+Id3gdJAFeAAMNzbOClnZtmB2UB359zfXue5GhcZS8DtF+dconOuCnANUNXMKnqdyd+Lfh9Q7LzvrwF+9yjLVXPO/Z78+SAwl3NTU4Hsz+S51f/MsR70OE+qOef+TP4FTQLeI0D2TfI88GxgmnNuTvLigNwvFxtLoO4XAOfcMeBLoAEe7xN/L/q1QBkzu9bMIoCHgfkeZ0oVM8uW/CYTZpYNqAfEXH4tvzcfaJv8dVtgnodZrsp/fgmTNScA9k3yG38Tge3OuRHn/Sjg9sulxhJo+8XMIs0sd/LXWYC7gR/xeJ/49Vk3AMmnU40EQoFJzrlXPY6UKmZWinNH8QBhwEeBNBYzmw7cybnbrf4JvAx8CnwMFAd+BVo65/z+Tc5LjOVOzk0POGAP8NR/5lT9lZndDnwNbAGSkhe/wLm57YDaL5cZSysCaL+YWSXOvdkayrkD6Y+dc4PMLB8e7hO/L3oREbk6/j51IyIiV0lFLyIS5FT0IiJBTkUvIhLkVPQiIkFORS8iEuRU9CIiQe7/AbyYtaBwQD1jAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(norm_rpd_calculation_points[0:31],\n",
    "    normalised_transverse_rpd_s6[0:31])\n",
    "zdisk_plots.plot_fitted_model(\n",
    "    norm_rpd_calculation_points[0:31],\n",
    "    31.,\n",
    "    params_optimised,\n",
    "    params_covar,\n",
    "    trans_model_with_info,\n",
    "    plot_95ci=True\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(params_optimised,\n",
    " params_covar,\n",
    " params_1sd_error) = zdisk_modelling.fitmodel_to_hist(\n",
    "    norm_rpd_calculation_points[0:31],\n",
    "    normalised_transverse_rpd[0:31],\n",
    "    trans_model_with_info.model_rpd,\n",
    "    trans_model_with_info.initial_params,\n",
    "    trans_model_with_info.param_bounds\n",
    "    )\n",
    "print('')\n",
    "print('Initial parameter guesses:')\n",
    "print(trans_model_with_info.initial_params)\n",
    "print('')\n",
    "print('Parameter bounds:')\n",
    "print(trans_model_with_info.param_bounds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(norm_rpd_calculation_points[0:31],\n",
    "    normalised_transverse_rpd[0:31])\n",
    "zdisk_plots.plot_fitted_model(\n",
    "    norm_rpd_calculation_points[0:31],\n",
    "    31.,\n",
    "    params_optimised,\n",
    "    params_covar,\n",
    "    trans_model_with_info,\n",
    "    plot_95ci=True\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "normalised_transverse_rpd_s8 = np.load('./normalised_transverse_rpd_smoothed_Churchman-8.npy')\n",
    "calculation_points = np.arange(fitlength + 1.)\n",
    "norm_rpd_calculation_points = calculation_points[calculation_points > 0.]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(params_optimised,\n",
    " params_covar,\n",
    " params_1sd_error) = zdisk_modelling.fitmodel_to_hist(\n",
    "    norm_rpd_calculation_points[0:31],\n",
    "    normalised_transverse_rpd_s8[0:31],\n",
    "    trans_model_with_info.model_rpd,\n",
    "    trans_model_with_info.initial_params,\n",
    "    trans_model_with_info.param_bounds\n",
    "    )\n",
    "print('')\n",
    "print('Initial parameter guesses:')\n",
    "print(trans_model_with_info.initial_params)\n",
    "print('')\n",
    "print('Parameter bounds:')\n",
    "print(trans_model_with_info.param_bounds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(norm_rpd_calculation_points[0:31],\n",
    "    normalised_transverse_rpd_s8[0:31])\n",
    "zdisk_plots.plot_fitted_model(\n",
    "    norm_rpd_calculation_points[0:31],\n",
    "    31.,\n",
    "    params_optimised,\n",
    "    params_covar,\n",
    "    trans_model_with_info,\n",
    "    plot_95ci=True\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.11603027e+01 1.30711725e-02]\n",
      " [2.40419243e+01 3.66960752e-02]\n",
      " [5.80236918e+00 1.55535221e-02]\n",
      " [1.28784895e+02 1.16939837e+00]\n",
      " [1.99900011e+02 2.28218321e+00]\n",
      " [1.05942593e+01 5.64415618e-03]]\n",
      "SSR = 5.0241686788911083e-05\n",
      "AIC = -399.31223265322444\n",
      "AICcorr = -394.44266743583313\n",
      "\n",
      "Initial parameter guesses:\n",
      "[20, 40.0, 5, 1, 1, 20]\n",
      "\n",
      "Parameter bounds:\n",
      "([0, 0, 0, 0, 0, 0], [50, 50.0, 20, 1000, 1000, 100])\n"
     ]
    }
   ],
   "source": [
    "(params_optimised,\n",
    " params_covar,\n",
    " params_1sd_error) = zdisk_modelling.fitmodel_to_hist(\n",
    "    norm_rpd_calculation_points[0:31],\n",
    "    normalised_transverse_rpd_s5[0:31],\n",
    "    trans_model_with_info.model_rpd,\n",
    "    trans_model_with_info.initial_params,\n",
    "    trans_model_with_info.param_bounds\n",
    "    )\n",
    "print('')\n",
    "print('Initial parameter guesses:')\n",
    "print(trans_model_with_info.initial_params)\n",
    "print('')\n",
    "print('Parameter bounds:')\n",
    "print(trans_model_with_info.param_bounds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x at zero bg is -0.05299779236206779.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Code\\perpl-python3\\models_2d_distances_normalised.py:101: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  / x_values # Normalisation\n",
      "C:\\Code\\perpl-python3\\models_2d_distances_normalised.py:104: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  / x_values # Normalisation\n",
      "C:\\Users\\fbsacu\\Miniconda3\\envs\\PERPL\\lib\\site-packages\\numdifftools\\extrapolation.py:489: RuntimeWarning: invalid value encountered in less_equal\n",
      "  converged = err <= tol\n",
      "C:\\Users\\fbsacu\\Miniconda3\\envs\\PERPL\\lib\\site-packages\\numdifftools\\limits.py:173: RuntimeWarning: invalid value encountered in less\n",
      "  outliers = (((abs(der) < (a_median / trim_fact)) +\n",
      "C:\\Users\\fbsacu\\Miniconda3\\envs\\PERPL\\lib\\site-packages\\numdifftools\\limits.py:174: RuntimeWarning: invalid value encountered in greater\n",
      "  (abs(der) > (a_median * trim_fact))) * (a_median > 1e-8) +\n",
      "C:\\Users\\fbsacu\\Miniconda3\\envs\\PERPL\\lib\\site-packages\\numdifftools\\limits.py:175: RuntimeWarning: invalid value encountered in less\n",
      "  ((der < p25 - 1.5 * iqr) + (p75 + 1.5 * iqr < der)))\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3QUZfvG8e+dRg0lEEpooVepC9KL9CIdlaZSpEhTOgLSRKUjIihNQJpgpb1UwSCCkChIICAdQkvoHVKe3x/E9/DjDSSEbGZ3c3/O4SSZ7GSucczF8OzMPGKMQSmllOtyszqAUkop+9KiV0opF6dFr5RSLk6LXimlXJwWvVJKuTgPqwPEJnPmzMbf39/qGEop5TSCgoIuG2N8Y/ueQxa9v78/gYGBVsdQSimnISKnn/Y9HbpRSikXp0WvlFIuToteKaVcnBa9Ukq5OC16pZRycVr0Sinl4rTolVLKxcV5Hb2ILACaAGHGmBIxy9oAo4GiQAVjTKwXvYvIKeAWEAVEGmNsiRNbWeVBZBRHLt7i79AbREZFUy6PD0Wze+PhrucMSjmq+NwwtRCYCSx+bFkw0BL4Kh7r1zLGXH7+aMpq/5b6gXM3CD53g79Db3D5+Gn8L5ygQPhJ3KOj+TlHEc74F6NY/mxU8PfB5u9DmdwZSOnpbnV8pVSMOIveGBMgIv5PLAsBEBH7pFKWMMbwy+EwtoSEcejMZe4Hh+B/6QQFwk5S6cppmkTfRzJlxqtEcbLUqY6bpycXNm8n+qfvuXP/Ifsz5WNOjqIcylOcHIXzUt7f59GfvD6kT+Vp9e4plWzZ+xEIBtgkIgb4yhgz52kvFJFuQDeA3Llz2zmWetLFG/cZsyqILIu/otapP2nk5U6kfz5SliqJX+fO5HilKh7p0//Pejk7dQDAPHxIqW0BnF6zmQe/zyXq58sc9c7KmhzF+Khgedq2rkrnqnnx1CEepZKcvYu+ijHmvIhkATaLyGFjTEBsL4z5S2AOgM1m0/kNk0h0tGHpnjNs/fwb2v+6lIjWb1Bry3w8vZ7vDFy8vMhUvw6Z6tcBHv3roHTwIcr+uJ4bK2aye+8aWjfrwoi2L1Pe38ceu6KUegq7Fr0x5nzMxzAR+RGoAMRa9CrpHQu7xYSvNlFzyVSa+Gam4C8b8C+UJ1F+toiQ9qXilHqpOGbkQLw+/YIKs/oz9/fGrHqzA0MbFcMnjVeibEsp9Wx2+3e0iKQREe9/Pwfq8ehNXGWxB5FRzFgfzNLXetNq/hh8Bg+i1Y6fEq3knyQiVBzWG9v+XbTwvEKVwZ3oOmg+K/eeJTpa//GmlL3FWfQishzYBRQWkVAR6SIiLUQkFKgErBORjTGv9ROR9TGrZgV+E5H9wB5gnTFmg312Q8VX0OlrDO05hSLvtCanf3bK7d9Nwy7Nk+SN9bQ+GWj4w9cUWTyX7juWcLFHT96etpEjF2/ZfdtKJWdijOOdUdlsNqPPo09ct+5HMGvxL+SYOo4o73T4fzaJGpWLW5YnOtqwZdJc3D+fxpqSdcnWsxt96xchtZdDTpGglMMTkaCn3aukl0AkA6HX7jKt7UBeGjeAG2/3oNVvqy0teQA3N6HekG6U3v8HldNE8FKvtvR870sOnr9haS6lXJGePrm489fvsaJtX3KGheK7fRuN8sc605hlMmVKx2urZhO0Yz9tOnfl67PHaTb1A6oVdKycSjkzPaN3YReu32NVmx6kvXKRcj8uo7yDlfzjylUrRdldWykXdoRdb/Xih8AzVkdSymVo0buoSzfu8WPLrrjfukmF7xZTKo/jX7vulzkdjbb/RIasmbjdsQOzNwTjiO8hKeVstOhd0KWb91jz6lvcexBBpRXzKOsEJf+v9Km8eGvVTG7VbUzO7u35eOGvROklmEq9EC16FxN+8x4bGrYlHE+qfTObcv6ZrI703FJ4uNNz+mAu9R1C+eE9GPHJt9x7GGV1LKWclha9C7ly6z6b6rXhtFcGan49gwr5nK/k/+XmJnQd0I5rM+ZQY/ZYxvSfwdU7D62OpZRT0qJ3EVdv3mPrK80JSZud2nMnUym/85b8415vXR23VT9SYcMyZrw5jNNX7lgdSSmno0XvAq7fvMf2Wq8S5JOXerM+pkqBzFZHSlT1KhUmz4a15DwVws/NOrH/9FWrIynlVLTondz1G3cIqN6Q37IWpcGMMVQv5LiXUL6IcgWyUus/33IvXUaCGzbn94PnrI6klNPQondid+9HsKNmEzbnLEPjKSOoVTiL1ZHsKn8Wbzp9+xn/lKlGaMs2/BFy3upISjkFLXonZYxhXbte7Emfh8afDqF20axWR0oSWbxT0vvLDzlQtianm7cm8MgFqyMp5fC06J3Uhk+/4l7IYTINH0aDEtmsjpOkfL1T0OerD9lftiYnm7Xir6MXrY6klEPTondC+7fuRmZ9zs6+Y+lbu5DVcSyRJV1K+nz5IX+VrsmxV1vy9/FLVkdSymFp0TuZ8PPhXOnUmXlthzHhrcq4uSXfCdqzpU9Jry9HElSqBv80bsHBk2FWR1LKIWnRO5GIiEiCGrZiQZW2jO73KhlS61R8OTKkotfskewpWYOQRs05fDrc6khKORwteiey+c0+7MyYj2YD3qZEjvRWx3EYuXxS02v2CPaUqM7BBs04euay1ZGUciha9E7i92nzubzvIFF9+9PGlsvqOA4nT6Y09Jg1gl3Fq/N3/WYcD71idSSlHIYWvRM48Vsgt6dOZW3XEYxqau3MUI4sn29aeswazs7i1dhfrymnzukdtEqBFr3DuxN+hbPt32Jaq0FM61KNFB7uVkdyaAWyeNN95nB+K1aNfXWacPbCNasjKWW5OIteRBaISJiIBD+2rI2IHBSRaBGJdTLax17rLiJ/icjaxAicnERHRbGnfktmV3yd4b2akCNDKqsjOYXC2bzpOuMDthetRmCdVzkXpvPQquQtPmf0C4EGTywLBloCAfFYvx8Q8nyxFMBvnfqyNV0+qvVoR9WCrvWgMnsr5peOLjM+IKDgy/xR+1UuXL5ldSSlLBNn0RtjAoCrTywLMcYciWtdEckJNAbmJThhMnVg1kJOBQZz8e136Vkjv9VxnNJLOdPz5owR7Mxbjt9rN+XS1dtWR1LKEvYeo58ODAai7bwdl3LpwBHCPp3E0vZDmPJ66WR9U9SLKpM7I21njGR3zpfYUbsZ4df1efYq+bFb0YtIEyDMGBMUz9d3E5FAEQkMD0++N71ER0RwpE17JjboxfTOVUmfytPqSE7P5u9D689HsSdbUX6t3ZwrN+9ZHUmpJGXPM/oqQFMROQWsAF4RkSVPe7ExZo4xxmaMsfn6uuYz1ePjt3f6syZnOdq93Yii2dNZHcdlvJwvE80+H01g5gL88kpzrt3SslfJh92K3hgzzBiT0xjjD7wB/GKM6WCv7bmCo6s3EbZzD5fbdebNSnmsjuNyqhTITMPPx/JXxtxsrtOS63ceWB1JqSQRn8srlwO7gMIiEioiXUSkhYiEApWAdSKyMea1fiKy3r6RXdO9q9e42Kcfs1sOYGKb0ojouLw91CjkS+0ZH/F3Wj821W3Nzbta9sr1iTHG6gz/w2azmcDAQKtjJKlfazTmqyzl6TK2Z7KZRMRKGw9e5PfeQykWcZ3GW1aRNqW+F6Kcm4gEGWNiva9J74x1AH9N+oJ9t93I3761lnwSqV88GxU++4QQ93SsbfAGdx9EWh1JKbvRorfY5YNHuDJzNhte68MHjYpaHSdZaVwyO2U++4QjUalY90pLbty5b3UkpexCi95C0RERHG7dnkn1ejK548uk8tLn2CS1ZqVzUGbmRPam8WNblYZcuqyPS1CuR4veQrt6DGR1jrK07tiQ4n76fHmrNC3lR+OZ49jiX569Vepx+rROS6hcixa9RU6u3czFgN2ce70TXarmtTpOslejkC9dZgzjJ1sTDtesx6Hgk1ZHUirRaNFb4MG165zr/R4zm7/PlNfK6CMOHESZ3BkZOu09ltfrxLmGjQn8bb/VkZRKFFr0FtjTsiOzbK0Z/GYNsqVPaXUc9ZgCWbwZP6E7S9r05/rrr7N9dXwe0KqUY9OiT2J/T5nFvpuGXG+0ouFL2a2Oo2KRI0Mqpo17k2VdRvOwZzfWLVxtdSSlXogWfRK6dvgo4TNmsbpVbz58tZjVcdQz+KTx4vPhrfnx/cnIyGGsmrzI6khKJZgWfRKJjozkYKv2TKzTg0kdK5ImhYfVkVQc0qTwYPp7Ddn64Sy8Zk7jm6HTcMQ7yZWKixZ9Etn97mDWZCtF0w4NKJ0rg9VxVDyl8HBnYpfq7J84F6/vv2Vx+/7cuh9hdSylnosWfRI4sWEbl7b9xtk2neihs0U5HXc3YWQbG1dmL+T86QtsKV+bQ4dOWx1LqXjTorez+9dvcK5HL2Y2e5+pbcvirpdSOiUR4d06Raj1zRdsKNeQc3Xr8/3cn3QoRzkFLXo729X6beaWbcGQTrX0UkoXUDFfJj75YgBrB0zG/eOxzO8wgOt3HlodS6ln0qK3o8DP5nEk/C553mxLveLZrI6jEolPGi+mv9+YKwtWcOXUOTZXrEvQ36esjqXUU2nR20nYPye4MWUq69q+x/DG+lRKVyMidKlViAYrviKgfAOuNGzA0tk/EB2tQznK8WjR20FUZBQHW7Zjep3uTHq7Cik99amUrqpUrgyMnzWQbUOmknbCWGa3H0jYTX3csXIsWvR2sK3PcLZlKkybd5pRJJtO8O3q0qX05OM+Dbm77Afunj7Ltsr1Cdh7zOpYSv2XFn0iO7xxB3c2buLC2+/SsaJO8J1ciAhtK+ej6Q/z2VOxIfebNWFmh8EcOnPV6mhKadEnpjs3bnGxew/mthrAhNfK6ATfyVDhbN6Mnz2Is4t+4PrFy5ypXJUZA6dx9sodq6OpZEyLPhH9+lpXlpZszLDu9fFJ42V1HGWRlJ7udKlbjD5rvuLw5HlIwHYCy1dn5pRvuXL7gdXxVDIUZ9GLyAIRCROR4MeWtRGRgyISLSKxzjouIilFZI+I7I957ZjEDO5odny+iIvnwsnfszNVCmS2Oo5yAOlTefLeG5Vps/lbDr8/igzzZrK2ciPmLt3G3Yc6GblKOvE5o18INHhiWTDQEnjWw7ofAK8YY0oBpYEGIlIxISEdXeiRkzyc8Ckb3xrMgHqFrI6jHEy29CkZ1qcpFX9dz5lX25J1+HvMr92e5Vv+JiIq2up4KhmIs+iNMQHA1SeWhRhjjsSxnjHG3I750jPmj8tdZHzv3gNCmr7GrPrdmdCpGp7uOhqmYlcgizcjx3clx/Zt3Chcgkxvvca0Jt35cuVOzly5a3U85cLs2koi4i4i+4AwYLMx5o9nvLabiASKSGB4eLg9YyUaYwybWnZibd6X6d7/Dfwzp7E6knIC5fx9GD53BO7rNyNZs+E7rB+7Kr7CuK7jWL7jKDfu6tMxVeKya9EbY6KMMaWBnEAFESnxjNfOMcbYjDE2X19fe8ZKNGtHTuPSuXDyD+yrjzhQz0VEqF0qF4MWjqNi0A7ujJ1AutBT+LzRlCXVWjBu1NdsOXhRh3ZUokiS2S+MMddFZDuPxvqD43i5U/h99XbcFn9N8Jgv+ax2QavjKCeWI0MqOretgXmjOgdCr3N86Vqyf7eCqHkTmVisMh5vtKVe3XKUyplBJ5JXCWK3ohcRXyAipuRTAXWACfbaXlI6fvwcd3u9y/JOo5jT/mW9Xl4lChGhZK6MlBzakYhB7Qn4+wzp5iwh24ThnBwfybp8pYisXI28datTuVAWCmRJq//vqXiRuJ6nLSLLgZpAZuASMIpHb85+DvgC14F9xpj6IuIHzDPGNBKRksAiwJ1HQ0QrjTFj4xPKZrOZwMDAhO2Rnd24+4BtFeqwpmwjxnzWj5wZU1sdSbm4G/ci+HX3YS6s2YDHrp34nT3KzVTeHClQGreq1clfuzKVC/qS2ye1Fn8yJiJBxpjYL3d3xIkTHLXoo6INi1r15OzVuzT4ehov58tkdSSVDJ29epe9uw9xae0GUuzdRY7zJwj3zsTxQmXwrFGTAtXLUyFvJvJk0uJPTrToE8nCsXNIseRrohYuo0PlvFbHUQpjDCcu3yHot/1cXreR1EG7yRYeyrXU6TmR7yV4uSI5a1WhfKHsFMnmrWP8LkyLPhGsW7cLz26d2DVhPqM6VLE6jlKxio42HA+/zb49IYRt3IrHn3vJGXqUCDcPjuYqwsNyFchSpyZlX8pDqZwZ8ND7PlyGFv0L+vvYRc69Upf/vDWEyaPa4eWhvxzKeYReu8ve4DOc2/gL7P4DvxPBpIx4wOG8LxHVsAnlmtemSkFfUnjovAnOTIv+BVy6cY+Nr7QguGBZhswbRaa0KayOpNQLuXz7AX8cucDRFatJs3UDeS6e4Eiuotyt24CSrRtTo1g2UnslyZXXKhFp0SfQg8goZnYYgteJY9T+aTHF/HQSEeVaHkRG8fs/YexfuR6PjespdPogp7L6c71mfQq+0YxXSuYifSpPq2OqeNCiT4DIqGimjF5A4W++wOuHn2lUNreleZSyt8ioaPacvMLen37BrF1N4WP7CEvvy4WWHWjd53UKZPG2OqJ6Bi365/QgMooJYxZS/uupnP/ia7o0q2BZFqWsEB1t2Bd6nYDVAaSdP5uMNy9ztF0PLXwHpkX/HO49jOKjMQuptmgK4bMX0fHV8pbkUMpRXL3zkOUrtpHi86lkvHmFo+170Kr36xTMqoXvSLTo4+n2g0hGj15EnW+mcGfuYlo1jPW/mVLJ0tU7D1mxYhspZk4lw80r/NO+J617vaaF7yC06OPhxt0IRoxZRJOlU4j8eimN65ZN0u0r5Sz+W/ifTyXDrUeF36rXaxTSwrfUs4peLwjn0eVmg0YtoumSyXgs0pJX6ll80njxbpf6tNy5hksfjCfL6pUcqFyLiZO+1SkSHVSyL/qLN+4z8MNvaLN8CqmWraB2bS15peIj42OFf37IOIrO/Jip7YZw5MJNq6OpJyTroj979S79Ry+lw8opZFyxgqo1SlsdSSmnkzGNF326NSTThv/gc/4Uu+s1Z+X2EBxxWDi5SrZFfzz8Nv3HLqfrqkn4frscW9VSVkdSyqlVLupHm43LOFe1Lunbt2T8pFXceaBDOY4gWRZ9yIWb9B/3LT2/m4Tfqm8pVUXP5JVKDL7eKRjyxVDOfDKTMjPHMbHjBxy+qEM5VktWRX8/IorPthyl3+jlvPfjJPJ89y1FKumZvFKJyd1NeOfN2mTcsIEcZ46ws34rVgYc0aEcCyWbot92JIyGk7cRNv5jRm/7kiI/rSR/RS15peylcrEctNi8kvCXa5CuXQvGTf6e2zqUYwmXf0Rd6LW7jF1ziBPbdzFi21wyN29CyS3zEHd9JKtS9ubrnYJBXw5n4eIK2Eb259OgINpNHaYPCExiLnvD1IPIKObtOMmXmw7RYcdy6t4+SZElC0hdolgipVRKPY/dIec4+GZ3bnukoPqKuZTJ42N1JJeS7G6Y2nE0nIbTd/CfBT8x79thtG9UhjJ/7dKSV8pCFYvmoPGW7/BM6cUvHd4l+NwNqyMlGy41dHPhxj0+WhvC1qCTDPpjCVXlOvm3rsUjr87vqpQjyJY+JU2/n8d/mnRgzZvv4bVshj46IQnEeUYvIgtEJExEgh9b1kZEDopItIjE/hAdkVwisk1EQmJe2y8xgz/pxt0I6k0L4Or6jaxa/SGtOjag0I5ftOSVcjA5fdJQ9+dvyH7jEiveHsSJ8NtWR3J58Rm6WQg0eGJZMNASCHjGepHAAGNMUaAi0EtE7DZ2kj61J7OuBTDx2u8U27GFDD26I6Iz3ivliPx901J1zTLynT/K1++M5OzVu1ZHcmlxFr0xJgC4+sSyEGPMkTjWu2CM+TPm81tACJDjBbLGqeqnI8i1fg3u2bPbczNKqURQIHsGyq5ZRYl/ApnVcxwXbtyzOpLLSpI3Y0XEHygD/PGM13QTkUARCQwPD0/YdlKlStB6SilrFMvtQ5HV31H+r1+Y1nsCYbfuWx3JJdm96EUkLfA98J4x5qn3Qhtj5hhjbMYYm6+vr71jKaUcRKkCWcm9+geq/L6Gie9P5+qdh1ZHcjl2LXoR8eRRyS81xvxgz20ppZyXrUgOsnz/PTW3LGf84C+4cS/C6kguxW5FL4/eCZ0PhBhjptprO0op11CptD/eK1dRZ818xgyfp49LSETxubxyObALKCwioSLSRURaiEgoUAlYJyIbY17rJyLrY1atAnQEXhGRfTF/GtlpP5RSLqBGhUJ4LF1B/ZWfM2r0Yh5GRlsdySW47CMQlFLOa+3GQLw6dWD/hDkM6ljd6jhOIdk9AkEp5dya1LdxrPcw/McMZNOB81bHcXpa9Eoph9RlcAcuFitHUN9hekPVC9KiV0o5pBQe7jSbP4mC5//hs9HzdLz+BWjRK6UcVm7ftKSfP496381i+rIdVsdxWlr0SimHVrdqcY73HUG+MQPYsD/U6jhOSYteKeXwugxsx/mXKrCv3zDOXNHx+uelRa+UcnheHm60mD+J/JdOMOPDOTyIjLI6klPRoldKOYVcmdKQcd5c6v34JdOW6Hj989CiV0o5jTpVinHi/ZEUGjuA9ft0vD6+tOiVUk6l8/ttOVu6Egf6DuH0lTtWx3EKWvRKKafi5eFGq3kTyHv5NDNHfsn9CB2vj4sWvVLK6eT0SUOm+fOo+9Ncpn2z3eo4Dk+LXinllGpXKsLxAWMo9NEQth26aHUch6ZFr5RyWl36teFigRJsGzGRO/r8+qfSoldKOS0vDzde/vwTKu/ZwOxvd1odx2Fp0SulnJqtcHaO9hxM1k9GciD0htVxHJIWvVLK6b3Zvz0PvdOzfPQXREbpUy6fpEWvlHJ66VN5knP6ZGquX8w3mw5YHcfhaNErpVxCvUqF+LtVZ+6NHKETlTxBi14p5RJEhHbj+uBz9zpzJi7BEefDtooWvVLKZfhlSEX0J5OosmQ6a/aetDqOw4iz6EVkgYiEiUjwY8vaiMhBEYkWkVhnHX/aukopZU+vv/oyB6o04sjQ0Vy/+9DqOA4hPmf0C4EGTywLBloCAQlYVyml7MbdTWgydQSFTh5g9lfrrI7jEOIsemNMAHD1iWUhxpgjCVlXKaXsrWiODJwfOo6in41l19Fwq+NYzmHG6EWkm4gEikhgeLgeGKXUi+nUuSGnCpZiy5CPk/0TLh2m6I0xc4wxNmOMzdfX1+o4Siknl9LTHdv0j7Ht2cSC75L34xEcpuiVUiqxVS3uR3C3QWQY+wH/XLpldRzLaNErpVxalwHtuOXtw/KRnxMdnTyvrY/P5ZXLgV1AYREJFZEuItJCREKBSsA6EdkY81o/EVn/rHXtsxtKKRU7nzReZJs8gUprF/NdwGGr41hCHPHuMZvNZgIDA62OoZRyEcYYPuk8htQnjtBl42K8U3paHSnRiUiQMSbW+5p06EYp5fJEhCZj+5H3dAiLl2+3Ok6S06JXSiULJXNl5PA7g8gwcWyye+iZFr1SKtl4q3crHnilZMmUb6yOkqS06JVSyUa29CmJ+GAUpRbPYM+x5HNjpha9UipZebNFJYKLVWTbBxOSzeWWWvRKqWQllZc7RcZ+QJmda/kp4JDVcZKEFr1SKtl5tUI+fm/UkQvDR3HnQaTVcexOi14pleyICE1H9CT7hZMsWbLZ6jh2p0WvlEqWyubx4VD3wWScOJbQa659uaUWvVIq2ercvSlX0/qw8pP5VkexKy16pVSy5ZchFVFDR1Js2WyCjl2yOo7daNErpZK1Ts1sBJaszq9DPnbZyy216JVSyVpqLw+KjxpM8d2bWPdrsNVx7EKLXimV7DUrn5dfm3Ti0rDh3H3oepdbatErpZI9Nzeh5bCueF+5xIqF6+Newclo0SulFGDz9+FQt8FknDiO89dd63JLLXqllIrxTteGnM3oxw8fz7U6SqLSoldKqRg5M6YmcuBQCiz7ikNnrlodJ9Fo0Sul1GO6NLURWKQSW0dMtDpKotGiV0qpx6RP5Ynf0AEU3vITv+07aXWcRBFn0YvIAhEJE5Hgx5a1EZGDIhItIrFORhvzugYickREjonI0MQKrZRS9tS+ZmG2VW3OgSFjXeImqvic0S8EGjyxLBhoCQQ8bSURcQe+ABoCxYC2IlIsYTGVUirppPBwp+LQXuQ5tJd1W/6yOs4Li7PojTEBwNUnloUYY47EsWoF4Jgx5oQx5iGwAmiW4KRKKZWEXi2dk1+aduHK8JHcj4iyOs4LsecYfQ7g7GNfh8Ysi5WIdBORQBEJDA9PPnM5KqUck5ub0PL9jqS6fpkfljn3M+vtWfQSy7KnDnYZY+YYY2zGGJuvr68dYymlVPxULpCZPzu+h9cnY7l+96HVcRLMnkUfCuR67OucwHk7bk8ppRJdp26NuZLCm++nL7U6SoLZs+j3AgVFJK+IeAFvAKvtuD2llEp0RbKl4/y7A8j+5VTOXr5tdZwEic/llcuBXUBhEQkVkS4i0kJEQoFKwDoR2RjzWj8RWQ9gjIkEegMbgRBgpTHmoL12RCml7KX769X4O89LrB813eooCSLGON41ojabzQQGBlodQyml/mvKD0G81Ls9fjt+pUT+rFbH+R8iEmSMifW+Jr0zViml4uGdRiXZXrYevw8ea3WU56ZFr5RS8ZAupScFB/Yhz+5f+G3PP1bHeS5a9EopFU/tq+ZnQ532HBs0nCgnejSCFr1SSsWTl4cbtd7vRIYzx1m/ZqfVceJNi14ppZ5Dk5J+bG3Vg3sjnefRCFr0Sin1HNzchLY9W3I30vDj3B+tjhMvWvRKKfWcKufPzIG338N78niu3H5gdZw4adErpVQC9Hy7Nkd88/HzuFlWR4mTFr1SSiVAgSze3Ondn5xL53D8/DWr4zyTFr1SSiXQuy1sBBSvwS+DPrI6yjNp0SulVAJlTpsCvwF9ybN9HX/sd9z5ZbXolVLqBXSuVZj/1Hid4P7DHXZ+WS16pZR6ASk93ak+sCuZjx9i/X/+sDpOrLTolVLqBTUrnZOtLXtwa/hwh7yJSoteKaVekJub0K53GyIeRPDdnJ+sjvM/tOiVUioRVMyXieDO/Uk/ZXgO67oAAAdkSURBVDxhN+9bHef/0aJXSqlE0uPtOhzz9efn8V9aHeX/0aJXSqlEks83Lff7DiT30q84GnrV6jj/pUWvlFKJqHtzGzuLV2fLkI+tjvJfWvRKKZWIfNJ4kWtgP/JvX8vOP09YHQfQoldKqUTXsWZhNtV8jUODRjrETFRxFr2ILBCRMBEJfmyZj4hsFpGjMR8zPmXdfiISLCIHReS9xAyulFKOKqWnO9UHdSfbsWDWrtttdZx4ndEvBBo8sWwosNUYUxDYGvP1/yMiJYB3gApAKaCJiBR8obRKKeUkXi3lx7ZWPbg9ciR3H0ZamiXOojfGBABPvn3cDFgU8/kioHksqxYFdhtj7hpjIoFfgRYvkFUppZyGiNCu7+tE33/Aqrk/W5oloWP0WY0xFwBiPmaJ5TXBQHURySQiqYFGQK6n/UAR6SYigSISGB4ensBYSinlOGz+PoR0G4jP5HGcvXLHshx2ezPWGBMCTAA2AxuA/cBT//1ijJljjLEZY2y+vr72iqWUUkmqd+d6hOQowo/DJluWIaFFf0lEsgPEfAyL7UXGmPnGmLLGmOo8Gv45msDtKaWUU/LLkIoMw4dTZO1ytgVZc7llQot+NfBWzOdvAbEOQIlIlpiPuYGWwPIEbk8ppZzW23WL8+srbTg0wJqnW8bn8srlwC6gsIiEikgX4FOgrogcBerGfI2I+InI+sdW/15EDgFrgF7GGMeeWFEppezAy8ONxqP7kuPUYRYv2ZLk2xdjrL+Y/0k2m80EBgZaHUMppRLVRx8vo+Cimby8czP+mdMk6s8WkSBjjC227+mdsUoplUS69W7JdW8flo79iqQ8ydaiV0qpJJIlXUpSjR6D7fu5bNp3Jsm2q0WvlFJJqF3DMux7uT5/fjA+ye6Y1aJXSqkk5OHuRu2JH1AqeCdzV+5Mkm1q0SulVBKzFcjCkU79SDtxHMfCbtl9e1r0SillgY6D3yJVdCRzJy+z+xuzWvRKKWWBzGlT4DH2I6qtmMnqv0Ltui0teqWUskir5lU4W7w8f4ydxq37EXbbjha9UkpZxN1NqPTZeGrtWc8XP9rvJlEteqWUslDpQtk51bYrKadOIOTCTbtsQ4teKaUs1npsH/JeP8cXs1YTbYc5ZrXolVLKYhnTpMB97HjKXT/Nw6joRP/5Hon+E5VSSj23Jm3rQdt6dvnZekavlFIuToteKaVcnBa9Ukq5OC16pZRycVr0Sinl4rTolVLKxWnRK6WUi9OiV0opF6dFr5RSLk6Sciby+BKRcOD0E4szA5ctiJPYXGU/QPfFUbnKvrjKfkDS7EseY4xvbN9wyKKPjYgEGmNsVud4Ua6yH6D74qhcZV9cZT/A+n3RoRullHJxWvRKKeXinKno51gdIJG4yn6A7oujcpV9cZX9AIv3xWnG6JVSSiWMM53RK6WUSgAteqWUcnEOX/Qi0kBEjojIMREZanWeFyEip0TkgIjsExH7TfluByKyQETCRCT4sWU+IrJZRI7GfMxoZcb4esq+jBaRczHHZp+INLIyY3yISC4R2SYiISJyUET6xSx3uuPyjH1xquMiIilFZI+I7I/ZjzExyy09Jg49Ri8i7sA/QF0gFNgLtDXGHLI0WAKJyCnAZoxxuptARKQ6cBtYbIwpEbNsInDVGPNpzF/CGY0xQ6zMGR9P2ZfRwG1jzGQrsz0PEckOZDfG/Cki3kAQ0Bx4Gyc7Ls/Yl9dwouMiIgKkMcbcFhFP4DegH9ASC4+Jo5/RVwCOGWNOGGMeAiuAZhZnSpaMMQHA1ScWNwMWxXy+iEe/mA7vKfvidIwxF4wxf8Z8fgsIAXLghMflGfviVMwjt2O+9Iz5Y7D4mDh60ecAzj72dShOePAfY4BNIhIkIt2sDpMIshpjLsCjX1Qgi8V5XlRvEfk7ZmjH4Yc7Hici/kAZ4A+c/Lg8sS/gZMdFRNxFZB8QBmw2xlh+TBy96CWWZY471hS3KsaYskBDoFfMEIJyDLOB/EBp4AIwxdo48SciaYHvgfeMMTetzvMiYtkXpzsuxpgoY0xpICdQQURKWJ3J0Ys+FMj12Nc5gfMWZXlhxpjzMR/DgB95NDTlzC7FjK3+O8YaZnGeBDPGXIr5BY0G5uIkxyZmHPh7YKkx5oeYxU55XGLbF2c9LgDGmOvAdqABFh8TRy/6vUBBEckrIl7AG8BqizMliIikiXmTCRFJA9QDgp+9lsNbDbwV8/lbwM8WZnkh//4SxmiBExybmDf+5gMhxpipj33L6Y7L0/bF2Y6LiPiKSIaYz1MBdYDDWHxMHPqqG4CYy6mmA+7AAmPMeIsjJYiI5OPRWTyAB7DMmfZFRJYDNXn0uNVLwCjgJ2AlkBs4A7Qxxjj8m5xP2ZeaPBoeMMApoPu/Y6qOSkSqAjuAA0B0zOIPeDS27VTH5Rn70hYnOi4iUpJHb7a68+hEeqUxZqyIZMLCY+LwRa+UUurFOPrQjVJKqRekRa+UUi5Oi14ppVycFr1SSrk4LXqllHJxWvRKKeXitOiVUsrF/R+OfnODgnsTtAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(norm_rpd_calculation_points[0:31],\n",
    "    normalised_transverse_rpd_s5[0:31])\n",
    "zdisk_plots.plot_fitted_model(\n",
    "    norm_rpd_calculation_points[0:31],\n",
    "    31.,\n",
    "    params_optimised,\n",
    "    params_covar,\n",
    "    trans_model_with_info,\n",
    "    plot_95ci=True\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "file_extension": ".py",
  "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.10"
  },
  "mimetype": "text/x-python",
  "name": "python",
  "npconvert_exporter": "python",
  "pygments_lexer": "ipython3",
  "version": 3
 },
 "nbformat": 4,
 "nbformat_minor": 2
}