Common mean direction test.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## D. Heslop, J. L. Scealy, A. T. A. Wood, L. Tauxe and A. P. Roberts (2023). A bootstrap common mean direction test. Journal of Geophysical Research.\n",
    "### Example common mean direction test (see Section 3 of manuscript) performed on the data of Ao et al. (2013).\n",
    "\n",
    "Please note that because this is a numerical procedure the results may differ slightly from those in the paper."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Import required libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-W- cartopy is not installed\n",
      "    If you want to make maps, install using conda:\n",
      "    conda install cartopy\n",
      "-W- You haven't installed a module for plotting maps (cartopy or Basemap)\n",
      "    Recommended: install cartopy.  With conda:\n",
      "    conda install cartopy\n",
      "    For more information, see http://earthref.org/PmagPy/Cookbook#getting_python\n"
     ]
    }
   ],
   "source": [
    "import pmagpy.ipmag as ipmag\n",
    "import pmagpy.pmag as pmag\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function ```form_Mhat``` estimates $\\hat{M}$ according to equation (4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def form_Mhat(mhat):\n",
    "    \n",
    "    #input - mhat, estimated sample mean direction\n",
    "    #output - Mhat matrix according to equation 4\n",
    "    \n",
    "    b = np.matrix(mhat[0:2][:,np.newaxis])\n",
    "    c = mhat[2]\n",
    "    \n",
    "    if c==0:\n",
    "        A1 = np.eye(2)-(b*b.getH())\n",
    "    else:\n",
    "        A1 = c/np.abs(c)*np.eye(2)-c/(np.abs(c)+np.abs(c)**2)*(b*b.getH())\n",
    "        \n",
    "    A2 = -b\n",
    "    Mhat = np.hstack((A1,A2))\n",
    "    return Mhat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function ```form_Ghat``` estimates $\\hat{G}$ according to equation (5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def form_Ghat(X,Mhat):\n",
    "    \n",
    "    #input - X, collection of directions (one per column)\n",
    "    #input - Mhat, Mhat matrix for mean direction\n",
    "    #output - Ghat matrix according to equation 5\n",
    "\n",
    "    n = np.shape(X)[1]\n",
    "    term1 = np.power(np.linalg.norm(np.sum(X,axis=1)/n),-2)/n\n",
    "    X = np.matrix(X)\n",
    "    Mhat_T = Mhat.getT()\n",
    "    Ghat = np.matrix(np.zeros((2,2)))\n",
    "    \n",
    "    for u in range(2):\n",
    "        for v in range(2):\n",
    "            for i in range(n):\n",
    "                Ghat[u,v] += Mhat_T[:,u].getT()*X[:,i]*X[:,i].getT()*Mhat_T[:,v]\n",
    "            Ghat[u,v] *= term1\n",
    "    \n",
    "    return Ghat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function ```form_Q``` creates rotation matrix to rotation ```b``` to ```a``` according to equations (9) and (10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def form_Q(a,b):\n",
    "    \n",
    "    #input - a, destination direction (unit vector)\n",
    "    #input - b, starting direction (unit vector)\n",
    "    #output - Q, rotation matrix so Qb = a\n",
    "    \n",
    "    a = np.matrix(a)\n",
    "    a = np.reshape(a,(3,1))\n",
    "    b = np.matrix(b)\n",
    "    b = np.reshape(b,(3,1))\n",
    "\n",
    "    c = b-a*(a.getT()*b)\n",
    "    c /= np.linalg.norm(c)\n",
    "\n",
    "    alpha = np.arccos(a.getT()*b)\n",
    "    A = a*c.getT()-c*a.getT()\n",
    "\n",
    "    Q = np.eye(3)+np.multiply(np.sin(alpha),A)+np.multiply(np.cos(alpha)-1,a*a.getT()+c*c.getT())\n",
    "    \n",
    "    return Q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function ```find_CMDT_CR``` provides a collection of points on the boundary of the 1-$\\alpha$ confidence region for the common mean direction according to the procedure in Appendix B.\n",
    "\n",
    "N.b. ```find_CMDT_CR``` should only be used if null hypothesis of common mean direction cannot be rejected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_CMDT_CR(Ahat,Tc,mhat12):\n",
    "    \n",
    "    #Ahat - combined covariance matrix\n",
    "    #input - Tc, T value on the boundary of the confidence region\n",
    "    #input - mhat12, estimated common mean direction\n",
    "    #output - mCI, sequence of points along confidence region\n",
    "    \n",
    "    [D,V] = np.linalg.eig(Ahat)\n",
    "    \n",
    "    idx=np.flip(np.argsort(D))\n",
    "    D = D[idx]\n",
    "    V = V[:,idx]\n",
    "        \n",
    "    mCI = np.zeros((3,201))\n",
    "    y = np.matrix(np.zeros((3,1)))\n",
    "    for i in range(201):\n",
    "            theta = i*np.pi/100\n",
    "            \n",
    "            ylen = np.zeros(201);\n",
    "            phi = np.linspace(0,np.pi/2,201)\n",
    "            for j in range(201):\n",
    "                y[0] = np.sin(phi[j])*np.cos(theta)*np.sqrt(Tc)/np.sqrt(D[0])\n",
    "                y[1] = np.sin(phi[j])*np.sin(theta)*np.sqrt(Tc)/np.sqrt(D[1])\n",
    "                y[2] = np.cos(phi[j])*np.sqrt(Tc)/np.sqrt(D[2])\n",
    "                ylen[j] = np.linalg.norm(y)\n",
    "                \n",
    "            idx = np.argsort(ylen)\n",
    "            phi0 = np.interp(1.0,ylen[idx],phi[idx]);\n",
    "            y[0] = np.sin(phi0)*np.cos(theta)*np.sqrt(Tc)/np.sqrt(D[0])\n",
    "            y[1] = np.sin(phi0)*np.sin(theta)*np.sqrt(Tc)/np.sqrt(D[1])\n",
    "            y[2] = np.cos(phi0)*np.sqrt(Tc)/np.sqrt(D[2])\n",
    "                \n",
    "            mCI[:,i] = np.ndarray.flatten(V*y)\n",
    "    \n",
    "    #Check if points are in the correct hemisphere\n",
    "    mCIbar = np.mean(mCI,axis=1)/np.linalg.norm(np.mean(mCI,axis=1))\n",
    "    if np.arctan2(np.linalg.norm(np.cross(mhat12.T,mCIbar)),np.dot(mhat12.T,mCIbar))>np.pi/2:\n",
    "        mCI *= -1\n",
    "        \n",
    "    return mCI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example processing\n",
    "#### Import the example data sets and transform to Cartesian coordinates (1 direction per column). Estimate the mean direction via equation (3)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "file = \"Ao_normal.csv\" #normal polarity directions from Ao et al. (2013)\n",
    "F1 = pd.read_csv(file,header=None).to_numpy()\n",
    "X1 = np.transpose(pmag.dir2cart(F1)) # normal directions in Cartesian coordinates (one direction per column)\n",
    "\n",
    "file = \"Ao_reversed.csv\" #reversed polarity directions from Ao et al. (2013)\n",
    "F2 = pd.read_csv(file,header=None).to_numpy()\n",
    "X2 = -np.transpose(pmag.dir2cart(F2)) # inverted reversed directions in Cartesian coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Estimate mean directions of data sets and pooled data set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "n1 = np.shape(X1)[1] #number of observations in first data set\n",
    "n2 = np.shape(X2)[1] #number of observations in second data set\n",
    "n = n1 + n2 #total number of observations\n",
    "\n",
    "X12 = np.hstack((X1,X2)) #form pooled data set\n",
    "\n",
    "mhat1 = np.mean(X1,axis=1)\n",
    "mhat1 /= np.linalg.norm(mhat1) #mean of first data set\n",
    "\n",
    "mhat2 = np.mean(X2,axis=1)\n",
    "mhat2 /= np.linalg.norm(mhat2) #mean of second data set\n",
    "\n",
    "mhat12 = np.mean(X12,axis=1)\n",
    "mhat12 /= np.linalg.norm(mhat12) #mean of pooled data set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Form $\\hat{M}$ and $\\hat{G}$ for each data set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "Mhat1 = form_Mhat(mhat1) #Mhat of first data set\n",
    "Ghat1 = form_Ghat(X1,Mhat1) #Ghat of first data set\n",
    "\n",
    "Mhat2 = form_Mhat(mhat2) #Mhat of second data set\n",
    "Ghat2 = form_Ghat(X2,Mhat2) #Ghat of second data set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Form $\\hat{A}_0$ (equation 8) and find minimum eigenvalue ($\\lambda_{\\textrm{min}}$) and corresponding eigenvector ($\\hat{m}_0$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ahat = Mhat1.getH()*np.linalg.inv(Ghat1)*Mhat1\n",
    "Ahat += Mhat2.getH()*np.linalg.inv(Ghat2)*Mhat2\n",
    "Ahat *= n\n",
    "\n",
    "D,V = np.linalg.eig(Ahat)\n",
    "idx = np.argmin(D)\n",
    "Lmin = D[idx] #minimum eigenvalue\n",
    "mhat0 = V[:,idx] #eigenvector corresponding to pooled sample mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Transform data sets to estimated common mean direction $\\hat{m}_0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q1 = form_Q(mhat0,mhat1) #rotation matrix for first data set\n",
    "X10 = np.matmul(Q1,X1) #rotated version of first data set\n",
    "\n",
    "Q2 = form_Q(mhat0,mhat2) #rotation matrix for second data set\n",
    "X20 = np.matmul(Q2,X2) #rotated version of second data set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Main bootstrap sequence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "B = 10000 #number of bootstrap iterations\n",
    "Lmin_b = np.zeros(B) #predefine output array for minimum eigenvalues\n",
    "T_b = np.zeros(B) ##predefine output array for Tb (equation 11)\n",
    "\n",
    "for i in range(B): #loop through boostrap iterations\n",
    "    idx1 = np.random.randint(0,n1,n1) #select observation indicies with replacement\n",
    "    X10_b = np.asarray(X10[:,idx1]) #form bootstrap sample from rotated version of first data set\n",
    "    \n",
    "    mhat10_b = np.mean(X10_b,axis=1) #mean direction of boostrap sample\n",
    "    mhat10_b /= np.linalg.norm(mhat10_b)\n",
    "    mhat10_b = (np.asarray(mhat10_b)).flatten()\n",
    "    Mhat10_b = form_Mhat(mhat10_b) #\\hat{M} for bootstrap sample\n",
    "    Ghat10_b = form_Ghat(X10_b,Mhat10_b) #\\hat{G} for bootstrap sample\n",
    "    \n",
    "    idx2 = np.random.randint(0,n2,n2) #select observation indicies with replacement\n",
    "    X20_b = np.asarray(X20[:,idx2]) #form bootstrap sample from rotated version of second data set\n",
    "    \n",
    "    mhat20_b = np.mean(X20_b,axis=1) #mean direction of boostrap sample\n",
    "    mhat20_b /= np.linalg.norm(mhat20_b)\n",
    "    mhat20_b = (np.asarray(mhat20_b)).flatten()\n",
    "    Mhat20_b = form_Mhat(mhat20_b) #\\hat{M} for bootstrap sample\n",
    "    Ghat20_b = form_Ghat(X20_b,Mhat20_b) #\\hat{G} for bootstrap sample\n",
    "    \n",
    "    Ahat_b = Mhat10_b.getH()*np.linalg.inv(Ghat10_b)*Mhat10_b #bootstrap estimate of \\hat{A}_0 (equation 8)\n",
    "    Ahat_b += Mhat20_b.getH()*np.linalg.inv(Ghat20_b)*Mhat20_b\n",
    "    Ahat_b *= n\n",
    "    \n",
    "    D_b,V_b = np.linalg.eig(Ahat_b) #Eigenvalues and eigenvectors\n",
    "    Lmin_b[i] = np.min(D_b) #minimum eigenvalue for boostrap sample\n",
    "    T_b[i] = np.matmul(np.matmul(np.transpose(mhat0),Ahat_b),mhat0) #Bootstrap T for pooled data (equation 11)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Main CMDT results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CMDT value = 8.22\n",
      "CMDT critical value = 12.27\n",
      "CMDT p-value = 0.14\n"
     ]
    }
   ],
   "source": [
    "p = (1+np.sum(Lmin_b>=Lmin))/(B+1) # p-value (step 8 of CMDT, Section 3)\n",
    "# (n.b., if p > 0.05 cannot reject null of common means at alpha = 0.05)\n",
    "\n",
    "alpha = 0.05 #set significance level\n",
    "Lmin_c = np.quantile(Lmin_b,1-alpha) #test critical value\n",
    "# (n.b., if Lmin > Lmin_c reject null of common means at alpha significance level)\n",
    "\n",
    "print(\"CMDT value = {:.2f}\".format(Lmin))\n",
    "print(\"CMDT critical value = {:.2f}\".format(Lmin_c))\n",
    "print(\"CMDT p-value = {:.2f}\".format(p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot histogram of bootstrap values and compare to test statistic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEOCAYAAAB4nTvgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de3RU9bn/8ffDpQYICAEEJClQD9pqW1EiVWwrraJA5aIFRAs/tBaoYpVf6amAN44Kh9ZLsb8Wf1LPOSKKiNqKsaVVWKK2oIg98agIQiuQCIJgERC5+pw/ZocOIcnMhOzZs5PPa62smf3dl3l2ZsGT72V/v+buiIiI1KRR1AGIiEjuU7IQEZGUlCxERCQlJQsREUlJyUJERFJSshARkZSaRB1AbZjZQGBgy5Ytx5x88slRh9MgrFmzBoBTTjkl4kjSF8eYRbLh9ddf3+bu7TM5x+L8nEVxcbGvXLky6jAahD59+gCwdOnSSOPIRBxjFskGM3vd3YszOSeWNQvJvptvvjnqEEQkQkoWkpYLLrgg6hBEJELq4Ja0lJaWUlpaGnUYIhKRWNYsKjq4/+Vf/iXqUBqMCRMmAGr/F2moYpks3L0EKCkuLh4TdSwicXTgwAHKy8vZu3dv1KFISPLy8igsLKRp06Z1cr1Qk4WZ/V/gB4ADbwJXAc2Bx4GuwHpguLv/Izh+MnA1cAi43t3/FGZ8Ig1VeXk5LVu2pGvXrphZ1OFIHXN3tm/fTnl5Od26dauTa4bWZ2FmnYHrgWJ3/zLQGBgBTAKWuHt3YEmwjZmdGuw/DegHzDKzxmHFJ9KQ7d27l7Zt2ypR1FNmRtu2beu05hh2B3cToJmZNSFRo9gEDAbmBPvnAEOC94OB+e6+z93fA9YBvUKOT6TBUqKo3+r6+w2tGcrd3zezu4GNwKfAc+7+nJl1cPfNwTGbzeyE4JTOwCtJlygPyqq1fft2Zs+efXi7WbNmjBo1qi5vQwLTp0+POgQRiVCYzVBtSNQWugEnAi3MbGRNp1RRdtTj5WY21sxWmtnKnTt3HrHv008/PYaIpSa9e/emd+/eUYchEiv9+/dnzpw5qQ+s5OWXX865aWrCbIa6AHjP3T909wPAb4HewBYz6wQQvG4Nji8HipLOLyTRbHUEd5/t7sXuXpyfnx9i+JJs2bJlLFu2LOowpJ7Lz88//NOoUSOaNWt2ePvRRx/N+Hp9+vThwQcfrFUsmZ47depURo488u/hRYsWMXr06JTnmhnr1q07vP2Nb3zj8NxmuSLM0VAbgbPNrDmJZqjzgZXAJ8BoYEbwujA4/hlgnpndS6Im0h1YEWJ8koEpU6YAes5CwrV79+7D77t27cqDDz6o2QNyRGg1C3d/FXgS+CuJYbONgNkkkkRfM1sL9A22cfe3gQXAKuCPwHh3PxRWfCISH5999hkzZszgpJNOom3btgwfPpyPPvoISIzsGjlyJG3btqV169acddZZbNmyhZtuuomXX36Z6667jvz8fK677rqjrpvpuTfccANFRUW0atWKnj178vLLLwPwxz/+kenTp/P444+Tn5/P6aefDhxZO1m3bh3nnXcexx9/PO3ateOyyy4D4Jvf/CYAp59+Ovn5+Tz++OMsXbqUwsLCw3GWlZVx6aWX0r59e9q2bVvlvYQt1Ocs3P024LZKxftI1DKqOn4aMC3MmESkahWz9CYbPnw41157LXv27GHAgAFH7b/yyiu58sor2bZtG0OHDj1q/zXXXMNll11GWVkZRUVFR+1P1y9/+UuefvppXnzxRdq3b8/111/P+PHjeeyxx5gzZw4ff/wxZWVlHHfccZSWltKsWTOmTZvGX/7yF0aOHMkPfvCDKq+b6blnnXUWt956K8cffzz33Xcfw4YNY/369fTr148pU6awbt06HnnkkSo/65ZbbuHCCy/khRdeYP/+/VTMmP3SSy9hZrzxxhtUzEqRXIM/dOgQF198Md/+9reZO3cujRs3JorZtmM5N5SZDTSz2erQFmkYHnjgAaZNm0ZhYSHHHXccU6dO5cknn+TgwYM0bdqU7du3s27dOho3bkzPnj1p1apVWtfN9NyKWkiTJk2YOHEi+/btS7tvoWnTpmzYsIFNmzaRl5fH17/+9bTOW7FiBZs2beKuu+6iRYsWGZ1bl2I93UeXLl003YdIHampP6p58+Y17m/Xrl2N+4+lVgGwYcMGLrnkEho1+ufft40bN2bLli2MGjWKsrIyRowYwY4dOxg5ciTTpk1La5qLTM+95557ePDBB9m0aRNmxs6dO9m2bVta9/Dzn/+cW265hV69etGmTRsmTpzI97///ZTnlZWV0aVLF5o0ifa/61jWLCT7Zs6cycyZM6MOQxqooqIiFi1axI4dOw7/7N27l86dO9O0aVNuu+02Vq1axbJly3j22Wd5+OGHgdQPpmVy7ssvv8zPfvYzFixYwD/+8Q927NjB8ccfT8UCcqk+q2PHjvzmN79h06ZNPPDAA1x77bVHjICq6d43btzIwYMHUx4bJiULSUuPHj3o0aNH1GFIA/XDH/6Qm266iQ0bNgDw4YcfsnBhYiDlCy+8wJtvvsmhQ4do1aoVTZs2pXHjxExBHTp04O9//3u1183k3F27dtGkSRPat2/PwYMHuf3220l+1qtDhw6sX7+ezz77rMrPeuKJJygvLwegTZs2mFlacfbq1YtOnToxadIkPvnkE/bu3ctf/vKXtH5vdUnJQtKyePFiFi9eHHUY0kDdcMMNDBo0iAsvvJCWLVty9tln8+qrrwLwwQcfMHToUFq1asWXvvQlzjvvvMPPO9xwww08+eSTtGnThuuvv/6o62Zy7kUXXUT//v05+eST6dKlC3l5eUc0rw0bNgyAtm3bcuaZZx71Wa+99hpf+9rXyM/PZ9CgQdx3332HJ/mbOnUqo0ePpnXr1ixYsOCI8xo3bkxJSQnr1q3j85//PIWFhTz++ON18FvNTCzX4K5Yz6J9+/Zj7rzzziP2jR07Npqg6rk4rmcdx5iz5Z133uFLX/pS1GFIyKr7nmuzBncsaxbuXuLuY5s1axZ1KCIiDUIsR0PVRBMLiojUvVjWLNKl5zBEROpGvU4WIiJSN2LZDJXUwR11KA3GAw88EHUIIhKhWNYs1MGdfaecckrOza8vItkTy2Qh2VdSUkJJSUnUYYhIRGLZDCXZd8899wAwcODAiCMRkSioZiEisTd9+vRqpyEHePTRR7nwwguP+XMqr2hXF9avX4+ZRT73UyqxTBaaolykfps3bx7FxcXk5+fTqVMn+vfvz5///Odqj58yZcrhRYaq+s/3e9/7Hs8991zocddnsUwW6uAWqb/uvfdeJkyYwJQpU9iyZQsbN27k2muvPTxxYGW5/hd5fRFasjCzU8ysNOlnp5lNMLMCM3vezNYGr22SzplsZuvMbI2ZXRRWbCKSmz7++GNuvfVWfv3rX3PppZfSokULmjZtysCBA7nrrruAxKR7Q4cOZeTIkbRq1YqHHnqIqVOnHp4AsGKZ0tatW5Ofn8/y5ct56KGHjlgw6O2336Zv374UFBTQoUMHpk+fDiQWGjrnnHNo3bo1nTp14rrrrmP//v0p454/fz7FxUdOtfSLX/yCQYMGAfD73/+eM844g1atWlFUVMTUqVOrvVbXrl2PmLQz+d4AXnnlFXr37k3r1q05/fTTszb3WWgd3O6+BugBYGaNgfeB3wGTgCXuPsPMJgXbN5rZqcAI4DTgRGCxmZ2sdbhzw9y5c6MOQUI0YcIESktLQ/2MHj16pFwTZfny5ezdu5dLLrmkxuMWLlzIE088wcMPP8y+ffv42c9+dnjfSy+9RLdu3dixY8fhBYOSV7PbtWsXF1xwAT/5yU8oKSnhwIEDrFq1CkjM8PqLX/yC4uJiysvL6d+/P7NmzWLChAk1xjNo0CDGjBnD2rVr6d69O5BoSps4cSIALVq04OGHH+a0007jrbfeom/fvvTo0YMhQ4bUeN3K3n//fb7zne8wd+5c+vXrx5IlS/jud7/L6tWrCfu5s2w1Q50P/M3dNwCDgTlB+Ryg4rc1GJjv7vvc/T1gHdArS/FJCkVFRce82plIKtu3b6ddu3YpV4U755xzGDJkCI0aNSLT5uhnn32Wjh07MnHiRPLy8mjZsiVf+9rXAOjZsydnn302TZo0oWvXrowbN44XX3wx5TWbN2/O4MGDeeyxxwBYu3Ytq1evPlyz6NOnD1/5yldo1KgRX/3qV7n88svTum5ljzzyCAMGDGDAgAE0atSIvn37UlxczB/+8IeMr5WpbA2dHQE8Frzv4O6bAdx9s5mdEJR3Bl5JOqc8KDuCmY0FxgIUFBSEFrAcqWL+/MsuuyziSCQMubIKYtu2bdm2bRsHDx6sMWEcyx8uZWVlnHTSSVXue/fdd/nxj3/MypUr2bNnDwcPHqRnz55pXfeKK65g4sSJ3HrrrcybN48hQ4bQvHlzAF599VUmTZrEW2+9xf79+9m3b9/h9S8ysWHDBp544okjnnk6cOAA3/rWtzK+VqZCr1mY2eeAQcATqQ6touyoxTbcfba7F7t7cX5+fl2EKGm4//77uf/++6MOQ+q5c845h7y8PJ5++ukaj6tpCdNUy5sWFRXxt7/9rcp911xzDV/84hdZu3YtO3fuZPr06aS75s+FF17Itm3bKC0t5bHHHuOKK644vO+KK65g0KBBlJWV8fHHH/PDH/6w2uu2aNGCPXv2HN7+4IMPjoh91KhRRywv+8knnzBp0qS0YjwW2WiG6g/81d23BNtbzKwTQPC6NSgvB5L/XCgENmUhPhHJEccffzy3334748eP5+mnn2bPnj0cOHCARYsW8dOf/jSta7Rv355GjRpVu0zpxRdfzAcffMDMmTPZt28fu3btOrzq3q5du2jVqhX5+fmsXr06oz+QmjRpwtChQ/nXf/1XPvroI/r27Xt4365duygoKCAvL48VK1Ywb968aq/To0cP5s+fz4EDB1i5ciVPPvnk4X0jR46kpKSEP/3pTxw6dIi9e/eydOnSw8u1hikbyeJy/tkEBfAMMDp4PxpYmFQ+wsyOM7NuQHdgRRbiE5Ec8uMf/5h7772XO++8k/bt21NUVMSvfvWrtDuDmzdvzk033cS5555L69ateeWVV47Y37JlS55//nlKSkro2LEj3bt354UXXgDg7rvvZt68ebRs2ZIxY8Zk3Ox6xRVXsHjxYoYNG3ZEM9qsWbO49dZbadmyJbfffjvDhw+v9hp33HEHf/vb32jTpg233XbbETWUoqIiFi5cyPTp0w//bu66665q1/2uS6Euq2pmzYEy4Avu/nFQ1hZYAHwe2AgMc/ePgn03Ad8HDgIT3H1RTdfv0qWL33TTTTXGoGVW60YclyiNY8zZomVVG4a6XFY11A5ud98DtK1Utp3E6Kiqjp8GTAszJhERyVwsJxLUehbZl9xuKiINTyyThbuXACVdunQZk+pYrcldN9q1axd1CCISoVjODVVbmniw9h566CEeeuihqMMQkYg0qGQhtadkUf+EObhFolfX36+ShUgDlJeXx/bt25Uw6il3Z/v27eTl5dXZNWPZZ6EObpFjU1hYSHl5OR9++GHUoUhI8vLyKCwsrLPrxTJZZNLBLSJHa9q0Kd26dYs6DIkRNUOJiEhKsaxZSPZlYwpkEcldShaSloqplkWkYYplM5SZDTSz2XpuIntmzZrFrFmzog5DRCISy2Th7iXuPjbTFbKk9hYsWMCCBQuiDkNEIhLLZCEiItmlZCEiIikpWYiISEqxTBbq4BYRya5YDp3VE9zZp9XmRBq2UGsWZtbazJ40s9Vm9o6ZnWNmBWb2vJmtDV7bJB0/2czWmdkaM7sozNhERCR9YTdD3Qf80d2/CJwOvANMApa4e3dgSbCNmZ0KjABOA/oBs8ysccjxSZruvvtu7r777qjDEJGIhJYszKwV8E3gPwDcfb+77wAGA3OCw+YAQ4L3g4H57r7P3d8D1gG9wopPMvPss8/y7LPPRh2GiEQkzJrFF4APgf8ys/82swfNrAXQwd03AwSvJwTHdwbKks4vD8qOYGZjzWylma3cvXt3iOGLiEiFMJNFE+BM4H53PwP4hKDJqRpWRdlRK7O4+2x3L3b34vz8/LqJVEREahRmsigHyt391WD7SRLJY4uZdQIIXrcmHV+UdH4hsCnE+EREJE2hDZ119w/MrMzMTnH3NcD5wKrgZzQwI3hdGJzyDDDPzO4FTgS6AyvqOq7Zs2cfft+sWTNGjRpV1x9RL2keLpGGLeznLH4EPGpmnwP+DlxFojazwMyuBjYCwwDc/W0zW0AimRwExrv7oTCD00N96Vu0aFHUIYhIhEJNFu5eChRXsev8ao6fBkxLdV2twS0ikl2xnO5DU5Rn3x133MEdd9wRdRgiEpFYJgvJviVLlrBkyZKowxCRiChZiIhISkoWIiKSUixnna1LGkorIpJaLJNFWKOhNJS2em3bto06BBGJUCyThdazyL6nnnoq6hBEJELqsxARkZSULCQtkydPZvLkyVGHISIRiWUzlGTf8uXLow5BRCIUy5qFmQ00s9nqkBYRyY5YJgtN9yEikl2xTBYiIpJd6rOQtBQWFkYdgohESMlC0vLII49EHYKIRCiWySLM9Sw0/YeIyNFi2WeRrQ5ujbb6pwkTJjBhwoSowxCRiKSVLMzsy7W5uJmtN7M3zazUzFYGZQVm9ryZrQ1e2yQdP9nM1pnZGjO7qDafKeEoLS2ltLQ06jBEJCLp1iz+v5mtMLNrzax1hp/xLXfv4e4Vy6tOApa4e3dgSbCNmZ0KjABOA/oBs8yscYafJSIiIUgrWbj714HvAUXASjObZ2Z9a/mZg4E5wfs5wJCk8vnuvs/d3wPWAb1q+RkiIlKH0u6zcPe1wM3AjcB5wC/NbLWZXVrTacBzZva6mY0Nyjq4++bgmpuBE4LyzkBZ0rnlQdkRzGysma00s5W7d+9ON3wRETkGaY2GMrOvAlcB3wGeBwa6+1/N7ERgOfDbak491903mdkJwPNmtrqmj6mizI8qcJ8NzAbo0qXLUfslHCeffHLUIYhIhNIdOvsr4DfAFHc/PEQoSAQ3V3eSu28KXrea2e9INCttMbNO7r7ZzDoBW4PDy0k0c1UoBDalfysSpuQhxSLS8KSbLAYAn7r7IQAzawTkufsed59b1Qlm1gJo5O67gvcXArcDzwCjgRnB68LglGeAeWZ2L3Ai0B1YUbvbqjsV/0nqmQsRacjSTRaLgQuAik6C5sBzQO8azukA/M7MKj5nnrv/0cxeAxaY2dXARmAYgLu/bWYLgFXAQWB8RXLKBQ39mYuxYxNdTqphiDRM6SaLPHc/3Jvs7rvNrHlNJ7j734HTqyjfDpxfzTnTgGmpggnzCW6p2rvvvht1CCISoXRHQ31iZmdWbJhZTyCyP7U1RbmISHalW7OYADxhZhUdzp2Ay8IJSUREck1aycLdXzOzLwKnkBjiutrdD4QamYiI5IxMZp09C+ganHOGmeHuD4cSleScHj16RB2CiEQo3Yfy5gInAaVAxQglByJJFurgzr6ZM2dGHYKIRCjdmkUxcKq758QT0+5eApR06dJlTNSxiIg0BOmOhnoL6BhmIJLbRo4cyciRI6MOQ0Qikm7Noh2wysxWAPsqCt19UChRSc4pLy+POgQRiVC6yWJqmEGIiEhuS3fo7Itm1gXo7u6Lg6e3I1uYSB3cIiLZle6yqmOAJ4EHgqLOwNNhBZWKnuAWEcmudJuhxpOYXvxVSCyEFKxRIQ3EOeecE3UIIhKhdJPFPnffH8wgi5k1oYqFiaT++vd///eoQxCRCKU7dPZFM5sCNAvW3n4CKAkvLBERySXp1iwmAVcDbwLjgD8AD4YVVBzMnTv3iDUu6vviSN/97ncBeOqppyKORESikO5oqM9ILKv6m3DDSU8ujIaqvBhSfV8cafv27VGHICIRSnduqPeooo/C3b+QxrmNgZXA++5+sZkVAI+TmJRwPTDc3f8RHDuZRA3mEHC9u/+pqmtGNd2HVokTkYYqk7mhKuSRWAq1IM1zbwDeAVoF25OAJe4+w8wmBds3mtmpwAjgNBJrcC82s5NzaWlVEZGGKq0ObnffnvTzvrvPBL6d6jwzKwS+w5H9G4OBOcH7OcCQpPL57r7P3d8D1pEYrisiIhFLtxnqzKTNRiRqGi3TOHUm8NNKx3Zw980A7r456XmNzsArSceVB2WSA84/v8pl00WkgUi3GeqepPcHCfoaajrBzC4Gtrr762bWJ43PsCrKjuonMbOxwFiAgoJ0W8LkWN1yyy1RhyAiEUp3NNS3anHtc4FBZjaARD9HKzN7BNhiZp2CWkUnYGtwfDlQlHR+IbCJStx9NjAboEuXLnowUEQkC9JthvpxTfvd/d4qyiYDk4Pz+wA/cfeRZnYXMBqYEbwuDE55BphnZveS6ODuDqxI7zYkbP379wdg0aJFEUciIlHIZDTUWST+QwcYCLwElNXiM2cAC8zsamAjiZFVuPvbZrYAWEWiqWu8RkLljvr+HImI1CyTxY/OdPddAGY2FXjC3X+QzsnuvhRYGrzfDlTZW+ru04BpacYkIiJZku7cUJ8H9idt7yfxUF0kzGygmc3WX7siItmRbs1iLrDCzH5HYoTSJcDDoUWVQlRPcKdS0xPe9X3uKBGp39IdDTXNzBYB3wiKrnL3/w4vrPon7rWgiy++OOoQRCRC6dYsAJoDO939v8ysvZl1C560lgbgJz/5SdQhiEiE0l1W9TbgRoKhsEBT4JGwghIRkdySbgf3JcAg4BMAd99EetN9hEId3NnXp08f+vTpE3UYIhKRdJuh9ru7m5kDmFmLEGNKKVc7uFNJ7gBXh7eIxEm6NYsFZvYA0NrMxgCLyZGFkOJKtSIRiZOUNQszMxKLFX0R2AmcAtzq7s+HHJuIiOSIlMkiaH562t17AkoQIiINULp9Fq+Y2Vnu/lqo0aQpF9bgbmiGD69xRnoRqefSTRbfAn5oZutJjIgyEpWOr4YVWE3i2sEdZ9dee23UIYhIhGpMFmb2eXffCPTPUjySo/bs2QNA8+bNI45ERKKQqmbxNInZZjeY2VPu/t1sBCW5Z8CAAQAsXbo02kBEJBKphs4mL3X6hTADERGR3JWqZuHVvI9Ufeng1kN6IhIXqWoWp5vZTjPbBXw1eL/TzHaZ2c5sBFgVdy9x97HNmjWLKoQ6p4f0RCSX1Zgs3L2xu7dy95bu3iR4X7HdqqZzzSzPzFaY2Rtm9raZ/VtQXmBmz5vZ2uC1TdI5k81snZmtMbOL6uYWRUTkWGUyRXmm9gHfdvfdZtYU+HOwJsalwBJ3n2Fmk4BJwI1mdiowAjgNOBFYbGYnax3u3HDllVdGHYKIRCi0ZOHuDuwONpsGPw4MBvoE5XNIrM19Y1A+3933Ae+Z2TqgF7A8rBglfUoWIg1buhMJ1oqZNTazUmAr8Ly7vwp0cPfNAMHrCcHhnYGypNPLgzLJAdu2bWPbtm1RhyEiEQmzGYqgCamHmbUGfmdmX67hcKui7KgRWGY2FhgLUFBQUCdxSmpDhw4F9JyFSEMVas2igrvvINHc1A/YYmadAILXrcFh5UBR0mmFwKYqrjXb3YvdvTg/Pz/UuEVEJCG0ZBGs0906eN8MuABYDTwDjA4OGw0sDN4/A4wws+PMrBvQHVgRVnwiIpK+MJuhOgFzzKwxiaS0wN2fNbPlJBZTuhrYCAwDcPe3zWwBsAo4CIzXSCgRkdwQ5mio/wHOqKJ8O3B+NedMA6aFFZOIiNROqB3cYakv033EyTXXXBN1CCISoVgmC61nkX2XXXZZ1CGISISyMhpK4q+srIyysrLUB4pIvRTLmoVkX8WMuHrOQqRhUs1CRERSimXNor52cCevb5FMa12ISNRiWbOoj+tZ1ERrXYhI1GKZLEREJLti2Qwl2Tdx4sSoQxCRCClZSFoGDhwYdQgiEqFYNkOZ2UAzm622/OxZs2YNa9asiToMEYlILJNFQ+vgzgXjxo1j3LhxUYchIhGJZbIQEZHsUrIQEZGUlCxERCSlWCYLdXCLiGRXLIfOaory7Lv55pujDkFEIhRasjCzIuBhoCPwGTDb3e8zswLgcaArsB4Y7u7/CM6ZDFwNHAKud/c/hRVffTJ37twjpgQJYy6pCy64oE6vJyLxEmYz1EFgort/CTgbGG9mpwKTgCXu3h1YEmwT7BsBnAb0A2YF63dLCpWb48JonistLaW0tLTOrysi8RDmGtybgc3B+11m9g7QGRgM9AkOmwMsBW4Myue7+z7gPTNbB/QClocVY1xVrklkw4QJEwCtZyHSUGWlz8LMugJnAK8CHYJEgrtvNrMTgsM6A68knVYelFW+1lhgLEBBQUF4QeeY6qYvFxHJhtBHQ5lZPvAUMMHdd9Z0aBVlflSB+2x3L3b34vz8/LoKU0REahBqsjCzpiQSxaPu/tugeIuZdQr2dwK2BuXlQFHS6YXApjDjExGR9ISWLMzMgP8A3nH3e5N2PQOMDt6PBhYmlY8ws+PMrBvQHVgRVnwiIpK+MPsszgVGAW+aWcUwminADGCBmV0NbASGAbj722a2AFhFYiTVeHc/FGJ8koHp06dHHYKIRCjM0VB/pup+CIDzqzlnGjAtrJik9nr37h11CCISIU33IWlZtmwZy5YtizoMEYmIpvuop5KH2tbFE91TpkwB9JyFSEMVy5qFZEY1MBE5VrGsWUjm6rqmISINi2oWDZBqGiKSqVgmC3Vwi4hkVyybodTBfewybZaaOXNm2CGJSA6LZbKQupVODa1Hjx5ZiEREclUsm6Ek+xYvXszixYujDkNEIqKahaTlzjvvBLRinkhDFcuahTq4RUSyK5bJwt1L3H1ss2bNog5FRKRBUDOUAP8cHVV5ZFTFEq6bN28+vK0H+kQanljWLCQ8lZv2Um2LSMOgmoUcpar1vr/3ve9FEImI5IpYJgszGwgMbN++fdShNBgdO3aMOgQRiVAsm6HUwZ19b7zxBm+88UbUYYhIRMJcg/s/zWyrmb2VVFZgZs+b2drgtU3Svslmts7M1pjZRWHFJbWjh/JEGrYwaxYPAf0qlU0Clrh7d2BJsC7XEzEAAAjZSURBVI2ZnQqMAE4LzpllZo1DjE1ERDIQWrJw95eAjyoVDwbmBO/nAEOSyue7+z53fw9YB/QKKzYREclMtju4O7j7ZgB332xmJwTlnYFXko4rD8qOYmZjgbEABQUFIYYq1alqtFQFLawkUj/lSge3VVHmVR3o7rPdvdjdi/Pz80MOSzKl5zBE6qds1yy2mFmnoFbRCdgalJcDRUnHFQKbshyb1OCqq65K+1gt4SpS/2Q7WTwDjAZmBK8Lk8rnmdm9wIlAd2BFlmOTGtS2ya9yTaNi+hBQIhGJk9CShZk9BvQB2plZOXAbiSSxwMyuBjYCwwDc/W0zWwCsAg4C4939UFixSeZee+01AM4666xjuk5y8vj0009VCxGJidCShbtfXs2u86s5fhowLZ1r6wnu7HvppZeAY08WNVF/h0juypUO7ozoCW4RkeyKZbIQEZHsiuVEghIvNT2XISLxoJqFiIikFMuahTq4s2/cuHFRhyAiEYplsnD3EqCkS5cuY6KOpaHI1tPyGkorkptimSwk+5YtWwZA7969s/aZNT3QB0omItmkZCFpWb58OZDdZAE1d47roT6R7FGykHqjcvLYsmULHTp0iDAikfpDyULqrc8+++yIbc1LJVJ7sUwWGg0l6dq8eXOVTVmaWkQkM7FMFhoNJXVB/R0i6YtlspDs+9GPfhR1CKHSyCuRmilZSFo+97nPRR1CVlVOHmq2koZOyULSsnTpUgD69OkTaRxhOpY5rCrXRJKpViL1QSyThTq4s+/1118H6neySKW2ySTV8yBq8pI4iGWyUAe3xFnl5JHp/sqUXCQbci5ZmFk/4D6gMfCgu8+IOCSRnFZTzSXTWotqOVKdnEoWZtYY+DXQFygHXjOzZ9x9VbSRicRHTTWTVKO+Ul0rVRNaskyOTed8iVZOJQugF7DO3f8OYGbzgcGAkoVIHTmWjvxMmsgybU6r6/OTZVqjyuRax1obS/XZuZI0zd2jjuEwMxsK9HP3HwTbo4Cvuft1SceMBcYGm18G3krz8scDH2cQTibHh3VsO2BbxDGEee1cuL8wfxdxu79Mfxf1+f7CurdMjw/r2FPcvWWaxya4e878AMNI9FNUbI8C/l8Nx6/M4NqzM4wl7eNDPDaU+8uh30Xk9xfy7yJW91eL30W9vT/933L0T64tq1oOFCVtFwKb6ujaJSEeH9axmQgzhvp8f2H+LsK6bi4cm6lciDkXvrtMj8+F+wNyrxmqCfAucD7wPvAacIW7v13N8SvdvTiLIWaV7i/edH/xVZ/vDWp3fznVwe3uB83sOuBPJIbO/md1iSJQN71fuUv3F2+6v/iqz/cGtbi/nKpZiIhIbsq1PgsREclBShYiIpJSbJOFmfUzszVmts7MJkUdT10zs/Vm9qaZlZrZyqjjOVZm9p9mttXM3koqKzCz581sbfDaJsoYa6uae5tqZu8H31+pmQ2IMsZjYWZFZvaCmb1jZm+b2Q1BeX35/qq7v9h/h2aWZ2YrzOyN4N7+LSjP+LuLZZ9FMC3IuyRNCwJc7vVoWhAzWw8Uu3u6DwblNDP7JrAbeNjdvxyU/Rz4yN1nBAm/jbvfGGWctVHNvU0Fdrv73VHGVhfMrBPQyd3/amYtgdeBIcCV1I/vr7r7G07Mv0MzM6CFu+82s6bAn4EbgEvJ8LuLa83i8LQg7r4fqJgWRHKUu78EfFSpeDAwJ3g/h8Q/0Nip5t7qDXff7O5/Dd7vAt4BOlN/vr/q7i/2PGF3sNk0+HFq8d3FNVl0BsqStsupJ19uEgeeM7PXgylO6qMO7r4ZEv9ggRMijqeuXWdm/xM0U8WyiaYyM+sKnAG8Sj38/irdH9SD79DMGptZKbAVeN7da/XdxTVZWBVl8WtPq9m57n4m0B8YHzR1SHzcD5wE9AA2A/dEG86xM7N84ClggrvvjDqeulbF/dWL79DdD7l7DxIzYvQysy/X5jpxTRZhTguSE9x9U/C6Ffgdiaa3+mZL0F5c0W68NeJ46oy7bwn+kX4G/IaYf39Be/dTwKPu/tuguN58f1XdX337Dt19B7AU6Ectvru4JovXgO5m1s3MPgeMAJ6JOKY6Y2Ytgo42zKwFcCHpz64bJ88Ao4P3o4GFEcZSpyr+IQYuIcbfX9BJ+h/AO+5+b9KuevH9VXd/9eE7NLP2ZtY6eN8MuABYTS2+u1iOhgIIhrHN5J/TgkyLOKQ6Y2ZfIFGbgMSULPPifn9m9hjQh8TUz1uA24CngQXA54GNwDB3j11HcTX31odE84UD64FxFW3EcWNmXwdeBt4EPguKp5Bo168P319193c5Mf8OzeyrJDqwG5OoHCxw99vNrC0ZfnexTRYiIpI9cW2GEhGRLFKyEBGRlJQsREQkJSULERFJSclCRERSUrIQEZGUlCxE6piZ3Wdmzc3swUrlPc1sXFRxiRwLJQuROmRmBSQe4nLgi8GaCPPNzNz9deAb0UYoUjtKFiK1YGZfMbMNZnZNpV09gFUkZi590t2nAh8Dxwf795pZh+xFKlI3lCxEasHd3yQxJ9n/qbSrANhBYtK5/wnKmgeTuAH8A2iVlSBF6lCTqAMQibGtwGmVyt4lMatnd6CjmQ0nMUldhc4k5uIRiRUlC5HamwEcZ2Zd3H1DUPYmcK27j6l8cLBewk5335fNIEXqgpqhRGrBzPoBLYDfk1S78MTMnI+aWfMqTjsRuCs7EYrULc06K5IhM8sDVgCDgKuAT9z959FGJRIu1SxEMncz8LC7ryfR7JRymUozG2RmGjYrsaU+C5EMmNkpQF/g3KDoTRIL5WBmVwLfAj4lsWZzUxKJZDiJUVKfBcecB/wd+Czui1pJw6FmKJE6EiSC/e4+z8yWuPv5ZjYFWAScDmwjsZrep+7+uJk95u6XRxiySNrUDCVSt3YGrx8Gr/uB4yod80nwalmJSKQOKFmIiEhKaoYSEZGUVLMQEZGUlCxERCQlJQsREUlJyUJERFJSshARkZSULEREJCUlCxERSUnJQkREUlKyEBGRlP4XODUZuMD+9VQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.hist(Lmin_b,bins=int(np.sqrt(B)),color = \"0.6\", ec=\"0.6\");\n",
    "axes = plt.gca()\n",
    "y_min, y_max = axes.get_ylim()\n",
    "plt.plot([Lmin,Lmin],[y_min,y_max],'--k',label='Test statistic')\n",
    "plt.plot([Lmin_c,Lmin_c],[y_min,y_max],'-k',label='Critical value')\n",
    "plt.ylim([y_min,y_max])\n",
    "x_min, x_max = axes.get_xlim()\n",
    "plt.xlim([0,30])\n",
    "plt.xlabel(r'$\\lambda_{\\rm{min}}^{(b)}$')\n",
    "plt.ylabel('Frequency')\n",
    "plt.minorticks_on()\n",
    "plt.rcParams.update({'font.size': 12})\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test direction of interest (in this case GAD) against common mean of two groups\n",
    "\n",
    "#### N.b. should only be calculated if null hypothesis of common mean direction cannot be rejected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "95% confidence region threshold T = 20.43\n",
      "Tested direction T = 156.01\n"
     ]
    }
   ],
   "source": [
    "T_c = np.quantile(T_b,1-alpha) #boundary value for confidence region on common mean direction\n",
    "m = np.transpose(pmag.dir2cart((0.0,ipmag.inc_from_lat(40.15)))) #GAD direction in Cartesian coordinates\n",
    "T_m = np.matmul(np.matmul(np.transpose(m),Ahat),m) #T-value for m-direction\n",
    "print(\"95% confidence region threshold T = {:.2f}\".format(T_c))\n",
    "print(\"Tested direction T = {:.2f}\".format(np.squeeze(np.asarray(T_m))))\n",
    "#If T_m <= T_c then direction of interest is consistent with the common mean direction at alpha significance level"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot histogram of bootstrap values and compare to test statistic (see Figure 6(a) in manuscript)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAETCAYAAAD6R0vDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de3wV1bn/8c8DoXLZhFsQFWOwqODBCi14bVV6RFsoUARUqnjgtBUVOZUj1p+iIFWgRUWxtSiptgh44VIvJwi94AGlXktbUGlRwSNyKVRQIgHC9fn9MRPcbBKSHXZm7+x836/XvJi91sysJ4TkYc2aWcvcHRERkVSql+4AREQk+yi5iIhIyim5iIhIyim5iIhIyim5iIhIyim5iIhIyuWkO4CakpeX5+3atUt3GJV67733AOjQoUOaIxGRbFeV3zd/+ctftrh766NtK+uSi5n1AfqccsopLFu2LN3hVKp79+4ALFmyJK1xiEj2q8rvGzNbm4q2si65uHsRUNStW7dr0x2LiEgmufPOOyNrK+uSi4iIlK9Hjx6RtaUBfRGROmL58uUsX748krbUcxERqSNGjhwJRDPGm3U9FzPrY2aFxcXF6Q5FRKTOiiy5mFk7M1tgZp+Z2SYze9jMcsK6i81slZntNLPFZlYQd56Z2SQz2xpu95qZVdSOuxe5+7BmzZpF8WWJiEg5ouy5TAX+BRwPdAEuAoabWR7wLDAGaAksA2bHnTcM6Ad0Bs4EegPXRRe2iIgkK8rkcjIwx91L3X0T8DugE9AfWOnuc929FBgHdDazjuF5Q4DJ7r7e3TcAk4GhEcYtIiJJijK5PAQMMrPGZtYW6MkXCWZF2UHuvgNYE5aTWB/ud6ISW7duZebMmSkKXUSk9ps4cSITJ06MpK0onxZ7GbgW+ByoDzwBPE9wm+uThGOLgabhfiz8HF8XMzPzhGU0zWwYwW00WrZsya5du1L9NYiI1Frnn39+ZG1F0nMxs3rA7wnGVpoAeUALYBJQAuQmnJILbA/3E+tzgZLExALg7oXu3s3du8VisdR+ESIitdxrr73Ga6+9FklbUd0WawnkAw+7+2533wr8BugFrCQYrAfAzJoA7cNyEuvD/ZWIiEhSRo8ezejRoyNpK5Lk4u5bgP8DbjCzHDNrTjBQvwJ4DjjDzAaYWUNgLPC2u68KT58B3Gxmbc3sBGAUMD2KuEVEpHqiHNDvD3ybYHxlNbAP+G93/wQYAEwAPgPOAQbFnTcNKALeAd4FXgzLREQkQ0U2oO/uy4HuFdQtAjpWUOfAreFWqbIp91u3PurlCEREpJqybvqXsjf0GzVqlO5QRETqLE1cKSJSR0yZMiWytpRcRETqiC5dukTWVtbdFhMRkfItWrSIRYsWRdJW1vVcNKAvIlK+8ePHA9GsSJl1PRcN6IuIpF/WJZdEmrxSRCR6WZ9cNHmliEj0sj65iIhI9DSgLyJSR0ybFt3MWVnXc9GAvohI+Tp06ECHDh0iaSvrkouIiJSvqKiIoqKiSNrKuttiIiJSvsmTJwPQp0+fGm9LPRcREUk5JRcREUm5rEsuZtbHzAr1fouISPpEklzMrCRh229mv4irv9jMVpnZTjNbbGYFcXVmZpPMbGu43WtmVlFbelpMRCT9IhnQd/dY2b6ZNQE2A3PDz3nAs8APCZYzvgeYDZwbnjIM6Ad0Bhz4I/Ah8GgUsYuIZIsop8NKx22xgcC/gKXh5/7ASnef6+6lwDigs5mVLXs8BJjs7uvdfQMwGRgabcgiIrVffn4++fn5kbSVjuQyBJjh7h5+7gSsKKt09x3AmrD8sPpwvxMiIpKU2bNnM3v27EjaivQ9FzM7CbgI+EFccQz4JOHQYqBpXH1xQl3MzCwuQZVdfxjBbTRatmyZwshFRGq/Rx55BIArr7yyxtuKuufyH8Cf3P3/4spKgNyE43KB7RXU5wIliYkFwN0L3b2bu3eLxWKJ1SIiEpF0JJcnEspWEgzWAwcH/NuH5YfVh/srERGRjBVZcjGz84G2hE+JxXkOOMPMBphZQ2As8La7rwrrZwA3m1lbMzsBGAVMjyhsERGphih7LkOAZ919e3yhu38CDAAmAJ8B5wCD4g6ZRvCI8jvAu8CLYZmIiGSoyAb03f26I9QtAjpWUOfAreFWKa3nIiJSvnnz5kXWVtZN/6I39EVEypeXl0deXl4kbWVdchERkfJNnz6d6dOnR9KWkouISB2h5CIiIrVa1iUXTbkvIpJ+WZdcNKAvIpJ+WZdcREQk/SKduFJERNJnwYIFkbWl5CIiUkc0btw4sray7raYBvRFRMo3depUpk6dGklbWZdcNKAvIlK+OXPmMGfOnEjayrrkIiIi6afkIiIiKafkIiIiKZd1yUUD+iIi6Zd1jyK7exFQVFBQcG26YxERySRLliyJrK1Iey5mNsjM/mFmO8xsjZldEJZfbGarzGynmS02s4K4c8zMJpnZ1nC718wsyrhFRCQ5kSUXM7sEmAT8J9AUuBD40MzygGeBMUBLYBkwO+7UYUA/oDNwJtAbqHBVSxERKd/999/P/fffH0lbUfZcfgLc7e5vuPsBd9/g7huA/sBKd5/r7qXAOKCzmZUtezwEmOzu68PjJwNDI4xbRCQrzJ8/n/nz50fSViTJxczqA92A1ma22szWm9nDZtYI6ASsKDvW3XcAa8JyEuvD/U6IiEjGiqrn0gZoAAwELgC6AF8F7gRiQHHC8cUEt84op74YiJU37mJmw8xsmZktKykpSe1XICIiVRZVcil7LvgX7v5Pd98CPAD0AkqA3ITjc4Ht4X5ifS5Q4u6e2Ii7F7p7N3fvFovFUvoFiIhI1UWSXNz9M2A9cFhCAFYSDNYDYGZNgPZh+WH14f5KkjBz5sxkDhcRyUqNGjUiqnkXo3zP5TfAf5nZ74C9wEhgPvAccJ+ZDQBeBMYCb7v7qvC8GcDNZraAIDmNAn6RTMN6oVJEBBYuXBhZW1Eml3uAPOB9oBSYA0xw99IwsTwMzALeBAbFnTcN+DLwTvj5sbBMREQyVGTJxd33AsPDLbFuEdDxsJOCOgduDbdKmVkfoE/r1q2rH6yISBa65557ABgzZkyNt5V1c4tpPRcRkfK99NJLvPTSS5G0lXXJRURE0k/JRUREUk7JRUREUi7rkktF67noXRcRqetatWpFq1atImmrzqznonddRKSu++1vfxtZW1nXcxERkfRTchERqSNuv/12br/99kjayrrbYiIiUr7XX389srayrudS0YC+iIhEJ+uSi97QFxFJv6xLLiIikn5VHnMxs77AAnffV4PxiIhIDTnxxBMjayuZAf17gMfNbDYw093frKGYRESkBsyaNSuytqp8W8zdOwM9CJYs/q2ZvWdmd5pZuxqKrVo0oC8ikn5Jjbm4+wp3/zGQD9wIXA6sMbNXzOxqM0v7GI4G9EVEyjdy5EhGjhwZSVtJJwMza0+wFPEjQMNw/1fACGDeEc5bYmalZlYSbu/F1V1sZqvMbKeZLTazgrg6M7NJZrY13O41M0s2bhGRum758uUsX748kraqnFzM7EYze4NgGeI2wDXu3sHdJ7j7TOBi4NJKLjPC3WPh1iG8bh7wLDAGaAksA2bHnTMM6Ad0Bs4EegPXVTVuERGJXjID+j2BycAL7r4nsdLdd5pZ/2rE0B9Y6e5zAcxsHLDFzDq6+ypgCDDZ3deH9ZOBa4FHq9GWiIhEIJnbYgOB5+MTi5k1MLNjyj67+x8qucZPzWyLmb1qZt3Dsk7Airhr7ADWhOWH1Yf7nSiHmQ0zs2VmtqykpKSKX5aIiKRaMsnlD0DXhLKuwO+reP7/A74MtAUKgaJw/CYGFCccWww0DfcT64uBWHnjLu5e6O7d3L1bLBarYlgiInXDaaedxmmnnRZJW8ncFjuTYLwl3lsEYyGVSngv5gkz+x7QCygBchMOzwW2h/uJ9blAibt7FeMWERGgsLAwsraS6blsIxjIj9cG2FHNth0wYCVxCcrMmgDtw3IS68P9lVSDVqMUEYlGMsnlt8BTZnaGmTU2s68AM4A5lZ1oZs3N7Ftm1tDMcszsauBCgltqzwFnmNkAMyt7tPntcDCfsI2bzaytmZ0AjAKmJxH3QXqxUkTqsmHDhjFs2LBI2krmttgdBE+LvQUcA5QCvwFGV+HcBsB4oCOwH1gF9HP39wDMbADwMDCL4NbboLhzpxGM1bwTfn4sLBMRkSS8//77kbVV5eTi7qXAjWY2AsgDtlR13MPdPwHOOkL9IoLEU16dA7eGW6XMrA/Qp3Xr1lU5XEREakBSK1GaWTOgA8ETXJQ9sOXu/5vyyKrJ3YuAooKCgmvTHYuISF2VzJT7Q4FfEjy9tTOuygluW4mIiADJ9VwmAAPdfWFNBSMiIjWnS5cukbWVTHLJIXiRUkREaqEpU6ZE1lYyjyJPAu7MhGn1j0TruYiIpF8yieK/gTuB7Wb2cfxWQ7FVi9ZzEREp3+DBgxk8eHAkbSVzWyyaiEREpEasX78+sraSec/l5ZoMREREskcyi4UdY2YTzOxDMysOyy4NX6oUERE5KJkxlweBM4CrCd5tgWACyRtSHdTR0IC+iEj6JTPmchlwirvvMLMDAO6+wcza1kxo1aM39EVEynfeeedF1lYyyWVP4vFm1hrYmtKIRESkRvz0pz+NrK1kbovNJVjk62QAMzueYCbjZ2oiMBERqb2SSS6jgY8Ipr5vDnwAbAR+kvqwas7MmTO1aJiI1EkDBgxgwIABkbSVzKPIe4CRwMjwdliVp9yPUmVT7mugX0Tqqq1boxvFSOZR5C+XbUBT4OS4z1VmZqeaWamZzYoru9jMVpnZTjNbbGYFcXVmZpPMbGu43Wtlc/2XQ2/oi4ikXzID+qv5Yt37MmU9l/pJXOeXwJ/LPphZHvAs8EOgCLgHmA2cGx4yDOgHdA7b+yPwIfBoEm2KiEiEqtxzcfd67l4//LMecAJQCFxT1WuY2SBgG/BSXHF/YKW7zw1XuxwHdDazspUphwCT3X29u28gWGp5aFXbFBGR6CW1EmU8d99kZiOB94GnKjvezHKBu4GLgR/EVXUCVsRdd4eZrQnLVyXWh/udqhu3iEhddfHFF0fWVrWTS6gD0LiKx94DPO7u6xKGTGLAJwnHFhOM65TVFyfUxczMEh8oMLNhBLfRaNmyZRXDEhGpG8aMGRNZW8ksc7yUL8ZYIEgqnQh6I5Wd2wXoAXy1nOoSIDehLBfYXkF9LlBS3pNq7l5IcKuOgoKCjHuSTUSkrkim5/JYwucdwAp3/6AK53YH2gEfh72WGFDfzP6NYGB+SNmBZtYEaE8wbxnhn52Bt8LPnePqRESkinr27AnAwoU1v1p9Mu+5PHEU7RRy6Jv8txAkm7JJL+8zswHAi8BY4G13XxXWzQBuNrMFBD2nUcAvjiIWEZE6Kcr3/JK5LVbp7S8Adx9bTtlOYGfctUqAUnf/JPw8gGAqmVnAm8CguNOnAV8mmBkAgh7UtKrGLSIi0UvmttipwACCd1TWAicBZwO/BUrDY6o0zuHu4xI+LwI6VnCsA7eGW6Uqe0NfRERqXjLJxYDvuftvDxaY9Qcud/f/THlk1aQp90VE0i+ZiSt7As8nlL0A9EpdONHR5JUiUtf07t2b3r17R9JWstO/3Aj8PK5sOLAmpRFFRBNYikhdc8stt0TWVjLJ5YfAc2Z2K7ABaAvsI5i+RURE5KBkHkX+m5mdSjCh5AnAP4HX3X1vTQVXHRrQFxEpX/fu3QFYsmRJjbeVzJjLIdz9FeBL4UuPGUNT7ouIpF8y67l8hWCSyl8Bj4fFFwG/roG4IqFBfRGRmpFMz+URYKy7dwTKboW9DHwj5VFFRIP6IiI1I5nk0ongDXoIX5Z09x2A7j+JiMghknla7COgK7CsrMDMziZ4RFlERDLcFVdcEVlbySSXMcCLZvYowUD+7cD1QEa9Ca+nxUREyjd8+PDI2kpmmeP5BG/ptyYYaykA+rv7H2ootmrR02IiIuXbuXMnO3furPzAFKhSz8XM6hM8KfZv7h5d6hMRkZTp1SuYrStj3nNx9/3AfqBhzYYjIiLZIJmnxaYAc8zsIjNrb2ZfLttqKrgo6F0XEZHUq/S2mJkd5+6bCBbzAuhBMP1+GQfq10Bs1ZLsgL7edRERSb2q9FzeB3D3eu5eD/ifsv1wq1JiMbNZZvZPM/vczN43sx/G1V1sZqvMbKeZLTazgrg6M7NJZrY13O41Myu/FQ3oi4hkgqoM6Cf+Ir+omm39FPiBu+82s47AEjP7G8Gqls8SzLpcBNwDzCaYIBNgGNAP6EzQS/oj8CHwaDXjEBGpk4YOHRpZW1VJLolLF1fYazjiRdxXJlzTgfYEL2audPe5AGY2DthiZh3dfRUwBJjs7uvD+skE79YouYiIJCHTkkuOmX2TL5JK/YTPuPv/VqUxM5sKDCWYMuZvwAJgArAi7lo7zGwNwXQzq8I/V8RdZkVYJiIiSdiyZQsAeXl5Nd5WVZLLvzh05uOtCZ8dqNITY+4+3Mz+CzgP6A7sBmLAJwmHFgNNw/1Y+Dm+LmZm5u6H9KrMbBjBbTRatmxZlZBEROqMgQMHAtG851JpcnH3dqlsMHxn5k9mNhi4ASgBchMOywW2h/uJ9blASWJiCa9dCBQCFBQUHFYvIiLRqPZiYSmQQzDmspJgsB6AcPGxsnIS68P9+PEbERHJMJEkFzM71swGmVnMzOqb2beA7wH/CzwHnGFmA8ysITAWeDsczAeYAdxsZm3N7ARgFDA9irhFRKR6kpkV+Wg4wS2wRwkS2lpgpLu/AGBmAwhe0pwFvAkMijt3GsGYzjvh58fCMhERyVCRJBd3/4QjvB/j7ouAjhXUOXBruFVKU+6LiJTvhhtuiKytqHoukXH3IqCooKAgo9aZERFJtyuvvDKyttI5oC8iIhFat24d69ati6StrOu5iIhI+a655hogg9ZzERERSUbWJRcz62NmhclMpa81XUREUivrkkt1ptzXmi4iIqmVdclFRETSTwP6IiJ1xKhRoyJrS8lFRKSO6NOnT2RtZd1tseoM6IuI1AXvvfce7733XiRtZV1yqc6AvohIXXDddddx3XXXRdJW1iUXERFJPyUXERFJOSUXERFJOSUXERFJuax7FFnruYiIlO/OO++MrK2oljk+xsweN7O1ZrbdzP5mZj3j6i82s1VmttPMFptZQVydmdkkM9sabveamVXUVnWfFps5c6bmGBORrNajRw969OgRSVtR3RbLAdYRrEbZDBgDzDGzdmaWBzwblrUElgGz484dBvQDOgNnAr2BlD9Lt2vXLs0xJiJZbfny5SxfvjyStqJa5ngHMC6uaL6Z/R/QFWgFrHT3uQBmNg7YYmYd3X0VMASY7O7rw/rJwLXAo1HELiKSLUaOHAlk8XouZtYGOA1YCXQCVpTVhYloTVhOYn243wkREclYkScXM2sAPAk8EfZMYkBxwmHFQNNwP7G+GIiVN+5iZsPMbJmZLSspKUl98CIiUiWRJhczqwfMBPYAI8LiEiA34dBcYHsF9blAibt74vXdvdDdu7l7t1gsltLYRUSk6iJLLmFP43GgDTDA3feGVSsJBuvLjmsCtA/LD6sP91ciIiIZK8r3XB4BTgd6uHv8Y1nPAfeZ2QDgRWAs8HZ4ywxgBnCzmS0AHBgF/CK6sEVEssPEiRMjayuS5BK+t3IdsBvYFDdccp27PxkmloeBWcCbwKC406cBXwbeCT8/FpZJFjhw4ABbtmxh27Zt7N+/P93hSA2oX78+zZs3Jy8vj3r1NClIOp1//vmRtRXVo8hrgSO9+LgI6FhBnQO3hlul9IZ+7bJ+/XrMjHbt2tGgQQOO8H6s1ELuzt69e9m8eTPr16/npJNOSndIddprr70GRJNksu6/EVrPpXbZsWMHbdu25Utf+pISSxYyM770pS/Rtm1bduzYke5w6rzRo0czevToSNrKuuRytDQFTPR0qyT76Xtc9+g7nmDXrl1KMCIiR0nJpRyaY0yyUc+ePXniiSeSPm/p0qV06NChBiKSbJZ1ycXM+phZoRKEHK1YLHZwq1evHo0aNTr4+cknn0z6et27d+exxx6rVizJnjtu3DgGDx58SNnChQsZMmRIpeeaGatXrz74+YILLuC9996rerAiZOF6Lu5eBBQVFBRcm+5YpHaLn0KoXbt2PPbYY5FNVy5SE6ZMmRJZW1nXcxGpaQcOHOBnP/sZ7du3p1WrVlxxxRV8+umnAJSWljJ48GBatWpF8+bNOeuss9i8eTN33HEHS5cuZcSIEcRiMUaMGHHYdZM996abbiI/P5/c3Fy6du3K0qVLAfjd737HxIkTmT17NrFYjM6dgwku4ns/q1ev5qKLLqJZs2bk5eVx5ZVXAnDhhRcC0LlzZ2KxGLNnz2bJkiWceOKJB+Nct24d/fv3p3Xr1rRq1arcr0UyU5cuXejSpUskbWVdz0Vqv+7dux9WdsUVVzB8+HB27txJr169DqsfOnQoQ4cOZcuWLQwcOPCw+htuuIErr7ySdevWkZ+ff1Tx/fznP+f555/n5ZdfpnXr1vzoRz/ixhtv5Omnn+aJJ56guLiYdevWccwxx7B8+XIaNWrEhAkTePXVVxk8eDA//OEPy71usueeddZZjB07lmbNmvHQQw9x+eWX89FHH/Htb3+b0aNHs3r1ambNmlVuW2PGjOHSSy9l8eLF7Nmzh2XLlgHwyiuvYGasWLGCU045BTh0evb9+/fTu3dv/v3f/52ZM2dSv379g+dK5lu0aBFAJD1w9VxEkjRt2jQmTJjAiSeeyDHHHMO4ceOYN28e+/bto0GDBmzdupXVq1dTv359unbtSm5u4rys5Uv23LJeTk5ODqNGjWL37t1VHhtp0KABa9euZePGjTRs2JBvfOMbVTrvrbfeYuPGjdx33300adIkqXMl/caPH8/48eMjaSvrei56Q7/2O9JCRo0bNz5ifV5e3hHrj7bXArB27Vouu+yyQ97dqF+/Pps3b+aaa65h3bp1DBo0iG3btjF48GAmTJhAgwYNKr1usudOnjyZxx57jI0bN2JmfP7552zZsqVKX8O9997LmDFjOPvss2nRogWjRo3i+9//fqXnrVu3joKCAnJysu5Xh6RY1vVc9Ia+1LT8/HwWLlzItm3bDm6lpaW0bduWBg0acNddd/H3v/+d1157jfnz5zNjxgyASmcgSObcpUuXMmnSJObMmcNnn33Gtm3baNasGWUrUVTW1nHHHcevfvUrNm7cyLRp0xg+fPghT4gd6Wv/+OOP2bdvX6XHSt2WdcklVWbOnKl7yVKu66+/njvuuIO1a9cC8Mknn/DCCy8AsHjxYt555x32799Pbm4uDRo0oH79+gC0adOGDz/8sMLrJnPu9u3bycnJoXXr1uzbt4+7776bzz///GB9mzZt+Oijjzhw4EC5bc2dO5f169cD0KJFC8ysSnGeffbZHH/88dx2223s2LGD0tJSXn311Sr9vUndouRSgV27dvHXv/413WFIBrrpppvo27cvl156KU2bNuXcc8/lzTffBGDTpk0MHDiQ3NxcTj/9dC666KKD75vcdNNNzJs3jxYtWvCjH/3osOsmc+63vvUtevbsyWmnnUZBQQENGzY85Jbf5ZdfDkCrVq342te+dlhbf/7znznnnHOIxWL07duXhx56iJNPPhkI3pEZMmQIzZs3Z86cOYecV79+fYqKili9ejUnnXQSJ554IrNnz07B36pkGytnQcesUFBQ4HfcccdRX2fYsGEpiKZiZU9GHWmcIJv94x//4PTTT093GBIBfa/Tr+yBjyPNuGBmf3H3bkfbVtaNymlAX0SkfFFO45N1t8U0oC8iUr6ioiKKiooiaSuy5GJmI8xsmZntNrPpCXUXm9kqM9tpZovDlSvL6szMJpnZ1nC717Twh4hI0iZPnszkyZMjaSvKnstGYDzw6/hCM8sDngXGAC2BZUD8COEwoB/QGTgT6E2wZLKIiGSoyJKLuz/r7s8DWxOq+gMr3X2uu5cC44DOZla27PEQYLK7r3f3DcBkYGhEYYuISDVkwphLJ2BF2Qd33wGsCcsPqw/3OxGRwsJCLR4mIpKkTEguMaA4oawYaFpBfTEQK2/cxcyGheM6y+KnSz9aWhtGRCQ5mfAocgmQODtfLrC9gvpcoMTLeUHH3QuBQgjec0l9qCIitVeUd2EyoeeykmCwHgAzawK0D8sPqw/3VxIx3R6Tqpo4cWKF0+oDPPnkk1x66aVH3U7iipGp8NFHH2FmmjssS+Xn56dk8taqiPJR5BwzawjUB+qbWUMzywGeA84wswFh/VjgbXdfFZ46A7jZzNqa2QnAKGB6VHHH27VrlxJMHfTUU0/RrVs3YrEYxx9/PD179uRPf/pThcePHj364KJc5f2yvvrqq/nDH/5Q43GLJJo9e3Zk0/VE2XO5E9gF3AYMDvfvdPdPgAHABOAz4BxgUNx504Ai4B3gXeDFsCwtNP5StzzwwAOMHDmS0aNHs3nzZj7++GOGDx9+cKLKRPofv2SyRx55hEceeSSStqJ8FHmcu1vCNi6sW+TuHd29kbt3d/eP4s5zd7/V3VuG263ljbeUMbM+ZlaoJCBHq7i4mLFjx/LLX/6S/v3706RJExo0aECfPn247777gGCSx4EDBzJ48GByc3OZPn0648aNOzjhZNmywc2bNycWi/H6668zffr0QxbYWrlyJZdccgktW7akTZs2TJw4EQgW5jrvvPNo3rw5xx9/PCNGjGDPnj2Vxv3MM8/QrduhU0M9+OCD9O3bF4AXX3yRr371q+Tm5pKfn8+4ceMqvFa7du0Orl5Y9vWWfW0Ab7zxBueffz7Nmzenc+fOdXaOPDlcJgzop5S7FwFFBQUF16Y7FkneyJEjWb58eY220aVLF6ZMmVLpca+//jqlpaVcdtllRzzuhRdeYO7cucyYMYPdu3czadKkg3WvvPIKJ598Mtu2bTu4wFb8apHbt2+nR48e3HLLLRQVFbF3717+/ve/A8EMxA8++CDdunVj/fr19OzZk6lTpzJy5MgjxtO3b1+uvfZaPvjgA0499VQguLU3atQoAJo0acKMGTPo1KgLsYoAAA0YSURBVKkT7777LpdccgldunShX79+lf6dxNuwYQPf+c53mDlzJt/+9rd56aWXGDBgAKtWrUJz+0kmDOiLZKStW7eSl5dX6aqL5513Hv369aNevXokO6fd/PnzOe644xg1ahQNGzakadOmnHPOOQB07dqVc889l5ycHNq1a8d1113Hyy+/XOk1GzduzHe/+12efvppAD744ANWrVp1sOfSvXt3vvKVr1CvXj3OPPNMvve971XpuolmzZpFr1696NWrF/Xq1eOSSy6hW7duLFiwIOlrSfbJup6L1G5V6VFEpVWrVmzZsoV9+/YdMcEczdM369ato3379uXWvf/++9x8880sW7aMnTt3sm/fPrp27Vql61511VWMGjWKsWPH8tRTT9GvXz8aN24MwJtvvsltt93Gu+++y549e9i9e/fB9V+SsXbtWubOnXvIRIh79+7lm9/8ZtLXkuyjnotIBc477zwaNmzI888/f8TjjjSPamVzrObn57NmzZpy62644QY6duzIBx98wOeff87EiRM5wnDjIS699FK2bNnC8uXLefrpp7nqqqsO1l111VX07duXdevWUVxczPXXX1/hdZs0acLOnTsPft60adMhsV9zzTWHLPe8Y8cObrvttirFKNGbN28e8+bNi6StrEsuUQzo63HkuqFZs2bcfffd3HjjjTz//PPs3LmTvXv3snDhQm699dYqXaN169bUq1evwmWDe/fuzaZNm5gyZQq7d+9m+/btB1e13L59O7m5ucRiMVatWpXUUz45OTkMHDiQH//4x3z66adccsklB+u2b99Oy5YtadiwIW+99RZPPfVUhdfp0qULzzzzDHv37mXZsmWH/GIaPHgwRUVF/P73v2f//v2UlpayZMmSg8snS+bJy8sjLy8vkrayLrlEsZ6L3nepO26++WYeeOABxo8fT+vWrcnPz+fhhx+u8uB348aNueOOO/j6179O8+bNeeONNw6pb9q0KX/84x8pKiriuOOO49RTT2Xx4sUA3H///Tz11FM0bdqUa6+9liuvvDKp2K+66ioWLVrE5ZdffshtvalTpzJ27FiaNm3K3XffzRVXXFHhNe655x7WrFlDixYtuOuuuw7pAeXn5/PCCy8wceLEg3839913HwcOHEgqTonO9OnTmT59eiRtaZnjo9CoUSOuueaao7qGljnW0rd1hb7X6VeV3zepWuY463ouUdK7NCIi5VNyOUozZ87ULTIRkQRZ9yiymfUB+kT1Epd6LyIih8u6nksUA/rlUe9FROQLWZdc0kVPkFWfni7KfvoeZ4YFCxZENoOCkksK6RZZ8po0acKGDRvYs2dPlV8QlNrD3dmzZw8bNmygSZMm6Q6nzmvcuPHBmRpqWtaNuUjtcuKJJ7JlyxbWrl2r6eqzVE5ODs2aNYvs5T2p2NSpUwEYPnx4jbeVdckl6gH9RIWFhYd8TsW7MNmsXr16HHvssRx77LHpDkUk682ZMweIJrlk3W2xdA3oV6RsLObXv/51ukMREYlMrei5mFlL4HHgUmALcLu7VzwhUoYpG4uJH/BXb0ZEslmtSC7AL4E9QBugC/Cima1w95XpDSs58QP+hYWF5OTksHnzZg4cOKDEIyJZJeOTi5k1AQYAZ7h7CfAnM/sf4BqgVs/tvW/fvoOPaCYmHtB4jYjUXhk/caWZfRV4zd0bxZXdAlzk7n0Sjh0GDAs/ngG8W8VmmgHFaTgOII/gVl8qr6k4U3vNVMdZ1Rhrou1kjlWcqb1mbfi3CdDB3ZtW8diKuXtGb8AFwKaEsmuBJZWctyyJNgrTcVwycSZ5TcWZwXGm89+m4sz8OGvLz1BlW214WqwEyE0oywW2p7CNosoPqZHjkpHMNRVnaq9ZG+Ksia8nGYozdWrLv80jqg23xZoAnwGd3P2DsGwGsNHdKxxzMbNlnoI1CWqa4kyt2hBnbYgRFGeq1bU4M77n4u47gGeBu82siZl9HfguUNlEXoWV1GcKxZlatSHO2hAjKM5Uq1NxZnzPBQ6+5/Jr4BJgK3Cb16L3XERE6ppakVxERKR2yfjbYiIiUvsouYiISMplXXIxs5Zm9pyZ7TCztWZ2VQbEdIyZPR7Gs93M/mZmPePqLzazVWa208wWm1lBOuMNYzrVzErNbFZcWUbFaWaDzOwf4fd6jZldkGlxmlk7M1tgZp+Z2SYze9jMctIZp5mNMLNlZrbbzKYn1FUYkwUmmdnWcLvXzCzqOM3sXDP7o5l9amafmNlcMzs+0+JMOOYuM3Mz65GJcZpZYzObamZbzKzYzF452jizLrlw6DxkVwOPmFmn9IZEDrAOuIjgTdkxwJzwF08ewdNwY4CWwDJgdroCjfNL4M9lHzItTjO7BJgE/CfQFLgQ+DDT4gSmAv8CjieYF+8iYHia49wIjCd4SOagKsQ0DOgHdAbOBHoD10UdJ9CC4ImmdkABwTtvv8nAOAEws/bAQOCfCVWZFGchwff89PDP/z7qOFPxJmambEATgsRyWlzZTOBn6Y6tnFjfJpgzbRjB9DbxX8MuoGMaYxsEzAHGAbPCsoyKE3gN+EE55ZkW5z+AXnGf7wOmZUKc4S+a6VX9uwv/zofF1f8AeCPqOMup/xqwPeHfRsbECSwEegEfAT0yLU6gA/A5kFvB8dWKM9t6LqcB+939/biyFUC6ey6HMLM2BLGuJIhtRVmdB+/1rCFNMZtZLnA3MCqhKmPiNLP6QDegtZmtNrP14e2mRpkUZ+ghYFB426Et0BP4HZkXJ1Qe0yH1ZM7P1oUEP0tlMiZOM7sc2OPu5S1cnylxngOsBX4S3hZ7x8wGxNVXK85sSy4xDp+crZjgtklGMLMGwJPAE+6+isyL+R7gcXdfl1CeSXG2ARoQ3Gq4gOB201eBO8msOAFeJvhB/BxYT3Cr6XkyL06oPKbE+mIgVpPjBJUxszOBscCP44ozIk4ziwETgZEVHJIRcQInEkz0WwycAIwAnjCz08P6asWZbcklinnIqs3M6hHcpttD8A2EDIrZzLoAPYAHy6nOmDgJbtUA/MLd/+nuW4AHCG49ZEyc4ff79wTjGE0IZsVtQTBWlDFxxqkspsT6XKDEw3slUTOzUwhuOd3k7kvjqjIlzp8AM939/yqoz5Q4dwF7gfHuvsfdXwYWEyzOCNWMM9uSy/tAjpmdGlfWmUO7zGkRZvnHCf7XPcDd94ZVKwliLDuuCdCe9MTcnWCQ9GMz2wTcAgwws7+SQXG6+2cEvYDy/nFnTJwEA6P5wMPuvtvdtxIMPPfKsDjLVBbTIfWk8WcrfIptEXCPuydOBZUpcV4M/MiCpwQ3EfxbmGNm/y+sz5Q4366kvnpx1vTgUdQb8AzwNMH/FL9O0IXrlAFxPQq8AcQSyluHMQ4AGhL8r7bGB/UqiLExcFzcdj8wL4wxY+IMY72b4Gm2Ywl6A0sJbullWpwfEixqlwM0B54juC2atjjDWBoCPyXoSTcMy44YE3A9wQMKbQlun6wErk9DnG0JxoJ+XMF5mRJnq4Sfp3XA5WW/AzIozgbAaoKnBHPC35vb+eJBjmrFmZYfuJrcCP63+DywA/gYuCoDYiog+F92KUEXs2y7OqzvAawi6J4uAdqlO+YwrnGET4tlWpzhD8RUYBuwCfg50DAD4+wSxvAZwUJRc4Fj0xln+H31hG1cZTEBBtwLfBpu9xJOIRVlnMBd4X78z1JJpsVZznEfcejTYhkTJ8G44OsEvzf/Dlx2tHFqbjEREUm5bBtzERGRDKDkIiIiKafkIiIiKafkIiIiKafkIiIiKafkIiIiKafkIiIiKafkIhIxM/uOmX0/3L/PzM6Kq8s3s/LmdhOpVZRcRKI3DHgq3D+DYGoNADyYjfo4M2uejsBEUkXJRSRCYdKo5+6lYdEJwP1m9lczGx6WLQW+lZYARVJEyUUkhcxsiJmVhFupme2P+7yNYNW/teGxeQRz4Y0GzidY4Q+CyS7/LR3xi6SKkotICrn7E+4ec/eyhaLml3129+Z8MWEgBOuRz3L3T8OezL/C8rQtviWSKkouIjWnC4cuDwvBmkMnh/tnAvsBzOy7BLMQE9b/A5FaTMlFpOZ0AZbHF7j7NuCAmTUEvgI0M7M5wHf5YgXQCwlWsBSptXLSHYBINjKzXIJVPZeXUz2NYC2fHyRWmFk+sMmD1TZFai2t5yJSA8zsAmA+0NzjfsjMbAJBz2Qz8B/uvjNNIYrUKN0WE6kZXYC3ExLLGUB7d7+AYP3376crOJGapuQiUjMOG28BLgAWhvsLgW9EGpFIhJRcRGpGZw5PLi2A4nC/mOAdF5GspAF9kRrg7t3KKf4MaBbuNwM+jS4ikWip5yISnT/xxbQu3wJeTWMsIjVKPReRiLj7O2a21syWEryN/x/pjkmkpuhRZBERSTndFhMRkZRTchERkZRTchERkZRTchERkZRTchERkZRTchERkZRTchERkZRTchERkZRTchERkZT7/zdcoPOlXcGyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.hist(Lmin_b,bins=int(np.sqrt(B)),color = \"0.6\", ec=\"0.6\");\n",
    "axes = plt.gca()\n",
    "y_min, y_max = axes.get_ylim()\n",
    "plt.plot(np.squeeze(np.asarray(T_m))*np.ones(2),[y_min,y_max],'--k',label='Test statistic')\n",
    "plt.plot(T_c*np.ones(2),[y_min,y_max],'-k',label='Critical value')\n",
    "plt.ylim([y_min,y_max])\n",
    "x_min, x_max = axes.get_xlim()\n",
    "plt.xlim([0,160])\n",
    "plt.xlabel(r'$T_0^{(b)}$')\n",
    "plt.ylabel('Frequency')\n",
    "plt.minorticks_on()\n",
    "plt.rcParams.update({'font.size': 12})\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Generate points along boundary of 95% confidence region and plot (see Figure 6(b) in manuscript)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAEUCAYAAABkhkJAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOyddbikZfnHP98tll26Y+lGOpVGJCQWpEsxABFBQEr8iSygUiIqjRICAoIKKClIl3R3Nyy5zdb398f9DGd29vSZOe+cmedzXe91zsxb99R7v3fLNplMJpPJFE2/ogXIZDKZTAayQspkMplMnZAVUiaTyWTqgqyQMplMJlMXZIWUyWQymbogK6RMJpPJ1AVZIWUymUymLsgKKZPJZDJ1QVZImUwmk6kLskLKZDKZTF2QFVImk8lk6oIBRQuQydQTkmYEFgDmBeZJf+cGZknLzOnvYOL3MxBYN+1+PzAJmAxMBMYCo4DR6e8nwIdp+QB4Hxjp3FAykwGyQso0IZJmBZYBlgOWBRYHFgEWJRRQa0ygRbmMBsYTiueLsm3GEr+pGQjFtTAtSmxmWvdIjJf0BvAG8DrwIvA88Bzwhu2p3XuVmUzfQ/nmLNOoSBKhFFYFVgM2IqyZcsUwGXiNUAglpfAWYcGULJmRtie2c547AGxv1IEsszKt5TU/LYpwEWAxYM6y3aYCHwPXAI8CjwFP2h7fwUvPZPokWSFlGgZJAwnls27ZMl9aPRV4gbCKrgf+RFghr9me1MPz3gHtK6QuHGsuwmpbDjgLeJewrmZPm0wmlNO9pcX2+z09byZTD2SXXabPkqyO5YBvAJsSFtBMafXrwH+JuM7DwFO2x0kqxW2uraIoj1frQLY/Au6R9AHx+zwWuJCw9FYD1iYU7f7AIQCSngFuSctdtsdUS55MpjfJFlKmTyFpBkLxDAe2ARZKq14GbgVuB+6x/W4b+98CzG57jdpL230kbQ/8A1jL9kOtrJ+BUFAbEAp5fSJ2NQm4E7gW+LftN3pN6Eymh2SFlKl7UubblsDOwDcJF9Y44GbgRuAW26938lgnABsC69Zzdpuk7wAnAkvZHtuJ7WckLKfNCEW9bFr1BKHYrrD9Uo3EzWSqQlZImbpE0mDgZ8DSwNaEEvqQCPD/C7gtB/fbRtLShBW5LbBeevpR4AHgotasrkymaLJCytQVkpYBvgfsTWScjQMuT8udticXKF6fRNIwwrrcFVgzPX0TcD7h1vuirX0zmd4kK6RM4UgaBOxIBOrXBaYQmXBfBZ6wvVmB4nUKSSNsjyhajvaQtCFwB3A1oZiGEWnlFwNn2X65OOkymdw6KFMgkhaQdCzwJvBXojbnSGAh29sC5wGbSFqgQDE7yzFFC9AJvg2MAfYkap+2AG4DDgRelHS9pC0k5etCphCyhZTpdSQtDxwB7AH0B24AzgD+U96ZIMVBXgAOs31qEbJ2Fkm2raLlaIuU9PA+cI3tvSrWLQDsC+xH3BQ8D5wM/LW9guBMptrkO6FMryHpa5KuBZ4BdiIKP5e2vbXtmyrb5Nh+Efg+cEXvS9sxkkZIsiSnx07LiIJFa41JRAzpd5UrbL+b3I0LE9bTBOAC4FVJh0qauTcFzTQv2ULK1BxJXyMKPDclGoyeDpyRikAbgnq3kLpCKjjejHCfbkzEmU4GzuxMCnom012yQsrUDElrAMcRtUMjgZOAcxrxotZICqkcSWsDI4h404fEZ3h2TrnP1ILssstUHUmLSboCeIhodfMzYHHbpzaiMkocW/5A0hBJS6bYTZ/F9v9sf5PIfnwSOJVIgPh2Tn7IVJtsIWWqhqTFgEOBfYjU7d8Cv7U9qlDBehlJ2xBxr6mAgB1s31ysVNVB0kbAKcAaRPfxU23/tVChMg1DVkiZHpPulPcjMuVEBMR/afudQgUrAElzE41dh5Q9PRZY0PbnhQhVZdLnvSsRV1qQaE+0be6bl+kp2eTO9AhJqxMdtc8E3gFMJCw0nTJKLEVktJUzhRgC2BDYnmr7MsKF9wXxmp+T9PPU9DWT6RZZIWW6haSZJZ1BxIlK6cIrEhlZf0iZWs3Im0TX7XIGAW8XIEvNkLQFkazyC1pmTP0aeDJ1hMhkukxWSJkuky5GTxOtfk4HlrX9V9ufEReo9YGtChSxMGy/DRxOjDj/PP09wPbIQgWrIulm4yRi5Mfptt+0vROwOTHD6Q5JZ0uapUg5M32PHEPKdBpJswO/B74DPAvsbfv+im36E/GFvzVzI1RJSwBLAi90djRGX0LSssRcqcrPfwiR6n8IMe12X9s3FiBipg+SFVKmU0jaBLiIGAl+AvDr3CW6hb7QXLU3kbQW0U18BeAcov1To6b8Z6pEVkiZdklziX5D3PE+D+xp+5Fipao/GrUwtiekBIdfEaUALxHfnTyHKdMmOYaUaZPklnmQUEZnAqtnZZTpLLa/sH04sAkwI3C/pCNzQW2mLfIXI9MqknYHHgbmB7ayfYDtcQWLVVf0seaqhWH7dmAl4J/EWPZ/S5qzWKky9Uh22WWmIbno/kCMI7gH2LWJa4o6TXbZdUzKzvsRcBrwAbBLZVJEprnJFlLmSyQtCNxFKKMTgY2zMspUCwdnAesAk4E7Je1TsFiZOmJA0QJk6gNJ6wD/AGYCvmX7moJF6msc2/EmnSNZElsThcYvAX93A7kybD+SOnxcDpwnaVXg4DwMMJNddhkkXUiMt36d6En2TLESNTeSTiMa1A4mhuX9G9i9kZQSfFmz9htievBIYH3bLxQrVaZIskJqYtKd+K+An6endrNdl9NZmwVJ8xE3BuXth8YBX7X9VCFC1ZA0nuNJooj4LWAT2y8VK1WmKHIMqUlJNSKXEsroAuBR4E+SVipUsMxsTN+cdRIwewGy1JR0Q3Q+oYyOoiU1fN1CBcsURlZITYikWYH/ALsTCmlvYDgwCviXpKEFitfsvEp8DlPTY6f/nyhMotpxMLAb8HPbJwJfAz4F/itp+0IlyxRCdtk1GZLmAW4GvgJ8N40RKK1bHVjR9kUFiZcBJC0JXEl00X6dSI9+slChakByT34fOKEUH5M0F/AvYtLwPrYvKFDETC+TFVITIWlh4BZgIWB72zcVLFJhpEF6AkY2WrJAXydZ6P8ENgMOtf27gkXK9BLZZdckpLvue4B5gU2bVRlJGiTp38R8ojeBm1NgvaP9hkjaTNI3UvFw5foR1Ze2OUlNWIcDVwGnSjq2iedrNRXZQmoCkjK6g8jc2tT248VKVBySjieafZaU0HjgLNuHtbPPvERPv9kIq+p9Iuvtk7JtcqeGKpPSws8j3HrHA8dka7axyYWxDU6FMtqkEWMRrSFpN+B7wGjg+DIlvB4tyoj0/3odHO5UoqffwPR4JuAxSRvYfqMKsi5INLCdgyiCvaGnx2wEbE9JnRwMHA0gKSulBiYrpAYmDYm7g+ZTRvsS/dKGEBezzSWtZftZ4AWidc2gtPmk9Fx7LE2LMoKwkhYCnpY0U9l5SxfKf9reoZOyzgs8TlhfA4BdJP0EeAXYlsg6O7uRJs52BdtT0+cJoZQMHFOgSJkakl12DYqkBYiY0axET7qmUEYAkl4FFit7ysDvbB8maQ7gf0QszcAnwNq2P2zneKcSTUErY02jgJ2A14AX0/FEFLIeYfvMTsh6OFGcPKjs6c/S4yGEwvyYyH78qKPjNSppZMWfCas3Jzo0KDmpoQFJrf1vAeYGtmgmZZRoLZbTDyDFfVYEtgd2BL7SnjJK/IJQYpUMITLB9qo47xDgZ52UdUagf8Vzs6RjQFhmswPf7eTxGhLbU4l2Sn8nEh2+X7BImRqQFVKDkVxINwBLAMObdELnqYSVUmIccGHpge0Jtm8F/gvsI+kOSVdKWqq1g9keD3wduKbiuAMIy2mLVnbrbILD1UD5KPhxQGWT0YFE3KqpsT0F2JMo6v6TpE65RTN9h+yyayAkDQCuJS6Q29u+tmCRCiGlCH+XcO+MAUbYfrCV7U4ADgSGEt0QRhMW0ztp/WzAfIQy+CewIDE2oT/TWjVTiCaopQ4XY4HjbJ/cSXk3JGJesxKpzrMTzW7LMwHXz9N6g1SndAuwGhEbvbdgkTJVIiukBiFdhM8C9gN+ZPucgkWqeySNZlrLYwJwpO0/pkD6HwgFVFI0bVk9k4GNCdferMDFwDndzQaTNIiw8nYg4lSH2L6xO8dqVFJHh/uIzMSv5YasjUFWSA2CpHOAHwIn2z6yaHn6ApLG0KJsICkkorXSo7TEcSpxWvqlfe603ZrbrlpyDgN+TzQhvZtQmk0/Tj6VNNxPWJCb5tEVfZ8cQ2oAUprwD4H3iPkyfRJJC0i6S9IoSc+kwW215ExaYkJTiVjOP4DlCaunLSYADxHD8y4hLJmaIGlmoih3W2BlohFuHp4I2H6ZsCQXAu5OlmWmD5MtpD6OpOWJDLCPgWHAc8A2tl8vUq6uktJ6nwMWJ5IFTLirlqxVunM650+IjLsPgaNsv5RGcNzP9BbS2PT3Ytv7d+E8g4FFgQ/Luzt0ct8ticmqs5Q9PRGYv6vHaiSSi/oo4NfEjcFSRL1Wpz+XTP2RLaQ+jKTZiSSGsUS3gc0JpfSApFna27cOmZ+40y0Va5fiNWt29gCShkm6JGXNHZlazyCpn6QVJK1WcRfdj3CBHQzsWopDpDT5kwhX0OfE+3swoby27qIyWgN4l7Co3pV0SGf3TUxh+tiV0vPNzM8JZXQZsBJwCvCj1Nkh01exnZc+uBAX0xuIu+V1y55fBvh+0fJ14/XMml6Ly5bRwDqd3H8O4AOikNSEEjmP6FJxO5FtN4q4m543ne+JdI7RRLeEWSqOuSSRrDBfJ84/opXnlGQqf01jgVW68L4MJjpJTCjb/4qiP6+iFyL78WBavDz9gZvSd+hrRcuXl+4t2WXXR5H0f0SFf8Nk1Ek6Cfgx4SobB9wLfNNRFNnRvnsCZzNt1txk4FjibrqUQj0J+DehKL5Py6jwycD5tvfr4Dwi4jgbE+19TrE9qrXmqslK/ZhpW3SNAfa3fUlHr6nsOLMR7XKWBu4kuk60F+NqSpLH4GGibmtV2x8XLFKmi2SF1AdJdSu3AX8D9nADfYiShgOrE+14Lm3vwpuUw3bEILu5CUVRqZCuZfqkg5eAkURPu3LGAHParixMLT/nmURnhqFEEsRrwKrA+FYUkohedLOWPT2OqJ15oK1zNAKpQHuy7Qm9fN7ViXTwW4jC8A5vZjL1Q44h9THSxNfLgZeBHzaSMgKw/S/bxzim1s4h6XJJT0j6cytxsT8TWW7HEW1lBtASWxlH1AM9TMSCSkwkXHWtFZnOQMThWiUlJ+xLS6r4DMCypeNLclpGpNdiImFiDBGLGg+c3sjKSNJQSTcTiniMpLNT8kiv4Cge/imwFdDmSJFMfZItpD5EuuO+nnAXre0G7lGXLv7PEIkOAwlr5AkiPjA11aA8ybQNT78g7oznJOJrJxI3XdcT1tBUYjDfBmn7yh52YwglfxmtkFKwP2FaF9xooivEPyotpLL95iZGxr9r+8WOXntfRtJ5RJeJ0hDDscBhvelWTr+TvxE3A19zc7bP6pNkhdSHkPQjohvDgbbPKFqeWiJpPUKpzFz29Hiitc9rktYEbmXadOhRwEa2H6s4loj4yyDgeduT0vP3AWvQMlpiDOH+G0d0C3/TFWMfJN1KZDTOQCQZfJqO/VFbCqmZkPQikYJdzpW2d+llOWYjbljGA6s5ptBm6pzssusjSFqGKAK8mSjobHQ6Snd+lrCISndUJjLRnq88kIMXbD9VUkaJrYnMrE/S8TYhFNRbRIzuDUl7SFpC0m6SNiYSSfqn800lFNIEInkiA28Q70uJUpytV7H9GRHrW4pICc/0AbKF1AdITVPvJ4pGV7T9bsEi1RxJA4kOBcsS7p/xRIbZlqW4WSoKvpLobP4ysJPt6RRSF845G/AO0xbETiSU4GRaFGJ5ksIE4Gjbv+3BeYcAFwDDCevssBRD63OkjukP0GJ1vg+safvzguT5LTGyfkvnfoB1T1ZIfQBJRxLxkJ1tX1W0PL1FytT6JbACcZE7sb0MuCqcbyViqOHM7WxWGsJXzvtE1t6fiaSFLv2oJP0F2JmWuMs4YCvbd3TlOPVCmse1MaHMb3GM7yhKlhmIBJZZCXfvqKJkyXRMVkh1jqSlCV/49e7kWOzM9CQrczHgc1cM5Esxpv2JoXrDOjhUKQ19QCvrxhIFsl2yliSNBOYqe8pEk9zODvnLtIOktQgPw7nOrYXqmhxDqmPUMrZ5PHBAweL0WSQtSrj0HgPelPT7pIRKnAiczrTKaCrxvn9Qcbh+RCr5VKZv3zOUGNjXVSrdWROBph1XXm0cs7B+T7QW2qCj7TPFkRVSfbM3sD7wU9vvFS1MPZE6g/9F0p2Sjk4WUPn6fpIOl3QP4bIZRiiMGYj3ddvSdkTdSqUbbjLwTaKlUOXz1xPp6JXrYPppr51hf1omxY4jXIB/6sZxMm1zNPAqMWl2ho42zhRDa26HTB2Q/PAnEIH8i4qVpr6QNCuhZOYivsNrENlU3ynb7FSiiLW1mUYzEqMcriEy5tpK136PcPOVM4ho+Ppjot6pnHF0I9vO9n8kfZUoyh0DXJZjHdXF9rhUNnEzcQNyQsEiZVohW0j1y2+IQOyPG60bQxXYjLB2SjdUQ4DdKzp570fbA/bGEy48Uhr49UybqgyRRPESYQWVWz3jiBjPrLRkkkH0yNvR9hVdfTFJjqds/9b2OVkZ1Qbb/yHmXR0taeGi5clMT1ZIdYikvYhWOL+3/UzR8tQhbVk05c9XKnETimg8oXzOlXR/ujDtCpxP1CN9RIwu3zTdCPyYsJS+IFK87yF6pfUvnfBq4FAYYDgf6aiev7zqIWmIpOGSdpA0R9Hy1AE/TX//VagUmVbJCqnOSMH23xLXutUlrVCwSPXIf4iMtlLG2zhiJMMXZducTcs02ClE4sBPiPd1JsLCWpMYTTHR9r6257Q9t+2Dy9LLDyVcgyKU0EjC1Tfp7XTgbYFTY/38wG+QjHRcDV53l5C0HGHlXUbUOb0gqdIF2TSk2NEuxHVvZUl7FCxSpoKc9l1nSNqJKPb8K7Al4Ro6Ezgou+5akLQgMURvEaJ/3W/KO4OrZRrsdoSF82uiK8OCFYcaR9SnvN7KORYhOj8Mrlh13CQ4on/Z8/0AhwVVvu0O2P9s5bhDgXG2nayWQ5JcN9j+eydefodI2oZwT5W7FacAN9veqhrn6EtI+hrRbHdJoiXVKkQG5Rq5I3gdUa3BSnnp+UJkgL1C1B31J4LmZwDnFC1bX18IV02p1VD58gUwR9l2XydaB40nEifGVu6zJ0yaCj5m+mMZ8NEwdSp4crgGXwc2SMdemsj0mkwkL+xOtNopyTWWGKPe09c6oDW50/J00Z9FQZ//skRz3s3T493T+/HtomXLS8uSLaQ6Io23/h2whe2by56X8wfVIySdDBzeyqpTbB+RtlkMeIqW8RKTiBuDaVzbY2DqEOgnmIQ9KO1r25K06dxw0wdpn7TjF0S3iVuBhWmJdU0klFN58sV4YGhPPm9J8xKKsNKymwj8yXZT1rSV/46SBf0gMT14aRfYTSLTQo4h1QlptMHPgVvLlRF8OVcn0zMeZtpC1snAXSVllFiPaZMhyt1dXzKALzVKeaNWJH0FuO430/+uZiBCTfMzbeLFVFom1rZ7zi4ykrCQKvkfcEQrzzcF5b8jh5vuMKI+rd0pwZneIyuk+uFAInj+i6IFaTRSMLuy7mQ8sE3Fc5+2svsUph3wx7UwJV3ZhiCV4jHHAsO/DQN+kJ6oaPHwEdOnlkPK1ivjk57egKSL7ZbE6xlLWGj7297A9rh2d24iHL0C/wv8LMX1MgWTXXZ1QCr0fA24z/bWRcvTaEham0h8KG+aOpYIaD9ftl2pq/ryRAHsRCIAvh3TFpGPnwKD+7VYOwamOm7wVHqiQtMsA6wFnEuLYnoFWJFpbwzH2a7KxTEp4oWAD2yPrsYxGw1J6wD3AkfYzmMqCiZbSPXBQcDsRGfrTPWZxPTf9X7AxNSC6EpJjxPFsl8hdMlHxOeyJdMqIwMfP9iisKAlJVwQ2QoVymgSYflcSiilHxMuvD9SYX0x/RTbbmP7C9svZ2XUNrbvI7Ivj0xu80yBZAupYNKIhTeBu21vW7Q8jYik/sBdwKpE26Bx6fFORFr3vEzfRmsiMbRvRaZ3q71MNFG9/z2YPF8MgFsJuHkYXPROxKtK8aJJwMW2925FroFEHdTKhKLrR8ztuatsm6GE8hqSZJ6NmIz7Qtpke2JO1uO2b+naO5OBLy3oB4BDbf+uaHmamdzLrnj2JqyjE4sWpFGxPUXSN4iA/krAQ0Svu00IN15rv4NBwAKEe61SIS1J9ETz/BFjeoFI6V6BSCiYsWzbAUw7WqJcrkmSNiJ62M0G3GP7jdL6NDDwEWCeJMNgwtXYj3AlTiGm3g4irL3TbB/d8TuSKcf2/yTdAfxU0hmu4cytTPtkC6lAUu+1l4HXbG9YtDzNhqRNieLR1lw1k4iOEOsTnR06696eSCihyu3ntP1JF+U7BjiK6TPxIFx9/QllVH7u+bt6ngxI2pxw3X3f9oVFy9Os5BhSsexKBJ1PKlqQJuVuWvrUQVhDk4HRhBv1B0Qn8cvo/FiJgbT+u1q7/IGkjSVdJukCSSu2caxhtK6MIKywyrvJSYS1nek6/yEa6R6RapQyBZAtpIJIPeseIS44K+Rao9oiaRNgHUIBXVxyy0iaHRhBuNzuJopXRVyc1ibceh8RSus02u4g3h4GdrB9dTrnVkR7qCFp3ThgHdtPVsi8I3AprSulcUnOkntwKvAOsLjLWihlOk/qbXcpFYXpmV6k6FYRzboQF0cD+xUtS6MvwMFE7GUKkQR3PzCgg32+l/aZSlz8nyPife8SPetabRvUzvIUsDHRBeLhinVTgfMrzj8HMYJkUivHmghcQmTsvUooyyeAJYp+r/vyQrg/3weuK1qWZl2yhVQQki4nJpIOsz2maHnqCUlLALsRF98rbL/Sg2P1JxRKeaxlDLCb7eva2e8zorFtiVJ38Zkpd8lddx0MHQoXXQTf/S6MHQtbt1lKNopwB44ClqtYd7nt3dO5ZwaeBuarkBviPfkDMUU4/3irjKRjiemyS/Xke5fpHtlXWgCS5gd2BC7Mymha0riNxwk32gjgsdSSpzP7Dkxzf/ZNoxcgLuiVWXIQWW3tUemaG0i4zqb9zQxNNax/+cu0j1tnFiLFvDSqvMQ44M9lj7cjGutWKiMIN903sjKqGecSlvT+RQvSjGSFVAx7EZlYZxctSB1yPOHW6k+8RzOl59ol1fSUxr2fBjwsaWtH08xHmLbvnIhBe+1xMy3zliDcatVQAgOI9O1DCTfew8Autm8r22YwbQ8hBFg6xb4aAkkLS1pB0hqSfiBp0xRj7XVsv0vMu9ordbrI9CZF+wybbSEuNC8SjT0Ll6fg92JuYizADGXP3cn0MZM7OnGsPQlXXPl+H5Sd55a0/lXSOIiK/RcgCk/HEHVF1xIKqXSsz4EbaRnrMJG2Rzy0t0wG/tLBa1mYcOtNbeMYk4AZi/78qvD5i7AMyyf5jiXcmpeRkq4KkGuL9D7vWPR71GxLtpB6n3WBpYiR2U1LqrF5iyhSfbPMLXcZ03aqHgtc3olDzsv0nbJnB7A90vamtmeyvbjLOiEkWUQ02fwaYZ0tDQxnWlffzMRQt+toGUsxiOuvh9tvjy1uvz1iSu0zlhgc2Ca23wQ2JLpzv04U25ZaDI0F/uDGGJewM1H6MJgWq3AIYRVvQ/xWiuAW4G3g+wWdv2nJnRp6n7OIu+tnihakKCRtSHRNmCEtQwk3yVLAeUT85OC0+WnpuY64i7A+SnGXSUQ2XWeYl2i/097vQUSSwY7pfwH92GqrqYTC3IONN55K+27wL4CzbH/ekUC2HyMUJJJmJFoVLUa4Gq/saP8+wldomT1VyVTi/S6CwYQrdTtJa9p+qCA5mo6cZdeLSBpMuCNKjTgfI8Yq/8FN9EFIOgA4mWlb7JhIxe72OGlJ3yEU/oyEdbGt7ZGd2G8o4aqrZBKdm090A2Hp7UC0D6pkcjrWY0RCwnhJGxOxxHHE5/9CK/s1NJL2BM6hdaU0jqjPe60X5VmZiO1tXybTWbZ/3FsyNDvZZde7lDpH70JYAFMJP3XTKKPEy0w/G+jdzigjSYMknSbpFUn/S40xAbB9MeFaG2R7nc4oo8SAVuSB+H1UFplOZdoEibFEyvYIoglq+bqJRKHtqkTHh/WTMhpOuP72IobDPSRpVUlflbRMJ2VuBC4DrieUzxgiu20qET/buTeVUWJBwlV7OeEyfZiKDhuZ2pItpF5E0pXEF31Bp2p6STO5yVK/U8zmHCIRYRJhLW5u+4FO7HshodBL1tVYYBXbL/dAnv7EBbFy5DdJvsnpfJOJi+WFwHeIC+iJtv9Qdqzlidc2jEjQOLDy85X0BNHktYQJd94XhEX2N+AHzXCjkr4LyxE1X08R7/EXRbz2NA9rgO0J6fFPiSa8S9t+qbflaUayQuol0piJD4ELbB9QtDz1QLp4zw085U42BJU0jmldfV8AR5YrhW7KsiuRMl6Z6juamF+0BqGMzrH9Tg/PVeoO3hZjgW87tRrKFIOkYUTizS9td1h6kOk52WXXe3yTuJA2SkC6x9h+1vadnVVGiUkVj6fQ0hy1J7JcAWzAtAWrJtoEXWX7INtH91QZJa7pYP0MRDp8pkBsv01Mk92+aFmahayQeo/hwMfAfUUL0sc5lhalMYmoD6qKkrf9IOFSfZOIZbwC/BBYKrn1qsVKtF9k+wVNnIVZZ1wLrCJp4aIFaQayy64XSL7pD4F/2f5uweL0eSRtTyj494FTu5C80JVzzALcQbjWTBTLblSNeJ+kl4ghf+VMJmqNBhI1agc2Qwyp3pG0NPHZH2D7zKLlaXSyQuoF0lTQ24kRBP8sWJymIGWrLQ+8avuJbux/JjEPqRRTmgCcZ/ugKsj2JDEavZx3Cbfup7bf6uk5MtVD0vPAm7Y3K1qWRie77HqHrWiZQJqpMZL2Jmp+LgLuSx2c29p2bknnS7pb0gll/Re5uTgAACAASURBVMtWZdoEh8HA7qkTd/n+M0jaWtLOkubtpIiPtPLcSNtPZmVUl/wb2CglJmVqSLaQegFJjwGf2d64aFkaHUmzAh8wrTIZT6SGv1ix7YxErGZBosPDeOB221tJOgPYh2k7bk8BbrK9ddp/JqIbxCLpfOOB9Ww/3YGMqxIdF0odxccB37V9VddfcabWpFH3/wG2tH1j0fI0MtlCqjGS5iF6oN1StCxNwnxMn4k3kRgVX8l6TDvmYUbgG5LmBo4iYlTl9Ae2kHSYpM2JnnRLkopxifESf+pIwNQWaAPgKqJAdpesjOqae4hEk02LFqTRyb3sas8m6e+thUrRPLzJ9N0VBgDPduUgtkdL+iUxIqS87qk/8GtC6b3CtMW0IgpiO3P8R4jmopk6J3XXuIeskGpOtpBqzybAZ7QeN8hUmdQFe0vgU8KFNg7Y1fZ7rWx+D/ARYUGRtr+1LGvvKkLBlUaWlxhE9DqrTEwAGCZpRA9fRqb+uAVYoQtxwkw3yDGkGpMydF4uxR0yvUOqG5oH+Mh2pQuvfLu5gJOITuP3AMfa/qJs/UzA3sDqwLeYvhHoFFrGVNwJDLc9qlqvI1MfSFqHKJLNmbI1JCukGpIudiOBn9s+oWh5Mt1H0mzAa7Q++tzEb6kqU05T3HENwsp7INcjFU/KvvwcONP2oUXL06hkl11tWSf9vbdQKTI9xvZnhPv1VabvsiCgGi2FkLQm0Q39MiKz6xpJ+XdaMMlqfpjihgY2BfmLXlvWIYLfecBXA2D7UdtLAPsz7VTbccDpVTrNFUTW3qzE5NSvk5Mf6oV7gdVSuUCmBmSFVFvWAJ5skHHTmRbOJcYSjCbGVpwNnFKlYy9Q8XgwMSk2Uzz/I1o7tZbMkqkCWSHViDTnZVXg0aJlyVQXB8fYnsX2zLYP68mk2wqeJhIlSnxBztCsF0q/5VULlaKByQqpdiwMzEGTKCRJAyQtLmn2omXp4+wIvEGkoE8imsfmllP1wRtEoslqRQvSqGSFVDv2Sn9XlLSVpIWT1dRwSFoKeB14Enhf0tHFStR3sf0GkYK+NDCX7WMKFqnpkTRU0lrA9wkX7bcKFqlhyWnfNULS34EdKp7+vu0Li5Cnlkh6hhgoV7rBGQtsY/v2bh5vZuBoYrT1fcBv26slymRqiaR3gfnTw8nE93xgFd20mURuHVQ7JhFWwyrACkQg9I4C5aklyzCttT2AeN1dVkiSBhHZTEsRAf2vA2sD2/VczEymWxxJJLA8RaT+n0u45F8vUKaGJCuk2rEs8Jztz4kLbCPXIr3HtD3cJhP1Ot1hbWBRWnrEDSEams5nu7LZaSZTc2xfUvpfUslSWo6skKpOjiHVgFTIuAzwfNGy9BI7E3eQnxPuumuAf3XzWAOYvvDU1PnNU2f610laUNJ2ktZt1HhiE1D6TS9bqBQNSl3/yPsww4gO0S8ULUhvYPt+SUsQbrqPgMd70O7mf0Qz2hmJmo8JRKZiVToh1JBjgBFtrUxTg6+jpffdTZJ2ym2B+ha2P5L0CXHDmaky2UKqDYumv68VKURvYnuk7VtsP9aTi6ztccBXgWsJn/0lwBYNcOG+nGjMOkv6uzkwvFCJMt3lNWIoY6bKZIVUG0pf1jcKlaKPYvs92zvZXsn2vrZHFy1Ta0gaIcmSnB47LSMqthPRebycAURgPNP3eIOskGpCVki1ofRlfbNQKTI1xfYI2192+S79b3tExXYmBgSWpwlPJXdg6Ku8Diya44DVJyuk2rAI8EHuYZcpY1tahv1NBP7P9n3FipTpJm8QMc65ihak0chJDbVhPiIVOtM8HNveStuvpsSPeYDP881Kn6b0256PmHeWqRLZQqoN8wIfFC1EpveodNO1sc1U2+9nZdTnKf228zjzKpMVUm2YB/iwaCEymUxNKP22KxNVMj0kK6QqU5ZRlS2kTKYxKf22s0KqMlkhVZ8ZiIDnJ0ULkslkasJnRPeQOYsWpNHICqn6zJL+jipUikwmUxNSGv9oYtR8popkhVR9Sl/SuizmzGQyVSErpBqQFVL1yRZSJtP4jKLlt56pElkhVZ+h6e/YQqXIZDK1ZAwwU9FCNBpZIVWfQenvxEKlyGQytWQS0Y0+U0WyQqo+pe4Xkzu7Q2dm6fRk++7uU8R+zbhvkecuYt/e+v7WeJ/J5E43VUd9v6t/fSFpS+B6YG3bD3ZyH5cadNZi++7uU8R+zbhvkecuYt/e+v7Wch9JtwKDba/XleNn2idr+OqzX/p7lqQx6f/HbR9clECZTKbnSPo9MYQSYCXgiwLFaUiyy676lEzOdu+yOjtLp7vbd3efIvZrxn37qty9+d2o530S2b1UZbLLrspI2gK4Efia7Qc6uU922TXxvkWeO7vsurePpNuAAbY36MrxM+2TLaTqU0pmyBk4mUzjMpAuJC5lOkdWSNVnUvrblfhcu7N0qrB9d/cpYr9m3LfIcxexb299f2u5zwBafuuZKpFddlVG0leB+4GtbN9QtDyZTKb6SHoMeMv28KJlaSSyhVR9Si2Dcp+rTKZxmYXcHqzqZIVUfUpNVbNCymQal5nJDZSrTq5Dqj6lu6bceDFTMyT1A/YAlgeeAy61PbVYqZqKWcgKqepkC6n6lIphZy1UikzDIknAecABRBPfHwEXpOczNUbSDMQgzuyyqzLZQqoytqdI+hiYu2hZMg3LYsDWwBK2x0r6HfAysBTwYqGSNQel3/aHhUrRgGQLqTZ8AMxbtBCZhmUW4GPbYwFsjwNGkt3EvUXpt/1BoVI0IFkh1YYPgXmKFiLTsDwPDJR0hKRFJB1KzOF6pmC5moXSbztbSFUmK6TakC2kJqOz/euqge0JwObAN4B7gS2AzW2P7y0ZmpzSbzsrpCqTFVJteAdYKAeZm4pjevNktl+zvZntYbY3tf1Kb56/yVko/X2nUCkakKyQasPrwGCy2y6TaUQWBd5PlmqmimSFVBveSH8XKVSKTE3p6fiJTJ9lEeKmM1NlskKqDVkhNQG2R9hWaVxB6X/bIwoWLVNbFqHlN56pIlkh1YZX09+lCpWiACTNLGklSXMWLUsmU20kDSRcdq92sGmmG2SFVANsjyYCnssULUtvImkT4BXgcuBlST8oWKTepKfjJzJ9g8WJhgLPFy1II5IVUu14Hliu/IlGzrpL7VQuB3ax/RVgLeAESUsWK1nvkN10jUvF77b0m36uCFkandw6qHa8Bewq6WRgxbT8ELi+UKlqx/zAF7ZvB7D9kqRHCSvx5UIly2R6xn3JVfcULW2DXihQnoYlW0i1Y0Yi9fsg4mJ9O/BxoRLVlveBIWlAIZIWAlYlK6NM3+d24BOiAHkrYKrt3Fi1BuSJsTVC0paENbS97auLlqc3kLQ1cBHwEpHQcZztPxYqVCZTRSQ9B3xoe8OiZWlEssuudtwJTAVWBppCIdm+TtLywNLAm7bfLFqmSiQtChxKNCK9zvZVhQqU6TNIGkJ8t68sWpZGJSukGpHGAjwPrFa0LL2J7Q+p0x5fkhYE7gMuAB4HfiNpbttnFStZ7yOpP+F+mhu413bOGuuYlYgwx6NFC9KoZIVUWx4Dvi5Jzr7RemBP4F+2fwEg6THgKqBNhZQyrIYA48o/Q0lrAsOJAXkXJEXcJ0jK6BpgPuBZIhtyb9v/Klayumf19PexQqVoYHJSQ215gEhoWLhoQTIADKRloi+EMmnzpkzS2kSLmI+BNyWtm57/Ji3ZkksCD0qqq+7uktaSdIKkoyUtULH6W4Rl9DXbewHbAWf3upB9j3WA94gM2kwNyAqpttyb/q5bqBSZEv8A9pL0c0n/Bu4BHm5tQ0kzEVbEQbYHA/sB/5A0G1EE+z3bR9veG7gRuKJXXkEnkLQFcB3wBTEq4UFJw8o2mQ94zPbk9PgRYJ5GrpOrEusS7s3s7agRWSHVlqeIO/KskOoA288B+xCjIlYEbgNWlvTDVjZfChhp+5r0+CZiztWaREJEecLGG8BG5TtLWkbSjyTtIWmwpPkkXSnpWUnXSKql1fxLYJ/Ua+8Awi35o7L19wHbSVoxue+OBu7JF9q2SfHHRWi5yczUgKyQaki6A30AWK9oWTJfsgpwpu1Fbe9CxJUOamW7D4BhkuZPmXnPACsQGZMfAadJWlrShsCB5TumFkr3EAkt3yMyLm8m+p/tQlhlt0iasfovD4jpseWzet4BZio9sP0ocBhwFzAB2BDYo0ayNAqlm8qskGpITmqoPXcCx0may/ZHRQvTbCQX2zFED7LHgEHERbjEBKB/5X6235V0EvBg2mcQMAI4F7gbeDqtm7XsXCUL40Pgu7avT26w/6bzH5WskKckbU+UBDxQtr+ATYi41JO27+vmy76GUJg/ImJFhwDfrnh9f5V0GTDQ9sRunqeZ+DowmsjOzNSIbCHVnluB0oUm04uk/nr/JSyGvxCdI9YC9pG0d4q1XARcmLaXpIMkPZky8D4EtgdmJnr0HZ+y6a4FHrE9W2vjJwADT6TnTCivIUTnDiQNAmYHKkeOnwacQWRzXS7piG6+9OOIG6F/Ar8HDrF9W+VGDrIy6hybAnfYnlS0II1M7tRQYyQNIFw8V9nep2h5mglJGxAX5NVtOymC94DvAPsSbqyrCReeJe0L/ATYm7CI/kK4tn4BnGz7ckmDiYv9abavKDuXS4opWR4TgAOIUQX/BT5Nm/4Z2Dyt/5btqWmf5YFbgOVsj0qZcc8Bi9tu5JZTdY+kxYku9j+xfXrR8jQy2WVXY2xPlnQbsGmuR+p1BEwpe2yie8bjtrdtZfudgCNtPwAxERbYkVBQ10nahwhsP8D01frl4yd+RFhenxKKpx9wA7A+0SXit8AZJWWUmBd4pdQjLbkMRwJz0dg9EPsC30h/bylUiiYgu+x6h/8QF7JlixakyXiCiPFckHoLXgo8BLzbxvbjiJToEvMSBbGPEAkNJwK7EokQc0m6WtI7kh4mMtkAsP257W8RDXbfJyyhHxCdz58GRrXi+nkKWEbS1pL6S/oeccP4evdffqZKfJOoPcodvmtMtpB6h+uIwsNtyXNUegVJyxCp2uMJJbIJMa/p2Has1JOAa1On8kFEivjGad0YQiktQWTP/RhYDPgfUfj8mKQlbH9ZNGl7qqR5SIHw5BZ8glB002D7o5TocAlx8/IssLXtL7r/LmR6SsqE3Ay4MHs3ak+2kHoB228TxYfDi5alLyNpNUnDJS0naX9JP5O0qmJs+vclHSipNDb+fOBU28sTCuAT4EHb49o6fspq247IcpsDWM/206lW51oi0+oFIv60PDDC9jdsL52Of2Irh72DyLKcQdIKRHr1nW2c/17biwODbK9o++kuvkWZ6rMJkZCS2yr1AtlC6j2uBY6VNK/tD4oWpq8h6VQixvMccZH4H3A/Ud8znnDPvQf8UtK3CIXxN4CUJHBDeu4f7ZxjdcL19g5h9UyUdDCwBmG1rJhigtcTc57eTPuVfkeztXLYHxCuwjFE2vChttutZbE9pb31mV5lOPG53VGwHE1BtpB6j2uJIPt2RQvS15C0DtF/bUUiOeB2YpT0EcC/gVltD7f9Q2B/wvX2PJGQgKSZiThAR+7SS4manR8Q8aQDgM8Iy+jTslY7rwKTgfMkHU60DhpMK3fRtj+xvSXhrn0b+JWkC1Jrokwdk240tgVuzOnxvUO2kHqPpwh3z65EcWWm8ywKPGz7c0mzE9bRV4kWPlMJRV/iaaIYdDhwU2oLNB9Rk9OeddSPmHUzlMimOiKdYyzwU2CcpCOJjg3fJiyyqYSbbipwOnBexTEXA3YA5iGU3O6EojwJOIdIjiide4jt8savmeLZmPjs/la0IM1CVki9RApoXw4ck2pMPiTSSeexfXGx0tU9jxOdB5YlMhavBz4nrJJliJrW5QiX3fHArbafT7U9ywOf2e5olPraRFr4AYTrrdTlYG8imeFK4OeE6+5zIp18eJJtSpn1BCHQikT90d+JJIgZgJdsvyHpAKKuheRePB+YUdIrwA62czZXAaR+dbsBl9l+l7h5HE1Y5ZnewHZeemkh0r5NuJzeS/+/SCpQzku7791ewCgijfpDIg33Y6Le5xBgJOFmu4SwNrp6/AeJIti3CWUziXDFzZrOtW/6rGZL23+LUDBtHe8q4MD0/z6EZXRuerx6OuaSSe7VCStvf8KtmL8PxXzHvpd+k1MIK3ks8Nei5WqmJXdq6GUkjSfulq8lLp7Xu85Te1OW2QxuJ0Othudel2heOpAY8XAf8L6rHPiX9BZRuPohoZiOIJIWBhPuvleAVW3vm7bvD0wkMuKmSBphe0TZ8W4Ffmv7phQveoaI2V5IuO9+lvbfyfaOZfuNAha1/Uk1X1+mc6RygT2JG5B5gGNsH1esVM1DTmroff5I3A3/0vY/+4AyOoxwW3wi6RZJc3Sw/UyS5k2NQnt67q8SjUJvJ1r8/A5YpwbKaABhmfye+E1cRKRx/xvYw/ZhRGxqU7UM4tsdeL5MlmMqDnsjMCK1nRlGZNndndaNIOJTuwJrShqS5FghrR/VCZnzb7cG2H7B9tFEnPJj4OSCRWoq8pe69zmFcAd9v2hBOiJ1N9iPiNMMJVxWrU4WTY1JjyPGNjwLPKCeT1H9AdFwdE2iEPVFKkY99JSkDG4l0rw3JpTB88AfbB9p+34A23cTsZ7nJT0H/IZQSm1xGhFDui/9/RuRyPBa2vct4CXiLvy51P/uv8B+rohHVci7XGr8OknSy8mCzFQRSfMDWwJ/sj2ho+0z1SMrpF7GMYLiWuDbqdlnPbMOcIntt5IlcDJtDxvclkizXpTov3YnPc8m7Ee4TiYAhxPKYtUqv2+HE0p0eduzEn3m7gO+KukqSZuXNrT9KyLdfGciI+9bkqw0dqL0f3LfTbX9f8AChJvuKOJzP4NwC+5K1FVNIdyRNwEb276sLUElDSS6fpxDuBIPAf4pae7qvR0Zovluf1IX+EzvkRVSMZwPzEkExuuZ94C1ytxDa9N2H7g1gStsj3QEJs8kCkp7wn8Jd9dj6e9qRMPSlXp43HKWI34Hp0oaThTFbkBkVl0HXCTpm2Xbf0iMo1iP6Pg9zfgJ4kI2TtJdkq4j3ocDgIuJxqoDCGtsFHHBu5XoBNDf9rMdyLoI0M/2ubYn2f43UUqwt6QtVLuBf01D+q7vTUzQfbFoeZqOorMqmnEhLlqvAHcXLUsHcpZGLdxPuJw+BL7axrb7EhfXAenxXsB9PTz/UkQW2j8I62Kn9L6t1MXjLEPEn84g2gGVnh9KKN13iCSDl4hRIbeXbbMncE36fxCRcv5cel/eBJZJ65z+Hk90A/8WUSg7lYhFLJfWX52e+4JIlviMsADvJdxE/dt5HbMTimyB9Hhxwv37SNr/KWDOor83fXkhCqgN7Fa0LM24ZAupABzurzOB9SStUrQ8beHwn38DOIG4uK7uNJqhFS4kssYek3QjUTC6fw9FeJm40M5AXMh3S891usdbql26m0gqeB34h6TtU+xoe6KO6F7CepmZ6GF3VdkhytNQ9yEsnBVtb0i4MM+W9GPg9pSEsRcR5zqFUHTjiQLeP6REj1K6/2QifnRmOubyxGC965Nrbjpsfwr8CrhP0jmEInrK9uqExXYn8MvOvjeZVjmAKC1os4g6Uzty2ndBpI4DbwOX2967aHmqQUqFXo+4AD9ge2QVjjmI6JTwFcI99Vt3IdAs6XRgpO3jUvr1HcAqhGVxP/AGkWCyIFEQ+yiRYXcEYcmcDOxt+4bUT+8D2yenY69M1C/dkGTbi8igPB9Y0vYuku5M5xtKKMVZCKU0kFCQKxOtiN4hBvfdAlxs+6J0DhGdIbZN+5+S5Fw5yf0L2zembXcAvm07t6fqBpKWJBJnjnNZCn+m98gWUkGku91LgT0UIwr6PLan2L7T9r+roYzSMSfaPtH2t23/qivKKDGYUDAApxLf+eeImNRsRCLGbkQB7NHEnfF3iJ6DOxDxhFGpBdEEYBdJsyVFcTLhjtve9s8IN91Awn05l6TjCbfj8UTywvWEwj6M6PZwHqEYJxNuwSnEvKYFyuQ/gOgQcSXhkruNiDnNRxQC/zR1Ep8xyfpQF98fJC0vaS9J36hGun4f5ifE55RbexVF0T7DZl6Izg1TgeOLlqVRF8Ll+C6RQj6RcKG9Q1zY9ydqje4mLJyzgaEV+/8fYcn8mXAXPkTEcd4jYkjnlm07FxETOohQXjcDWxMW1DkVx90/nXMs4WrrR2Qovk5k25W2ex5Yo+zxX9I+JxIZgRNo6SR+GTCwi+/PrkSW4V+J4t0LacJOEemzG0fMPSpcnmZdCheg2RcisP0pMHPRsjTiQsSF7qFlfPnDRFLJ1en/49rZd/702cybHs9OJHasTmTKbU4kQWxMWDV/JbWaIVLm7wCeBM4CtqGVZAxiVtNdSVGOI7UbKlv/AtEhovT4WeCmsscHE9bT7N14b/oTltqK6fGMhPW4cVeP1dcXYgS9SckneSlmyc1Vi+ckwtWzD5EJlqkuZxMxmtmJ2M4phIKah4jn/LZ8Y0kbEOnqbxAW0btO86tsfyrpNULJHUMonUFEwsd4wiLaP217H7BRqmO6lHDdLSfpCkfnByRtSrgLnyWSKh719F0ozgMukfR/wELExNpzyta/Q7Qv+rQb781QIqHi6STzeMVE2wXa3avBSLHFA4BrbeeJzgWSkxrqAEm3E6nJS9geX7Q8jYSkkUTd0omEq2siYfmsCJxh+6yybU8n0rwfJLLtXiOUzhFEH73hhDK4lBhfviuRxHAF8KztX1ScW4RFtb3tuyXNSiRN7EV0Ej+TiC/NRiRubGT7mVaOsTcR0xpDWDA7E3EuEfVNJ9r+SzfeGxEW3J+JllarEgW667uJOo6n9linEG2p7i9anmYmK6Q6QNJGRL+2Q2z/vmBxGgpJzxNWy+NELdOahGvmL8D+tqem7Q4hkh7OI2J7E4iYzknAkYSF8zoxgvyXRGr4N9M2HxMFq192dUjHnJmIzwx1+qGlFkE3EUrmFEdxK5KOBua2/ZP0ePF0/AnA321/np4XEXzfj7BuzgdOcjd/yCmz7B/EDdE4YB/bTZPynKyj1wjrdPOOts/UlpxlVwfYvoMIsh8laWjB4vQpUg+9HSWdJumIVt6/wwkL5hgiGeFFoph1vzJlNAD4NZGQ8CNiRPrcRHzlE9vLAoNtL+Gow3qLaMT6NJHAsCiwfCsZamMIl1ppEN8SRLzpCSIbb0zFtgPTdmsSzT1XJZIiHpQ0J3xZfXsVoaiGEC2JTuludpztl22vTHQOmbOZlFHiQCKhIddv1QFZIdUPvyTiGj8uUojUqXvzei7YreAXREHpO0Ts53ZJg0srkwWyKZER90+i08SbFccYTLi/3iasocGEtbAskb79PaJOqcR/iISA4cSk2AnEBX2+8oMm5bEDcLxivMVjxDiDJwhX25mSNpO0M9EpotTH7iTgUNt72/4WkYX3k7JDn0ukkA8jlOE3iC4W3cb22O5aWX0VSbMQNyzX2/5f0fJkssuurpB0E+FSWrKbQeqenn9D4u77KcJF9S8i66suvyTJshkNLG77vWQl3En0mLu6i8e6h3jdqxBZdP2IpIaRRILD5sB2tu9NLtbfE9YRRBzmE2Au2x+XHXMmor3PSKLbxEdOY8qTrPsSHcMnAqfavimtexL4ru1H0+MDieyv/dPjt4F1bb+RHv+CcAseVfGa5iWU61slazDTQqoT+wWRVv9I0fJksoVUbxxJZIP9vPSEpAUlrdVL57+YqPTfhOiMsAlx912vDCC+wx/BlxbJ+4QrCwBJq0t6UtJYSfelmElr7Eg0L12eUERnEr3tNrC9JxGz+WPa9l7CgvohMfX16nT+cmW0bjrO5UQiwpYlZVSS1dEkdUPbm5aUUeJWYpbS7JKWItxKt5atfxXYIp1nIPE5vVp27n6S/kzUMN0P3FNy+TUrkgZI2qbk2pQ0jMhsvDwrozqi6LzzvEy7ABcQjTeXJYoyxxIFizUtViQu7lOI4HzpuT8T83mKei/WJi7EjxAxnumKPonC1osIRfI9Iomg1Hx0DqKAdVcixfsgIobUYfEo0TPumLLHiwDvlD0utSFyK8uxRDHulmnbRQlFuUInX/fg9D0YRyRM/LRi/VcI9+Kd6fVcTWpqm9b/kEhtn4lQ2KcT7YgK/34X+F3aI302dxGxuQvT72zRomXLS9nnVLQAean4QCJW8QUR5Dbwd2CxXjr3k0SWFUTNy1vA2gW9D0sTrq7vJMV0C/DHVrabhcg0eyEpiPIi0q8Dd1Vs/xqwVCfOvzlhdSxGJBucS9xNt7W9y/6fB/i4Yv3VxLjyrrwHJZf67kR86azSd4FIFd8EWIuym4i07hzgx2WPVyGasBb+/S5qIWJ++6bv1NT025ru+5SXgj+nogXISysfStz5Gjiol8+7PDHe4U0iNnNwge/BYcDpZY8XJGIwXTnGykT8Z0h6XMqcm6uT+x+UbgwmEu1/2uyGUKGQBqQL34bp8XxE0sUq3XgfDk7K9ru0WF4LdrDPEUSKe//0+GekERrNviRF/nZSSosULU9epl1yUkMdImk2om5mNLCa7Um9eO6BhHX0sVPtSxFIOoDIiCulTK8A3Gx7wfb3nOYY/QgXzTJEN4QFiNY+I7p4jAG2J3aw3Yjy46YuDJcRcaQliS7lJ3X2vOkYCxNxLAir+QBiYu/Ltk9pZ7/BxHDBBYjWR/MBm9h+vSvnb0Qk7QlcQtSgnV20PJlpyQqpTpG0HeHmacpiWUlzEY1MryXiJIcQLpbTu3CM8wll9BxxIR8PrOXp2/PUhPQaliViT691cp8FiZY+rxKW8gpEOvsAolbtNuBp2yd0cJz+hDtvRuBBlyVUNCupU8YLRIHzOs6Zh3VHVkh1SsoGKo0rWM72OwWL1OtImo9wWc1BNBT9Zxf2nZ+wiobZHptSxJ8G9nIdnYbflgAAFcpJREFU1pykz/ssoi3QKMI6/gpRk/R1Ij35GCIetJrtl9o4VKYNJP2RsDLXdM6sq0ty2ned4rhTOJAIqJ/VjHNqbL9v+2e29+2KMkrMSGSpLZ7cfSZqhWastpxVYk+i/ulUQnFOIWYl3UAkthxNWEu/LCkjSetIelrSaEm3JRdfphUkrUMoozOzMqpfskKqY2y/QlyIhhN3zpnO8zHRlftu4EYiFjMfMXKiHlmJiBNtSWQN3k3cjFwDrEbLNNsZJe0taZG07mgi5ncb8K/USmlIK8dvWlJM7Xwia/TnHWyeKZCskOqf3xOxlNNTTCLTOY4hXJ7XEZl1A4H/9EYsRTFRdtHkJuwsrxCdxXeyfQ0RQ/qI6EJ9Y1qWIwqndyAatD5i+2rbnxF1WosR9U6fSnpJMWI9E/V8ywL72h5dtDCZtskKqc6xPZmYdjob0fus6Vx33WQ5IqNuT9srEO6aYbU+qaSfE3fidwNPp4aqrW3XX9IvJD2QWkY9lVbdL+kBIuX8EeBD25cQad9b2z6csKLGACtJGpT2W4moyfoeUVg7grCYZqjBy+wzSFqDSHu/2PbNRcuTaZ+skPoAtp8iLjA7ExXnmY55DtgxtdHpT7QG6tbwNUkzSfqVpMsl/Sylxre23SbEzcNSthciRllc3MZhf0MU3x5OtBe6hpg4+wmRWXgZMbPputTBfHYi27AUX3yI6Epxp6RTCYvpFds3OPgrUWvTtHGl5Lq8lLAaDy5YnEwnyFl2fYR0Ub2TuEit5NRYs9lQTGBdg0jdvaKtFO7Uyfl6YhhfqZP3Vl112SXlc3s6303EDcFoYBeX/XgkjSCslgVtH5Kem5mwcKZLpJD0DjEI79X0+DSimHYiMaLiA6IQdjciTlRqJLsx0Vn8SqKL+XKE0hlJuO2Wt/1ZijE9BSycXHpNh6QziVlYm9i+rWh5Mp2g6MrcvHR+IWIEowjF1L9oeQp4/f9HJCecSDQ4vYp2evwR7WJWIpR4t94vIq7zFKk9D+EO+4iKbglEFt/2hJttxvTcDkTNUGvHfY2yzg1E38BDK7bZlIgtLUBYSM8QmXcvEe67ymOeTBTiXkJ0hjig6M+swO/KlukzObVoWfLS+SVbSH0MSd8hpp0eZ/uYouXpLdIoh/eJ0Rzvp9jJk8APbN9bw/NuQFzU1kyP+xPW1lddZqVKMuECvxDYgFA4KwDD3Urdk6T9ibHlpxAjKvYgxiC8X7bNUUS7oiPS4zkJt9xs7ci7PrAE8KTT+IpmI3Xyfpxos7SW7QkFi5TpJF3JAsrUAbYvlvR14GhJ99i+pWiZeolZiLqiDwBsT5T0OpHsUUseAoZIOonIdNuLiEW9mdx05TcFpcr/84hhgI/aHll+MElycJakD4g7+c+Br5Uro8QbwDaSBjlaF21E9BlsE9t3EwkVTUlysV5BzJ/aKSujvkW2kPogKcj9INEsdBXb7/biuVcH/kTc1T9GDJKreTwr9ZR7lLjQn0V0uv4jsKLtD2t0zvmJZIhZiEatCxB33kcT3cgHAQ/bHp/0TJsZkJI2IxTVMGI8+e4dvW/JGruSaHr7GlE4u61jjHqmFdKNwxHE+3t50fJkukZWSH0UScsTd++PAV93WfNPSUNtj63BOeck4hg/BW4m2vnvAazsXugPJ2khYk5QKalhX9sP1ehcixOWxn+Irglb02KhXA/MS8yqmolQju+0pZAkLUp8VrsQsa+fEhmTq7mDH2BSxOsQ7ZMebMWKaloqv+dl/R/Ptb1fcZJluktWSH0YSbsQ7olzbe+X7qhPJS6c61ZbKaW7/J/Z/np6LKLmZt3esJIqZFmBNB+IUMo/tP1eFY//Z+BN28elx4cRg92eJxIldrU9RVKpIPVFt9FFPH1OO9veIT0W4aZbxAWMqm8EJH2XyELcxvaTkr4CPED0L9wwu+r6JrkOqQ9j+29ExtkPJR1MTE89iEhTrsUP8lNgkdSKBcJlOCuR+ddrSJqdsNAuAdYner9dl6yJajEHoXxKPE+kWy8J3FhmEV5PJFqMaOdYI4Fly4pUlyAyAJu+A3cPeIqIgd8raTeidmsMsH1WRn2XbCH1cZJVdAtRnzKFmPNyXo3OJaJ4cwliOut2xBTVERXb9SP67y0MPGT7/irLsRlwlO2Ny+QqZVS9VaVz7Ee4JHcg3te/kdLMiYLW4UTN0FlESvk+7RxLRPHrkoTrbjhwbK0+p2ZB0gLETdhqxGe0frW/a5neJWfZ9XGS2+h7RObXZGIgXa3OZcWAs12IpIZDbN9Qvk26+F5KzCF6ADhc0gm2z6qiKKOB+SUNdAwvnJ2I5VTTRXku0Yz1QcKTcB7RV7AfEcN6i2iG+jqwTXsHSu/b7oQiGgZc1FoqeKbLvAc8QSikv2Rl1PfJFlKDIGkxQgGM/f/2zj3qyrrK45+vSqVCGqnloDJqeWlklNQyIBGvSBrULNFloxPqmHgvMstsLSgdaaFrCNG8DOo4Ol5mqYhTmSULmWDMwgtaiKFJXhISFYmAwdzzx96H9/jy3t/znvOcc/ZnrWcdznOe8zy/c3jX8z3799v7u/EU4hU1Gseh+M17fzPbEON6BhhoZhsqdI0t8MXrAbjL9Ql4N9lvVOL8ZdcRbLLqab1/EJ5l96Jlo7eaEL6BlwOXm9mltR5P0ntyDalBMO9Iejz+q/6BSA2vBTsCS8vE50U8kvhgpS4QAvAP+DTYtsBlwMWVOj+ApItwX7m1km4qNymNOqKXzeyFzsRI0lBJF0g6uT0PvKT7RMR5Oe75950aDyepEBkhNRiR+noPnq48tjwdvErX3xWvFzoZT5u+gJbU8Lr4Y5N0AvA9vGj1DTx54jkzm9TN84wHrsbXnvbDC2ePiWnGpIdIOhaYg6fQH1OpyDupPRkhNRjmvXTOBEYDt0XSQzWv/xK+xnQ9vtYzFrfPqQsxCo4AZkYE9BaeXnx4D84zAzd0PTfevxUe2SU9JKyR7sGz7MamGDUWmdTQgJjZrEiNngaslnRmNQXB3Fl5j5JNTrWuW0H+hEc0JYbghqrtImlyebZhrDOVCokxs3cl/QbIJos9RNJQvOHiS8BoM1td4yElFSan7BqYKNq8BJgJnF+n4lB1wpFiIZ65uApfmzvWzBaF0OwIrC7/dd6WdZCknwBLcZfyIXitzFFmtrg6n6RxkPT3wMPAOmCEmXXo6ZfUJzll19hcijs3nAvMKGWNJR1jZquAT+EC8mu8vmmRpI/jEc+zwCpJZ3RyqlNxV4c38HWksyotRpK2isSJoepey/S6oUyMNuA2WSlGDUpGSA1OiNA0YBIZKfUKSU8AN5vZjBCnx/H6p9ZMie0q3B18LdAP72G0qILjGYA7kO+EJ0y8SYNNZbUSo8PMbFmNh5T0ISlITUArUboRmFgNM9RGQt5/6S/A+0qp3pJuBhaa2Y2tp+wkfQHP1NuAN9XbDve8269SN1VJ03AxmhC7bgD+bGYN0a5b0qdwwV1HilFTkFN2TUBERBfhdRv/DNwRN9ik62zEkx2GA0jaGp/WGyHpvtj3obLj/w6PWJYAn8Gb9r2L/zCoFPsA95nZuyGS98W+ukfe8+th4C3g0BSj5iAFqUmIYs5Lga/jzgZzWhfPxnpEW1NQTU+I+gTgXkmzccuarfGsudvx1uVzy4xnnwP2BR6N947D25HvVsFh/QYYL2nLSO8/MfZ1iqT3S7pK0tOSHpE0vILj6jKSNiuYjlq6n+BF1SPM7IVqjyupDTll14RIOh2f3lmE2/eviP3T6aPWFY2CpMG4l916XIg+ambrY1r0l7jp68PxfD7wadw6aQe8TcYrZnZ2hcayLV4g+jHACF89M1vThfdej/vqXYpHVTPwm//SSoytK8Tf4eV4O/gXY99EfK3zV8AYM3ujWuNJak9DZuUkHRN1Sq/j1jv/K2kMMBJ3VZieYtQ+5n2flkv6CD4F927sN0nv4G7gpecjgX/HI9KNeGJDxSyOzGytpKPw7rXCLZu66qt3IrB3/Bh5IiKkz+Fp6tXiF3ir8QckjcDF8et4S4+TzCzbczQZGSE1MbFo/AA+9bQNbjd0fCY8dE5EQPfjSQuzgCNxq6GDzOwvrY7dFugXrg+FQNKrwNFm9kw8vxuYa2bXVXkcRwIPAiuBnfF2HheY2TvVHEdSDHINqYkxs8eAQ/BIeUvgoRSjrhHrQifi60KTcPPYUa3FKI5dWyQxCi7D1xG/GtN3n8R7PlWbpcAruBhdB5ybYtS8ZISUlNZFbsab/N0AnFdtU9a+JNpV7AC8mcamLUg6HjgaL9y92sw6tEfqg+uPwH3ptga+bWZXV/P6SfFIQUqATZ1nLwO+ibsojzezV2s7qq4R7s/749HKPeXrKJL2A2bjLcm3BL5iZnfWZKAJsGm68yzgB3gixlgzW1LTQSWFIKfsEsA7z5rZt4CTgKHAk7FgXmgkfRf4V9zI9GLg5pJFUjzOxhu4DQRG4BZKe/XBOCZX+pyNSLhL/Ce+VvRz3JYpxSgBMkJK2kDSvrj32ifwqGlKEdeWJO2AR0V7mtnrUay6BBhnZk+GSerzZrZ92XvuAe4ys7srPJbNzFWT9yJpf/zvak88o+773cgKTJqAjJCSzYhfrJ/GU5a/gxd8Dq7tqNpke3xd6HUAM1sHLMen5wBWA0g6IB4HAAcAGyR9WdL4ELFNyJkg6UZJk9sq3Ey6h6QtJJ0HPIp7/40ysytSjJLWpCAlbRKZYRNwc9ChwGJJpxTMMXw5sF7SJEkfkvQl4ON4ASqRrXUG8FDY+zwFPIb7+R0NfAVYEEJV4vu4O/oiYA9gXmvRKhGCZZIsnltsk/vgs9Ylkv4Gd12YgVsBHWBm82s7qqSo5JRd0imSdsfbeA8H7gbOqXZGVntI2gO4BY98ngdON7PHWx3zsXj9ZWAqcJOZ3Rriejuw2Mymhr/f28AgM1sVrz8CXGlmczoZR07ZlRHf3QnAD/Esuq8B16fTfNIR6dSQdIqZ/T5cB75BtPOOKZi72rrBSNXrFBs+Z4d2cswyYFmM7aN49FNyU1gE7BqH9sMteNaUvf4m7iZQM0IoJ9GSSXhFEVwM2vt/llQqcB2HWwCdUk1LoqR+ySm7pEtEFt4VwIHA73HbodmSBpUfFxHLM+EC0SskDZL0Y0mvSVoo743TW+YDl4S56C64+/l88GlK4KfALZIOlnQ+7ls3rwvnnVKBsW1GRBq3AZ/FXTV2w6cg+/XF9bo6Jkk34GnbrfefhieWjMZ/wAxLMUq6SgpS0i3M7GlgGO45djSwJKr9SzfIK4HBePV9j4m6qB/hHVsPBP4NeFDSwA7f2DlfAz6AT80tBW4xs3vLXv8S3jbieuAo4Agz+1NnJzWzyb0cV3vsDBwOfMHMbsfX9PoDB/fR9ToloqL/A86OOq9Svdc83EbpKWB/M5uWrgtJd8g1pKTHxNrMDOBY3NF6Fl4T9G0z+5denns3PCtrUGlaSNJcYKqZPdSrgfu5+gF/LXqmV3wPj+Hfw18jYnoMmFTL5IBIqf8d8GRs5+NZjd8EZhX9e02KSUZISY+JtZnP4WsFA3AxWovbwfSWNXgkMBA2raPsQqRy9xYz21gnN82X8B5HsySNwrMAt8HXZmrJW3j23CjgQjyC3cvMbqyT7zUpIClISa8w5348i20+nhiwWNKVem8H1e6e903gGuCRSKP+GfA0tb8RV5WIDsfhU4xTgB2Bw6PmqurEOtFoPCo6Ge+i+09mdpaZrarFmJLGIafskooSSQ7fA76M/4qeClzTkx5LMT01DneifgG4tYiOEc2CpEPw/9sj8Wy/i4F7M5U7qRQpSEmfEDYxU/Fsq5Xx7+sq9cs+nCMG4k3pNmv5UG0kTe7DxIaaIukgPDobA7yOd3m9tpEc4ZNikIKU9CnyTqRTgCOAPwJXATdYF9pst3M+4Zl8pwKv4WtXY8zst5UZcc9oxMJYScPwJIXj8RYV04CZRaiBShqTFKSkKkRh7WTgMHwq7xpghpmt7OR97ym+lHQcvrA/3MzeknQmcIaZtVn3JG/PfhgepV3XVzfTehGkzoqWQ/DH4EI0AheiH+Ct7d+uziiTZiWTGpKqYGaPmNko3LR1LnAJsFzSLIX5aWskfQB4TtKpZbs/ATxoLR1Y74p9bb3/PGAmXld0MDBf0jYV+UDUl5edpP6SFuDp2e29PhFP3/9vvJbsQmA3M/tuilFSDTJCSmqCpL2BrwKn4GnMC4CrgdlmtiGO+SKeQj7azH4a+z6Pr2EMN7O3JZ0OTDSzg1qdX7gQHWRmy+L5g3hixO198HkKHyGFTRJmdmDZvr2As4EJeBv2x/GI6A7L7rpJlUkvu6QmhJ3MWZK+hd8MzwHuBFZJ+g/gJlysXsNdoks8gK9H/U7SK3hr8jFtXEK4qeercT2L4/v3zSeqC24FpkeSwhDgNHxabiPep2gm8GhmzSW1IiOkpBBI2gK3IjoNT/UuGZ3OA/7RWrVTl7QnnmX32/ZSyiX9F25xMwU3Jr0WOMTMnu+D8Rc6y07SVsAXcdHfCLwPeA5317jVzF6r4fCSBEhBSgqIvBPsRcB5eJRjeBuIO/EpvRVdPM8AYDruJrASt9tZ0CeDLiAhQsOB8XgriB1xMXoWmAgszGgoKRIpSEmhibWmk2LbBxenXwJzgPuBJd25qUZx54eBHxf5Zhyu6SNx89fufL4BeKQ5Frd1Ggisw7+vO/GEkPWVH3GS9J4UpKQuiKSEIcDn8ZttKYnhD7it0M+AuZ05c0u6DRhpZrt2dFytiYy3a4HBZvaHDo7bEv8ujsTdyYfh051v4G7pc3ARytqhpPCkICV1SVgUHQccg7dn2C5eWoxn7JW25a3qmJ4CXjGzthIhejqW6QBmdmEFzzkC+B/gODP7Udn+bfAU9uGxDQO2j5efAH6OC9GCbP2Q1BspSEndE2slB+IRwqHAZ2jJpvsj3iH2cbxPzx14kefFFbz+PAAzO6yC59wOLyC+Fq8N+mRsQ/AICHwt6Bd4FuLDXenblCRFJgUpaThiGmsIHkEcAgwF9qWlEHw9fpNfgt/UlwHLgReBld1dW+qNIEnaGi9CHQzsDuwdY90n9pVYhUdAi/DIb2G6ayeNRgpS0hTEVNdI3ANvHbArftPfpdWh64GXgRV4Zt4KvMXCarwFxJp4XA+8g2etzcSTLc7Bo5et4rE/Xmw6IB4HAh8BdorHnePf5azDRfLZOMfTwC3AS0VOwkiSSpCClDQ1kZW2Oy1Ryt/iIrUTLcLxYbzQtres471CtwKPypaXbS9ng7ukWUlBSpJOiCnA/nikU4p23k9LNDQRnw78IR41lSKnP9MSVa3JdOsk6ZgUpCRJkqQQpNt3kiRJUghSkJIkSZJCkIKUJEmSFIIUpCRJkqQQpCAlSZIkhSAFKUmSJCkEKUhJkiRJIUhBSpIkSQpBClKSJElSCFKQkiRJkkKQgpQkSZIUgv8HotuorPPY0IgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mCI = find_CMDT_CR(Ahat,T_c,mhat12) #generate 201 points along boundary\n",
    "DIc = pmag.cart2dir(np.transpose(mCI)) #convert to dec and inc for plotting \n",
    "X1c = pmag.cart2dir(np.transpose(X1)) #convert to dec and inc for plotting\n",
    "X2c = pmag.cart2dir(np.transpose(-X2)) #convert to dec and inc for plotting\n",
    "\n",
    "ipmag.plot_net(1)\n",
    "ipmag.plot_di(di_block=X1c,color='black',edge='black')\n",
    "ipmag.plot_di(di_block=X2c,color='black',edge='black')\n",
    "ipmag.plot_di(di_block=DIc,color='red',edge='red',marker='.',markersize=2)\n",
    "ipmag.plot_di(dec = 0, inc = ipmag.inc_from_lat(40.15),marker='s',color='c')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}