{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[91mA new version of galpy (1.8.1) is available, please upgrade using pip/conda/... to get the latest features and bug fixes!\u001b[0m\n",
      "/Users/webbjere/Codes/clustertools/clustertools/analysis/profiles.py:37: FutureWarning: all profiles are setup such that the returned radial bins and profile values are in linear space and not normalized by the effective radius. Previously select profiles had unique returns.\n",
      "  warnings.warn('all profiles are setup such that the returned radial bins and profile values are in linear space and not normalized by the effective radius. Previously select profiles had unique returns.',FutureWarning)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import clustertools as cts\n",
    "from clustertools.analysis.functions import meq_func\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load a snapshot of a cluster in file 00000.dat, which has position units of pc and velocity units of km/s in clustercentric coordinates. Stellar masses are in solar units and were generated using a Salpeter IMF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster=cts.load_cluster('snapshot',filename='00000.dat',units='pckms',origin='cluster',ofilename='orbit.dat',ounits='kpckms')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once initialized, any function within ``clustertools`` can be called internally via ``cluster.function_name()`` or externally via ``output=function_name(cluster)``. An internal call sets variables within ``cluster`` while an external call will return the calculated values but change nothing within ``cluster``."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the event that the snapshot is not centered, the centre of the cluster can be found via:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.14946119556654258,\n",
       " -0.22087743738380078,\n",
       " -0.038780116570714403,\n",
       " 0.014544487149392757,\n",
       " -0.097072929729564258,\n",
       " -0.028231116712277449)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cluster.find_centre()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which not only returns the position and velocity of the centre of mass, but also sets ``cluster.xc``, ``cluster.yc``, ``cluster.zc``, ``cluster.vxc``, ``cluster.vyc``, ``cluster.vzc`` Alternatively, if you don't want to set any ``StarCluster`` variables just call the function externally:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.149461195567 -0.220877437384 -0.0387801165707\n"
     ]
    }
   ],
   "source": [
    "xc,yc,zc,vxc,vyc,vzc=cts.find_centre(cluster)\n",
    "print(xc,yc,zc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is important to note that, by default, the centre is calculated to be the centre of density. Similar to NBODY6 and phigrape, once the centre of density is found for the entire cluster population, a centralized subset of stars within an ever decreasing radius are used to find the true centre of density. The parameters ``rmin`` and ``nmax`` set the minimum radius that can encompass the subset of stars and the maximum number of stars within the subset. \n",
    "\n",
    "Setting ``density=False`` will instead find a central subset of stars to find the cluster's of centre by removing stars beyond ``nsigma`` standard deviations of the previously calculated centre. For systems with a large number of escaped stars, which make finding the cluster's centre difficult, it may also help to tell the function where to start. For example, it is possible to tell the cluster to start looking for the centre 1 pc away from the origin along the x-axis and remove stars beyond two standard deviations via:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.47255599444213536,\n",
       " -0.046784905321136626,\n",
       " 0.03368417753435262,\n",
       " 0.0073060547641532469,\n",
       " -0.05167863301943116,\n",
       " -0.0025985741474006487)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cluster.find_centre(xstart=1.,density=False, nsigma=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other functions that can be called include:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Half-Mass Relaxation Time:  50.5231000207\n",
      "Core Relaxation Time:  19.6957365444\n",
      "Lagrange Radii:  [0.73326598977617097, 1.0765583789354674, 1.3646287125197105, 1.7124821982817706, 1.9617869129434311, 2.3496933983866723, 3.1257216188480137, 3.9773141578345004, 5.0211933936198365, 8.3242591110873079]\n",
      "Virial Radius:  2.51937714706\n"
     ]
    }
   ],
   "source": [
    "print('Half-Mass Relaxation Time: ',cluster.half_mass_relaxation_time())\n",
    "print('Core Relaxation Time: ',cluster.core_relaxation_time())\n",
    "print('Lagrange Radii: ',cluster.rlagrange())\n",
    "print('Virial Radius: ',cluster.virial_radius())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All functions allow for projected values to be used instead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Half-Mass Relaxation Time:  33.8815120123\n",
      "Core Relaxation Time:  20.6789810497\n",
      "Lagrange Radii:  [0.49097752209716733, 0.7013820946414554, 0.95595297621692432, 1.2127055798991366, 1.5030278513713957, 1.8280543898615964, 2.3000812878962309, 3.068226381254116, 4.1294671402214904, 8.1792344411375115]\n",
      "Virial Radius:  1.57999162536\n"
     ]
    }
   ],
   "source": [
    "print('Half-Mass Relaxation Time: ',cluster.half_mass_relaxation_time(projected=True))\n",
    "print('Core Relaxation Time: ',cluster.core_relaxation_time(projected=True))\n",
    "print('Lagrange Radii: ',cluster.rlagrange(projected=True))\n",
    "print('Virial Radius: ',cluster.virial_radius(projected=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The core radius of the cluster can be calculated by finding where the density drops to 1/3 the central value. The parameter ``mfrac`` sets the mass fraction of the cluster that defintes the central region (default 0.1). When ``projected=True`` the core radius is where the density drops to 1/2 the central value. ``rcore`` also has a plotting feature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.865233682914\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAD6CAYAAAC/KwBlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvVklEQVR4nO3deXyU9bX48c9Jwr5kEnbZE1pRUSEEkGqrKCjan3Upi9baapXE1nu9XRS0v/5u7W2vCGqr3rYaQKttrSJgrdxWlOBSW2VJgihWFElcQAUkCSjImvP743kmDJOZyfbMPLOc9+s1ryTzbAeFOflu5yuqijHGGBNNlt8BGGOMSW6WKIwxxsRkicIYY0xMliiMMcbEZInCGGNMTJYojDHGxJTjdwAm9YjIEmByyFu1QL37fUHI+9Xu1wCQH/q+qo4NuV8AKFDVKq9jNca0n9g6CtNaIrIFJzHMCv9wF5HZwDxggaqWhh0rAMqAYlXNC3m/EigCpqhqeZzDT0vuf/dCYAvQC0BV5yTyfnZNGlNVe9mrVS+gDqcFEOnYNECBeVGOB5y/dse8Ny/WPe3V7P+PlUBJ2HuTgcpE3c+uSe+XjVGYtqhQ1ermT2tKVeuBpW53U/C9Oaqa19Z7ZjIRmQegqgtC31enZVbhdhPG9X52TQbwO1PZK/VeRGktuMditijcc2ZjrQev/l8oMC3KsSL3eCCe97Nr0v9lLQrTKu44w5Z23qYapwvKtIOITHO/jTYJINhCK4nX/eyazGCJwrRWLVDRzntUcXSWlGm7KQAapctOnW6+xvPidD+7JgPY9FjTKu4/lHZNYw39R+j29xa4r1mqujTk2Er3/fzgMREpwWmN9MLpBigLXiMiRRydtlvoPuuYmVeh3HGSeThJa5d7TaWG9U0nsYLmT6G6hee19X52TQawRGH8NguYgTNtNtx0nCb+PGgcZJwb/K0u2A0mIsHf7gKqOj94sYgsEZFKDVmzEXKsCFiCMyU3NHGtFJGxsRJMEslv/pS4n2fXZADrejK+cj/0H49xLNjCKMVpPdSHHK/Gad3MAYpCWyOuMqDITSjhVuEMuId3MZQCJSIyOcI1ySZAy7rwAnG8n12TASxRmGRX636tj9JvXI3T3RSeJILHwOmiauS2TAKRupjCkk+LichsEVEPXoFWPLYlv9XWx/l+dk0GsERhUkW0NRa1EH3wMYqSGPcLPqu4FfdDVeerqnjwqm/Nc1sgkAT3s2tSnI1RmFQRa0pufUtv4v7GHgBq3YHxSILnJLtaWjaoWh/H+9k1GcAShck0wQ+B8hizm1Jl1lN9C87J52j3XTzuZ9dkAOt6MhlFjxYxTIdZKxXQ2EqKJkDLpzO35X52TQawRGEykefz4H0azK50v0b8s4Tca2Uc72fXZABLFCYTLSVsJlS4GOMXEfk0mB2cVhxt4D04xbelpdvbcj+7JgNYojBeS3SXTqufp86eAtXBKqHhQga8k5qbVObjrP2IZCbOviBNZni5Cw7bfT+7JkP4XZXQXun1wqkMqziL41p6TYF7TUmMY7OjXLsSqItyrKiZ+24hrEooToKIWvk2GV843SXheydMc/98gdb+P2rt/eya9H/ZDnem3dydwGbifPgGQg7V44wHVGiUkhju7nah11UDY1W1PsqxUlUtd59ZytG+5HqcmUzT3RbBqgjXTtemO/LNc88L/oa4S0PKgKSK1uzG5q46XwLM0Sgzv5J5F7lkviZdWaIwxhgTk41RGGOMickShTHGmJgsURhjjInJEoUxxpiYLFEYY4yJyYoCJlDv3r112LBhfodhjDFNVFZWfqKqfSIds0SRQMOGDaOiosLvMIwxpgkReS/aMet6MsYYE5O1KOJIRC4ELgz+PGLECB+jMcaYtrFEEUequhxYHvy5uLh4lo/hGGNMm1iiSHa/+2rT9066GMbPgoP74JHpTY+P/gaMuQL27oLHv9X0+LjvwKivw+6t8ESEEkxf+jc4/nz4ZDMs/37T41+5EQonwUevwYpbmh4/5z9hyAR4fw2s+q+mx6fOhQGnwJbn4e93Nj1+4d3Q+wvw1tPw8q+bHr+0DHIHwcZlsO7Bpsdn/B669YL1j8Crf2p6/Iol0LErrF0IbzzZ9PjVf3W+/vNeePuZY4916AzfXOZ8/+J8qH7x2ONd82DmH53vy2+FD9Yde7zncfD1hc73T98MH79+7PFehfC1e53vn7oBdoXtANv/ZDj/duf7ZbNgz4fHHh88Dibf6ny/+Juwr+7Y4wVnwpmzne//+HU4tP/Y4188D06/wfne/u41PZ4qf/c8ZmMUzRCRAhEpEZGM2qjEGGOCrChgDG7N/mKcCqQzVXVse+5XXFysrZn19PxbO3h+0w7GD89n/PB8+vbo3J7HG2NMVCJSqaoRN2yyrqcY3JLUVSIyzY/nV+/cy9LKrfz+FWfW2vDe3Rg/LL8xcQzO7+pHWC129YqrAfjd1N/5HIkxpj3SLlG4exHMA1aq6tIY55XgtBTq3a/Vsc73wzVnDOdbE4fyxod7WFuzi7U1tTy98SMWV3wAwMBAl8akMX54PgW9uyEiPkdtjEk3aZMoRKQMZ1vMdUAJMTY/dzerOWaDGhGZJyIFybZpTYfsLEYPDjB6cICSrxTS0KC8tf1T1tbUsramlpc2f8Kf128DoHf3jowLaXGM7N+T7CxLHMaY9kmbRBG6g1q0vZDdYwU4WxzmhV0/R0TqRGSBtm6D+4TKyhJOGNCTEwb05NtfGoaqUvPJ3sbEsaamlqc3fgxAj845xySOkwfm0iHb5i8YY1onbRJFK5QC5VGOVeO0RpKqVRGLiFDQpzsFfbpz2fghAGyt28e6d48mjuc27QCgS4dsioYGGD+sF+OH5zNmSIDOHbL9DN8YkwIyMVFMJnaimEIKJYpIBuV1ZVBeVy4ZMwiAnZ8eOCZx3L3qbVShQ7Zw6qBAY4tj7NA8enTu4Fkc5w07z7N7GWP8k4mJoghYHOVYsEWRVvr06MQFJw/ggpMHALB73yEq3juaOMr+Xs1vX9hClsBJx+U2Jo5xw/LJ79axzc+9bORlXv0RjDE+ysRE0ZxAhPfyEx1EPOV27cA5J/TjnBP6AbD3wGHWv1/P2ppdrKmp5Q+r3+OBf9QA8MV+3d3E0YsJw/Pp17Plazk+P/w5AF1yunj/hzDGJExGJQp36iw4U2Jbcn4BzpjGZKDIHSTfoqoL4hKgT7p1yuGML/TmjC/0BuDA4SO8tnV3Y4vjz1Xb+OPq9wEY2qtr41qOCcN7MTi/S9Qpud8r/x5g6yiMSXVtThQi0hPnA7ReVZ/zLqTkoarVwJzWXueu0WjShTVkyBAvwoq7TjnZjBvmdD1dPwkOH2ngXx/taUwcK9/czpLKrQD079m5satqwvB8RvTtbms5jEkzzSYKNyEsBcbirE0oAQSowe2mEZFK4BxV3RO3SD2gqvXuh1ggzs9ZADRpdRQXF6dkvZSc7CxOGRTglEEBrv1yAQ0NyuYdnzmLAN+tY3X1Lp7a4BSny+/WkXHD8hg/vBd7Dxyma6eMarQak5Za8q94Ec4g70pgHLAEqML5TbscyANm4KyG/m58wkyYAC3slmqJdN2PIitLOL5/D47v34MrJzprOd6v3ccady3H2ppannljO12G7CYnS7jqd2sbWxwnDwzQMcfWchiTSlqSKKpV9ebgDyLydWC4qrq1kqnBqYd0UzwCjINyoDDKsXzAs71KM2U/ChFhaK9uDO3VjRnFgwH4aPfnlKx8hE/3H2bbjs+Z/9ZbAHTKyaJoSF5j4hgzJI8uHW0thzHJrCWJ4p3QH1R1mYhE+sBLlW6VKpwpspEUAGVePShdWxQtMSC3C9eMngHAxSPOZNdnB1j3bp3T4nh3F//z3GbuUcjJEk4ZlNs4q2rssDx6eriWwxjTfs2WGXeTwjpgsqreGXasJ5Cvqu+KyE2qekf8Qm05EVFgeqQif+5MpsrwEh7ujKg6IC9eJTxaW2Y8ne3Zf4jK9+oau6pe21rPoSOKCJw4oGdji2PcsHx6de/kd7jGpL1YZcZbkihyccYfhqvqeWHHbgduwlnJfLuq7vYm5PZxE0VptGmswVpQqjon5L0ynKmvcVuVnWmJom6/s7taXue8Zs6Ezw8eYf0HRxNH1ft17D/UAMCIvs5ajrFD8sjr1oGO2dl0zMlyXtnO105hP3fMySInS2wGljEt1K5E0cyNc4ECVV3f5pt4RERm4wy2F+F0IdXjjDdUA3PCWwlhZcYLgXVelxmP0PU0a/PmzV4+Iqm1Zz+Kg4cbeH3bbjdx7KLi3To+PXC4VfcQ4dhEEpJEjk0q2XTMjpxsOuZkMaJPdyaf0I/crtYlZtJX3BKFe/OeOC2K6cBtqnpX2PG5wNxknzqbCJnWovBy46IjDUrNJ5/x2YEjHDzc4LyOON8faPy54eixkJ8PxDh28HADBxq/P9LkvAOHGzjcoORkCRMLe3H+qAGce1I/elt3mEkz8d7h7hZVvU5ECoHLgLvCjt8O3OK+jGmT7CxhRN8eCX+uqrJh626e3vgRKzZ+zI///Do/efJ1ioflc/6o/kwd1Z8BuVaixKQ3LxJFLYCqTol0UFV3i8jjInJ2uq7gjiaTZz2lCxFp3Djq5qkj2fTxpzy98WNWbPyIny3/Fz9b/i9OHRzg/FH9OX9Uf4b26uZ3yMZ4zouup/tUtdmFdiJyY/isqUxjXU/pZcvOz1ix8WNWbPyY17c58zhOGNCTqSf15/yT+/MFK2diUki8u55ERCap6vPNnefBs0wKmXn8TL9DiKvCPt25ftIIrp80gq11+xqTxt2r3uZX5W9T0KebkzRGDWDUwJ6WNEzK8qJFkQtUAveFD2SHnZdxLYpMn/WUqXbs2c8z/9rOio0fsbq6liMNyqC8Lkw9yRnTKBqSR5btZW6STFxnPbkPKMKZiloBlKnqAxHOuV9Vr2v3w1JYpnU9fbzX2bu7f7f+Pkfin9q9Byl/czsrNn7MPzZ/wsEjDfR1N5IqPbPABsJN0oh7onAfUgQsBMbglPOo4mjdpBk4K6UzajA7XKYlinQfo2itT/cf4rlNO1ix8WNWvbkDBK760jC+e2Yhee3YSdAYL8R7jAIAVa0CxroL2abjlCUfi5MwZmR6kjCmR+cOXDR6IBeNHsjWun3cXb6ZRS9V8+ia9yn5SgHfOWM43awsu0lCntd7VtUFqjpFVfNVNUtVi1V1ldfPMSaVDcrryp3TT2XF97/CxMJe3LXybc6843kefvldDh5u8Ds8Y47h2a8vIjKakNIZthLb1lGY5n2xXw8WfKuYqvfrmPf0Jn761Bss+kc1P5zyRb526kCybdDbJIF2tyhEJFdEKnC6mJbibHBUJyJvi8gl7b1/KlPV5apaEnzl5ub6HZJJUkVD8nis5DQe/s54enbuwA8Wb+Cr977Eqje349U4ojFt5cX02MdxypCXu28Vu6/pOEX3KlR1fLsekiYybTD7hQ9eAOCswWf5GUbKaWhQ/vr6R9z17Fu8u2sfxUPzmHP+SMYNy/c7NJPG4l0U8PbQHfDCjhXhlChvCC9RnokyLVGY9jl0pIHHKz7gnvLN7Pj0AGeP7MtN5x3PCQN6+h2aSUOxEkVcNy9W1Sq3BlS5W0XWZJCa3TXU7K7xO4yU1SE7iysmDOXFmyZx8/kjqXi3lgvufYnvP7ae93ft8zs8k0G8aFGcA+Sq6hPNnBe15ZGuMn1ltq2j8NbufYco+/sWHvxnDYePKJePH8K/nzOCvj06+x2aSQNxbVG4U18LROQ2EUl8HegkZoPZxku5XTswe+pI/n7TJC4bP5hH177POXe9SOV7tX6HZtKcF7Oegluh3gzUi8hmEblPRC5xNzVCRC4F3mnvs4wx0LdnZ35x8ck8+4Ov0Lt7J658YC0vv/OJ32GZNObFGEWh+xoLfBdYD8wEluFMkz2CM6CNiAwLvVBEnvHg+cZkpII+3VlcehqD87py1UPrWPXmdr9DMmnKi0RRD+xS1fXuquwZqpqPkzy+CyzCKTG+ANgiIkdE5Bl3cDtif5gxpmX69ujMYyWnMbJ/D0r/UMn/vvah3yGZNORV9dibgHxgXbRBbbcc+RT3NRYoAlRVs9sdQIrItOmxr3z4CgATj5vocyTp79P9h/jOQ+uofK+O279+CjOKB/sdkkkxCake6z5ouKq2aD6kiARwFuNlTF2LTEsUJrH2HTxM6R8qeWnzJ/zsayfx7S8N8zskk0ISto6ipUnCPbced+zCpKdNtZvYVLvJ7zAyRteOOSz6djHnntiPnz71Br953uaPGG+0qCigiJxNHAr9qepCL++XbDK9KOC8tc7vAbaOInE65WTzmyuKuHHJBu545i32HTzMjeceb9uwmnZpafXYc4EFIlIJLG5ucZ1xqOpyYHnw5+Li4lk+hmMyRIfsLH45YzRdO2bzm+e3sPfAEf7z/5xo26+aNmtR15Oq3uyOJSwAzhWRXe5aibPjG54xpi2ys4TbLjmZa84YzkMvv8vNT7zGkQarQmvaplVjFKq6SlWvU9VeONVir3MX2M1196MwxiQJEeEnXz2BG875Ao9XbOU/HlvPoSO2KZJpvTZvXKSqy4Bl7rTXGcB8ERkOlAFLVfVdb0L0l4jMxlkrAhBQ1fk+hmNMq4gIP5zyRbp1zGbu05vYf+gIv/5GEZ07ZMysdOMBr6fH5gIlOCuzFSdpPJ6qu90Fk4SqLnB/ngxMV9XSttwv06bHvrrjVQBG9x3taxzG8YfV7/H/ntzI6SN6sfBbxXTtaPtzm6MSto4i7KHDgVJgGrAFKEu1QXARqQOGu1N5g++pqrZpVDDTEoVJPssqt3LT0g0UDcnjwavH0bNzB79DMknCl/0oVLUmZBD8ZpxB8FoRWRzPQXARCYhImYhMa+a8EhGZHfJ1WtjxApyupvoI1072Nur09OqOVxtbFSY5fH3sIH7zjSI2bK3nGwtXU7v3oN8hmRQQ142Lgtw6UNe5NaAW4AyCB2dOjfbiGW5yWILT9VXSzLnzcMcb3PpU84FxbldTUEGUy+txtng1zbin6h7uqbrH7zBMmPNPHsCCbxWzeftnzCx7he179vsdkklyCUkUodyZUzNCZk61qb8/wn1LVXV6c4PNbkuhJPw8VZ0D3OKWFomlFqeulTEpa9LxfXno6vF8WP85M8peYWud7Zhnokt4ogilqstU9bsJfmwpToKKpJpmWiNYkjBpYmJhL/547QTq9h5k+v2vsOGDer9DMkkqrokifP+JJDEZJyFEUo1T3ZYY5wRiHDMmpYwZksdjJRPJEmHa/S+z6KVq4jXBxaSueLcoCkXkfneHu2RRBOyKcqwad48MVa3G2bGvyViFqkZrkRiTck48rid/u+HLTDq+L7/465vM+n0F9ftskNsc5UmiEJGzReRGEbk0tBURXMkN7BaRa714VgIEQr6fi9MCAcCdGbUg0QGlqjnj5zBn/By/wzAtkNu1A2VXjuXWC0/k729/wgX3vETFu7YXt3F4sWf27Th9/vOBpTi72O0Skd+KyKngJAySYKZQyEB1fUvOdwe8A+4U2hJgXEsW27nnV4S/du7c2ebYU9HI/JGMzB/pdximhUSEq04fzrLvfomc7CxmLljNb194hwarEZXxvFiaWQDkqepud5FdEU612cuAUre8cbCr5k4PnpdQbSnZ4a7kbtLyKC4uzqh/cbbDXWo6eVAu/3vDGdzyxOvMX/EWq6tr+eWMU+ndvZPfoRmfeJEo1qnqbmjcuKgGWIaTJCbjDA7nkgSbFKlqvZu4Aol4XqbvR7HgNSdXWqJIPT07d+DXl4/h9MLe3Lr8DS645yXuuWwMEwt7+R2a8YEXYxRR/+aoarmqznEX27V49zsfBWhht1RLqOpyVS0JvnJzc726tTFxJyJ8Y8IQ/nL96XTvnMMVi1Zzd/nbVq48A3mRKBaLyI0e3CdRyoHCKMfyASvGZEyIEwb0ZPm/ncHFowdyd/lmvrloDTtsNXdGaXeiUNX1wHoRuc+DeBKhiujlOQqAJV49SEQuFJEFwdfu3bu9urUxCdWtUw6/nDmaO6adwqsf1HP+PS/x97cza3JGJvNi1tMsYCXOmMQREXlGRH6UxBsZleGulQjlzogqAh736kHW9WTSzfTiwTz1b6fTu3snvvXgWuav2MRh2wwp7bW7zLiIPI7z4VuA8wE8FucDN3jjcuBZYFmiNjMSEQVKg/tIRDg+DxrrOwXfKwO2eLkxUYTB7FmbN2/26vZJr2a3Myw1PHe4z5EYr31+8Ag/W/4Gj637gOKhedx7+RiOC3TxOyzTDnHdj0JEblfVm8PeywXG4cx4OgcncdS5hQDjwq38Os59VgHOoHQFzmrrOeHlwt11EQH3vEKc2VtL4xUf2H4UJv385dVt/PiJ1+mQk8Vd00/lnBP6+R2SaaN4J4pc4BbgE2JsgSoiw1Nk5lPcZFqieOGDFwA4a/BZfoZh4qzmk71c/0gV//poD9eeMZzZU0fSMcfXeqOmDRK2w52IjHEHtw3W9XT1iqsB+N3U3/kciYm3/YeOMPdvb/LwK+9x6uAAv758DIPzu/odlmmFhO1wF54kRGRWEg9qx50NZptM0blDNj+7aBT3XVFE9c7PuODel1ix8SO/wzIe8WLW06XRkoGqLgR6icjt8dz+1BiTHM4/eQB/u+HLFPTpznV/rOKnf9nI/kNH/A7LtJMXLYofA5Xu1Nh1InKbiEwKHnQryN6MM8hsjElzg/O7sqR0IrO+PJyHX3mPr977EoteqmbHp7ZIL1W1u9aTqhaH1HQaA9wM3OxOUa3GWeBWTfRFbmkr02s9mczVMSeL//vVE5lY2ItfrdzML/76Jrf97U3O+EIfLhlzHOee2J9unbwoNWcSwdPB7MabihTh7ONQ7H7dBUxX1Vc9f1gKybRZTx/v/RiA/t36+xyJ8ds7Oz7lyfUf8uf129hW/zldOmRz3kn9uHjMQM4Y0ZucbJsl5beEzXqKEcDtwG2quifuD0timZYojAnX0KBUvl/HE1Xb+OtrH7Jn/2F6d+/E1049jkvGDGTUwJ64FZ5NgvmeKNwg5qrqLQl5WJLKtESxomYFAFOHT/U5EpOMDhw+wvObdvLk+m08t2kHB480UNinG5eMGchFowfa9NoEi5Uo2t1JKCI9cRbc7QLKY3QvRdunOm1l+hjF4rcWA5YoTGSdcrKZOqo/U0f1Z/e+Q/z19Y94cv027nz2be589m3GD8vn4jED+erJA8jt2sHvcDOaFyuzn8UZqC7Aqe9UjVOBtRyoUNU9bjJZqKoz2xlvSsu0FoUtuDNt8UHtPp7a8CFPVG1ly869dMzOYtLIPlwyZiCTRvalU0623yGmpbi2KICVqnqHW8pjCjATuA6Y4z48eF6ze00bY8zg/K5cP2kE3zurkI3b9vDn9dt4asOHPPPGdnp2zuGrpzjjGcVD88jKsvGMRPBkfpqIDHNrPC11X4jIOTitjABOl5SV9jDGtJiIcPKgXE4elMuPLxjJP7fs4sn123hy/TYeXfs+AwNduHjMcVwyZhAj+nb3O9y05slgtrsnRS6xxygynnU9GdN+ew8c5tl/fcyf13/IPzbvpEFhRvEgbj7/BPK7dfQ7vJRlRQGTRKYlirr9dQDkdc7zORKTrnZ8up8HXqrhgX/U0L1zDnOmjmRm8WDrkmqDpJgem4kyvXqsMYny9vZP+cmTG1lbU8uYIQF+cfEoTjrOinC2Rrz3o+iJU+/pHGBWsOtJRIYBJGpXu1SQaS2KJ995EoCLR1zsaxwmM6gqf16/jf/+65vU7TvIt780jB9O+SI9OtvU2paId5nxW4C5QCVuPScRuR/YAmxxCwWe6sFzTIr5yzt/4S/v/MXvMEyGEBEuLRrEcz86iysmDOWhl9/lnLte5KkNH2I9J+3jRaIQVd2tqtep6hPubKcSnN3usnGmy97htjyMMSaucrt24OcXj+LJ751Ov56dueHR9XzzgTVs2fmZ36GlLC8SRXiqnuK+NwdAVauB6TgtD2OMSYhTBwd48vrT+flFJ/Ha1t1Mvfvv3PnMW3x+0PbHaC0vEkVeWGthGlAfOjahqrsBm4ZgjEmo7CzhyonDeO5HZ3HhKcfx6+ffYcqvXmTVm9v9Di2leJEo5gDPiciNInIfzjjF4xHO+8SDZxljTKv16dGJX84czaOzTqNzh2yuebiCkt9XsLVun9+hpQSvFtzlAvNw9p6oUtUZ7vs9g6XFrXps5s16+vzw5wB0yenicyTGHHXwcAMP/rOGe8qdqeo3nPMFrjljOB1zMntPDN/WUYhIBdATOBdn46I74vawJGTrKIxJXtvqP+e/lr/BM29sZ0Tf7vz8olFMLOzld1i+8TNRTMOZAfWsqt4ZtweliExrUTy26TEALht5mc+RGBPdc5u289On3uCD2s+5ZMxAbrlgJH17dPY7rISLd/XYqFS1sUigyTzPvPsMYInCJLezR/ZjYkFvfvvCO9z/4hbK39zOTecdzxUThpJtpUCAVgxmB1daG2NMuunSMZsfnXs8K77/FU4ZlMt//uUNLvrNP9jwQb3foSWF1ozelMUtiiQmIgUiUiIiK/2OxRgTX4V9uvPHaybwP5ePYceeA1z823/ykydfZ/e+Q36H5qvWJIqIfVfpTESKcGZyBYB8f6MxxiSCiHDhqcex6kdncvWXhvOnNe9z9l0vsKxya8aWAmlNosgTkbdF5DYROTsTSnKoapWqLsDZ3tUYk0F6dO7Af154Isv//QyG9urKj5ZsYGbZat76+FO/Q0u4Fs96EpE6nNIcAY6W7agCKoCVOJsW7Ylx/ejWbmokIgGc9Rkr3YHxaOeVuHHVu1+rY53fWu7srVtUdWx77pNps56MSRcNDcqSyg+Y+/Qm9h08woIrx3LW8X39DstTXlWPLcKpEjsdZ0/sJ4BCnL2wlwJ1IrJZRH4rIpdGGPye14qAy0RkCc7U2pJmzp0HBFR1vqouUNX5wDgRmd3S5xljTCxZWcLMcUNY9cMzGdGnOyV/qOQfmzOn2ESLE4Wq1rgL5uqBXFWdrqr5QB4wA1iEU8/pOmAJTonxXSLymIhci1uCvIXPKnXvPz/WeSJSAJSEn6eqc4Bb3BaJ8clDGx/ioY0P+R2GMZ7p1b0Tj1w7gYLe3bjm4XW8vCUzkkWr11Go6ipglbtP9i5VfQKnRbEUGst5TAbGuV9nuK94jAKVAuVRjlXjtEbmu3FNwyl53pzFXnZbZbIXt74IwFWjrvI3EGM8lNetI49cO4HLF67mmocqeOjqcUwoSO8V3W1ecKeqC0VkuIjcBCwJVot1K8Uuc19A44d0PKbXTiZ2opiCmyhs8Z8xxitOy+I0Ll+4mqsfWsfvvzOe4mHpOzGyXVWwQrqjxrrdS9HOWwrUtOdZURQBu6IcqyYDp/QaYxKjT49O/GnWBPrndubbD66l8r06v0OKG0/KJarqMlVdJCKzROTsKKct9uJZrRTw6D7p+6uCMabN+vbozKOzTqNPj05c9eBaXk3Tldye1tVV1YVAjbs3xbCwY55Wjg0ZqK738r5hzyhwZ1WVAkUiMs+ditvcdSUiUhH+2rlzZ7xCTUqdcjrRKaeT32EYE1f9enbm0ZLTyOvWkSsfWMPrW3f7HZLn4lY9VkS+DuSp6qJ23kdxSpQvDXs/ANQBpe6iuPDr5gGzVTVpqnrZOgpj0te2+s+ZWfYKn+4/zCPXTmDUwFy/Q2oVr9ZRtIqqLsOZHXWfiEyKw/3r3W8DXt/bKyJyoYgsCL52706/3zSMMY6BgS48Ous0unfK4ZsPrOHNj6KuP0457U4UIjLMLelxrYjMFZFn3IV3R4B3cKewikiixygCxLFbqiVUdbmqlgRfubmp9RtGe92/4X7u33C/32EYkzCD87vyp1kT6NIhmysWrUmbch8tnh7rzmoqxFk4F3wFQk8J+b4KZ3psNbDF/RqPeknlbkyR5OOUF/FNhB3ufIwm8dZ8tAaA6069zudIjEmcob268adZp3HZgle4YtFqHis5jRF9e/gdVru0Zh3FApxFc4LzoV9J00RQ7a6jSJQqnCmykRTgc2l0VV0OLA/+XFxcPMvHcIwxCTK8dzBZrObyhWt4rOQ0Cvt09zusNmtt19MqnK6kKap6rqpep6p3uNNj1yc4SYCTCJoMvrgD3UXA4wmOxxhjAGdvi0dnTUBVuXzBamo+2et3SG3WmkRRparn4nT3THYHqe9zp8KObu7idpYlj7iOQVWrgQXuDKdQ84A5IQPevrDBbGMy24i+PXjk2tM43OAki/d2pWayaE2Z8ftVtUlnc1htp1ycrqgqVX0u7Lx1qjquhc+a7d6vCKcLqR5nvKGaCAkgrMx4IbAuGes1Zdr02B88/wMAfjXpVz5HYoy/3vxoD5cvXE3XDtksLp3I4PyufofURKzpsZ6vo4iQOOpxxjVmqWp6V85qRqYlCmPMURu37eaKRWvo0TmHxaUTGRjo4ndIx0hooojw8Fzgx8CNqpod14clmQiznmZt3rzZx4iMMX56fetuvrFoNXldO7K49DQG5CZPsvA1UYQE8Y6qZtb80DCZ1qK4u/JuAL4/9vu+xmFMMnn1g3quXLSGXt07srh0Iv16dvY7JMCnldkRJN2YgYmvDTs3sGHnBr/DMCapjB4c4OFrxvPJZwe5fMFqduzZ73dIzUpYolDVmxP1LGOMSWZFQ/J46OpxfLxnP99YtIadnx7wO6SYPEkUItLTnSa72H3NbcmU2XRn02ONMdEUD8vnd1eNY1vd51yxaDW7PkveZOFFracxODObfoyzo9x0YA5Q6e6Xndpr19sh02s9GWNim1DQiweuKub92n1csWgNtXsP+h1SRF60KG4GClU1X1XzgTychHEnzqrpehE51YPnmBTTr1s/+nXr53cYxiS1LxX2ZtG3xlHzyV6+uWgN9fuSL1m0e9aTiNwUa1Mid7/sucBYVU2furttkGmznowxLffi2zuZ9XAFx/fvwR+vnUBulw4JfX68Zz3F3CbUXSF9Hk5ZjYxiYxTGmJY684t9KLtyLJs+3sO3HlzLnv2H/A6pkReJokJE7ot1gluTKeM+JTN9jGLe2nnMW5txvx8Y02aTRvblvivG8q8Pd/PtB9fyaZIki3YnCncnuywRWScil8Q4dVd7n2VSy6baTWyq3eR3GMaklMkn9uN/Li/ita27ufp369h74LDfIXkzPVZVS4ElwDIR2SUivxWRS4MVY0XkUpxigcYYY5oxdVR/7r1sDOs/qOfqh9ax76C/ycKzBXeqOh8YASwCLsNZiV3nbol6C7S71LgxxmSMr54ygF/NHE3Fu7Vc81AFnx884lssnq7MVtVqVZ3jTpMtBL4LPOF+H0wcm919LC6xxGGMMdF97dTj+OWM0ayu2cWs31ew/5A/yaI1W6G2iqrW4GyfugBARIbjrK+YDMzE2SmvDkjb0uOZvmf20J5D/Q7BmJR38ZiBHG5Qblq6gdI/VFJ25Vg6d0hsIe6EVY9t8mAncQRUdb0vAfjA1lEYY9rq8XUfMHvZa5w9si/3fbOITjneJotkqR57DFWtyaQkYYwx7TFj3GBuu+Rkntu0g+sfWc/Bww0Je7ZvicKkv1tfvpVbX77V7zCMSRvfmDCEn190EuVvbuffH63i0JHEJAtLFCZu3tvzHu/tec/vMIxJK1dOHMZPLzyRZ97Yzn88tp7DCUgWcRvMNsYYEx9Xnz6cIw3KL/76JtlZG/jVjFPJyY7f7/2WKIwxJgVd++UCjjQoc5/eRLbAXTNGk50lcXmWJQpjjElRpWcWcrhBueOZt8jOyuKOaaeQFYdkYYkijjJ9HcXI/JF+h2BM2rt+0giONCi/XPk2OVnC3EtP9jxZ+LaOIhPZOgpjTLz88tm3uPe5d7h8/BD+++JRrU4WsdZRWIvCGGPSwA+mfJHDDcpvX9hCp5wsbv3aSZ7d2xKFMcakARHhpvOOB+CUQQFP722Jwhhj0oSIMHuq92ODtuDOGGNMTJYojDHGxGSJwhhjTEyWKIwxxsRkicIYY0xMtuAugURkJ+BlOdVcYLeP92jrta29rjfwSRuek+m8+Pvhh2SIOxExxOMZ7bnnUFXtE+mAJYoUJiILVLXEr3u09drWXiciFdFWjJrovPj74YdkiDsRMcTjGfGK27qejDHGxGSJwhhjTEyWKFLbcp/v0dZrvYjbNC9V/zsnQ9yJiCEez4hL3DZGYZKejVEY4y9rUZhUsMDvAIzJZNaiMMYYE5O1KIwxxsRkZcZNShOR2UC9+2NAVef7GI4xaclaFCZlBZOEqi5Q1QVAlYiU+R2XMenGxihMyhKROmC4qtaHvKeq6u3O8sZkOOt6MgklIgFgHrBSVZfGOK8ECOB0KwWA6tDzRaQAp6upPsK1k1W13Mu4jclklihMQrhdQvnAOqAEWBnj3HnArtDxBhGZJyIFIe8VRLm8HiexGGM8YonCJISqlga/dxNBRG5LoURV88KunyMidW7Rs/oYj6rFSUjGGI/YYLZJNqVAtG6japzWSCyWJIzxmCUKk2wm4ySESKqBKSHfRxKIccwY0waWKEyyKQJ2RTlWDRQDqGo1UO92VR3DBrKN8ZYlCpNqAiHfz8VpgQAgItOwulDGeM4Gs03ScKfOwtGV1jGp6nwRme1OpQUoDB00N8Z4wxKFSWlWssOY+LOuJ5M0Qqa9BnwMwxgTxhKFSSUBWtgtZYzxjiUKk2zKgcIox/KBigTGYozBEoVJPlVEL89RACxJYCzGGGww2ySfMqAy/E13RlQR8HiiAzIm01mLwvglYqkNdyHdggj1oOYBc5qp82SMiQPbj8IkhLvJ0DicVkEBzqB0Bc5q6yYJIKzMeCGwLlZZcmNM/FiiMMYYE5N1PRljjInJEoUxxpiYLFEYY4yJyRKFMcaYmCxRGGOMickShTHGmJgsURhjjInJEoUxKUBEpoVs7BSvZ5TE+xkmNdmCO2PiRETKcPb4LnLfqgJqw06rBua5pUui3Wc2UK6qVXEJ9NhnldkugSacFQU0Jk6CH7giokC1qo4NP8dNAltEpFRVm+z3LSJFQK9EJAnXPBFZoqrTE/Q8kwKs68mYOHI/6MHZZ6MJdyvXaqAsSrfPPFWdE6fwIsVTDVSHxG2MJQpj4myy+3VljHOC3U6TQ98UkWk43VWJVgYs9OG5JklZojAmvqa4XyO2KFzBjZrCk8JMYLHnETUjOF4iItE2kDIZxhKFMfE1GWd8oj7SQRGZjJMoyiMMaBclcGwi3GJgmk/PNknGBrONiZPmxifcMYkynCQxJezY5GjXucfn4ezXsUVV57vnB+9RAMyKlJzcZwY3hdoC9ALmRji33D1vfrQYTOawRGFM/EQcn3C7dKbhdC3NizTbCedDf0ukm7ozpebi7BK4RUTAabXMCTm+Chgbdl0Rzp7jpapaHvJ+GXDMlFhVrXKTjzGWKIyJo+Bv+ONC+vsLcRJIVaTpsiEKgHVRjvVS1fqQe/ZyZ08FleNMc23sunKTRCXOboKhSWIlzlqPSGsn6kUkYNvPGksUxsRPcHyiyfRWEVkpIltUtTDKtQGcbWDDrwtNIMXu17lR7lHA0QHyhTjJKbwrKT/G9bXu8SZxmMxiicKYOGhufAKnC6hMRCaH/oYfooCmq7gb1zm4P04h8kB5MIHUu7GU4KwOb9JqaKZVU4+TsEyGs1lPxsRHc+snAu7X9kxBjTbgHfzwr3C/TgeIMhZiTLMsURgTH82tnwgej1bjqRqn2ycitwsqQORENINjWxrFtG3hXgDrdjJY15Mx8RJz/QRHu4eiJYp6Ynf7BFssxyQit8srAMwKeTvA0dZFa+THKlZoMoe1KIzxWAvGJ+BoEmgch3DHEoKqid0tNRYgQiKahzNovTTkvaitCREJuGsyYsVoMpwlCmM8IiIF7tqD4AdvvogURSn2d8xv6m6SCE0s63Cm0kYzGad4X+PqaXf9RKQB6jk43VHh8RYBtxBh1pN7LFaiMxnE9qMwxiPuh30gwqGl4V047hhDGU4X0zqcVkDo+oYAUBlp+qx7rI6js5gaz4lWadZNYNM5dhFfVZQZV8GkU28D4AYsURiTtNzFcHPC6z25rYglQGG8xhBEpBI4xxbbGbCuJ2OSWRlOmY9wU+BolVevBVd8W5IwQZYojElS7oB0pA2Eionv+EEpx86aMhnOEoUxyW1O6Kwkd3yiiDhtaOS2Jgp8LG9ukpAlCmOSmPuBvcudPVWCU9gPoMQdw/DaHNsv24SzwWxjUoA7gF0ez3GDRDzDpCZLFMYYY2KyridjjDExWaIwxhgTkyUKY4wxMVmiMMYYE5MlCmOMMTFZojDGGBPT/wdPFmxDVyOUtAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cluster.rcore(plot=True)\n",
    "print(cluster.rc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible to easily measure the mass function and its slope within a given mass range:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD6CAYAAABEUDf/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAApk0lEQVR4nO3de3xdVZn/8c/TS1p6TZMGUaG0JyCgopA2DIIMDCTFC4hA0jLqULE0wSuDIw38fs7PcbyUqIN3MSkdQUVtElAug0BSFBUcbBMKqIAlp7QCw5jmQm/Qlnb9/lg7bZKeJOe+T06+79frvOBc9t4PPeXJyrPXepY55xARkfw0IewAREQkc5TkRUTymJK8iEgeU5IXEcljSvIiInlMSV5EJI9NCjuAbJg7d66bP39+2GGIiGREe3v7NudcSaz3xkWSnz9/Phs2bAg7DBGRjDCzLcO9p3KNiEgeU5IXEcljSvIiInlszCR5M2sNOwYRkbFmTCR5M6sCKrJ6Uedg59+yekkRkXTL+SRvZoXA0qxf+KGvwPffCS8/n/VLi4ikS84neaAGWJv1q775Itj3Cvz0Mti7i63du6m88SFKr7+XyhsfYmv37qyHJCKSqJxO8mZWBrQBfVm/+JEnQtV/wv/+Ce6o4cpbHqWzayf7naOzayfLb12f9ZBERBKV64uhKpxzXzGz7Nbj+x1fCed/Ge67jve/NpGvOF81OuAg2rUrlJDkkD179tDT08OOHTvYv39/2OGIpMXEiROZOXMmRUVFTJkyJeXz5WySN7MaoDGMa2/t3s3yW9cT7dpFZO5J3PGWD/GxP/2YZ90buGP/WUwwiJRMDyM0CezZs4etW7cyZ84c5s+fz+TJkzGzsMMSSYlzjn379rF9+3a2bt3KvHnzUk70OVmuMbMI0OOc6wvj+stvXX+oNLNtF9VbL+HVN57BDZNXUz5hE6UlM1izrDyM0CTQ09PDnDlzmDt3LgUFBUrwkhfMjIKCAubOncucOXPo6elJ+Zw5meSBKudcS6IHmVmNmW0Y+ujq6kroPNGuXRwItr494GDTtj1M/eBtFBTNo7nw27R+ZAHziqclGp6k0Y4dO5g1a1bYYYhkzKxZs9ixY0fK58m5JB/MiU84wQM45xqdc4uGPkpKYjZnG1akZDoTgoHhwdLMtCL4QBPs3wc/uQz2pP6HL8nbv38/kydPDjsMkYyZPHlyWu415VSSD+bEFznnomHGsWZZOaUlM5hoNrg0M/d4WHILdD0Nt18JB3SzL0wq0Ug+S9ff71y78VoDlJvZwiGvRwDMrCF43uyca8tUEPOKp9H66bNjv1l6Lry7Hu79DLR9DhZ/MVNhiIikLKeSvHPuK7Fe729r4JyrzXJIsZ22ArqegUe+DXNPgLJ/CjsiEZGYcqpcM6a86waInAP3XAPPPRx2NCIiMeXUSH4EReBr9mFNqzzMxElQfQvcXAlrPwQrHoSiBWFHJZLzWlpa6Onpob29nWg0Sm1tLVVVVWk7rqWlhWjU39aLRCIUFhZSUVGR0HkGnqOzs5OFCxdSU1MzbGx9fX3U1dXR0NAw6PWOjg4aGhooLS2lu7ubpUuXUlZWNup/a1o553L2AVQBzUAv4ID24HlhIudZuHChy5htzzq3ap5z3znNuVf6MncdGeTPf/5z2CFIEpqbm11nZ+fB5729vS4Sibj6+vq0HFdfX+8aGhoGHReJRBI6z9DPOOdcWVnZiDFWVVW5qqqqQa+1tra6ioqKw84z9NwjiffvObDBDZdHh3sjnx4ZTfLOuZc2PuD2fW6O+/W/nuUWf22d27JtV0avJ0ryY9XABNyvubnZ+fFmasd1dna6srKyQZ/p7e11ra2tCZ2nrKzMrVy5ctBnVq5cOeiHxUCtra1u5cqVhyX5SCQy6Nr91xr6uZGkI8mrJp8GH1pXwL/uu4KzJzzO0t5GNS8TiaG/NNLR0THo9f7yxdDXEz2utraW2trBczMGlmriPU99fT1Llw7ubt7X10ckEjkstr6+PgBKS0sPez0ajbJo0aLDrtXSktQyoKSNlZp8Tot27WKTO5fj7AWWT/olnd1vBIaZgimSZn19faxatYri4mIAuru7qa+vDzmqw0UiEVauXHlYsuxPlLGSaCLHtbW1UVdXd/DfI5HIoGPiPc/A+n3/+01NTaxbt+6w2JqamqipqaGxcXCbrf56/lBFRUUHz1lYWBjzM+mWcpI3s1n4G6OFwUt9zrnnUj3vWBIpmU5n106+/NoHKLUX+fzkH0D0fRBRog/FL6+Dl54MO4qRHXUyvPuGlE/Tn9jWrVt3MGnU1dVRW1t78CZgY2PjiDcNsynWD5+1a9dSVlY2YtIb7bj+RA3+pmlFRcXBkXt9ff3Bcydy/ba2Njo6Ouju7qa9vf2wHw5tbW0sWbIkZrz9vx309PQMOm9/L5poNJq1G7AJJ3kzm4+/IboYvyVfb/DoCz5SaGZFwGx8L/hmoMk5tz0N8eakNcvKD3at/HrhdZwx+XPQdLmfcVNcOvoJRJLQ19dHdXX1oAQPUFlZSWVlJQ0NDUSj0WFHyP2qq6sHJcnRNDQ0jHrOePX19dHY2Eh7e3tKx/WPnPv6+g7OlCkrK6OyspIVK1bQ3Nyc8PUrKiqoqKigpaWF+vr6QT8s+vr6Rh2N19TU0NLSwsqVKw++1l8SSkfjsbgNV6wf+gBOBTYA64FrgVPjPObaAcedEu/10vnI9I3Xw3RHnbthvnPfKnNud092rz1O6Marn9Ex9Eajc/4GJOA6OztHnbUStoqKisNuTiZzXP9/c3t7+6DP9fb2OsD19vamdP36+vpBf9ZDb+A2NDTEvKFaU1Pj2tvbD94Abm1tPfjdxCNrN17NrAlYAlQ758qdc191zj0Wxw+Qx4LPLgqOv8zMsr+VX7YVLYDLboPeLdD8Yd/UTCTNWlpaDrtBCIfqvi0tLWkbcSdi4cKFlJaWDnrEUldXR11d3WE18NHEOq7/v3PoyLr/+YYNG1K6fkVFBR0dHQcf8cbc0NBAYWEhGzZsYNGiRQe/m2x+LyOWa8xsNlAP1DnnNqdyoeD468wsYmY3BefM2xIOx54BF34T7vwY3Hc9vPdrYUckeaS/PDFSXbezs3NQqWA41dXVCV27vr5+xCQVT+mlsbGRysrKhBP8SMcNrc0PNDTe4c4TjUZZuHDhYTX4/h8W/X/u69cPnkHX1tZ2cEFUeXn5oMVVA28AZ7MWf9BwQ3z/GwArRno/lQdwaabOPfSR9XLNQPd/1rnPzXLu0cbwYshD471c01+GiPVrf/97Q0sXuaK5ufmwEkl7e/uoJYzRjlu5cuVhZZT29nZXWFgY93na29tdRUXFYeWd0cos9fX1MRdDDf18ouWpjJdrnHOrM/jD5fZMnTunVPwbvOnd8Ms6ePbwKVgiyeif/93WNrgZazQapampiUgkMmgmR65oa2ujp6eHSCRCNBolGo0eXPrfP9rtHxEnetz1119/2OyZVatWDXpttPP036wdqr6+npqammF/g+nu7j7stbq6ukFz8vundSb620uqzP8QyG+LFi1ysWpyWbNnB6w5H15+Hq5sg5I3hRdLnnjqqac46aSTwg4jVP3JcOHChRQVFdHT08OiRYsoKyujra2N5uZmSktL4yrZZENfXx9z5syJ+V4kEqGzsxPwM1DOO++8gyWTeI8D/wOtoaGB4uJiuru7B5VOEjlPY2MjnZ2dFBcX09nZOeyfY//1+vvhLFmyhNraWsrKyujo6Dh4L6B/Jk6iU1nj/XtuZu3O3/s8/L3RkryZrXXOHX53J06pHp8OoSd54IXNzzDth5VsPzCVa2b+B9+44jxtIZgCJXkZD9KR5OOZXVNmZsea2axkHgQbfox3H/75S6zYcw1H0c21279M7S2PhB2SiIwD8SyGKgVyp6g3RkW7drHfvYm6fTV8o+B7LOv7LrjzINjia2v37oMLqiIl01mzrFwjfRFJWTxJ/nZ8m99WINFlWsX4xVDjXn/rg18ceCfHv/YCH590Jzz6fTj9owAsv3U9nV07OeCgs2sny29dP/wWhCIicRo1yTvnqoP58hXAHPxUnY3xXsDMsnsrOUcNbH1w55wruKJkL9Pu/z9QfBwcX0m0axcHgtsjB5wf+YuIpCqu3jXOuZfxI3rMbIGZrcCP7tvc6M3I6kZ5f1w4bHPwvafBf54PzVfAla0HR/oHHEwwP/IXEUlVwv3knXObnXOrnXM3A3PMbIWZXRI0Lov5+VSDzEsF0+EffwYF0+AnS/nBkgWUlsxgohmlJTNYs6w87AhFJA+k1GrY+f41jwGY2XlmVgl040f4+duyIF1mHw2X/RRueQ9HP1BD66fuhElTwo5KRPJI2naGcs6tC0b4dwCLzOxKM7s4XefPW0cvhIu+C1t/D/d8GsbB4rR0GQ8L+WT8Stff77TvDGVm5wLVwFJ8b/nmsBdD5byTq2DbJnjoBig5Ac78VNgR5byJEyeyb98+CgoKwg5FJCP27dvHxIkTUz5PWkbyZnaKma0ys278RiGVwCpgjhJ8nM6ug7dcDK3/D56+N+xoct7MmTPZvl0VQclf27dvZ+bMmSmfJ+kkb2bzzewzZrYJaMeP3lcDpc6545zvI/9yyhGOFxMmwEXfgzecArdfCS/9MeyIclpRURG9vb1s27aNvXv3qnQjecE5x969e9m2bRu9vb0H+8+nIqEGZUGbgiVALVAGvAw0AQ0ujk1EwpILvWvitv1/YPU/wIRJfvvAGUeGHVHO2rNnDz09PezYsYP9+/eHHY5IWkycOJGZM2dSVFTElCnxTcRIqUFZcIJL8Im9AjCgBZ/YR+2da2aXBDdjQzOmkjzAi4/Bf74bCo+BDzT5naZERIaRUoOyoM7evwvuEufcBOfckngSfODw7dFlZG84FT7UAjv/BjdXwF//EHZEIjJGxdNq+ADQCPQ3Wx54gMV4rd9coApY4JxL/RZxCsbcSL7ftk1wWzVsfxEu/j689ZKwIxKRHDTSSD6eKZQdzrmrkrzwDcAYzK45Yu7xcOU6WPtBaLkCeqJw1r8c7FwpIjKaeGbXrE325M65PqAh2eMFmF4Ml98JJ1fDg1+AOz8Br+0NOyoRGSPi6UL51VQukOrxgm91cMlqKCr1C6b6tsDSH8ERsbcyExHpl7a2BsMJVsBKqszgH66Hixtg63/DmsXQo95vIjKyjCd5VK5Jr7df5ss3u7rg5vNg66NhRyQiOWzUco2ZrU/y3IVAUfBPSaf5Z/obsrdVwa0Xwvu/5/vfiIgMEc/smjn4xU+dMd6rxu//Guu9cnyCb002OBlBcalP9D/7INy+HHo3w1mf0cwbERkkniQfdc5dN/RFMzsPqBupnUHwme4U4pORTCuCy38Bd30SHvwidEfhwm/CJHVmFBFv1Jq8c27xMG+dOlq/mmBVbMwJ+pImk6b4m7HnXA+P/wR+fAnsTnS/dRHJV6nceI23w6QyTqaZwTnX+WmWf30U1lT6hVMiMu6lkuRnx/m5SArXkES8bYmfebO7B1af56daisi4lkqS32xmq0b6gJl9JoXzSzKOPQOubPMLpW69EJ5sCTsiEQlR0tv/OeduN7PFwaYh/bNvevAj91J8W+Koc+78tEQq8Ssu9Yl+7Yf8zJueKPz9tZp5IzIOpbTHq3Ou1syqgBsYXJaJAvXOudWpnF9SMK0I/unncNen4Fdf8on+wm/6G7UiMm6kvJG3c64FP5LHzBY457TWPldMmuJbFBeX+kTftxWW/tj/ABCRcSGtbQ2U4MOxtXs3lTc+ROn191J540Ns7d596E0zOHslXLoGnl/vNyHpjrV2TUTy0YhJ3sweCDbrnp+leCQJy29dT2fXTvY7R2fXTpbfGqMTxclVcPld8Eqv73mz5ZHsByoiWTfaSL4aPx++0czWK+HnpmjXLg4Ee3MdcP55TMe+w9+QnVYMP7wIHk96qwARGSNGTPLOuZedc6uDVa8V+ITfMiDhz8pKlDKiSMl0JgQTZyaYfz6s4lJY3gpHnwY/r4FfrYI4NnMXkbEp7pr8gIS/CJ/wDXgwSPhXKuGHZ82yckpLZjDRjNKSGaxZVj7yAf0zb97+j34Tkp/Xwmt7shOsiGTVqBt5j3oCswX4DbuX4puRNQNNzrntqYeXHmN2I+9Mcw5+8zX41Rdh3hlw2W2aeSMyBo20kXfKs2ucc5udc18NLnAVvjVxh5ndb2aXpHp+ySAzOPtaP/PmhXZ/Q3bbs2FHJSJplPYplEHCPw6f8E8zs2fNbK0Sfg47uQqW3QWvvgxrKuC5h8OOSETSJGPb/wUJ/7og4d+AEn5um3d6MPNmbjDz5mdhRyQiaZCNPV5xzj02JOEvNrP7s3FtSUBRBK5s9Qn/57Vw58f96F5ExqysJPmBgoR/lRqX5agj5sCH7oB3fho2/gS+9w7Y1BZ2VCKSpKwneRkDJhVAxedgeRsUzIDbLtWoXmSMGjXJm9kpWYhDctHRC9lafR8/LbiU/R238bf6Mv7W8V9hRyUiCYhnJK92wePY8tue5P/uuJRL9n6el/dP4ci7PgB3fkKjepExIp5WwwvNbC3QGuc5DRi4wqrHOXdHwpFJTujvi/M4x3HB3i9xzaQ7uGrjbdD5ILzvW3BcRdghisgI4knyfcBaYF2c5ywCGvCtD3qBFUlFJjkhUjKdzq6dHHCwzwq4vehKrlpyNfzio/DjS6Hsclj8JZiqrhYiuSieck2bc+6OoHfNiA98Qn8Wn+C/4pwr1ih+bIvZF+foRVD7Wzjzanjsx34GzrPxjgFEJJtS7l0DYGbn4kfvpfhdolYEST8nqHdNBj2/wY/qt/0FypbB4i9qVC+SZRnrXWNm84NFTa34Wnylc25JLiV4ybBBo/of+VF954NhRyUigaSTvJmtAjqBcuAq59xxzjn9zj4eTZ4Klf8OH3kACqbBjy6Gu6+GV3OmEanIuJVwkjezS8ysG6jDT69c4JzTNEuBY8qh9jdwxqeg44dw0xka1YuELO4kb2anmNl6fM19M7AwaE+g0owcMvkIWPwF+Mj9MGmqRvUiIYtnxessM7sJ6MDfWK12zi1yzj0WzwXUcXKcOuY0uOq3cMYnof3WYFT/q7CjEhl34hnJPwfUAPXOuSLn3O0JXuP6hKOS/DD5CD/bZvkDMGkK/Oj9cPc/w54dYUcmMm7EsxiqEF+i6TGzzwSvDV3VGstc/Hz5U5OOTvLDMafBVb+DB78Iv/+un1N/0bchck7YkYnkvXiSfJtzbkkyJzezQiCazLGSZyYfAed/CU56H9z5Mb8xyaKP+Fk5U2aGHZ1I3oqnXNOR7Mmdc31AY7LHSx6a93d+VP+OT8CGH8D3zoDoQ2FHJZK3Rk3yzrnrUrlAqsdLHuof1X/kPpg4GX74PrjnGtXqRTIg45uGBC0PRA437/TBo/qbNKoXSbds7AzVkIVryFhVMO3QqH5C/6j+07BnZ9iRieSFUW+8BgugklGIbztcmOTxMp70j+of/CL89/fgmV/CuZ+Ft18GEyaGHZ3ImBXP7Jo5+CmUnTHeq8bPnon1Xjk+wce72YiMdwXT4F1fhjdfBPdd52fh/P67fgbOceeBWdgRiow5o7YaNrMHnHOLY7x+Hn7Xp2FXvgaf6XbObUw10FSo1fAY5Bz86eew7vPQ+xwsONsn+zecEnZkIjknpVbDsRJ84NTRWhsEXSljXlhkRGbw1kvg4+vhXfXw0pPQeDbcvgL6toYdnciYkcqN13gbk/WkcA0Z7yYVwOlXwdUb4Z3XwFN3wbcXwgOfhVd6w45OJOelkuRnx/m5SArXEPGmzoaKf4NPtsPJ1fDId+Cbp8Aj34Z9r4YdnUjOSiXJbw42DhnWgF43Iinb2r2byjWdlD56ATXTv8ErR57qR/TfKYcnmuDAgbBDFMk5Ke3xamYNwLkcmn3Tgx+5l+Kbk0Wdc+enIc6U6MZrfqi88SE6u3ZywMEEg9KSGbRedABa/x+89AS8/u1Q+QWInB12qCJZNdKN13imUA7LOVdrZlXADQwuy0TxrYm1Y5SkTbRrFweCMckB559T+h5Y8BD8sQXW/btfTHVcJVR+Hl73lnADFskBKa94dc61BPu7TgBKnXMTgudK8JJWkZLpTAimyk8w/9w/mQBvWwKf2OD71z//B7jpTPjFx+HlF8ILWCQHjJjkzWxWIidzzm2O97MjndvMCs1spZk1mFm9mTUH/yxMJB7JL2uWlVNaMoOJZpSWzGDNsvLBH5g81e9E9amN8I6Pw5NN8O0yaPs8vKpdKmV8GrEmb2aXAp3pXsxkZqcAEefcHTHeKwSud87VDXm9GSjD7y3bl8j1VJMfp3q3+DYJTzbBtGI4uw4WXuGnZYrkkaQXQwVb/VWONosmwWBuACpiJfhADVARY9S+Cl/3T2oDExmH5hwLl66Gml/7+vwvV8J3T/MraVOYcCAylsSz4vWrQJOZPWBmq4JReELM7BQzu8HM7gd+5pz72ggfj6K59ZJObzgVLr8LPng7TJ4GzR+Gm8+D5x4OOzKRjEtoCmVQvlmK37e1DT9tMgr0cWhla3/nyf6plJXB5xpGGL3Hc+0afNviUudcQlsKqlwjBx3YD4//zJdxdrwIJ7zHL7IqOSHsyESSNlK5Jul58kHzsUp8Mi/k0Oi7P+lHgfX4PWJTvutlZp3BuWoTPVZJXg6zdzc8ehP89uuwbxeUXQ7nXA8zjwo7MpGEZSTJZ4uZRYB6/MKqutE+H4uSvAxr1zb4zVdh/c0wscDPzjnjk9pcXMaUMZnkgxuvNUAx0I1P8i2jHFMTHDPIvHnzFm7ZsiUTYUq+6O6EB7/gb8pOPxLO/BSULYOpCc0iFglFVpK8mZ2Syb7xQQuFiHOuMtFjNZKXuD2/wfew3/wbmDIbypfD6R+FGUeGHZnIsFLqJ5+AWjPL5O+4dfiplfUZvIaMd0cvgmV3w4oHofQc+N3X4etvhbuv9qN9kTFm1CRvZtea2bNmtmnA434zO3foR4G+IZ/bFONzSQkWQHUQoxwjknZvXAhLfuhbG5/yj7Dxp76PfdPl8EJH2NGJxC3eefJL8LXxx4DFzrnznXMPxvj4ZvxmIi8Hn68b5nPDMrPOEUbrPUCh2htI1hSXwoXfhH9+0m9a0vlrWP0PcMsF8GybFlVJzou3XHMuUO2cWzJKf5pK59yi4FEELE6k/02QvEdaCBUB+hJtayCSspmvg4rPwTV/9O2Mu5+FH18KDWfBky2w/7WwIxSJKZ5yzbnAy8F+rSOK8QOgjgTKK0Hybow1VXLAD4CmeM8nknZTZ/mZN1c/ARd9F17bA7cvh2+fCo82+vn3IjkknpF8bZxtg1uHvhAsgipOMKb6YCbNUKuBjmQWQ4mk3aQCOPVD8LFH4bKfwszXwy+vha+/BX59A+zW1saSG+LZNMTiOVHQzCyWwrij8eeJmlndkLp8BFjvnKtO5FwiGTdhApz4Hv/Y8nt4+Jvw61X+n2WX+5bHhfPCjlLGsXiSfLwbdg8n4WZjQdkmqdWtIqE59h3+8ben4OFv+VW0f1gNb70UzrwajnprXKfZ2r2b5beuJ9q1i0jJdNYsK2de8bQMBy/5Kp5yTaLllnQfLzK2HHkSXHwTXP24X0j1zL3w/TPhx1Xw3O9GnZGz/Nb1dHbtZL9zdHbtZPmt67MUuOSjeJJ8m5ldnMzJzewSfLdKkfFn9tFw/pf8jJxzPwsvPga3vNe3Of7zXb4jZgwx97IVSVI8Sb4RuDnR1axmNht/s/T7yQQmkjeOmAN/f61P9u+90d+Ubfonv4FJ+y1+hs4Aw+5lK5KEeBZDRfFdIJ8zs7fHc1Izmw9sAFY7555LJUCRvDH5CN8L55PtUPUDKJju2yV842TfPiHYh3bUvWxFEhB3g7JgtstngGZgLbDOObd9wPuzgEVALVAFtDjnlqY94iSoQZnkJOdg80Pwu29A9FdQMBMWXQGnfwxmvT7s6GQMSVsXSjOrwpdvBs646WPwNMmX8e0M4plbnxVK8pLzXtwIj3zLtzq2CXDie/2m4wvO9tM0RUaQ9lbDQbJfCizAT5Hswe8E1Qw0pWMnqHRSkpcxo2ezn3q58SfwSg/MWQALPwynfBBmlIQdneSoMblpSDopycuYs+9VeOpuaP8BbHkYJkyGky705Zz5Z4HFtUZRxomRknw8i6FEJNsmT4W3VftH1zN+Fs7Gn8Cf7oDi4/zo/u0fgOlahiIjS7jYZ2bzzWxV0Ct+/4DHJjP7npkdm4lARcatkhPgXavgX56Gixtg2lx44LNw44lw+5Xw3MNqeSzDSvTG6034rpKGv+EaDd4q5FD7Agc0OOc+lrYoU6RyjeSd//2zH90//jPY8zLMPSEY3V8G04rCjk6yLC01eTPbgE/m9Qxzc9XMFgDVwHXANmDRwGmWYVGSl7y1d7cv4Wz4AbywASZNhTe/39fuj/k71e7HiZSTvJmtBTYEu0TFe9EGYL5z7vy4I80QJXkZF1560if7J5pg7w448s1+dP+2pXBEYdjRSQaltJG3mZ0HbE4kwQMEfd83p2uPVxEZxVEnwwU3+tr9hd+CSVPglyvhP06EX3wM/rpetftxaNSRvJnd5Jz7aNIXSPH4dNBIXsatFzf6aZhPtsDenfC6twaj+yUwNdUu4pIrUhrJ41ewpiKnFkaJjCtvOMVvRP4vT8MFX/erae/9jB/d3/kJeKFdo/s8F888+W0pXiPV40UkVVNmwqKP+FYJL3b42v0fb4fHfgRHvc3fqD252n9O8ko2Ng0pTfF4EUkXM3jjQrjoO350/56vgTsA91zjR/d3Xw3Pa3SfT+La49XMZiUzFTLoKa85XCK5aOpsOG0FlF8Jz2/wtfvH1/r590UROHmJH93PPS7sSCUF8dx4PRW4IZmpkMHUy1XOuY3JhZceuvEqEqdXX/a7Vj3ZBJt/Czh4w6k+4b/1Upj5urAjlBjSMU++ASiMtz980Fu+GYiGPbMGlORFkrL9RV+3f6IJXnrC37RdcLafmXPiBTB1VtgRSiBdK16bgXOBBvyK141D3p8FVOBbEFcBbbmwEAqU5EVS1vWMT/ZPNkPfFr+y9oR3+xH+cRUwqSDsCMe1dG4aUgPcwPCbhvTX3+sSXTyVSUryImniHPz1D76c88c7fM/7I+b4VgpvWwLHnK5NTkKQiU1DaoBK4FSgCG0aIjL+7N8HnQ/6Ef4z98K+3TD7GDi5yo/wX/fmsCMcN7RpiJK8SGbt2QlP/5cf4Xf+Ctx+v7r25Gqf9GcfHXaEeU1JXkleJHt2dvnOmE82w/PrAYNjz/QboLz5Il/ekbTKSpIP5sTXBE/LgW58vb7TOXdzWi6SJCV5kZB0d/q+OU82QfezMLEAjl/sR/hvepffAUtSFspIPqjbfx844JwLdZtBJXmRkDkH/7MRnmiGP7bAzv+FKbPgpPf5Ef78s2DCxLCjHLNCK9eYWQVwv3Mu1G9PSV4khxzYD5t/48s5f77L976fcZRfbPW2anj9KdrsJEGh1uTN7FnnXKjropXkRXLUvlfgL/f5Ef6mB+DAPiic50f4J14Ax5ymEX4cRkry2SijdGThGiIyFk0+At5ysX/s7oGn74Gn7oY/NMLvvwPTS+DE98KJF8KCv9eiqyRkI8n3ZOEaIjLWTSuCssv949XtfmT/9D3+xm37LTBlNrxpMZx0oV9lWzA97IjHhFGTvJmtjbdnzTDyf46miKTX1FnBoqoq2PcqRH8NT98NT9/ra/mTpkLpeXDSBX6WzrSisCPOWfGM5CMpXkN/+iKSvMlT4YR3+ccFr8HW3/uSztP3wDP/BTYR5r/Tj/BPvABmvT7siHNKPK2Ge4D1KVyjQrNrRCTtnPO7XD0V1PG7N/nXjy73yf6kC6F4fOxZlNLsGjM7gF/UlExtvQiYrSQvIhnX9Qw8dZdP+v+z0b925JsPjfCPOjlvp2ammuQfcM4tTuHiG4a7eLYoyYuMM31bfS+dp+6BrY/4LQ4Lj/UJ/6QL4ejT8qpbZqpTKFOdAtmW4vEiIokpnAenf9Q/dm3zXTIHTc080k/NPOkCmJ/fUzPVoExExo+BUzM3tcLencHUzPN9wh+jUzPDXgwlIpIbRpya2QQTp8D8M30TteMX58WNW43kRUT2B1Mz/3KfH+Fve8a/PmfBoYQ//0y/QjcHqZ+8kryIJKL3OZ/sN7X6ZmqvvQKTjoAFZwVJvxLmzA87yoNUrhERScSc+XDaCv/Y9yps+V2Q9B/wD4C5b4LjKn3CP/YMmDQl1JCHo5G8iEgiujsPJfznfgf798Dk6RA5B46v8Im/8JishqSRvIhIuhSX+sfpV8HeXbD5t/BsK/zlAd9mAfwirOMrfcKfdzpMnBxauBrJi4ikg3Ow7S9BSacVtjzi++NPmRWM8hf7KZoZ6K2jkbyISKaZQckJ/nHGJ2HPDog+5JP+s22+5QL49gr9M3beuAgmZjYNayQvIpJpzsHf/nxolL/1v8Hth6mzfcvk/lH+jJKkTq+RvIhImMzgdW/xj3deA6/0+YVYm1p9Pf9PdwDw+X3L+F3xpaxZVs684mlpuXT+dOgRERkrjiiEt7wf3v9d+PTTXDX963zttSU8euAEOrt2svzWVLq7D6aRvIhImCZMoLXnKO5z7z/4UrRrV/pOn7YziYhIUiIl05kQtLqfYP55uijJi4iEbM2yckpLZjDRjNKSGaxZVp62c6tcIyISsnnF02j99NkZObdG8iIieUxJXkQkjynJi4jkMSV5EZE8Ni7aGphZF7Al7DjGgLnAtrCDkKzT9z72Heuci9kTYVwkeYmPmW0Yrv+F5C997/lN5RoRkTymJC8ikseU5EVE8piSvIhIHlOSl4Eaww5AQqHvPY9pdo2ISB7TSF5EJI8pyYuI5DEleRGRPKYkLyKSx7RpSJ4zsxqgEOgL/hl1zrUkcZ6VQDFQFrzUkMx5JPPS9Z3HOG+rc64y1fNIdml2TR4zs3qg2zn3lZFei/M8q5xzfcHzMmAd0Oacq05v1JKKdH3nMc5bBTQ75ywNYUoWKcnnKTOLAO3OuTkx3usFFvQn7VHOU4NP5tEhr1cArUC1RvS5IV3feYxjC4HVQJWS/Nijmnz+qgXahnkvCtTEeZ7SoQkewDnXhi8H1CYVnWRCur7zoWqAtUkeKyFTks9fFfj/sWOJAvHWVleaWfMw720A1KI2d6TrOz8oKM31/0CXMUhJPn+VAd3DvBcl/uQ83MgQoAjoSSQoyah0fecDVTjnOpIPScKm2TXjV2E8HxplNkUZ6nsylhQm8uHgfoy+3zFOI/k8FNwogwz+ih0kAID6TF1D4pfu7zy4iduTzI1ayS1K8pKwIKHU42fWDFcDlrGtSrOm8oOSfB4aMPoqzNAlmoE6JYHckc7vPJgTr+82TyjJj0+FJPlrfbCwptk5p1rt2FJIHN958FtakX5Dyx+68Zq/2oDSYd4rwk9/TEjQ2qBTCT5npeM7rwHKzWzhkNcjAGbWEDxvDtZKSI5Tks9fHRzqMzNUBGgY5r2Ygl/h+4YmeDOr0P/sOSPl73y41gfB91/hnNPitzFG5Zr81UCMedHBr+NlQFOM9yKxThS0MGCYEbwaVuWOtH3nkj+U5PNUUFNtDGroA9Xjb5r2DXzRzNqBzmCF48DXywhaF5hZ1YBHTVC+KczUf4MkJl3f+TCKgmMK0xCqZJEalOW5IW1nS4H1sWbFBImhCqgceNMtaGxVOMIl6lLpbijpl+p3PuQzVcBSfMuEQnxJKAqs0Bz6sUFJXkQkj6lcIyKSx5TkRUTymJK8iEgeU5IXEcljSvIiInlMSV5EJI8pyYuI5DEleRGRPKYkLxInM1tpZu1m5oKVwPEe1xwc0znCpugiGaEkLxIn59xXnHML8fueFgZL/kc0pNdLpXOuOlPxicSiJC+SuE587/Z42u5WAK3Bv/dkLCKRYSjJiySnAaiIoytjEUruEiIleZEkDOjquGS4zwS92hPegUsknZTkRZLXwsglmwrnXEe2ghGJRUleJHkNQJl2V5JcpiQvkqRgb9s+/MYbgwS7LWnvWwmdkrxIapqIXbJZNNxuSyLZpCQvkpoGIBJjn1TNqJGcoCQvkoLgxmoUvw8qAGZWgd8LVSR0SvIiqWsBagY8j6hUI7lCSV4kdQ34NgcVwXOVaiRnTAo7AJGxzjkXNbMoUG1moFk1kkM0khdJXHGM1xrwq18jzrm+7IYjMjwleZEEBL1qYnWfbAQKiV2qKcpgSCIjUpIXiZOZrQTa8VMmOwe2Gg5G7y0DetpgZpGgf3xd8NI69ZOXbDPnXNgxiIhIhmgkLyKSx5TkRUTymJK8iEgeU5IXEcljSvIiInlMSV5EJI8pyYuI5DEleRGRPKYkLyKSx5TkRUTy2P8HKR1AIZWli0EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "m_mean, m_hist, dm, alpha, ealpha, yalpha, eyalpha=cts.mass_function(cluster,mmin=0.1,mmax=0.8,plot=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calling the mass function function internally simply sets the value of alpha and saves the mass range over which alpha was measured with cluster.mmin and cluster.mmax "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.2623485836819826"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cluster.mass_function(mmin=0.1,mmax=0.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at how velocity dispersion changes as a function of stellar mass is also a proxy for a cluster's dynamical state, with the power-law slope eta evolving towards (but never reaching) -0.5. An eta of -0.5 corresponds to a state of complete energy equipartition. The external call to ``eta_function`` (shown below) and the internal call behave similarly to ``mass_function`` above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAD6CAYAAAC/KwBlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjAElEQVR4nO3da3Bc533f8e+fFEWLNy0BQlcbphZ2LCmxXYNA7KRNVFtAaufasUgxTtIwDaVlMpm2M05KyJ30RTPTyqBTv+nNC4ZJ2CatRMDuNKmTGQNM4s6kzQRYKLEnkmNrVxbliyIYwPImi6SAf1+cs8BisXuw2D3AWSx+n5kd7J7znOc8WCzOf5/rMXdHRESklh1JF0BERFqbAoWIiERSoBARkUgKFCIiEkmBQkREIilQiIhIpNuSLoDE79ChQ3748OGkiyEiW0gul/uOu3dV26dA0YYOHz7M1NRU0sUQkS3EzF6utU9NTyIiEkmBQkREIilQiIhIJAUKERGJpM5sWeXS7OucPD9JYeY66a69nDvRT3fnnqSLJSIJUY1CVjl5fpL8zDUW3MnPXOPk+cmkiyQidbo0+zqDn/4iPZ/4IwY//UUuzb7edJ4KFLJKYeY6i+Hq84sevBaRrWEjvugpUMgq6a697LDg+Q4LXovI1rARX/QUKGSVcyf66enax04zerr2ce5Ef9JFEpE6bcQXPdMd7tpPX1+fa2a2yPbU6GAUM8u5e1+1fRr1JCLSRro79zD+8UdizVNNTyIiEkmBQkREIilQiIhIJAUKERGJpEAhIiKRNOqpCjPLACmgGP4suPtYDPmOu/tgjX2ngU6gN9yUjeOcIiLNUqCoYGbDwKy7nynfZmbp8m0N5HsUGIg459PuXgxf9wIXzey4ux9r9JwiInFQoChjZmkg4+4Hy7e7+5CZzZvZSOlivs58U8DxGvsyBLWHpXzdfdrMjgHjZnZUNQsRSZL6KFY6BUzU2FcAMg3mmwGerbGvx90LlRvdfYKg6etUg+cUEYmFAsVKAwQBoZoCULV/IUrYjFS66Fdz2sxGa+ybAqpOqRcR2SwKFCv1ArM19hVo7KI94O7TEftr1WAAOoC5Bs4pIhIb9VGsT2o9icP+h5GoNLVGQYV61zpeRGSjKVCEwg5nqN1EtN780sBcI53f4fGl/pDhNdKs6jfp7u5u5JQiIlUpUGyco40Opw2D1jBwrFpHd4m7j1ClxtHX16e140UkNuqjCJV98081m1c4Z6KZIa2jwJCGxYpIK2iLQGFmOTPzBh5H13GaFHU0S4W1gY6omsAaxw8Do2FtQUQkcW3R9OTuR2LKagLoqbGvg2C46loyQL+ZVZYpDWBm2fD1aDhXYkm4jEdeQUJE6uYOb1yGa6/BtVfDn38Ht+2G/idiOUVbBIoYTbO81lKlNJCtsW9JrX6J0hIe7l51Al24v1gZJMxsoDKgiMg28OZNuB5e9K+9BlfLgsCKx2vw5hurj+96UIFig2SBXOXGsDmpF7hQZV+60WamsjwGYKlzutIg0XMtRGSrcIc3inC17CJ/7dWy53+3vO+7NaZQ3dEB++6G/XdD9w8Ez0uP/aXnd8FbUrEVW4GijLsXzGzEzIbdfahs1zBB53KxPL2Z5YBeMzuyxqQ6CJquMLNUeT7hzO1TwLMVfSYdBP0iqQZ/HRHZTG/egKvfhivfhqvfgivfWv382quwcHP1sTt3L1/oO3vg7T+4fMHff0/wc9/dsPcuuO32Tf/VFCgqhAsAZsL+giJBn8V4jRFIE6zRyR1e/I+zvHLsRTMrAE+GAeNimEetjvWhGttjc2n2dU6en6Qwc510117Oneinu3PPRp+2ZcshsoI73LgSXOivfHM5GCw9DwPB699Zfextd8CB+4JH9wfCi/7dKy/+++6Gt9wJZpv/u9XJ3DXkvt309fX51FQ9/e6BwU9/kfzMNRYddhj0dO1j/OOPbGAJW7scso0sLsD1meBCX37RvxoGgivhtlvXVx+7pxP23wcH7g0CQen5/jAwHLg3aP5p4QBQzsxy7l51mSLVKITCzHUWw+8Lix683s7lkDby5g24/A24/AoUX6n4+XIQBBbfXHnMjttg/73B4+6H4Z2DwfNSzaC0b9dbkvmdEqBAIaS79q74Jp/u2rutyyFbyM3rKy/8lcHg6qtAWauJ7Qgu8ne+Dd72/uDngfvgwP3LtYG9XbCjLaaYxUaBQjh3on9V38B2Loe0CHf47vzq2kDx5eXnlSODduyCO98KqbdBz6PBzzvfBqnu4PmB+2HnrmR+ny1MfRRtaL19FCKJWbgFxUswVwges3mY/3oYCC7BzWsr0+/aE174ywNA9/K2ffeoNtAg9VGISHLevBnUAsqDQel58RL4wnLa2/fBwQeCxwM/XBEU3g57OrZM53A7UaAQkea9eQPmX4a5/OqAcPkV8MXltLsPQEca7nsffN9jwfPOnuDn3i4FghakQCEi9bn1RtAsNFeoCAhhMCjvNH7LndDRA2/th/ccXxkM9nQqGGwxChQistIbV2DmK/DaC8Fj5oWgdnD5G6wIBnccDC783R+Ajp9ZGQzuOKhg0EYUKES2q5vXYeZvl4PBay/Aa1+BK99YTrNrD3S9K1hSoiMMAp3poA9hT0dyZZdNpUAh0u5uvQHf+WpYS3g+CAavPR90JJdqCDt3Q9f3wOG/H6w6etdDwePObo0iEgUKkbbx5s2g76A8GMx8JehHKHUm77gNOt8J9x+B9/1cEAy6HoKOB2DHzmTLLy1LgUJip8X9NtjiYnDxf+35lc1Gsy8uL0dhO4KmorseCkYWdT0Idz0c9CFowpmskwKFxO7k+cmlpTjyM9c4eX5Si/s14+qr8M1c2eM5uHE53Glw8O1BEHjXjwY/73owqDVso7WIZGMpUEjstLhfE964At/+q7KgMB2sYgpBs9Hd3wvvfgzu6w2ed70LbteaWLKxFCgkdlrcr05v3oTX/mY5IHwzF4xCKnUwd6SD0Ub3Hwke97wbdt2RaJFle1KgkNhpcb8q3IN+hfImpG9/CRZuBPv3HIK39gX9Cff3BjUGDT+VFqFAIbHr7tyjPolrr1X0K0wH90qGYG7Cfe+D92eWawt3vk0T1KRlKVCINGtxEV79a3jp/ywHhcuvBPtsZ9DB/PBPLQeFrgdhp/71ZOvQp1WkEddeg/yfwIsXg5+l+yWn3h6sb/T+XwqCwr3vUWezbHkKFCL1WLgFr/wlvDgRPF79UrB9zyF4x6PBTXJ6Pgj77kq2nCIbQIFCpJb5ry/XGApfhJtXgyGqb3s/fOhfwzsG4J73aIkLaXsKFCIlN1+Hl/98udYw+2Kw/c5uePfRoObwwA8HS2iLbCMKFLJ9uQdrIb04EdQcXv6/wXDV294Ch/8B9D8R1Bo636ERSbKtKVDI9vLdeSj82XKTUmnWc9eD8P1PQs+HgklumtgmskSBQtrb4gJ866+CWkP+InxjMlhJdfed0PMPoWcoaFK6861Jl1SkZW1IoDCzDwEpoDS1tODuf7IR5xKp6lvPweRvwVc+H9QisGCS2w/9WhAY7u/TXAaROsX6n2JmnwEy4csiMEcQLFJm5sCIu/9ynOcUWXLru/A3/zMIEN/MBTOgH/oJeOePQPqDsLcz6RKKbEmxBQoz+zVg1N1/KSLNY2b2tLt/Iq7zijD3Ekz9Njz334Law6HvgY+cgff+tEYoicQgzhrFZXe/GJXA3T9rZlrpTJq3uABfGw9qDy9OBDfqefDHgpFKD/ywRimJxCjOQOExpxNZ7fp3gprD1G8H93zedzc8chqO/AIcuC/p0om0pTgDxUEzewK44O5XKnea2WFggKCTW6R+7vCNKZg8G/RBLNyEwz8Eg78BD/64bu0pssFiCxTu/ikz+5fA182sWsNwEXja3X8zrnNKm7t5Hb48FjQvvfoluH0/9J6A/pPBvaBFZFPEOurJ3T8FfMrMHgDSZbsK7v5SnOeSNvadF2HqHDz3+8G9oe96GH7s0/Cex2H3/qRLJ7LtbMhA8jAoKDBI/RbehK/+cVB7KPwZ7NgFD/9k0Dnd/QPqnBZJ0KbPODKzJ9z9tzb7vNKirv4dTP9XyP1OsJzGgfvhQ78O7/t52H930qUTEZJZwqMngXNKK3GHS/8vqD08/weweCuYEPeRM/A9H9aMaZEWE+eEuy8AB+tI2gtowt125A7P/R78xX+G154P1lv6/ieh7xfh0DuTLp2I1BDnV7dRgqGvhYg0nSwv8SHbya3vwh/8M/jyaHCzn5/8D/B9j+k2oSJbQJyBYgJ41N0/G5XIzNJR+6UNXfkWPPOz8K3p4M5wP/Sr6pwW2ULinEfxkpnl6kiajeucsgV8IwfP/AzcuArHfx8e+vGkSyQi6xT3PIrn6kijYbPbxZdG4X/9SrDMxskvwD3fl3SJRKQBDd0V3sweMLNn15H+o+E9KmQ7WFyEiX8Dn3sC7j8CmT9VkBDZwhqtUfQCn1wrUbiUR4Gg/2LezIbc/R81eE7ZCm5chc+dgr/9PPT+PPzov4fbbk+6VCLShEYDxTRwJyzdzY4ad7D7JGX3qDCzR83so+7+uQbPK61s/uvwPz4GM1+BDw/D+0+p01qkDTTU9BT2Mwya2YsEtYUJM5s1s/dWJB0APlN23EVWrgEl7eLrfw5nPxTMrv65z8IHfklBQqRNNNOZ3Q8MEdQu5sLXZ4DypqU0q+dVXG7inBKjS7Ovc/L8JIWZ66S79nLuRD/dnXvWn1Hud+HzvwoHH4CPPQOH3hF7WUUkOQ13ZgPPuvtn3f0ld7/s7hPAU5Wd1lXuTaEbF7WIk+cnyc9cY8Gd/Mw1Tp6fXF8GC2/CH52GP/wX8MAj8MSEgoRIG2ooUBDMwM5X2Z6nvmU8pAUUZq6zGIbtRQ9e1+278/D7j8FfZuEDvwI/cwHuSG1IOUUkWY32UTwHfMLMKm8OMAyMw1Kto1ojdaqRc0r80l172RH+hXZY8LouM1+Fs48G/RI/+R/hw/9OC/mJtLFGaxQQrO102cy+ZmaTZrYAdAB9ZvZYuH/EzD5aOiBsllIfRYs4d6Kfnq597DSjp2sf5070r33Q1ybgtwbgjctw4g+h959sfEFFJFHm3niXgZn1Ao+HL8fd/WJYk+gFpsNlPT5JMJS2CPRuhXkUZpYhqPkUw58Fdx+LId9xdx+MO22lvr4+n5qaauTQ2tyDVV+/8Otw1/fCx/47pLrjPYeIJMbMcu7eV21fU+0F7j5NMOqpfNuKu9u5+1NhDaPD3Vt+eXEzGwZm3f1M+TYzS5dvayDfowTDhWNNuynevAH/++PwV78HD/0E/OPPwO59SZdKRDZJLA3LZnaAYChs0d2/Xrl/rRVlW0W4sm3G3Vd0yLv7kJnNm9mIuxcbyDcFHI877aa4NgPP/hy88hfwyBA88hSX5t/g5H/6YvPDakVkS2imjwIze8LM5oB5IAfkw4l3J2Mp3eY7RTCBsJoCjd9LIwPUuzbWetJurFe/DGc/CN/+azj6O/DBfwU7djQ/rFZEtpSGaxRmdoGgFjHEykl1R4BfNrPHt0J/RIUBogPFIMGkwrqF/TgTBB39saXdcC/8IXwuA29JwS/+Mdz3vqVdTQ2rFZEtp9EJd08QTLjrc/ez7n6x7HEm7BAZCdNtJb3AbI19BaBqR88aBsK+nLjTbgx3+OKZoLnproeDlV/LggQ0MaxWRLakRpueDq7V7xDub7fJd6n1JA5HT43EnXbD3Hwdxv4p/Om/hff8NPzC52H/PauSNTSsVkS2rEabnop1pptvMP9NF3YiQ/2/21r5pYG5ejq/15O24rgMVfpNursbGLZ6+ZvwzMfg21+Cwd+AH/znNRf16+7cw/jHH1n/OURkS2o0UNQ7+aKnwfzbwdF1DKddT9ol7j5ClVpIX1/f+ibHvDIJz/5sUKP42DPwrg+vtygi0sYabXqyte5YF87Inmsw/01X9m0+1Wxe4TyIuiborSfthvjqF+B3fwx23QFPjCtIiMgqDdUo3P2smV0ws1PAMwQT7IoEF9k+4FiYblNGPZlZjqAjer2OrWPGdYo6mqXCJqwOd69cXr2ptBvm3vfCwz8FHxmGPckPthKR1tPw8Fh3f9zMTgOfZbkpqtSoPbyZs7Dd/UhMWU1Qu7msA6hnXYwM0G9mlWVKA5hZNnw9ShDc6kobLuMev/13w2NnNyRrEWkPzS7hcQY4U7a+UyFcWXarmqZ2zSQNZGvsW1Krr6G0LIe7nyrbXPXiXyOtiEgiYlnCo3J9pxIzO1xtSY8oYd9GP8G9LSbWe3yTsgQzzCvLlCIIIBeq7Esn2nQkIrLBmlrCow5rfgOv4qcJJr1dBA6a2WPlS5VvpPCCPxIuDFhuGBiqHL4a9o3kwxnVa+kIj0nFnFZEZENF1ijM7L/Q+HISKRqYyezuj1dsei4sy2NAbqNrGOECgJmw/6VI0GcxXqPTe4I1OrnDZqTjLK8Ge9HMCsCTVQJP3WlFRDZL5P0owgX/pli5llO9OoBH3b2zwbJVK89H3f1zceXXrjbkfhQi0taauR9Fwd1/pIkTr2rTr+OYJwkWGhwnGO3zJ2W7NX5TRGSTrdVHcazJ/IcaPG6QoBbzlJkthrdb/Rrbe6a3iEgiImsU4WimhjV4fAF4wN0/BXwKwMzeR3BTpKbKIyIi67fRo56qCu+IV5W7XwReKl8ixN2fU5AQEUlGIoECmDezv1drp7u/VNE3ISIiCUkqUBjqmBYR2RJimZndIAcoq1kU3P1KcsUREZFqkqpRQLBG1ALBkhk5guaoSTP7YIJlEhGRCkkGil7gKYLZ2z8C/DLBDOeLZvZMguUSEZEySQaKIXf/VDii6aK7j7j7IEHfxWUzezbBsomISCipQFGkxrIg7l4Ml9eeMLMnNrVUIiKySlKBYoo1Fgx097NrpRERkY2XVKB4Chgys/eukS6/GYUREZHaEgkU7j4NfAKYNrOTEUk110JEJGGJzaNw9zPhvRZGzGwEGCNYMXYuTHKKxm58JCIiMUpywh3hzYDGwpsEPc7yarV54Cnde0I2wqXZ1zl5fpLCzHXSXXs5d6Kf7s49SRdLpGUlOTx2ibufcfc+d98RPt7p7p9NulzSnk6enyQ/c40Fd/Iz1zh5fjLpIom0tJYIFCKbqTBzncXwxo6LHrwWkdoUKGTbSXftZYcFz3dY8FpEalOgkG3n3Il+err2sdOMnq59nDvRn3SRRFpaop3ZIkno7tzD+McfSboYIluGahQiIhJJgUJERCIpUIiISCQFChERiaRAISIikRQoREQkkgKFiIhEUqAQEZFIChQiIhJJgUJERCIpUIiISCSt9SSyhhs3bjA3N8fVq1dZWFhIujgia9q5cyf79++no6OD3bt3N52fAoVIhBs3bnDp0iUOHjzI4cOH2bVrF2aWdLFEanJ3bt26xZUrV7h06RLd3d1NBws1PYlEmJub4+DBgxw6dIjbb79dQUJanplx++23c+jQIQ4ePMjc3FzTeSpQiES4evUqBw4cSLoYIg05cOAAV69ebTofBQqRCAsLC+zatSvpYog0ZNeuXbH0qylQiKxBzU2yVcX12VWgEBGRSAoUIiISSYFCREQiKVCIiEgkTbgTkS1penqaiYkJ0uk0c3NzpNNpBgYGYj1uaGiI48eP09vbW3X/2NgYhUIBgHQ6TSqVWpHX2NgYc3Nz5HI5CoUCp06d4ujRo6vKk81m6enpYXZ2tur5Gv1dY+PuerTZ48iRIy7xeP7555MuglSRz+f96NGjK7YdPXrUc7lc08fl83nPZDJ++vRpT6VSPj4+XjWv4eFhz2azS69HR0c9nU6veJ3P55dez8/Pezqd9uHh4aVt4+PjPjAwsCLf3t7eFcc1+ruW1PsZBqa8xjUl8YuaHgoUrUyBojVlMplVF/BcLrfqotvscel0umqgyOfz3tvbu2Lb/Pz8irTlQaRkdHTUg+/ntfMfHR1dERga/V1LFCj0UKDYYAoUrSmVSq341u0eXKjLL8JxHFcrUAwMDFQNBCX5fN6BVd/6y7eXzjs/P181TaNlrhRHoFAfhUij/vgpePXLSZci2j3vho98sqksytvZs9ns0utiscjk5CRnz54llUrFU946FAoFisUi6XR6xfZSGaanp6v2KTR6XDUTExMMDQ0tPU+n0yvyTafTnD59etW5isXi0v5S30aljo6OpbSl9zmOMjdDo55EJNLc3ByZTIYLFy5w6tQp0uk0mUyG06dPA/D0009vanlKF9taai2C1+hxUfmMjY3R19dHsVjk1KlTK/YNDw+vCqDPPvssvb29pFKppQt85XlLr0uBLY4yN0s1CpFGNflNfSuYmJhgYGCAYrFIsVhkcHBw1TfYtS5mx44dWzNNuWw2u+obdCsp1QSKxeLSCKbe3l4GBwd58sknGR0drXpcsVhkZGSEXC63tC2TyTA2NrYUdCGoJUAQBEq1i6QpUIhITaUhmBMTEwBVh3aWmmBqqXXhbFaxWGyoyavR40pKx1YGs4GBgaWgWC3/Y8eOMTo6uuK4bDbLqVOnmJ6eJp1OMzU1tSL/UoBttszNUqCowswyQAoohj8L7j4WQ77j7j64xnl7yjbNuvuZZs8r0qzx8fFV4/aLxSKFQmFzx/OzfKGem5tbcfEsXVRrfQtv9LhKpQt95YW79HpqamrVezI0NMTQ0FDV9yqbzVIoFJiamqKvr2/FvIzS82bL3CwFigpmNkzFBdrMhs0s3cxF28yOAjX/o8xsFJh096HybWaWcfeRRs8rEoeJiQmOHz++alsqlVqzmejYsWPrOtfw8HBknqWJbZXNWaX2+lqdu40eV021fMrPU25kZITBwcHIgFreGV4oFJbKEmeZm6FAUcbM0kDG3Q+Wb3f3ITObN7MRdy82kG8KOB6xfxiYqxKIBoBn13s+kbhNT09z9uzZFduq1TKq2Yimp4GBgRUXVKCu2k2jx1XKZDJMTU2tyGd6enpV4BwbG1s1i7o8XbURU9lsluHh4djL3JRa42a34wMYBkZr7MsBpxvM9zRwNHi7V+1LAw6k4/o9NI8iPppHEcwepsqY/XQ6HTmXYCNVm/BWOVt5fn7eT58+ve7jytWaR1GaZV2ZT/n7MT4+7tls1vP5/NIjl8t5JpNZStPb2+ujo6Mrjinf30iZK8Uxj8KC/QJgZjlgwsuaf8r2jQIpj+hjqJFn6WtABzDu7laxfxg46u49qw5uUF9fn09NTcWV3bb2wgsv8NBDDyVdjESVRupks9mlbYVCgZ6eHvL5fGIjlCYmJigUCktt+dW+uT/66KPkcrkVZVzruGKxyNNPP02hUGBsbIze3l4GBgZWNR8VCgWy2SydnZ3Mzs7S39+/1NlfLBY5eHBFw8SSdDpNPp9fKmPpf7XUYZ3JZNb9u0ap9zNsZjl376u6T4FimZk5MORV+iLCC/qqZqk68jzt7mfMbIDqgSJH0Fl+LGyiGgCm3b36bJw6KFDER4GiupGREYaGhpifn0+6KLKGOAKFJtytT2o9icNRTGt1RPcCc2Wd3RNAKuzIbt3B5LKt1ds/Ie1BgSIUfpuHYEhsHPmlCTqo68mvAyi6+5i7F919GhgCcmXlqnaOjJlNVT5mZmbi+BVEapqYmGBwcF2tsLKFKVBsnKO+xtyLsiDQ6+4T5fvCpqcC8Ilax7v7iLv3VT66urqaLbvImh5//PGkiyCbRIEiVPbNP9VsXmEz0poT9MrOOV0jSYFgtJRIS5mfn090prBsrrYIFGaWMzNv4LGei3CKOpqlwlpCxzo6o4tArZW95giGz4qIJKYtJty5+5GYsppg5RIa5TqAeoYSZYB+M6ssUxrAzEpjDEfD5qZCmHctxTrOKSKyYdoiUMRommAUUjVpIFtj35JqQ2theQkPdz9VsetZoHJb+Tk1zjVh7o6ZrZ1QpMXENf2hLZqeYpQFVo0jDpuTeoELVfY12zQ0AqRr5DNAMFtcErJz505u3bqVdDFEGnLr1i127tzZdD4KFGXCfoWRcHJduWGCiXjF8o3hZLl82ezrKB3hMamKcxYJhsKuqK2ETVQjlaOhZHPt37+fK1euJF0MkYZcuXKF/fv3N52Pmp4qeLAAYMbMThP0D/QQzKiuNoppgjU6ucMmp+Msrxx70cwKwJOlwBPO3C6W9V8A5LwFVo29NPs6J89PUpi5TrprL+dO9NPduSfpYm2ajo4OLl26BMCBAwfYtWuXmqGkpbk7t27d4sqVK8zPz9Pd3d10nlrCow01uoRHtaBw8vwk+ZlrLDrsMOjp2sf4xx/ZgFK3rhs3bjA3N8fVq1dZWFhIujgia9q5cyf79++no6OD3bt313VM1BIeqlHIkvKgkJ+5thQ0FsPvEosOhZnryRYyAbt37+bee+/l3nvvTbooIolQH4UsqRYU0l172RG2tOwwSHftTa6AIpIIBQpZUi0onDvRT0/XPnaa0dO1j3Mn+pMtpIhsOvVRtKE4+yi2U8e1yHamPgqpS3fnnm3XUS0ia1PTk4iIRFKgEBGRSAoUIiISSYFCREQiadRTGzKzGeDlpMuxxRwCvpN0IdqA3sfmJfUevt3dq94eU4FCBDCzqVpDA6V+eh+b14rvoZqeREQkkgKFiIhEUqAQEZFIChQiIhJJgUIkkPhNotqE3sfmtdx7qFFPIiISSTUKERGJpEAhIiKRFChERCSSAoWIiETSjYukbZlZBkgBxfBnwd3H1pnHKFAAsu5eMLM0MAAccfdTsRa4xZhZChgGxtf7vpXl0fTfYKtr9n1shc+gRj1JWzKzYWDW3c9EbasjnxzQW7F5wt0H4ylp6zGzLNABTBJc4I41eIGL5W+wVcX4Pib+GVSgkLYTfuPKufvBKvvmgQfcvVhnXsPAOJAON025+3RcZW11ZuY0cIGL82/QDhp9H8NjE/8MqulJ2tEpYKLGvgKQAer9Rjvr7rXyktri/Btsd4l/BtWZLe1ogOBiVE0BaNtmoxaiv0EbUaCQdtQLzNbYVwDWvda/maXMbCBsUpG1xf432O6S/AwqUMh2lFpH2s5w5E4fMAWkzGw8HMkijUslXYAtJPHPoPoopK2U/fMUY8z2QlnH63TYuZgDemI8R9vYoL/BdpfoZ1A1CpEI7j5UOTon7FjsCL/liWyoVvgMKlBIWyn7h0pt8KkKwLENPseWtIl/g+1u0z6DChSy3aSIp0mkwPK4dlmfFGqWisOmfQYVKKQdTVC77baDoENwTWaWC9uCZf1i+Rtsd63yGVSgkHY0Te1vWmlgtM58UkA+Ip9tM0O7AXH9Dba7FC3wGVSgkHaUpco4/XA0Ti9wocq+ahe1rLuvui1lmLY3PM+2V+u9Y51/g+2ulT+DChTSdty9AIxUqbIPA6tGkISLruXNrHLhtWp5QPDPeSbpZRU2UUetHbXeu/X+DbaJdb+PtMhnUIsCStuqWOK6B5istihb+I94FBgML3Dl+1LAJ8KXKYJ/9mw7BwkzOw30E3xjTRO8f1MEnacrLvJR7124v66/QTuK631shc+gAoWIiERS05OIiERSoBARkUgKFCIiEkmBQkREIilQiIhIJAUKERGJpEAhIiKRFChERCSSAoVIGzCz0+FKo25m8+s4bjQ8Jm9mWqhPqlKgEGkD7n7G3Y8AIwT3VD661jEV91wedHfdiEmqUqAQaS95gntBnKoj7QAwHj6f27ASyZanQCHSfrLAQEWNoZoOFCCkDgoUIm2mbHXWx2ulCe9noLvMSV0UKETa0xjRzU8D7q479EldFChE2lMW6K1x1zSRdVGgEGlD4U1tigQ3w1khvIta2954SeKnQCHSvi5Qvfmpr9rd6ERqUaAQaV9ZIF3lPswa6STrokAh0qbCzuoCcLy0zcwGAHViy7ooUIi0tzEgU/Y6rWYnWS8FCpH2liVY0mMgfK1mJ1m325IugIhsHHcvmFkBOGZmoNFO0gDVKETaS2eVbVmCWdppdy9ubnGkHShQiLSJcG2naqvGjgApqjc7dWxgkaRNKFCItAEzOw3kCIbD5suXGQ9rEWNla0BhZunw/hND4aaLuh+F1GLunnQZRESkhalGISIikRQoREQkkgKFiIhEUqAQEZFIChQiIhJJgUJERCIpUIiISCQFChERiaRAISIikRQoREQk0v8HsZM+I7zF4e4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "m_mean, sigvm, eta, eeta, yeta, eyeta=cts.eta_function(cluster,plot=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, as per Bianchini et al. 2016, the state of energy equipartition can be probed by the meq parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAD6CAYAAABpjtoAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAApmUlEQVR4nO3deXxcZaH/8c/TUpZuTNMGLWApE0DwKmAWEAFBSRDB7UJCXdAqpYmACvLSht6fXvVetSQol8VtUqv2ikoWcENRZgqUXdNEQMXtZpB67++qpelAFyilfe4fz5k0y0wyy5mcWb7v1yuvNmfOPOfpJJ3vnGc11lpERET8MCPoCoiISPlQqIiIiG8UKiIi4huFioiI+EahIiIivlGoiIiIbw4IugISrEWLFtmlS5cGXQ0RKSEDAwPPWGurUz2mUKlwS5cuZdOmTUFXQ0RKiDHm6XSPqflLRER8o1ARERHfKFRERMQ3ChUREfGNOuolL5u37mLF+n7iW3YSrp7DuuUNLFk4O+hqiUhAdKcieVmxvp+hLTvYay1DW3awYn1/0FUSkQxt3rqLphs2UrP6ZzTdsJHNW3flXaZCRfIS37KTfd7uCfus+15ESkMhPhQqVCQv4eo5zDDu7zOM+15ESkMhPhSqT0Xysm55w4Q+lVKye/duhoeH2b59O3v37g26OiLTat07F/PsC3vp/5/niQ7tZP6cQ/IuU6EieVmycDbRa84Kuho52b17N5s3b2bBggUsXbqUWbNmYYwJuloi0+boPS/xl39sZ8ni7Zy+5BmOPyacd5kKFalYw8PDLFiwgEWLFgVdFZFAHDzrAI4/YgGwgGeqZrNnz/NAKK8y1aciFWv79u3Mnz8/6GqIFIX58+ezffv2vMtRqEjF2rt3L7NmzQq6GiJFYdasWb70KypUpKKpD0XE8ev/gkJFRER8o1ARERHfKFREpKy1tbURi8UCK7uQ1y9GGlIsImUtHo/T1NSU8rFwOMzQ0FBG5SQSCdrb24lEIlmVnen1+/r6GB4eZmBggHg8TltbG83NzWPOj8ViDA4OArB161ZqampobW0dU5+Wlhba2tqor68nHo8TjUZpa2ujtrY2o39nvhQqIlLWwuEwAwMDhEKhMcdjsRj19fUZl7Ny5cqcys7knL6+PmprawmH3eTDRCJBXV0d8XicVatWjZwPjHyfPNbZ2TnmGEB7ezuJRILGxkY6OjqmLVBAoSIiZa6uri7tm2qmb7axWIxwOEw8Hs+67EzOGR4eHgkUgFAoREdHBy0tLSOBEYlE6O3tHfP8xsbGMXdOAGvXrp3WEBlPfSoiUtZGNw8ldXZ2pjyeSiKRAKCmpiansqc6J9nUlWzWSkoGQ/J4PB4vib4Z3amIVLDR7fiRSGTk+0QiQTQaHflk3NXVBUB/fz/Lli2b0NZfSmKxGI2NjRmf39PTQ2tr68hrkG/Z488Jh8OsWrVqzJ0K7A+z5PHVq1fT1NREJBIZCaTOzk5Wr1494RqJRIJ4PE4oFJpQbqEpVERSueta+Ntvgq7F5F7+GnjLdXkVMTw8TGtrKwsWLAAY0znc399Pe3s7dXV1I00wg4ODnHPOOSUdKtFolI6OjozOjcViXHzxxb6WneqcVM/p7u6mtrZ2pC+mubmZSCRCW1sbkUiEZcuW0draOqGvpru7m6amppGO+paWFtauXTvhvEJRqIhUqOQn5kQiQSKRoKmpaUxbfCKRYHBwcEybffIuplBaWlqyKj8SiWT1Sbyvry/tSKzxkq9Lpm/GmZSd6fUTiQRdXV0MDAyMOZ7seO/u7qa9vZ1QKDSmeS0cDrN69eqROtfW1tLU1MTKlSsn9McUjLVWXxX8VVdXZyvVk08+GXQVikI0GrXurWCscDhsI5HImGMdHR22sbFxuqrmu9raWrtt27aMzh3/b49EIra5uTmvsjO9fmNjo41Go2OORaPRMXWKRCIWsK2trZOWNTQ0lPLnm0qm/yeATTbNe4o66kUqXDQandAPkGyTH3882bRSipJ3XpnceQwODmbV75JJ2Zlev729nfb29gnXH92XAm4AwNDQED09PSOd+an6faqqqgAmDAQoFDV/SUnavHXXhB0nlyycHXS1SlIsFmPZsmUTjo3v5E2+KRayGaWlpSWr8zs6OjJu/kr+mzIRj8fp7x+7X3ssFhuZANnQ0DCmXymTsjM5p6uri6ampgmBEovFaGiYuKtqsrkr2Snf1tZGY2NjytdkujrsFSopGGNacTvVJLw/49baPh/KjVprU37MM8asAhYCyUbtSKprGmN6gbj3eNwYEwYagTprbVu+dZxKsbyZr1jfz9CWHeyzMLRlByvW95fsDpRBGxwcZO3atWOOpbp76enpIRwOEw6HGRwcZHh4eOScrq4uqqqqCIVCVFVVjfTNdHZ2jryZjZ7Il04hA6u/vz/jN9bm5uYJgxE6Ozvp7+9P2ameSdlTndPX10c4HB7zuifvbMLhMNFoNOXzko+Hw+GUfUw9PT1jOvwLTc1f4xhjOoCQtbbTWttlre0EGrw3/XzKbca9+ae7Zpe1tt0LnXZgrRcg44WBVcCQMcYCQ0DLdAQK7H8z32vtyJt5EOJbdrLPur/vs+57yV5y3sP4yXKxWGxCM9fAwMDIG153d/fI39vb2wH3RtzY2DhSZnJpkObmZmpra+nu7i7ov2Uq8Xh8pClovOQdyGS2bt2aU9mZnBOLxUYmQMbjceLx+MggiWRgJI+Nr/fQ0NDIz6+qqmrMBM1EIkEkEpnwoaGQdKcyivepv9Vau2D0cWttuzFmmzGmy1qbyKHcELAszWOtuLuOkXKttYPGmBYgaoxpHnfHEsOFTvLjyCZr7fQ0llI8b+bh6jkjdyozjPteshePx9NOAhw/lLatrY01a9bQ19dHW5v7DJNIJOjs7CQajRKLxUbKS07US44cy7aPohCqqqrSvqnH43G6urpoa2ub8Ek/Ho+PmcPT1tY2YS2tycqe6pzkyLtURtelt7eXzs5Ouru7WbhwIbB/5n1Sc3MzfX19I011Q0ND9Pb2TutcFeM68gVG7hjC1toJDbvGmAGg27tzybbcVbgmq15rrRn3WIe1NuVHJGPMNlxoNI06tiqXOqRTX19vN23alPH5TTdsHPNmXlM9N5BmJz+a4X7/+99zwgknFKiGlSEZHOObrZLDYZOh0tbWRlNTU0nPb6kEmf6fMMYMWGtTLpym5q+xGnFv/qnEgayHvRhjanF3F4k0p6xK08wFsAnIfMW7abBueQM11XOZaQw11XNZt3xi5+F0WLJwNtFrzmJozflErzlLnfQBSfXJe3BwkKqqqjHLmvT09AR+pyLTQ81fY9UC6Rp+40BmiwWN1Wit7TTGpPsfNdliPlXAcKoHvCa1etwggnRB6Lvkm7kIuL6YhoYG+vr6RjqCGxsbqa2tpb29faRJLNmJL+VPoZKdUDYne/0lky4YlG40mKc2xfMXeuXGcXcyYWNMFNdZn8imfiJ+SDeiK9nW39fXp7uUCqJQ8Xif/CF9M1W25YWB4Vzf6L3gAEi1kFDPqHIHvb6gAWDiMqpjy5twp7VkyZJcqieSsWg0WrITJiV76lMpnPGjtjLmBVwH7u5jTNOWN+w4Me5YDKgaFUQTeMOj68d/VVdX51JFkYz09fXR09MzYSKhlC/dqXistQljDGTZxJWKNycln8mSvUB7lqEUB1qYorlNZDqlmkQo5a0s7lSMMQPGGJvDVza/7SEyaBrz7jKqcu0895qyeq212YZDnP1zV0REAlEWdyrW2jqfioqRvl+iCtcxPpVW3Az88XUKAxhjkuuI93rNViO8+SxD6QLFmysTSzevpRQVy7IvIuKPsggVHw2yf+2t8cJAJM1jI9JNTEwu05JuORXv8cT4QDHGNI4KnxBuWZZ09Zu2mfV+CXoNL2stXrOnSEXzayJ8WTR/+ShCismGXpNWLdCT4rG8m5ySc1jS3KGMHjYTSXWOV4daMgi9YhPksi8zZ85kz54903Y9kWK2Z88eZs6cmXc5CpVRvH6QLq9fY7QOXMd5YvRBrzlqyJs1P5Uq7zmhcWXUAm3e35tHfbV6zWGjz09VN3Bh0jm+Oa0UhKvnMMO7UZjuNbzmzZvHc889N23XEylmzz33HPPmzcu7HDV/jeMtHpl8Q0/g+liiaUZixZiiA99r1lrG/hWKNxhj4sBKL6Q2eGWkGzQw0n/ijVBbMypYQriw6ijFQAG37Mv4PpXpUlVVxebNmwGYP38+s2bNUlOYVBRrLXv27OG5555j27Ztvsxb04KSFS6rBSWthV1bYc6iwlZqGu3evZvh4WG2b9/O3r17g66OyLSbOXMm8+bNo6qqioMOOiij50y2oKTuVCRzf7wL+i6FhhVw+tUwt/QnTh500EEsXryYxYsXB10VkbKgPhXJ3GHHw6veAY9+FW46Ee7+JOzYEnStRKSIKFQkc1VhuDACV/bDCW+DR77ihcunYOczQddORIqAQkWyt+gYuLALrvwVHP9WePgWuPFEiH4adqbfclVEyp9CRXK36Fi4aK0XLufDQzfBja+B2GcULiIVSqEi+as+Di76Blz5S3jlefDgja5ZLPZZ2JVyjzERKVMKFfFP9Suh+ZtwxaNw7Lnw4H+4ZrEN/65wEakQChXx32HHQ8u34IpH4NhGeOCLLlzu+Rw8vy3o2olIASlUpHAOOwFavg2XPwLHnAP3X++Fy+cVLiJlSqEihfeyV8HF6+HyhyF8Ntzf6cLl3i/A84mgayciPlKoyPR52T/Bsu/Ahx6C8FmwscOFy33XKVxEyoRCRabfy18Ny26Ftgfg6DPhvjVutNh9HfDCs0HXTkTyoFCR4Cw+Ed71XWi7H5aeCfd9wc1z2dgJL2hJepFSpFCR4C0+aX+4HHU63Pt5L1yuV7iIlJiCLH1vjHkT+/f6AIhba+/x/UKSt6yWvp8u///XrinsT3fBIQvgtA/DqW1wUP4bCIlI/iZb+t7XUDHGfB1o9b5NAMO4YAkBFuiy1l7u2wUlb0UZKkn/M+g68//0cxcur/8InNKqcBEJ2GSh4lvzlzHm40CvtXaG91VlrT3G+3OGtXYmEDPGrPHrmlLmjqiF93TDynvgyFNgw7+50WIP3AC7dwRdOxFJwc8+lWettRsmO8FaezsQ9/GaUgmOqIP39sBl97i/b/isGy324I0KF5Ei42eoZNqOpv2LJTdH1sElfbAiBoe/FmKfduHy0E3w4s6gayci+BsqC4wxlxlj5qd60Biz1BhzGa5/RSR3r2iAS26HFVFYfDJE/9U1iz10s8JFJGC+hYq19npgAfAXY8ze8V/AABCy1n7Rr2tKhXvFKfC+O+DSu+Hlr4Hop+Cmk9ymYS/uCrp2IhWpUEOKjwbCow7FrbVP+X4hyVtRj/7K1uZH3ez8+H0w5zA442qo+yAcODvomomUlWkbUiylJ9NQ2bx1FyvW9xPfspNw9RzWLW9gycIifbN++mG3nthTG2Huy+D0q6H+gzDrkKBrJlIWpmVIcRaVuWy6ryn5W7G+n6EtO9hrLUNbdrBifX/QVUrvqNfD8h/DB34Gi46DX6x2zWKPfg32PB907UTKWhDLtNQEcE3JU3zLTvZ5N7X7rPu+6C09HT5wJ3zgpy5cfn4t3HQy/DICe14IunYiZekAvwoyxtyN66ifSi2w2q/ryvQIV89haMsO9lmYYdz3JWPpGS5cnnrA9bnctcptdXzGNVD7fph1cNA1FCkbvvWpGGNW4oYLTza5cSHQmq4tTqZfWfapTMZa+MsDcO8a2PwwzDsczvTC5YCDgq6dSEmYlo56b8TXOdbab0xx3nXW2mt9uajkraxGf2XDWteRf+8a+OujMP8IOONjCheRDExLR703ZHggg1Mjfl1Tit/mrbtoumEjNat/RtMNG9m8tUjmjxjjtja+9Ofwvh/CoUfCzz4ON9dC/zp4aXfQNRQpSb521Ftrf53BOZqvUkGKftSYMVDzRrj0F/C+H8D8w+Gn17hw2fRNeOnFoGsoUlJyChVjzNHGmO4szr/Q22NFKkzJjBozBmreBCvuhkvugPmL4c6PwS21sOlbCheRDOV6p1ILXDfVScaYQ40xW4FlwMXGmF/keD0pUeHqOcww7u8lMWrMGDjmHLeu2Htvd5Mn77wabqmDgW/D3j1B11CkqOUaKoN4qw0bY940yV3Idbg9VpZZaz8EdBpjLszxmlKC1i1voKZ6LjONoaZ6LuuWNwRdpcwYA8c2wmUxeG8fzK2Gn1zl7lwG/1PhIpJGzqO/jDGfANrYv8bXNuBN1trHR53zZ6DFWvvYqGMf16KSxaNiR39ly1r4cxTu+4Lb7jh0FLzhE3DSu2DmrKBrJzKtCjX6qwFox82QX4Br4uocd06YifNWns3jmiLBMAaOOxdW3gvv6YHZVfDjD8OX6+HXt8Lel4KuoUhRyLmjHui21t5urX3KWvustTYGXDu+Kcxa+9y4p2sFSyldxsBxb3bh8u5uOPhQ+NGVXrh8V+EiFS/XO5UQMJTi+BCZLdUiUtqMgVeeB60b4V3fh4PmwY+ucOHy2PcULlKxcgoVbz7KamPMvHEPdQBRGLmbMSmeHsrlmiJFyRg4/nxoux/e9T04aC788HL4yinw+G0KF6k4+fSp9ALPGmP+bIzp93Z3rALqjTEXeY93jR7t5TWNqU9Fyo8xcPwF0PYALPsuzJoNP2iDr54Kj3fDvr1B11BkWuS19pcxpha42Ps2aq3d4N2h1AKD1tqnjDHXAYcCCaDWWvvmPOssPtLorwLZtw/++FO3WdjffwsLj4Gz2uHVF8GMmUHXTiQvge/86N25VFlr1xb8YpIVhUqB7dsHf7jThcs/fgcLj/XC5UKFi5Ssgi8oaYyZb4w52RizNNXj3igxBYpUnhkz4FVvhw89CC3r3ZyWOy6Dr74OftOnZjEpO3mFijHmMmPMMG7i4wAwZIzZaoxZ4UvtRMrFjBnwT++EDz0ELd8GMxNuXwFfez389nZ3RyNSBvKZUd+Dm9wYYewExzpcP8tW9Z8UPzV/BWTfPnjyh7CxA7b8AaqPd81ir3qnCyCRIuZ7n4ox5jJgm7X29knOuQhYMNWmXRIshUrA9u114XJfBzzzR6g+Ac5uhxPeoXCRolWIPpUFkwUKuH4UNBFSZHIzZroRYVc8AhetA7sXej8AXz8DnvyRmsWk5OQaKokMz9uWY/kilWXGTHhNM1zxKFz4Ddj7IvS8HyJnwpM/VrhIycg1VDJtM6vJsXyRyjRjJpzYAlf+Ei5cCy+9AD3vg8gb4Pc/caslixSxXEPFTLWTozeTfjjH8kUq24yZcOLFcMUv4Z+7YM8u6L7E3bn84acKFyla+Y7+ssBtwFO4JrEQUA+0AGj0V/FTR32J2PsS/KYX7u+E4Ti8/EQ4ezW88i1uiRiRaVSQyY/W2otxc1Nu9/4cwu0I2YVbokWBIuKXmQfAye+GK/vhnV+D3dvhtndD11nwx7t05yJFw5dlWkat9xX3VjCWEqE7lRK19yV4otvduWz7Cxz+Wnfncuy5unORggts7S9jzFJr7V+yfM6FuF0lh4BYts+X7ChUStzePW6J/fuvh8TTcHitFy5NChcpmIKv/TWJSA7PeRewFdgALDDGXDR6+XwRGWXmLKh9H3xkAN5+C+x6Br7XAt84B/4cVbOYTLtJ71SMMV/D7ZGSixBQb61dmOPzx9flImBgOu5cjDGtuPonvD/j1to+H8qNWmubUhzvxS11E7HWxo0xYaARqLPWthWyfrpTmdzmrbtYsb6f+JadhKvnsG55A0sWzg66Wum99CI8/n24/4vw7GY4oh7euBpqztGdi/gm5+Yvb7HITYxd2ytTVcA5foWKV58LrbV3+FVemmt04NYt65zsWA7lNgO91toJ/7ONMQO4PqnRYmkCyNf6KVQm13TDRoa27GCfhRkGaqrnEr3mrKCrNbWXXoTHv+eFy1/hyFPg7Guh5k0KF8nbZKFywBTPjVtrz83jwj05PGcl0I7blrjXWnvPqIdzvWvK9NphoNVaO2Z5GWttuzFmmzGmy1qbyKHcELBsklNiuH9z2Pt+k7V2cLrqJ+nFt+xkn/e5a59135eEAw6Eug/ASe+Bx26F+78Et14IrzjVhUv4jQoXKYip+lRa8iy/PcfnNeHujq41xuzztiz+M4Wfod+Ge4NPJQ605lhuK9A9yeNbrbUxa22X9zUhUApcP0kjXD2HGd577wzjvi8pBxwI9ZfCRwfhghvg2f+G7/wzfPM8GLpXfS7iu0lDxVr7VD6F5/j8OHC0tfZ6a+251toZuKX0z7XWrs6nPhloJH1TXxwXdlnxtlyOkfl6aZPxvX4yuXXLG6ipnstMY6ipnsu65Q1BVyk3BxwEDSvgo7+GC77kmsS+80741lsgfp/CRXwTyNraxpj56R6z1m4Anhq9DIy19tf5BlyGanEjz1KJ41YLyFbjJHceYxhjQsaYRq+Za7rqJ5NYsnA20WvOYmjN+USvOau4O+kzccBB0HCZC5fzvwjbnob/fAd863x46v6gaydlIKgNG7YZY05O96C19qlxfSnFIpTNyd4ora4MTl3onVuPGxgRMsZEvb6YgtVPKtgBB8EpK124vOV62PYUrH8bfOsCeOqBoGsnJWyqjvpCMRS40z1bo97AEz6VFwaGs+g47xl17qA3omsArx8p3/p5oTWhz2XJkiW5FCflYtbBcGor1L4fBtfDAzfA+rfC0jNdh/7SM4KuoZSYILeWswDGmJO9r7RNYiWqOdO5I9ba9vHhY62NAVVeGOTNGwBQP/6rurraj+Kl1M06GE5tg6seg/Oug2f+BN++AL79Vnj64aBrJyUkyFDpNMbsxX0aH8A1ifUbY94YRGVGvamH8i3Lm5OS92RJXD9JcsXnhHcs5EO5IqnNOgRedzlc9Ti8eQ1s+aPrzF//Nnj6kaBrJyUgyFCpBa7F9SOcC1yOa9rZYIy5LZuCjDEDxhibw1dzFpcJkUHTk9dMVWWtzWXC6Hhx9s9dmUoIn5ruRJh1CJx2hRcuX4B//AG+dR6sfztsfjTo2kkRC6pPBaDdWvvFUd9vALq8N+UOY0y3tXayCYMjrLV1PtUpRvq5MFW4TvSptAINxpjxdQoDGGOS66H1Wmtj3mz6mLU2kzk9ftRPJHMHzobTroS6D8Kmb8JDN8I33+wmT569GpacGnQNpcgEFSoJ0sy38Jp52owxK40xl1lrvzGN9Rpk4nIpSWEyWCAz3VIp3l1RY4q1vEK4FZnTXXP0cOS86yeSkwNnw+s/DPUfhP518NBN8M1z3bIvZ6+GV5wSdA2lSATV/LWJKeZUWGvXTnVOAURSXdO7e6oFJiw7M8mckoyvaa2dMOzYK7eWsUGRdf1EfHXgHDj9o3D1E9D0b/C/j8O6JvjOhfDX/qBrJ0UgqFC5Fmg3xpw0xXnpPsEXhNcP0uUN5x2tA9dclxh90Gu6GvJmzU+lyntOaNzxVNcDFyCd3iiwnOonUjAHzoHTr4KrnoDGz8L/PgbrGuHWi+C/1QpbyQq6SdekFzZmFbAGt0DiujTnrJmGpVlSXXf00vI1QH+q4cHem3sz0JSuY95r9lqGW2IlhGvCigMrkyHgBU3y3xnCBVBkdKDkUr9MaJVi8cXuHdC/Fh66GZ4fhmOaXLPYkX51d0oxCWznx6l4b7hdwKG4IbhRYNh7uA33xlrQpe4rnUJFfLV7O/yqCx6+BZ7fBse+2U2iPCKTm3kpFUUbKiOVcHctF7O/E3oIuNZae3twtaoMChUpiN3b4ZcReOTLLlyOO8+Fy+GvDbpm4oOiDxUJjkJFCuqF5+BXEXj4y/BCAo57ixcuJwddM8lDkHvUi0glO3g+vOETbrTYGz8Jmx+GrrPg++9xI8ek7ChURKTwDj4UzvoEXP0bOPtf4C8PQuQNcNt74X+fCLp24iOFiohMn4MPhbPb3Z3L2avdMvuRM124/O03QddOfKBQEZHpd0jI9a1c/QScda3bIOzrZ0D3JfC33wZdO8mDQkVEgnNICN642guXdohvhK+fDj3vh7//LujaSQ4UKiISvEMWwBv/xa2K/IZV8F/3wNdeDz3L4e9PBl07yYJCRUSKx+wqeNP/c3cuZ34c/ivmwqX3A/CP3wddO8mAQkVEis/sKjjnU2602JnXwJ+j8NXToPeDbm8XKVoKFREpXrOr4Jx/deFyxsfgz3fDV18HfSvcrpRSdBQqIlL8ZldB46fdqshnXA1/vAu+circfhls+VPQtZNRFCoiUjrmLITGz7g+l9M/Cn/4GXz1VLh9JTzzX0HXTlCoiEgpmrPIbRJ29RNw2ofhD3fCVxrgjjaFS8AUKiJSuuYsgnP/3TWLnXYlPPkjFy4/+BBsndY9/sSjUBGR0je3Gs79nLtzed0V8Lsfwpcb4AeXK1ymmUJFRMrH3MPgzZ93kyhfdzn87g4XLj+8AoZTbs4qPlOoiEj5mfcyL1yegFPb4Le3wy318KMrYfipoGtX1hQqIlK+5r0Mzlvj7lxOaYUneuHL9fCjD8O2vwRdu7KkUBGR8jfv5fCW61y4NFwGT/TALXXw44/AtqeDrl1ZUaiISOWYvxje0gFXPQb1l8Ljt8EttfCTqyCxOejalQWFiohUnvmHw/nXw0cfg7oPwmPfg5tr4SdXQ+KvQdeupClURKRyHXoEXPBFL1yWw2PfhZtfC3d+DJ7976BrV5IUKiIihx4BF3wJPvprqH0/DH4HbjoZ7rwGnv2foGtXUhQqIiJJhx4Jb73BhctrL4HB/4SbT4afflzhkiFjrQ26DhKg+vp6u2nTJl/K2rx1FyvW9xPfspNw9RzWLW9gycLZvpQtEojEZnjgS/DrW8HMgLoPuCX45x8edM0CZYwZsNbWp3xMoVLZcg2VVAGyYn0/Q1t2sM/CDAM11XOJXnNWAWotMs22Pe3C5bHvgpk5KlwWB12zQEwWKmr+kpwkA2SvtQxt2TESMPu8zyj7LMS37Ay2kiJ+WXAUvP1m+MgAnLQMNq2Dm06Cu9ph+9+Crl1RUahITlIFSLh6DjOMOzbDQLh6TnAVFCmEBUvh7bfAhzfBiS3wq7UuXH6+WuHiUahITlIFyLrlDdRUz2WmMdRUz2Xd8oZgKylSKFVHwzu+4u5cXt0Mv4x44fIvsP3vQdcuUOpTqXB+9qmoU14q1tYh1+fy+G0w80BoWAGnX+VWTS5D6qiXtPwc/SVS8bYOwf1fhCdug5kHeeFytdvvpYyoo15EZDosrIF//prrc3nVO+DRr8JNJ8Ldn4QdW4Ku3bRQqIiI+G1hDVwYgSv74YS3wSNfceES/VfY+UzQtSsohYqISKEsOgYu7IIrfwXHvxUeuhluPBGin4adW4OuXUEoVERECm3RsXDRWi9czoeHboIbXwOxz5RduChURESmS/VxcNE34MpfwivPgwdvdM1isc/CruGga+cLhYqIyHSrfiU0fxOueBSOPRce/A9357Lh30o+XBQqIiJBOex4aPkWXPEIHNvk5rrceCLc87mSDReFiohI0A47AVq+DZc/AsecA/df72bo3/N5eH5b0LXLikJFRKRYvOxVcPF6uPxhCJ8N93e6O5d7vwDPJ4KuXUYUKiIixeZl/wTLvgMfegjCZ8HGDi9c1hR9uChURESK1ctfDctuhbYH4OgzYeN1brTYfR3wwrNB1y4lhYqISLFbfCK867vQdj8sPRPu+4IbLbaxE154LujajaFQEREpFYtP2h8uR50O937eC5friyZcFCoiIqVm8Unw7u9D632w5DS493MuXO6/HnZvD7RqChURkVJ1+GvhPbfBynthyevc/JYbX+PmuwQULgoVEZFSd0QtvKcbVt4DR57iZubfeCI8cAPs3jGtVVGoiIiUiyPq4L09cNk97u8bPuvuXB78j2kLF4WKiEi5ObIOLumDyza4u5jYZ9xQ5Idughd3FvTSChURkXJ1ZD1ccjusiMLik90mYTee6PZ1KVC4KFRERMrdK06B990Bl97t5rxEP+XWFnv4Fnhxl6+XUqiIiFSKJafC+34Al/7CLQVz9ydds9jDXwZrfbmEQiVgxpiQMSYcdD1EpIIseR28/0fwwZ/DYa+C+L1gjC9FH+BLKWXGGNMKhICE92fcWtvnQ7lRa23TuMONQIcxpg8Y8o6N30ghYa2NeWX0AnEgYq2Ne4HUCNRZa9vyraOIVJCjToPlP/a1CUyhMo4xpgPYaq3tHH3MGBMefSyHcptxb/7jNQBhYNUkT+8EYt7fw0AzsMrs/2QRSxFWIiKZOXC2b0UpVEbxPvW3WmsXjD5urW03xmwzxnRZaxM5lBsClqV5OAQsSFeuMabDWts+6lAMaMeFC8Ama+1gtnUSESkEhcpYbey/IxgvDrTi7hqy1Qp04+4wxktMEiitQGTc4a3JpjARkWKjjvqxGnHhkUocyLqJyRhTiwuqRJpT1qR5XhjAWpuuPiIiRUehMlYtsDXNY3GgPocyGydrnpqkOa3NWtuV7nneqLFGjRwTkWKiUMlOKJuTveartMEwxfPGN3slLfQerwc2ASFjTNTrtxERCZT6VDyj3pQTPpUXBoaz7dj36lE32V0K0DOq3EFvxNoAUDNJua24vp3xdhhj/phNHYVFwDNBV6IM6HXMX1Cv4VHpHlCoFE5zjkOQW4FougfHjQRLHosZY6qMMa3pwsg7nvVdk0xkjNlkrc2lKVRG0euYv2J8DdX85Rn1yT+Ub1nenJRcJ0uuJv0ItMnEgZYcryki4ouyCBVjzIAxxubwlWqIbzohMmga85qvqnIZtWWMaYRJO+8nE2f/3BURkUCURfOXtbbOp6JipO+XqMJ1jE+lFWgwxoyvUxjAGJPsgO9NMd+kifRDmjHGDOBmz09oAhMRKQZlESo+GsQNK04lTPoRWSPS9aMkl2mZYn2uVMu4jBZi//pgqeqnmfUiEqiyaP7yUYQUc1G8Jq1aoCfFY342OYWZuJjkaJFUHfFeHWrJIPTEFxrw4A+9jvkrutdQoTKK1w/S5Q3RHa0DaB/f1+E1Rw15s+anUuU9JzTJOZM9Rpq6gQuTTi3fMj2mGO4tGdLrmL9ifA2N9WljlnIybun7GqA/1dL33ht8M9CUrmPea/ZahmvaCuGaqOLAyhQhNQT0TdZn4oXSau/bEC6sIgoUESkGChUREfGNmr9ERMQ3ChUREfGNhhSL4M8W0pW81bPX19cBRHPdertQ23iXknxfx2L4HVSfilS8dFtIjz+WQTkDTJznVNZbPXuTeauAftybYUuOb4a+/AxKlY+vY+C/gwoVqWjeJ7mB8VtIe49tA47OdNkc700wSoVu9WyMseTwZujnz6Ac5Po6es8N/HdQzV9S6fzcQlpbPeemUNt4V6LAfwfVUS+VzvctpCVr+hmUEYWKVDrft5DWVs9ZK8Q23hUtyN9BhYrI5EJZnKutngsjFHQFSkjgv4PqU5GK5fcW0p6st3quZAX6GVS6QH8Hdaci4hNr7YRFR71O0yrv06NIQRXD76BCRSqWn1tIT0FbPacxjT+DSjdtv4MKFZH0QvjTLKOtnnMXQk1jfpi230GFilQ6P7aQxhgzkGavG5maLz+DSlcsv4MKFal0g6T/BBcGejMsJ4S2es6VXz+DSheiCH4HFSpS6fzaQlpbPWcg3WtHlj+DSlfMv4MKFaloPm4hra2enap0D6R77bL9GVSIrF9HiuR3UAtKiuDPFtKVuNWzMWYV0ID7JBzGvX6bcB3DYwJhqu23M/0ZlCO/Xsdi+B1UqIiIiG/U/CUiIr5RqIiIiG8UKiIi4huFioiI+EahIiIivlGoiIiIbxQqIiLiG4WKiIj4RqEiUmGMMau8FW2tMWZbFs/r9Z4zZIzRIo+SkkJFpMJYazuttXVAF24P8+apnjNuj/Mma602HZOUFCoilWsIt5dJWwbnNgJR7+/DBauRlDyFikhliwCN4+5EUqlCYSIZUKiIVLBRqwBfnO4cbz8O7b4oGVGoiEgfkzeBNVprtXOlZEShIiIRoDbNboIiWVGoiFQ4bwOnBG7jpzG83QXLdpMx8Z9CRUTA7QOfqgmsPtUujSLpKFREBFwTWDjFvuca8SVZUaiICF5HfBxYljxmjGkE1EEvWVGoiEhSH9A66vuwmr4kWwoVEUmK4JZtafS+V9OXZO2AoCsgIsXBWhs3xsSBFmMMaNSX5EB3KiKVa2GKYxHc7PqwtTYxvdWRcqBQEalA3lpfqVYn7gJCpG76qipglaRMKFREKowxZhUwgBtCPDR66Xvv7qRv1JpgGGPC3v4p7d6hDdpPRdIx1tqg6yAiImVCdyoiIuIbhYqIiPhGoSIiIr5RqIiIiG8UKiIi4huFioiI+EahIiIivlGoiIiIbxQqIiLiG4WKiIj45v8AQCVRQbFO+gEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "m_mean, sigvm, meq, emeq, ymeq, eymeq=cts.meq_function(cluster,plot=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, as per Bianchini et al. 2018, MNRAS, 475, 96, if you can measure meq you can also measure the kinematic concentration ck, which is the ratio of meq within the half-mass radius to meq at the half mass radius:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0000000000004985"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cluster.ckin()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the cluster's galactocentric position and velocity are known, then its tidal radius and limiting radius can also be determined. Note when calculating the limiting radius, it is possible to select ``plot=True`` to see how the cluster's density profile reaches the background density. \n",
    "\n",
    "Note that you need to specify a ``galpy`` potential in order to calclated these values. For a Milky Way-like external tidal field, MWPotential2014 from Bovy J., 2015, ApJS, 216, 29 is useful. Different potentials can be used by assigning a ``galpy`` potential to the ``pot`` variable.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tidal Radius:  10.6993452944\n",
      "Limiting Radius:  5.35685836913\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAD6CAYAAAC/KwBlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAuW0lEQVR4nO3deXyU9bX48c/Jwk4yhF32EAUVBJIgKmpdgtLNBVlstbZWAZfbXlfAe29/9d72FoO21S5WwLq0xSub2tq6FFDrDgmRxQUVwioimhCQTUJyfn88z8CQzAxZnplnJnPer9e8ZphnOyrm5Ludr6gqxhhjTCRpfgdgjDEmsVmiMMYYE5UlCmOMMVFZojDGGBOVJQpjjDFRWaIwxhgTVYbfAZjkIyILgaKQryqBKvdzbsj35e57AMgJ/V5VC0LuFwByVbXM61iNMc0nto7CNJaIbMBJDJPr/nAXkWlAMTBHVafWOZYLzAYKVbVTyPcrgXxgjKoujXH4LZL7730gsAHoDKCq0+N5P7umBVNVe9mrUS9gF04LINyx8YACxRGOB5y/dsd8VxztnvY67n+PJcCUOt8VASvjdT+7pmW/bIzCNEWpqpYf/7T6VLUKWOR2NwW/m66qnZp6z1QmIsUAqjon9Ht1WmalbjdhTO9n16QAvzOVvZLvRYTWgnssaovCPWca1nrw6r+FAuMjHMt3jwdieT+7puW/rEVhGsUdZ9jQzNuU43RBmWYQkfHux0iTAIIttCmxup9dkxosUZjGqgRKm3mPMo7OkjJNNwZAI3TZqdPNd+S8GN3PrkkBNj3WNIr7P0qzprGG/k/o9vfmuq/Jqroo5NgS9/uc4DERmYLTGumM0w0wO3iNiORzdNruQPdZx8y8CuWOkxTjJK0K95qVWqdvOoHlHv8Uyht4XlPvZ9ekAEsUxm+TgYk402brmoDTxC+GI4OMM4O/1QW7wUQk+NtdQFVnBS8WkYUislJD1myEHMsHFuJMyQ1NXEtEpCBagkkgOcc/Jebn2TUpwLqejK/cH/oLohwLtjCm4rQeqkKOl+O0bqYD+aGtEddsIN9NKHUtwxlwr9vFMBWYIiJFYa5JNAEa1oUXiOH97JoUYInCJLpK970qQr9xOU53U90kETwGThfVEW7LJBCui6lO8mkwEZkmIurBK9CIxzbkt9qqGN/PrkkBlihMsoi0xqISIg8+RjAlyv2CzypsxP1Q1VmqKh68qhrz3AYIJMD97JokZ2MUJllEm5Jb1dCbuL+xB4BKd2A8nOA5ia6Shg2qVsXwfnZNCrBEYVJN8IfA0iizm5Jl1lNVA87J4Wj3XSzuZ9ekAOt6MilFjxYxbAmzVkrhSCspkgANn87clPvZNSnAEoVJRZ7Pg/dpMHul+x72nyXkXktieD+7JgVYojCpaBF1ZkLVFWX8IiyfBrOD04ojDbwHp/g2tHR7U+5n16QASxTGa/Hu0mn089TZU6A8WCW0rpAB74TmJpVZOGs/wpmEsy9IvRle7oLDZt/PrkkRflcltFfLeuFUhlWcxXENvSbXvWZKlGPTIly7BNgV4Vj+ce67gTpVQnESRMTKt4n4wukuqbt3wnj3ny/Q2P9Gjb2fXdPyX7bDnWk2dyewSTg/fAMhh6pwxgNKNUJJDHd3u9DryoECVa2KcGyqqi51nzmVo33JVTgzmSa4LYJlYa6doPV35Ct2zwv+hlihIWVAkkVjdmNzV50vBKZrhJlfibyLXCJf01JZojDGGBOVjVEYY4yJyhKFMcaYqCxRGGOMicoShTHGmKgsURhjjInKigLGUZcuXbR///5+h2GMMfWsXLnyC1XtGu6YJYo46t+/P6WlpX6HYYwx9YjI5kjHrOvJGGNMVJYojDHGRGWJwhhjTFSWKBLY+p17eW7tp1TX1PodijEmhdlgdgJbuHIrs/9VTteOrZlU2IcrT+9D707t/A7LGNNM175wLQCPjn3U50gaxhJFApt28WBGDchh3ttb+P0r6/n9K+s5f1A3rhrVl/MGdSM9TfwO0RiTAixRJLD0NOGCwd25YHB3tu3az/ySrcwv2cp1j5dyQnYbvnN6XyaN7EO3rDZ+h2qMacEsUSSJ3p3acftFg/jxhSey7IPPmLd8C79c8hH3L/uYi07pzlWj+nHWwM6kWSvDGOMxSxRJJjM9jbFDejJ2SE82fbGP/1uxhQWlW3n+3R3079yO747qy/iCPuS0b+V3qMaYFsI2LoqjwsJCjcXK7IPVNbz43g7mvb2FFZsqaZWexjeG9uC7o/oxsn8nRKyVYUwieXLdkwBcOfhKnyM5SkRWqmph2GOWKOInVoki1EeffckTy7eweOU2vvzqMCd268BVo/pyeX5vsttmxvTZxpjkZYkiQcQjUQTtP3SYv6/+lHnLN7N6227aZKZxybATuGpUP07rnW2tDGN8dODwAQDaZrT1OZKjLFEkiHgmilBrt+3miRWbeead7RyormFIryyuGtWPS4adQPvWNkxlTLwl4jqKaInCVmangKG9s5k57jSW/+eF/OzSU6k+rNz11FpG/WIZP3nmXdbt2ON3iMaYBGa/TqaQrDaZfO/M/lx9Rj/Ktuxi3ttbmF+6lT+/vZmCfp24alRfvjG0J20y0/0O1RiTQJrcohCRLBEZJyIXeBmQiT0RoaBfDr+aNJzld13If33zZCr3HeK2Bas5Y+Yyfv739yn/fK/fYRpjEsRxWxQikgUsAgqAJcAUQICNQMA9ZyVwoaq2uD4MEckFioAJqjrG73i81ql9K64/J5frzh7AWxsqmLd8C4+9uYmHX9/IWQM7c9Wofow5pTutMqyX0phU1ZCup4eBcpwkMRJYCJQB04GlQCdgIlAM3BibMP0hIvlAIU5CzPE3mtgSEc7K68JZeV3Y+eVBFpZu44nlW7j5iTK6dGjNpJG9uXJkX/rkWFFCY5rr0rxL/Q6hUY4760lE7lHVGSF/vgIYoKr31TnvTlW9NzZh+ktExgN3qWpBc+7j16ynpqqpVV79+HPmvb2Fl9Z9hgLnndSVq0b14/zBVpTQmJYk2qynhrQo1of+QVUXi8jkMOclxDxbEQngtG6WqOqiKOdNwWkpVLnv5dHOT0XpacL5g7px/qBubK86wJMlW3lyxRau/5NTlHDSSKcoYY9sK0poTGPsOrgLgE5tOvkcScM0JFGIiAwHioKtCFWd6x7IAnJUdRPOuIVvRGQ2TvdQCc44ypIo5xYDFao6K/Q7EckN/c4cdUKgLbeNOYkfXZDHsg92Mm/5Zn699CN+89LHFJ3cjatG9ePsvC5WlNCYBrjtlduAxFpHEU1DEsUCnN/QBwD31Tn2H8CdIjILuMfj2BpFVacGP7uJICx3cHqKqh6TylV1uojsEpE5qloVu0iTm1OUsAdjh/RwihKWbGFh6TZefO8z+uY4RQknFPSmc4fWfodqjPHIcaeyqOpuVb1BVS8Oc3gmUKiqd6nqbu/Di4mpOIPw4ZTjtEZMA/Tv0p67vn4yb911AQ9cOZwe2W245/l1nDnzJX78f++wvLwCW/lvTPJr1pxHNzlsEJGHRKRCRG6ve46IzHS7qBJFEU5CCKccaHFTYGOtdUY6lw7vxYKpZ7Lk1nP57qi+vPzhTibNeZsxv36VR9/YyM4vD1rSMCZJebEy+y5VvUFEBgJXAr+sc/we4C73lQjygfkRjlmLoplO7N6Ruy85leljB/Psmu3MW76F/372ff772ffJbptJXrcODOzanrxuHdzPHejdqZ3NoDImgXmRKCoBIi1GU9XdIrJARC5Q1Zc8eF6sBcJ816LXUMRC21bpTCzsw8TCPrz7yW5WbKxk/ed72bBzLy+t28mC0m1Hzm2VkUZul/YM7NaBvK4djrzndm1v5URMizRp0CS/Q2gULxJF7vFOUNV3ROQOwNdE4U6dBWdKbEPOz8UZ0ygC8t1B8g2qOuc4100hTMukb9++jQm3xRjSK5shvbKP+a5q/yE2fL6X9Tv3suHzfazfuZe123bz3NpPCfZQiUDvTm3J63q09RF872Q7+JkkNnbAWL9DaBQvEoWIyPmq+vLxzvPgWXGlquU4K9Abe90coF4yKSwstE56V6BdKwr65VDQ79jG2sHqGjZ+sc9NIEcTyZsbKvjqcO2R8zq3b8VAN2nktM/kcI1yuFY5XFPrvivVtbXUuJ8P19a63yk1tbVU1zjn1tQq1TVKTa3Sr3M7fnBWf84c2Nn26zAxtWPfDgB6tO/hcyQN40WimA6sFJE/qGrd8YlQvv+QVNUq9wdAwOdQTARtMtM5uWcWJ/c8dv5DTa3yya4DIcnDeX/+3U/58uBhMtKEzPQ00tOEzHQhIy3kc3oaGWlChvt98HO7VhlHvktPg9JNu/jn+58xuEdHfnj2AC4dfgKtM6zry3jvrtecIduWtI4iKncMYiJQKiKTgNmq+scwp+Y191lxEKCB3VImvtLThL6d29G3czvOH9wtJs84WF3D31Zt55E3NjJt0RpmvbCOq8/ox9Vn9KOLrQsxKcyTkqCqWoZTPC8dmCsiNSJSIiJ/cF8VOAv3EsFSYGCEYzlA8hRjMp5qk5nOxJF9eP7fz2He9aM4rXeA+5d+zFn3vMS0RattgyeTsjzbuMhNFgXuQO4EnLLkBTiVZicm0IynMpwpsuHkArPjGItJQCLC6LwujM7rwobP9/LoGxtZvPITFpRu4+y8Llx39gC+dlJXK1diUobnO9xFGshNILOBlXW/dGdE5ZM4LR+TAAZ27cDPLxvKHRcN4okVW/jTm5u59rEScru259rRAxg3opftO25avOOWGW/wjZzCgbk4ffylfm5iJCIKTI00jTVYC0pVp4d8Nxtn6mvMigImW5lxU191TS3Prf2UR17fyOptu2nfKp1Lhp/ApJF9GdY722ZLmQZ5ZesrAJzX5zw/wzhGtDLjzU4UIpINLONod07whhuA6ar6dLMe0PA4puFsrJRPSMLCWW09vW6hvzplxgcCJbEuM26JouVQVcq2VDG/ZAvPrv6UA9U1DO7RkUkj+3D5iF4E2tk6D5NcYp0oFuCU9g4W2it0XxNwfhCXqurpzXpIC2GJomX68mA1z67+lPklW1i9bTetMtL4+pAeTBrZhzMGdLaxDFPPxt0bARiQPcDnSI6KdaI4Zge8OsfycUqU10aoPptSLFG0fO9v38P8ki08/c4n7Dl4mH6d2zFpZB+uHNmXHFtNblzXvnAtkFjrKKIlCk+mx0aiqmVuDailIjIzls8yJhGcckIW/33pEFb8ZxH3TxpOj6w2zHrhQ84ufolfPPcBO7886HeIxjSaF4liiYiMi3aCu5e2tb9NymiTmc5lI3oxf+qZ/PPWc7nolO48/Fo55xS/zE//+i7bqw74HaIxDdbsRKGqy4BcEfmFiHT0ICZjWpSTunfk/itH8NLt53HZ8F7MW76Fr937MjMWr2FzxT6/wzPmuJqdKETkTmAWMAOoEpGP3dXYlwc3LHJbHOub+yxjkln/Lu0pHn8ar9x5Ht85vS9PvfMJF/zyX9w2fxXrd37pd3jGROTFYPZDOAPWAZzpqUXuK8DRqbLl7jlLVXVTyLUvptIgtw1mm1Cf7TnI3FfLmbd8CwcP1/CNIT351mk96dyhNTntM+nUrhWBdq1sU6cW6K3tbwFw5gln+hzJUTGf9QT8ou4COxEZgLOtaAFwIc7ahuDDluKU0piiqp2bFUASsURhwqnY+xWPvLGRx9/czN6vDh9zTAQCbTPp1L4VOe1a0al9KwJtM8lITyNNIE2ENHHKjqSJMKBLO8bl97bV4qbRYpoo3AfciVNQr0RVn4pwTjZO4ggmj3xAVTVl6jhbojDR7P3qMJsr9rFrXzWV+w+xa98hKvY575X7D1G59xC79h+ian81NaqoKrUKtarU1jqf9351mKw2GVx1Rj9+cFZ/ume18fsfy4SxrnIdAINzBvscyVExTxQhDxqgqhsbeG4AZzFeMpQf94QlChNrZVt28fBr5bzw7g7S04RLhvVi8rkDGNwj6/gXm7hJtnUUnrZPG5ok3HOrgjWXjDHeyO/biQevKmBLxX4eeWMj80u2srhsG+ec2IUbvjaQ0Xld/A7RJKEGzXoSkQuCM5i8pKpzvb6nMQb6dm7H3Zecylt3XcCdFw9i3Y4vuerh5dz9t/eorqk9/g2MCdHQ6bEXAWUiMv94i+uMMYkj0K4VN5+fx+vTz+e6swfw2JubuOaPK6jY+5XfoZkk0qBEoaoz3LGEOcBFIlLhrpW4ILbhGWO80DojnZ986xR+PWkYZVt2ccnv3uDdT3b7HZZJEk0ezBaRK4BJwAhgETBfVVd5F1rLY4PZJhGs3babqX8upXL/IYqvOI1Lh/fyO6SUs2rnKgCGdxvuaxyh4rEfxUScsuIDcHaQWxS6sM44LFGYRPHF3q+4aV4ZKzZWMuXcXKZdPIiM9JjWCDUJLqbVY1V1t6rOVdWLcPahEGCRiJSIyPWxGAQ3xjRPlw6tmXf9KK45sx9zXi3n2sdKqNp/yO+wUsaqnauOtCqSgae/QrhJ4143K00E8nAGwV+0QXBjEktmehr/c+kQiq8YyvLySi753Rus2+HbDsYp5YGyB3ig7AG/w2iwmLU1VXVjyCD4DJxB8Ep35pQNghuTICaN7MuTU8/gYHUN4x58k+fWfup3SCbBxKVTUlXfUdUbVDUHZ+bUDSEzp4bHIwZjTGT5fTvx9x+dzeAeHblpXhl/XfWJ3yGZBBL30StVXaaqE91igEuBqfGOwRhTX7esNjwx+QzOyM3h9gWreXndTr9DMgnC12kOqrpYVW/0MwZjzFFtMtOZe00hJ/fM4oa/rKRkU6XfIZkEENNEISL9Y3l/Y4z3OrbJ5LFrR9KrU1t++FgJ72+3AW6vTT99OtNPn+53GA0W6xbFQBF5yGY8GZNcOndozZ+vG0WH1hlc88gKNn1hW7Z6aXDO4IQqMX48niQKt2jgHSIyLrQV4Y5H3ADsFpHrvXiWMSY+egXa8ufrRlGrytV/XM6O3Qf9DqnFeGv7W0d2uUsGXuyZfQ/OoPQsnFIeG9wZTQ+KyDBwEgbO1qjGmCSS160Dj107kl37DnHNI8ttUZ5H5qyZw5w1c/wOo8G8aFHkAp1UNQ0YiLPQbhFwJc5iuxoReRFnZztjTJI5rXeAud8vZFPFfn7waAn76mzXalo+LxJFiaruhiOL7Bar6lR3zcTFwH3ARuAGD55ljPHBWQO78NvvjGDNtiqueWQFO/dYN1Qq8SJRdI50QFWXqup0d7Fdg3e/M8YknotP7cFvv5PP+9v38M3fvs6KjTZ1NlV4kSjmi8gdHtzHGJPgvnlaT565eTQdWmfwnblv88fXN9LcCtQm8TW7zDiAiFwIjLfFc9FZmXHTUuw5WM3tC1az5P3P+NZpPSm+4jTat87wO6yksXG308EyIHuAz5EcFdMy4yIyGVgCTA0OXIvI7VbDyZiWK6tNJrOvLmDa2EE8t/ZTLvv9G2z4fK/fYSWNAdkDEipJHI8XXU9j3NdU4GGcMYt7gZV1Ekd/D55ljEkQaWnCTefl8efrRlGx7xCX/u4NXnh3h99hJYVXtr7CK1tf8TmKhvMiUZS7C+vmurOdCoFOHJ3xdCRxePAsY0yCGZ3XhWd/dDYDu7bnhr+sZObzH1BdU+t3WAnt8fce5/H3Hvc7jAbzIlHMFJF73JXZ/eHIBkbBGU+F7hqLsH1fxpjk1yvQlgU3nMl3R/Vl9r/K+dZvXmfl5l1+h2U84tVWqDNU9T6clkSk82x6rDEtWOuMdH5x+VDmXlPInoPVjH/oTf7rmbXsPlDtd2immbzeCvWd0D+LyGQb1DYmtYw5pTtLbvsaPzirP08s30LRr/7FP9Z8atNok5gXs57GRUoGqjoX6Ox2Tdn2p8akiA6tM/jpt0/lmZtH061ja25+oozrHi9l2679fodmmqDZ6yhEpBQY4f6xDGeq7BJVfbnOeXe43VMpy9ZRmFR0uKaWx97cxK+WfIQq3DbmJK4d3Z+MdF/3TfPVjn3O7LAe7Xv4HMlR0dZReLXgrghniuwIoMj9WoFynORRDuSq6qRmPyyJWaIwqWzbrv389K/vsWzdTk7pmcXMcUMZ1ifgd1jGFfNEEeaB+TgJo9B9rwAmqOoqzx+WRCxRmFSnqrzw7g5++rf3+GLvV1xzZn9uv+gkOrbJ9Du0uHph4wsAjB0w1udIjoqWKGKy5l5Vy3BaEsEA7sFpVRhjUpiI8PWhPRl9Yhfue/FDHn9rE39fs51bik7iypF9UqY7av6H84HEShTRxOW/iqrOAO6Kx7OMMYkvq00m/3PpEP5682hyu3Tgv555l68/8Bovr9tps6MSkBeznrJEZKa74G54lFMrmvssY0zLclrvAPOnnsFDVxdQXVPLtY+V8L0/ruD97Xv8Ds2E8KLraRHOLne5gIpIObAQZ3vUUlXdIyJZwEgPnmWMaWFEhLFDenDB4G785e3NPLDsY77529eYUNCb2y8aRPesNn6HmPK86Hpaoqp5OKuyJwGrcHazWwLsEpEaYJf7Z2OMCatVRho/PHsAr955PteNHsDT73zCefe+wv1LP2L/Idt+1U9erKO4E1ioqpvqfH8hTisjACytu2o7FdmsJ2MabnPFPopfWMdza3fQrWNrvn9Wf7p1bE1W20yy2mSS3TaTrLYZZLXNpEOrDNLSxO+QG2zXQacOVqc2EasexV081lFMBrJxEsKqZt+whbJEYUzjlW6q5Of/+IBVW6sinpMm0DEkeXRq14ofX3giI/vnxC/QJBe3dRQiMsJaDpFZojCmaVSVPQcOs+dgNbsPVLPnQDV7Dlaz58Bh588Hne+cz4d5f/se9h86zDM3jya3awe/w6/nmfXPAHBZ3mW+xhEqbusoLEkYY2JBRMhul0l2u0z6NOD8rZX7ufT3b/DDx0p4+qbRdGrfKuYxNsZf1/8VSKxEEY1X02PvEZGS0OmxItLfdrUzxvihT0475nyvgO1VB7lx3koOHbaNlJrDi1lPdwEzcXawywUQkYeADcAGN4EM8+A5xhjTYIX9cygeP5S3yyv5yTPv2kK+ZvCi60lUdTfOlNjgbKcpODOhJolILvCQiIxXVVtFY4yJm8tH9Kb883389qX15HXrwORzc/0OKSl50aKom6bHuN9NB1DVcmACVsLDGOODW4tO4htDe/CL5z9gyfuf+R1OUvIiUXRyV14HjQeqQtdVuC2O5JnkbIxpMdLShF9OGM7QXtn8+5Pv8N723X6HxINFD/Jg0YN+h9FgXiSK6cBLbq2nP+CMUywIc94XHjzLGGMarW2rdB6+ppDstplMfryUnXsO+htPRlvaZrT1NYbGaHaicFsLFwJ5ON1Oi1T1RnBmRIWc2rm5z/KDiOSKyBQRsRIkxiSxblltmHtNIbv2VzP5T6UcrK7xLZYn1z3Jk+ue9O35jeVJmXFV3a2qN6hqnqpODDn0koh85E6TrfTiWfEUsgFTALAlnsYkuSG9srn/yuGs+WQ3ty9YTW1t02ZCHa5p3nTbFze9yIubXmzWPeIp1vtR3ANsAsar6r0xfpbnVLVMVedgmy4Z02JcfGoPZowdzD/Wfsr9Sz9q8HXVNbW8+N4Orn+8hME/eYGbnyhrcquk9LNSSj9LnioNMdnhLkhVF+GUIW8SEQkAxTgVaiPeR0Sm4PzWX+W+l0c73xiT2qacm8uGz/fym5fWk9u1A5eN6BXx3A93fMnC0q08/c4nVOw7RLeOrfnG0J48u2Y726sO8PA1hXTu0DqO0cdfgxOFiPSvWyE2VkRkNk5XTwnOmoyI4wMiUgxUqOqs0O9EJDf0O2OMCRIRfn7ZUDZX7GfaojX0yWlLQb+jvcu7D1Tzt9XbWVS6ldXbdpOZLhSd3J2JhX0458QuZKSn8Y2hPbhl/iouf/BNHr12JAMTsKaUVxpcFFBEXlTVi2McT7jnKjAhXAvBXcy3UlXr1eoVkV3AAFWt8iCG8cBdqlrQnPtYUUBjEsuufYe47ME32HvwME/fNJotlftZULqVF9/bwVeHaxncoyMTC/tw2Yhe5ISpF7VqaxXXP15CdY0y+3sFnJHbsDk7Qx8fCsDa76/19J+nOTypHisiFaoa95lLx0kUxUCuqk4Ic2wlMD/YqnB/2E9qwCPn132WJQpjWq71O/cy7sE32HeohppaJatNBpeN6MWEgj4M6ZWFSPQlYFsr9/ODR1ewpXI/s8afxuUjescpcm95VT22k4h8hDPmcGSbUy8CbIYiN5ZwynGm686C5o+XGGNaprxuHZh7TSF/emszY4f0YMwp3WmTmd7g6/vktOOpG0dzw19Wcuv81WytPMCPLsg7boJJJo2Z9bQb6ALM4Og2pyUi8gcRGVdnzUQ9oZVlPZQPVEQ4Vg6EzY7GGBNqVG5nfn9VPt8edkKjkkRQdrtMHv/h6YzL78WvlnzEHQvXRK1Y+9Dqh3ho9UPNCTmuGtOiyMcpz1GOM9B8Ec5CuwJgKqAiUo6TRJYCZXUGv4uBeI9xBDy6j62hMMZE1SojjV9OGEa/nPb8eulHbK86wEPfKyC7bWa9c5d/uhyAG4bdEO8wm6TBLQpV3eiuhagCslV1gqrmAJ2AicDDOPWcbgAW4pQYrxCRJ0XketwS5F5xp87ixhMT7qrsYpxEmO/OpprSgOumiEhp3dfnn38eq1CNMQlARPj3ohP51cRhlG6u5Io/vMnWyv1+h9VsjV5HoarLgGXuPtkVqvoUTt//IgARycYZOxjpvk90X0lXDN6tfDu9CdfNAebU/b6wsDDp/h0YYxpvXH5vema3ZeqfS7n8wTd4+PsjGd4n4HdYTdbkldmqOhd4R0TuDN3Jzi3nsVhVZ6hqoaqm4SSKqmZHe+zzg/cLeHlfY4zxwpkDO/PUTaNp2yqdK+e8xQvv7vA7pCZrVgmPkO6oArd7KdJ5i4CNzXlWEwSIYbeUMcYcT163Djx902gG98jixnkrefi1clSVQOsAgdYBv8NrME9KeKjqYgC3O2qDqr4U5rT5XjyrjqXAwAjHcgBbtGCM8VWXDq15csoZ3Dp/FT//xwdsqdzPvd/6JRnpsS615x1Paz2p6lwRGSAid+CUG98UciwWRQHLcGZjhZMLzI7BM+Pr0W/W/+7Uy+D0yXBoP8yrt9YQhn8XRlwF+ypgwTX1j4/8IQy5AnZvg6em1j9+1r/BoK/DFx/Ds7fUP37uHTDwfPh0DbwQZuPCC/8f9B0FW5bDsv+pf3zsTOh5Gmx4GV69r/7xb98PXU6ED5+HN39X//i42ZDdG95dDCWP1D8+8U/QvjO8Mw9WPVH/+FULoVU7WDEX3num/vFr/+G8v/Eb+KhOhc/MNnD1Yufzv2ZB+b+OPd6uE0z6i/N56d2wteTY41knwBVznc/Pz4AddVbmdh4Il/zG+fy3H0PFhmOP9xgKX7/H+bx4MuzZfuzxPiOh6G7n8/yrYf+uY4/nfg2+Ns35/JcroLrOvgwnXQyjf+x8tr979Y838e9eG+BBlPtG/y+/f2Mzw3csYlyrkvrXe/V3z2OepzS3O+o+jtMd5ZHZhFkr4c6Iyif8BkrGGBN3gnDnxYP4t/Pz+Mv+UmZWJ89ebg0u4dGkm4sMAKYBC1T15SbeQ4Gp7kyicMeLAVR1esh3s3G6wBKqKKCV8DDGHK6p5czHx3PwUA3//M6T9MxOjJ3uopXwaHaLQkT6i8gFInK9iMwUkRdF5GMRqQHW46xBWCoiDR6jEJFpIrJQRILt7mIRWSIis0PWTwBHEsQG95opbuJYkmhJwhhjADLS08jr2gGFZm2eFE+NKTN+Pc7AcW7IKxB6SsjnMmAxziruDe57gzf/aewP+UitDWOMSURtMtPp17k9b66s4JE3NnL9OZ6uR/ZcYwaz5+AsmhOcH/orqZ8Iyt09tI0xxkTRrWNrep3SnVkvfMjovC6c3DNquTxfNXbW0zKc8hxLVTXe6yKMMaZF6N6+OwDTrjiNi+9/lVueXMVf/210kwoSxkNj9qMoVdVCd4C6iKPTUjfgJI5Vx7k+KwHKkvvKBrONMXW98uFOfvBoCdedPYCffOsU3+Lwaj+KUnCmvwJzQ26eDRSJyJVANk7iKAuz6G4ZTv0nY4wxrvMGdeP7Z/bjj69v5PxB3Tj7xC5+h1SP59Nj6xQFzMYpoyHAZD92yEsk1qIwxgAUrygGYPrpzqz+A4dq+PbvXufLg9W8eMu5BNrV33Y11mI6PbauOkUBbwTuwUkUAa+fZYwxyWhd5TrWVa478ue2rdK5f9JwKvcd4j+eXkss17c1RcyLjbiJYzrxLwpojDFJY0ivbG4bM4jn1u7gqbJP/A7nGPGsSmX7VRtjTBRTzs3l9AE5/PRv7yXUhkdxSxSqOiNezzLGmGSUnib8auIwBLh1/ioO10TedzuePEkUIpIlIneIyHz3NVNEhntxb2OMaWn6ZfWjX1a/sMd6d2rHzy4bQunmXTz0rw1hz4m3ZpcZF5EROKu0q9yvAu77NBFZiDPb6cvmPscYY1qKu8+6O+rxy0b04qV1O7l/6cecc2JXhvm8jaoXLYoZwEBVzVHVHKATMAa4D6cEeJWIDPPgOcYYkzJ+dukQunVsza3zV7H/0GFfY/EiUZSGlvNwZzktU9XpqpoHTAIWiUjiFjIxxpg4uvvNu7n7zbujnpPdLpP7Jg5jY8U+/vcfH8QnsAi8SBQ50Q66+2VfDBR78CxjjEl6m/dsZvOezcc976yBXZhyTi7zlm9h2QefxSGy8DxpUYjIH6KdoKrlgFWVNcaYRrrtopM4uWcW0xat4fMvv/IlhmYnClVdDKSJSImIXB7l1IrmPssYY1JN64x0HrhyOF9+dZgZi9f4smrbk+mxqjoVp/z4YhGpEJEHRWRccFxCRMbhFAs0xhjTSCd178hdXx/MsnU7eWLFlrg/37MFd+6udHnAw8CVOCuxd7lbot4FznoLr55njDHJanDOYAbnDG7UNd8/sz/nnNiFn/39fTZ8vjdGkYXnefXYIzd29q0Y474uxFlfoTg74S0F/gksS6U9Kqx6rDGmOT7bc5CL73+VvjntWHzjWWSme1dcI67VY4NUdaOqzlHVCe76ioHAjcA7OFNmF2OFAo0xpsG6Z7XhnnFDWbNtN79Z9nHcnhvPWk/BxDExJHEUxev5xhiTKGa8NoMZrzWt/N3YIT2ZUNCb37+8ntJNlR5HFl48q8cew00c7/j1fGOM8ctn+z7js31NXxfx00tOpXendtwyfxVfHqz2MLLwfEsUxhhjmqZD6wx+PWkY26sOcPff3o/58yxRGGNMEirol8O/nZ/H4rJt/GPNpzF9liUKY4xJUj+68ESG9c7mP55ey47dB2P2HEsUxhgTZ8O6DmNY1+YX1c5MT+PXk4Zz6HAtdyxcTW1tbJY7WKIwxpg4u6XgFm4puMWTe+V27cBPvnUKr6//gkff3OTJPeuyRGGMMUnuO6f3oejk7hS/sI51O7xfw2yJwhhj4uzWl2/l1pdv9ex+IkLxFUPJapPJLU+u4mB1jWf3BksUxhgTd1VfVVH1VZWn9+zcoTX3jj+NdTu+5Jf//NDTezd7z2xjjDGN8+jYR2Ny3/MHd+Om8wZyygne1l+1RGGMMS3ItLGNq0rbENb1ZIwxJipLFMYYY6KyRGGMMSYqSxTGGGOiskRhjDEmqphthWrqE5HPgc1NuDQb2O1xOF7cszn3aMq1jbmmoed2Ab5oZBwtTSz+fnklnrEl6v9nzblPY67rp6pdwx5RVXsl+AuYk4j3bM49mnJtY65p6LlAqd//ff1+xeLvVzLGlqj/nzXnPl4937qejDHGRGWJwhhjTFSWKJLDswl6z+bcoynXNuaaWPw7a6kS+d9VPGNL1P/PmnMfT55vg9kmpYlIqaoW+h2HMYnMWhQm1c3xOwBjEp21KIwxxkRlLQpjjDFRWZlxYyIQkWlAlfvHgKrO8jEcY3xjLQpjwggmCVWdo6pzgDIRme13XMb4wcYojAlDRHYBA1S1KuQ7VVXxLypj/GFdT6bFEJEAUAwsUdVFUc6bAgRwupUCQHno+SKSi9PVVBXm2iJVXepl3MYkOksUJum5XUI5QAkwBVgS5dxioCJ0vEFEikUkN+S73AiXV+EkFmNSiiUKk/RUdWrws5sIwnJbClNUtVOd66eLyC4RmROuFRGiEichGZNSbDDbpJKpQKRuo3Kc1kg0liRMSrJEYVJJEU5CCKccGBPyOZxAlGPGtFiWKEwqyQcqIhwrBwoBVLUcqHK7qo5hA9kmFVmiMOaoQMjnmTgtEABEZDxWF8qkKBvMNinBnToLR1daR6Wqs0RkmjuVFmBg6KC5ManEEoUxEVjJDmMc1vVkUkLItNeAj2EYk5QsURjjCNDAbiljUo0lCpNKlgIDIxzLAUrjGIsxScMShUklZUQuz5ELLIxjLMYkDRvMNqlkNrCy7pfujKh8YEG8AzImGViLwrREYUttuAvp5oSpB1UMTD9OnSdjUpbtR2GSnrvJ0EicVkEuzqB0Kc5q63oJoE6Z8YFASbSy5MakOksUxhhjorKuJ2OMMVFZojDGGBOVJQpjjDFRWaIwxhgTlSUKY4wxUVmiMMYYE5UlCmOMMVFZojAmCYjI+JDNl2L1jCmxfoZJTrbgzpgYEZHZOPtw57tflQGVdU4rB4rd8iKR7jMNWKqqZTEJ9Nhnzbad/ExdVhTQmBgJ/sAVEQXKVbWg7jluEtggIlNVtd6e3CKSD3SOR5JwFYvIQlWdEKfnmSRgXU/GxJD7gx6cvTDqcbdbLQdmR+j2KVbV6TEKL1w85UB5SNzGWKIwJsaK3PclUc4JdjsVhX4pIuNxuqvibTYw14fnmgRlicKY2BrjvodtUbiCmynVTQqTgPmeR3QcwfESEYm0yZNJMZYojImtIpzxiapwB0WkCCdRLA0zoJ0fx7GJuuYD4316tkkwNphtTIwcb3zCHZOYjZMkxtQ5VhTpOvd4Mc6eGhtUdZZ7fvAeucDkcMnJfWZw46YNQGdgZphzl7rnzYoUg0kdliiMiZ2w4xNul854nK6l4nCznXB+6G8Id1N3ptRMnJ38NogIOK2W6SHHlwEFda7Lx9kXfKqqLg35fjZwzJRYVS1zk48xliiMiaHgb/gjQ/r7B+IkkLJw02VD5AIlEY51VtWqkHt2dmdPBS3FmeZ6pOvKTRIrcXb8C00SS3DWeoRbO1ElIgHbItZYojAmdoLjE/Wmt4rIEhHZoKoDI1wbwNmqte51oQmk0H2fGeEeuRwdIJ+Lk5zqdiXlRLm+0j1eLw6TWixRGBMDxxufwOkCmi0iRaG/4YfIpf4q7iPrHNw/jiH8QHkwgVS5sUzBWR1er9VwnFZNFU7CMinOZj0ZExvHWz8RcN+bMwU10oB38Id/qfs+ASDCWIgxx2WJwpjYON76ieDxSDWeynG6fcJyu6AChE9EEzm2pVFI0xbuBbBuJ4N1PRkTK1HXT3C0eyhSoqgierdPsMVyTCJyu7wCwOSQrwMcbV00Rk60YoUmdViLwhiPNWB8Ao4mgSPjEO5YQlA50bulCgDCJKJinEHrRSHfRWxNiEjAXZMRLUaT4ixRGOMREcl11x4Ef/DmiEh+hGJ/x/ym7iaJ0MRSgjOVNpIinOJ9R1ZPu+snwg1QT8fpjqobbz5wF2FmPbnHoiU6k0JsPwpjPOL+sA+EObSobheOO8YwG6eLqQSnFRC6viEArAw3fdY9toujs5iOnBOp0qybwCZw7CK+sggzroJJp8oGwA1YojAmYbmL4abXrffktiIWAgNjNYYgIiuBC22xnQHrejImkc3GKfNR1xg4WuXVa8EV35YkTJAlCmMSlDsgHW4DoUJiO34wlWNnTZkUZ4nCmMQ2PXRWkjs+kU+MNjRyWxO5PpY3NwnIEoUxCcz9gV3hzp6aglPYD2CKO4bhtem2X7apywazjUkC7gD20liOG8TjGSY5WaIwxhgTlXU9GWOMicoShTHGmKgsURhjjInKEoUxxpioLFEYY4yJyhKFMcaYqP4/cj1wwm1FtsIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from galpy.potential import MWPotential2014\n",
    "\n",
    "print('Tidal Radius: ',cluster.rtidal(pot=MWPotential2014))\n",
    "print('Limiting Radius: ',cluster.rlimiting(pot=MWPotential2014,plot=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also worth noting that it is possible to calculate the tidal radius iteratively. For example, lets assume a snapshot also contains a large number of stars that are likely not members of the cluster. By setting ``rtiterate>0``, ``clustertools`` will first calculate the tidal radius using the mass of all stars in the snapshot (rt_old). However it will then calculate the tidal radius using only the mass of stars within the initial calculation (rt_new). This process is repeated ``rtiterate`` times or until the ratio of rt_new/rt_old is not less than ``rtconverge``, which by default is 0.9."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tidal Radius:  10.6993452944\n"
     ]
    }
   ],
   "source": [
    "print('Tidal Radius: ',cluster.rtidal(pot=MWPotential2014,rtiterate=10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since not all snapshot will contain the energy of each star, ``clustertools`` can calculate these values as well by simply calling:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "ek, pot=cts.energies(cluster)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If energies were instead called internally, such that ``cluster.ek``,``cluster.pot``,``cluster.etot`` are set, a subset of stars could be defined that are gravitationally bound (E<0):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster.energies()\n",
    "bound_cluster=cts.sub_cluster(cluster,emax=0.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calling ``energies`` internally also set the virial parameter Qvir:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Qvir =  -1.40814083302\n"
     ]
    }
   ],
   "source": [
    "print('Qvir = ',cluster.qvir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, in case one needs to know the distance to each stars nearest neighbour:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.032610081984\n"
     ]
    }
   ],
   "source": [
    "distance_to_neighbour=cts.closest_star(cluster)\n",
    "print(np.amin(distance_to_neighbour))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Edit 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.9.1"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}