{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Macroparasite model from May and Anderson (1978)\n",
    "\n",
    "*Author*: Federica Giardina @federigia\n",
    "\n",
    "*Date*: 2018-10-02"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(deSolve)\n",
    "library(reshape2)\n",
    "library(ggplot2)\n",
    "library(scales)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Anderson and May (1978) Model F eqn (19)-(21)\n",
    "HPL_ode <- function(times,init,parms){\n",
    "  with(as.list(c(parms,init)), {\n",
    "  # ODEs\n",
    "  dH=(a-b)*H - alpha*P\n",
    "  dP=(beta*L*H) - (mu+alpha+b)*P - (alpha*((P^2)/H)*((k+1)/k))\n",
    "  dL=(lambda*P) - (gamma*L) - (beta*L*H)\n",
    "  \n",
    "  list(c(dH,dP,dL))\n",
    "  })\n",
    "}\n",
    "\n",
    "# L= number of free-living (infective) stages (e.g. larvae)\n",
    "# H = number of hosts\n",
    "# P = number of parasites\n",
    "# a = host birth rate\n",
    "# b = host death rate\n",
    "# alpha = parasite impact on host mortality\n",
    "# beta*l = rate at which hosts pick up infecting stages\n",
    "# mu = parasite death rate\n",
    "# k = aggregation parameter\n",
    "# lambda*P = rate of production of new free-living stages (i.e. parasite fertility)\n",
    "# gamma = death rate of free-living stages\n",
    "\n",
    "parms <-c(b=1.05,a=1.4,alpha=0.0003,beta=0.01,mu=0.5,k=0.1,lambda=10.0,gamma=10.0)\n",
    "\n",
    "init <- c(H=100,P=10,L=10)\n",
    "times <- seq(0,30,length.out=2001)\n",
    "HPL_out <- lsoda(init,times,HPL_ode,parms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAG1BMVEUAAAAAujgzMzNNTU1h\nnP/r6+vy8vL4dm3////njUASAAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO2dgXba\nShIF7Ue8yf9/8cbB2Agk1N1zp9Uzqjq7ebFNXDRMRWLA5O0PADTzdvQVAJgBQgIQQEgAAggJ\nQAAhAQggJAABhAQggJAABDSF9PGSva+LSLJMNEzSLDWGUYWyByEZLRMNQ0gdICSjZaJhCKkD\nhGS0TDQMIXWAkIyWiYYhpA4QktEy0TCE1AFCMlomGoaQOkBIRstEwxBSBwjJaJloGELqACEZ\nLRMNQ0gdICSjZaJhCKkDhGS0TDQMIXWAkIyWiYYhpA4QktEy0TCE1AFCMlomGoaQOkBIRstE\nwxBSBwjJaJloGELqACEZLRMNQ0gdICSjZaJhCKkDhGS0TDQMIXWAkIyWiYYhpA4QktEy0TCE\n1AFCMlomGoaQOkBIRstEwxBSB3iDSAABHJGMlomG4YjUAUIyWiYahpA6QEhGy0TDEFIHCMlo\nmWgYQuoAIRktEw1DSB0gJKNlomEIqQOEZLRMNAwhdYCQ1nm/52PwYR4lhKSHkFZ5f2TkYZ4k\nhKSHkFZ4yuj7wNSbEmtP5kmREJIIvWWjo4yUSqw9mSdFQkgi5JbtjvrHVGLtyTwpEkISobbs\ndNS3pRJrT+ZJkRCSiANC6tdSibUn86RICEmE2PLUSmpLJdaezJMiISQRWst6JmktlVh7Mk+K\nhJBE9AnpyZLTUom1J/OkSAhJhNSymcfnMAkneSXWnsyTIiEkEUrLdhpfw/RuqcTak3lSJIQk\nokdIK5Y/j5fpElOJtSfzpEgISUR2SB89Wyqx9mSeFAkhiRBaXuTwOMx2S00xlVh7Mk+KhJBE\nHBPSx84zt8FrUGLtyTwpEkIScVhIH3uvgghcgxJrT+ZJkRCSCJ3lVQDbw0hjKrH2ZJ4UCSGJ\nODqkj/2X59mvQom1J/OkSAhJhDykdcueZi8mW00l1p7MkyIhJBEyy8sVbxpmN6bdnEqsPZkn\nRUJIIgqFdP9dYjmVWHsyT4qEkESoQ9qweDSmmNY3ByNX3AshdYCQbrw+WriHMcb0ICyx9mSe\nFAkhiSga0uJ7moMqsfZknhQJIYkQh7RliWtcNW1fAx2E1AFCutEtpPvvXqIoQuoAId3oG9K9\n4uimCKkDhPTFzpIVDhPPSRQVIXWAkL7IC2npa8avJqQOENIXO8uy1zCynhxpEVIHCOmLg0J6\n0PckbZZvT4qEkERoLHt/m2cNkxDUK7TDSL/bloSQRMwV0s9vj+ypjZVhOt5ihCRCGtK25cBh\nDq2iPIQkYv6Q7jl0zZaEkEScK6QHDl3CNSAkERLL7V7ZtowyzJFr+hAISYQypBeWgYZZ5dDF\n3hNCEkFIbsmG5dAewhCSCEJyS3paCCnA3gi6u+eVRvFNCKmHJ0VCSCKUP9/wwjLOMPsSQtJD\nSJ/sdjTSMPsSQtJDSJ8QUhdPioSQRBCSW0JIerwhXS6Xnw/aRhRBSG4JIelxhnT5+v+VthFF\nEJJbQkh6AiH90DaiCEJySwhJjzukGU/t9ne/BxrGICEkPd6QLjOe2u13NNAwBgkh6eEx0gch\n9fKkSAYK6VrO5XY0IqR+lFh7Mk+KZJyQro+KbhXdHiP995eO1yuV28sfAaLsh3RN5xbTN2cI\n6dc3B1wpGArrqd1TSJ+0HXRF6H6odPnZXys0q/YocTYk86RIxjm1O0tIi8+tZZRQU4m1J/Ok\nSAhJRJeQXnTUtaYSa0/mSZEQkogeIe121KumEmtP5kmREJKIDiEZO+qQU4m1J/OkSAYL6c/l\nuaNZQ/J1JK2pxNqTeVIko4W0/AGKK20jimi2PG7aPQeSl1OJtSfzpEgGCmmTthFFqEK6fbie\nxd9hMnIqsfZknhQJIYkQh7TRw20Yc02xnkqsPZknRUJIIjqF9GhZaBw5+XoqsfZknhQJIYnQ\nhrS1/leG8dRk7anE2pN5UiSEJKJPSE+WDY2vpt2iSqw9mSdFQkgiWi2LTbvN1f56GH9OG0WV\nWHsyT4qEkESIQrp+EAxp8Wdbkiqx9mSeFAkhiVCGtH2wsA8T72lDLYeQOkBI6pA+Ft+pYlaE\n1AFCWglpzRLTiHrShkVIHSCku5BeLNfGYdRBNeVFSB0gpJSQbvQNyhgZIXWAkJ5DWrWIh8kv\nyoZyRkKy0jaiCF1IrxZTx2GO7GYctm49QlLRaHk/PKR7jlyqg0JIIjQhff72655ZtxwwzKHr\ncxgISYQspNs9s26pMMyR67UshCTiRCGtcuQirgAhiVCHtGEZYphHDl3hSRCSCFVItztmwzLE\nMEbJluXIHsIQkghCckuGe0KWkBJos7w/hrRlGWEYq2S4kF5JCEkEIbklhKSHkEx7DWMMY5UQ\nkh5CIqRunhQJIYkQhfR6r2GMYawSQtJDSKaHSGMMY5UQkh5CIqRunhQJIYkgJLeEkPScPKTv\n3e+dh0hDDGOWEJIeQjIdkIYYxiwhJD2EREjdPCkSQhJBSG4JIekhJELq5kmREJIISUh7ew1D\nDGOWEJIeQroP6YVlgGHMEkLSc+6QHne/X1jqD2OXEJIeQiKkbp4UCSGJEIS0v9cwwjB2CSHp\nISTTXsMIw9glhKSHkExndiMMY5cQkh5CIqRunhQJIYkgJLeEkPQQEiF186RICEkEIbklhKSH\nkAipmydFQkgiGiy352P3d78HGMYhISQ9hPQT0ktL+WEcEkLSQ0imM7sBhnFICEkPIRFSN0+K\nhJBEEJJbQkh6CImQunlSJIQkoj0kw6bdAMM4JISkpymk0VmGdPS1gZE58xHJsftdfxiPhCOS\nHkIipG6eFAkhiWgOybLXUH8Yj4SQ9BASIXXzpEgISURzSJZNu/rDeCSEpIeQTA+R6g/jkRCS\nHkIipG6eFAkhiSAkt4SQ9Jw4JM/ud/lhXBJC0kNIhNTNkyIhJBGtIZk27coP45IQkh5CMnVU\nfhiXhJD0EBIhdfOkSAhJBCG5JYSkh5AIqZsnRUJIIgjJLSEkPYRESN08KRJCEkFIbgkh6SEk\nQurmSZEQkoioxfV8bPVhfBJC0kNIpo6qD+OTEJIeQiKkbp4UCSGJICS3hJD0EBIhdfOkSAhJ\nBCG5JYSkh5AIqZsnRUJIItpCMu5+Vx/GJyEkPacNyXdAKj6MU0JIegiJkLp5UiSEJIKQ3BJC\n0kNIhNTNkyIhJBGE5JYQkh5CIqRunhQJIYloCsm6+118GKeEkPQQkqmj4sM4JYSk56whOc/s\nag/jlRCSHkIipG6eFAkhiSAkt4SQ9BASIXXzpEgISQQhuSWEpOfcIZl3v2sP45UQkh5CMnVU\nexivhJD0EBIhdfOkSAhJBCG5JYSk56QhefcaSg/jlhCSHkIipG6eFAkhiSAkt4SQ9BASIXXz\npEgISQQhuSWEpIeQCKmbJ0VCSCIaQrK/sKH0MG4JIekhJFNHpYdxSwhJDyERUjdPimTMkC6f\nfH/UNqIIQnJLCEmPN6TFR20jiiAkt4SQ9JwzJPemXeVh/BJC0hM4tfv5qG1EEYTklhCSnsBm\nwwSPkQgpxZMiISQR8ZAcTyNVHsYvISQ9lpCu5fw7qbv8mS8km6XuMH4JIekxhHR9VHSr6PYY\n6b+/dLxefVmEdPSVgRnYD+mazi2mbwgJ4Afrqd1TSJ+0HXRFcGrnlnBqp+eUIfk37QoPE5AQ\nkh5CIqRunhQJIYkIh+TZ/S48TEBCSHoIiZC6eVIkg4V02/5e0DaiiNaQjJaywwQkhKTH+YTs\nkrYRRRCSW0JIek75g32ElEKJYVSh7EFIRkvZYQISQtJzxpACu991h4lICEkPIRktVYeJSAhJ\nDyEZLVWHiUgISQ8hGS1Vh4lICEnPeUNyPR9bd5iIhJD0EJLRUnWYiISQ9Jw+JKul6jARCSHp\nISSjpeowEQkh6TlhSJG9hrLDhCSEpIeQjJaiw4QkhKSHkIyWosOEJISkh5CMlqLDhCSEpOe0\nIfl2v8sOE5IQkh5CMlqKDhOSEJKes4dkthQdJiQhJD2EZLQUHSYkISQ9hGS0FB0mJCEkPYRk\ntBQdJiQhJD2EZLQUHSYkISQ95wsp9DRS1WFiEkLSc9aQnLvfVYeJSQhJDyEZLTWHiUkISc/J\nQ7Jbag4TkxCSHkIyWmoOE5MQkh5CMlpqDhOTEJIeQjJaag4TkxCSHkIyWmoOE5MQkp7ThRR7\nGqnoMEEJIek5aUje3e+iwwQlhKSHkIyWksMEJYSk59whOSwlhwlKCEkPIRktJYcJSghJDyEZ\nLSWHCUoISQ8hGS0lhwlKCEkPIRktJYcJSghJDyEZLSWHCUoISc/ZQgo+H1tzmKiEkPQQktFS\ncZiohJD0nDMk9/OxNYeJSghJDyEZLRWHiUoISc+pQ/JYKg4TlRCSHkIyWioOE5UQkh5CMloq\nDhOVEJIeQjJaKg4TlRCSnpOFFN39LjlMWEJIegjJaCk4TFhCSHpOGZJ/97vkMGEJIekhJKOl\n4DBhCSHpOXNILkvBYcISQtJDSEZLwWHCEkLSQ0hGS8FhwhJC0kNIRkvBYcISQtLTFNJ4LEI6\n+srARHBEMloKDhOWcETSc66QwrvfFYeJSwhJDyEZLfWGiUsISc+JQ/JZ6g0TlxCSHkIyWuoN\nE5cQkh5CMlrqDROXEJIeQjJa6g0TlxCSHkIyWuoNE5cQkh5CMlrqDROXEJKeU4UU3/0uOEyD\nhJD0EJLRUm6YBgkh6TlvSE5LuWEaJISkh5CMlnLDNEgISQ8hGS3lhmmQEJIeQjJayg3TICEk\nPYRktJQbpkFCSHoIyWgpN0yDhJD0nCmk8LtDfhQcpkVCSHoIyWipNkyLhJD0nC+k0POx9YZp\nkRCSHkIyWqoN0yIhJD2nDclrqTZMi4SQ9BCS0VJtmBYJIekhJKOl2jAtEkLSQ0hGS7VhWiSE\npOdEIbXsfpcbpklCSHoIyWgpNkyThJD0nC6k2O53uWGaJISkh5CMlmLDNEkISc9ZQ3Jbig3T\nJCEkPYRktBQbpklCSHoIyWgpNkyThJD0EJLRUmyYJgkh6SEko6XYME0SQtJDSEZLsWGaJISk\n5zwhNe1+VxumTUJIek4akt9Sa5g2CSHpISSjpdYwbRJC0kNIRkutYdokhKSHkIyWWsO0SQhJ\nDyEZLbWGaZMQkh5CMlpqDdMmISQ9JwspuvtdbZg2CSHpOU1IbU8jFRumUUJIes4ZUsBSaphG\nCSHpISSjpdQwjRJC0kNIRkupYRolhKSHkIyWUsM0SghJDyEZLaWGaZQQkh5CMlpKDdMoISQ9\nhGS0lBqmUUJIes4SUtO7Q34UG6ZVQkh6/CFdfn7bNqIIT0jh52NrDdMqISQ97pAuhNSREmtP\n5kmRDBrSZYYjUsRSaZhWCSHpISSjpdIwrRJC0uMM6TLFY6SIpdIwrRJC0uMN6ZPvj9pGFEFI\nbgkh6bGEdC3nVtCQR6TW3e9SwzRLCEmPIaRrQN8ndV8h/feXbtdKzjKko68NzMd+SNcj0S2m\nb0YM6RchgYW3t/UP3rZzsZ7aPYX0SdtBV4Tj1C5+ZldqmGYJp3Y7ENIWhHQvISQPhHQHId1L\nCMkDId1BSPcSQrrv4t+vb59cP7j+7+FTdxf9/oiQIpZCwzRLCOkxpLdrLfch3X3q7a6p28fR\nkB5e00BIfSix9mSeFIkopLuiNj/18/nnkzz/E7LjhhR/7XetYZolhHRX0ncSb8+Z3H3q+7IN\nIW3SNqIIi6X1hyhKDdMuIaSnkG7na48Pm9ZCWj23O2FIMUudYdolhPQY0tvz8ebhUw9HpOdv\nRUg2S51h2iWEdLf87w9BGw+bCOkGIS0khPS9/N/eDCE9bDZ8f0hIAUudYdolhHRX0i2ox82G\nne1vQiIkQnoK6Subx82G26dET8gSUqlh2iWEpOdEIbXsflcapl1CSHpOEVL7AanQMAIJIekh\nJKOlzDACCSHpISSjpcwwAgkh6SEko6XMMAIJIekhJKOlzDACCSHpISSjpcwwAgkh6SEko6XM\nMAIJIelZhvTiZ9LXaBtRxL6l+d0hPwoNo5AQkp7zhNT0fGydYRQSQtLT9Nq7thFFEJJbQkh6\nHo9IGz//t07biCKcIUUtVYZRSAjpyv+2IaQVCOlBQki9Q3LSNqIIQnJLCOkKIS00u5cgpAcJ\nIV3pGdL9+6ZMEpJi97vMMBIJIXUP6Y2QNixFhpFICKl3SG9f76oy1WaDYve7zDASCSER0opm\n7wKKA1KZYSQSQiKkFc3eBQjpUUJInpDeVkIxhjTX80iE9CghpN4hbb2xMSFVGUYiIaTuIW38\n4y+EVGUYiYSQbiH9XkcQkou2EUUQkltCSIS0otm7gGL3u8wwEgkhuUIynaU9fd1xZjdGSJID\nUpVhNBJCcoW0HsrrkCbcbCCkJwkh9Q7p+73CJ9r+JqQnCSH1D2n5X0L6ttQYRiMhJEJa0ex8\nnZCeJIR0C6nnKxsIadVSYxiNhJA8Idl42mxwdDRQSK2731WG0UgIqWdIbwumCUlzQCoyjEhC\nSIS0onn9ZUJ6lhBSz5ACtI0ogpDcEkK6QkgLzesvE9KzhJCudN1s4NRu3VJiGJGEkHqHxGOk\nTUuJYUQSQuofkusPt40ogpDcEkLS83hq5/rDbSOKeG2RvBfXR5FhVBJC0jP7v0YheS+ujyLD\nqCSEdIXNhoXm5VcJaUVCSL1Dmm+zQXRmV2MYlYSQ+ofk+sNtI4ogJLeEkHqHNOtmAyEtJITU\nPaS5NxtaLBWGUUkIqXdIPEbatFQYRiUhpGIhDcBi0+7oKwO16Hlq56Lt7woRLy2qA1KJYWQS\njkhXbCHZDiyEZLQUGEYmISRXSH8s29mTP0YipDUJId1Cel/n+c1PCImQniWE1DukV58kpARK\nrD2ZJ0VSOSRrSW0jiiAkt4SQXCE1bDbMdWrX/pLVGsPIJITkCsnUwtwhyQ5IFYbRSQipd0iz\nbTYQ0qqEkNJCMv1ZQpJTYu3JPCmSlOeRvCE5aRtRBCG5JYRESCuaF18jpFUJIXlCsjH3WxYT\n0qqEkAhpRbP9pUVHhPQtIaSeId19bpbNBt0BqcAwQgkhpYRk/2mkthFFEJJbQkgJIXl+qK9t\nRBGE5JYQUveQfD8c2zaiCEJySwipd0jOnzFvG1EEIbklhNQ3JPdbNbSNKIKQ3BJC0nOXjv8t\nT9pGFLFteSekdQkh6Zn5eSRC2pAQ0pU+p3azhiR5Pvb4YZQSQuoZUoC2EUUYQ2q1HD2MUkJI\nVwhpodn8SiCk3/9Ysxw9jFJCSFcIaaHZ/IorpN/PLCxHD6OUENIVQlpoNr/iCGklo2VMhw+j\nlBDSFUJaaDa/Yg5pM6OfmA4fRikhpCuEtNBsfsW4abeT0ZXDh1FKCOkKIS00W18wHpBMHT09\nZupEibUn86RIuod0+h81N4VkzyilphJrT+ZJkRCSiKaQVkI5NKYSa0/mSZE0hfRrHUK6xxDS\niz6OqanE2pN5UiSEJKIhpL0w8o9NJdaezJMiISQR4ZCMRaTWVGLtyTwpEkISsWXZe8mqJ4X9\nmEQ5lVh7Mk+KhJBEmEJ6/rIzgT9/LDE151Ri7ck8KRJCEhELybv2r8PYYmroqcTak3lSJIQk\nIhKSf83fDWOuKdBTibUn86RIeB5JRCCkwFJ/GMYRkyuoEmtP5kmR9A+p4V/smzuk0Apf0Thr\nsvhKrD2ZJ0XCa+1EuEMKnXptDxPIadtcYu3JPCkSQhKxYdl855NQR3vDxGp6vhol1p7MkyIh\nJBEvQ3p+GinWkW2Yppwi1yoGIX1BSAvN+qc3QoouWM8wkp46hkVIXxDSQrP+6fUzu/D6DAyj\n7UmYFyF9QUgLzfqn10JqWIgNw3QNKlIZIX1BSAvN+qdXQmroSDJMflF22qd7hpCstI0oYt2y\nsmnXtGy0wxzWy4Q0hKTkPCHd3fYRS69hDluAs0BIIl6FdLdpd3fThyw5wxy2HMeFkEQYQvpo\n2K67WY4Z5rj1OQyDhnS5XH4+2BuhaU1ZMYX0c7tHLQf/uO+C41ZtQcYM6fL1/yt7IxhXTxuW\nkH5u9rClUkivOW5NH8KYIX0yXkg/t3rcMk5I+5Iny3EZtDNqSCOe2n3f6A2WqUPq5EmRjBrS\nGEek5aZde0eEFPGkSAhJxH5Igo4IKeJJkQwU0rWcf/t1o2w2rIXUZiEkvydFMk5I10dFt4pu\nj5H++0vH69XISkhHXyWYmv2QruncYvpmsJCOvkIwOdZTu6eQPmk76Iowndq1Wzi183tSJOOc\n2g0Y0uOmncBCSH5PioSQROyGJPlhG0IKeFIkhCRiPySJhZD8nhQJIYnYDUljISS/J0UyWEi3\n7e8FbSOKICS3hJD0OJ+QXdI2oghCcksISc+UP9i3+S6rDRZC8ntSJIQk4tnym5BeSwhJz4Qh\n/f699XbFLRZC8ntSJIQk4tHy+yEkkYWQ/J4UCSGJeLD8JqRdCSHpmS2k34S0LyEkPZOFdNcR\nIW1KCEnPXCH9JiSLhJD0TBXS79/PZ3aE9CwhJD0ThfT9XkEdDkiEFPGkSAhJxM3ym5CsEkLS\nM01IvwnJLCEkPbOEdNdRh5esElLIkyIhJBH/LHcZ9XjtNyGFPCkSQhLxadnuiJBWJISkZ4qQ\nFh11OSARUsSTIiEkEX8e/9EWQtqTEJKeCUJ66IiQdiWEpGf8kB47IqRdCSHpGT6kp44IaVdC\nSHoGD+n3TkeEtCYhJD1jh7TSUZ8DEiFFPCkSQmpnrSNC2pcQkp6RQ1rtiJD2JYSkZ+CQ1jsi\npH0JIekZN6T1jAjJICEkPcOG9LojQnolISQ9o4b009GDhZD2JYSkZ9CQ7o5HqyGpn0YipIgn\nRUJIDdyf170KSWckpIAnRUJIYZbbdYTklhCSngFDetj2Xlo6PUQipIgnRUJIQR6fPiIkt4SQ\n9AwX0tPTsITklhCSntFCen45w1pI8k07Qop4UiSEFGHlZUEvQhKKCSngSZEQUoC1lzMQkltC\nSHqGCmn1ZUGE5JYQkp6RQlp/eR0huSWEpGegkDZeprqw9NprIKSIJ0VCSE62Xu69HZLSTkgB\nT4qEkFysbNfdNPcfEJJFQkh6BglpuyNC8ksISc8YIb3oiJD8EkLSM0RIrzpa3F29XiBESCFP\nioSQzLzsaCWkDpt2hBTxpEgIycrLjF6EpLHfLITk96RICMnITkeE5JcQkp7yIe11REh+CSHp\nqR7SbkeE5JcQkp7iIe13dH939du0I6SIJ0VCSAYMHRGSX0JIeiqH9Hrb+1vz89t+u9+EFPGk\nSAhpD1tHmyG1yZ8shOT3pEhmCKkrdx1Z/8gipJ7XDeCBskck4/Ho/u+9jg+ROCJFPCmSGY5I\nbSO+xNwRIfklhKSnaEjmjJ5D6rLXQEgRT4qEkF7g6IiQ/BJC0lMyJE9HWyFF3VsWQvJ7UiSE\ntImrI0LySwhJT8GQfB0Rkl9CSHrKhWTfrrtpbr/puWlHSBFPioSQVnF39BhSn70GQop4UiSE\ntIa/o42Q/OodCyH5PSkSQloh0BEh+SWEpKdUSJGOCMkvISQ9lUIKdfR9d3XdayCkiCdFQkiP\nhDJ6DKnTXgMhRTwpEkJ6INjReki+b2GxEJLfkyIhpCXRjgjJLyEkPVVCCndESH4JIekpElK8\no9vd1XevgZAinhQJIf0Q2667aa7/ISS7hJD0VAipqSNC8ksISU+BkNo6WobUa/ebkCKeFAkh\nfdHY0WpIoW/02kJIfk+KhJCutHb0dXd1PrMjpIgnRUJI/2jNiJACEkLSc3BI7R0Rkl9CSHqO\nDUnQ0SKkbnsNhBTxpEgISdPRWkgN323TQkh+T4qEkDQdEZJfQkh6DgxJ0xEh+SWEpOewkJq3\nvb81n7/03msgpIgnRXLykGQd3YfUb6+BkCKeFMm5Q9J1tBJS63dctRCS35MiOXVIwo4IyS8h\nJD2HhKTs6N/d1f0hEiFFPCmSE4ekzIiQAhJC0nNASNqO7kLquNdASBFPiuS0IYk7eg5J820f\nLYTk96RIzhqSuiNC8ksISU92SPKOPu+u/g+RCCniSZGcMiTpdt1NQ0heCSHpSQ2pR0eE5JcQ\nkp7MkLp09BNSz007Qop4UiTnC6lPR4TklxCSnryQOnX0FJL0m/9YCMnvSZGcLaROGX3eXQkP\nkQgp4kmRnCykbh0Rkl9CSHqSQurX0XdIXR8iEVLEkyI5VUgdO3oMSS+4WgjJ70mRnCmknh0R\nkl9CSHoyQura0cefjIdIhBTxpEhOE1Kvbe9vDSF5JYSkp3tIvTu6hdR3r4GQIp4UyaAhXS6X\nnw8MI3bv6CGkThJCinhSJGOGdPn6/5X9Eft3lPIsEiGFPCmSMUP6xBFSQkeE5JcQkp6uISVk\nlPKK1Q9CCnlSJKOG5HiMlNIRIfklhKTHvdlw/8HrEXI6Snk6lpBCnhTJQCFd4/m3X7fo6PUI\nqR0RkkdCSHoMIV1P5q4VXS63/e///vLyj3131H4lX7EMqa8LYIv9kK7l3GL6ZjekW0rx62aC\nkKAC1lO7p5A+2Tmo9j+v+8jaa+DULuJJkYxzatcUku4W24CQ/BJC0tM5JN0NtkXOXgMhRTwp\nkhOElHFDJm3aEVLEkyIhJAmEFJAQkh7z80hPTyL9KRRS94dIhBTxpEhGC2n5AxRX2kaUkHRA\nIqSIJ0UyUEibtI2oIOvMjpAinhQJISkgpIiEkPRMEVL/h0iEFPGkSAhJASFFJISkZ+yQ0s7s\nCCniSZEQkgBCCkkISQ8h2SCkgCdFQkgCCCkkISQ9Q4f0nrbXQEgRT4qEkNrJOyARUsSTIiGk\ndggpJiEkPYRkg5ACnhQJIbWT9xCJkCKeFAkhtUNIMQkh6Rk5pMQzO0KKeFIkhNQMIQUlhKRn\n/JBSzuwIKeJJkRBSM4QUlBCSnoFDyjyzI6SIJ0VCSK0QUlRCSHqGDynnzI6QIp4UCSE1kviK\n1Q9CCnlSJITUSOqZHSFFPCkSQmqEkMISQtJDSDYIKeBJkRBSG7kPkQgp4kmREFIbhBSXEJKe\nSULqafoHIQU8KRJCaoOQ4hJC0jN4SFlndoQU8aRICOXfyzYAAAkoSURBVKmJ5IdIhBTxpEgI\nqYnkMztCinhSJITUBCE1SAhJz9ghpZ3ZEVLEkyIhpBayHyIRUsSTIiGkFrLP7Agp4kmREFIL\nhNQiISQ9Q4eUd2ZHSBFPioSQGkh/iERIEU+KhJAaSD+zI6SIJ0VCSA0QUpOEkPQMGdI7ITVJ\nCEnPyCElPkQipIgnRUJIcQipTUJIeiYIqZdlASEFPCkSQoqTf0AipIgnRUJIYQ44syOkiCdF\nQkhhliElrXBC8ntSJIQUhpAaJYSkZ8CQHs7sCMktISQ9hGSDkAKeFAkhRXnY/CYkt4SQ9IwX\n0uOeHSG5JYSkh5BsEFLAkyKZIaRDeNizA6jA6EekpL/3OCJFPCmSGY5IbSPGePoRCkJySwhJ\nz6gh/SKkuISQ9BCSDUIKeFIkhBTj6UcoCMktISQ9o4X0/MpvQnJLCEkPIdkgpIAnRUJIIZ5/\nOJaQ3BJC0jNYSCvvH0RIbgkh6RkzpF+E1CIhJD1Dh3TV6C0rEFLAkyIhpABr7wxJSG4JIekZ\nMqRfhNQkISQ9hGSDkAKeFAkh+Vk7syMkv4SQ9BCSDUIKeFIkhORn7cyOkPwSQtJDSDYIKeBJ\nkRCSm9UzO0LySwhJz4Ah/SKkRgkh6SEkG4QU8KRICMnL+pkdIfklhKRnvJB+EVKrhJD0jBPS\nOyGpJISkZ9iQfjRSyxaEFPCkSAjJycYBiZD8EkLSQ0g2CCngSZEQko+tMztC8ksISQ8h2SCk\ngCdFQkgutvbsCCkgISQ9hGSDkAKeFAkhudjsiJD8EkLSM0hI2wckQvJLCEnPmCEtNDrLCwgp\n4EmREJKDFwckQvJLCEkPIdkgpIAnRUJIDl6c2RGSX0JIeoYKae2AREh+CSHpGSKkV2d2hOSX\nEJKeEUN60KgsLyGkgCdFQkhmXh6QCMkvISQ9hGSDkAKeFAkhWbl19LHaESH5JYSkZ6CQ1g9I\nhOSXEJKe+iE9HpAIqVlCSHrGCWmjI0LySwhJDyHZIKSAJ0VCSDZ2thoIKSAhJD3DhLR1QCIk\nv4SQ9FQP6Z2Q9BJC0jNKSJsdEZJfQkh6ioe0f0AiJL+EkPS4Q7rc/b5tRAt7TyJ9EFJAQkh6\nvCFdUkMyHJAIyS8hJD21j0iGAxIh+SWEpGeMkF50REh+CSHpKR2S5YBESH4JIekZKqR1TavF\nBCEFPCmSgUK6tnO57jM4QvpMoAXTAYmQ/BJC0mMI6TugRUj//WXnD/59bPPecNW+t+z+3EJq\n+GYAXdkP6Xokejoa7YX0tfYbSqIjGAfrqd1TSJ+8OqTeVn/47O7ddmbHqZ1fwqmdnvoh7XRE\nSH4JIenpFlJrSU8dEZJMQkh6+ocULMl6QCIkv4SQ9PQLqe2QZH2EREgBCSHpMT+PdHnuyBpS\noCR7R4TklxCSHvcTsvfsjNBQ0ndHhNRBQkh6ev5g369wSY4DEiH5JYSkp+tPyDaHtLtl90FI\nAQkh6UkJyVuS48SOkAISQtLT9z0bYiU9n9gRklRCSHo6v/nJT0j2kt5dByRC8ksISU/vdxHy\nl7TSESFpJYSkJyske0nOAxIh+SWEpKf7+9p5S/IekAjJLyEkPXkh2Upyd0RIfgkh6en/Tqv3\nIe2W9O49sSOkgISQ9CS8ZbGjpLWOCEkuISQ9Ge/9bS4p0hEh+SWEpCc1pF/vL1O678h6YkdI\nAQkh6Un51ygeSlpP6X21I0LqICEkPTn/rMtDSGslRTsiJL+EkPQk/ftIjyU9pvQe7oiQ/BJC\n0pP1D409lXTX0vv7Q0eE1FdCSHrS/sW+lZJWWVzWdEtaLtQMIQU8KZKzhfS04yDqiJD8EkLS\nk/hvyBpKWl7OeEvaLtYIIQU8KZLzhbRf0vJC1lvSeLk2CCngSZGcOqTVkh4uY74lrRdsgpAC\nnhTJCUO6L+kppcdL2G9J8yVbIKSAJ0VyxpAWJf16lREh9ZQQkp7ckB5SMn/ppcZ+0QYIKeBJ\nkZw0pGUud8FsfNpwSzouG4eQAp4UyVlDeixpDd8t6bp0FEIKeFIkpw1pvyTnLem7eBBCCnhS\nJOcNaS8l7y3pvHwMQgp4UiRnDulVSf5b0v0nIhBSwJMiOXVI2ym5b0hC8ksISc9RIa2n5LwR\nr5rIH/JbCMnvSZGcPqSPhi3vhSb455wWQvJ7UiSEJIKQ3BJC0kNIRstEwxBSBwjJaJloGELq\nACEZLRMNQ0gdICSjZaJhCKkDhGS0TDQMIXWAkIyWiYYhpA4QktEy0TCE1AFCMlomGoaQOkBI\nRstEwxBSBwjJaJloGELqACEZLRMNQ0gdICSjZaJhCKkDhGS0TDQMIXWAkIyWiYYhpA4QktEy\n0TCE1AFCMlomGoaQOkBIRstEwxBSBwjJaJloGELqACEZLRMNQ0gdICSjZaJhCKkDTSEBwBWO\nSEbLRMNwROoAIRktEw1DSB0gJKNlomEIqQOEZLRMNAwhdYCQjJaJhiGkDhCS0TLRMITUAUIy\nWiYahpA6QEhGy0TDEFIHCMlomWgYQuoAIRktEw1DSB3o+BKh//7r973TmWmYmWYpMwwh2Zhp\nmJlmKTMMIdmYaZiZZikzDCHZmGmYmWYpMww/RgEggJAABBASgABCAhBASAACOoZ0uVz6ffNU\nrnNMMc/XEFPMUmqYfiFdbgtweL7urhnmuVynmWKWWsN0C+ny/cvoXP/Cm2Oe6yRzzFJrGEIy\nME1Il1JrT0GZYQjJwDQhfXIps/bauVzK3DGEZGCmkAqtPQF1/lYgJAMThXRZ/jI8hDQS84R0\n+fl1/FlKPeBj+9vANNvfl5//Dj/LWba/izxPJmGaJ2Qv//gzxSx/zvKELMCJICQAAYQEIICQ\nAAQQEoAAQgIQQEgAAggJQAAhAQggpON4++bz90dfG2iC++84CGkiuP+OhYAmgfvxWAhpErgf\nj+UW0u307t/HX5+8nvTBEHBPHctDSP/a+Qro9vAJRoA76liWIf38cvsKJQ0C99OxPJ7a/fxC\nSEPB/XQsr0Li3G4guJ+OZe+IBIPAvXUshDQJ3FvH8vLU7g85DQP307Fsh/SHx0gjwf10LC9C\n4gnZkeCeAhBASAACCAlAACEBCCAkAAGEBCCAkAAEEBKAAEICEEBIAAIICUAAIQEIICQAAYQE\nIICQAAT8HzC8cFfWexVlAAAAAElFTkSuQmCC",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "HPL_out_long <- melt(as.data.frame(HPL_out),\"time\")\n",
    "\n",
    "#Visualisation\n",
    "library(ggplot2)\n",
    "\n",
    "ggplot(HPL_out_long,aes(x=time,y=value,colour=variable,group=variable))+\n",
    "  # Add line\n",
    "  geom_line(lwd=2)+\n",
    "  #Add labels\n",
    "  xlab(\"Time\")+ylab(\"Number\")+scale_y_log10(breaks = trans_breaks(\"log10\", function(x) 10^x),\n",
    "              labels = trans_format(\"log10\", math_format(10^.x)))\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}