{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Based on Hartfield & Otto 2011 which cite Haldane 1924:\n", "\n", "> Once established, the frequency of A1B1 can be modeled by the standard deterministic equation for haploid selection (Haldane, 1924): $p(t) = \\frac{p_0}{p_0 + (1 - p_0)(1-s)^t }$" ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [], "source": [ "p = lambda t: (p0)/(p0 + (1 - p0)*(1-s)**t)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Simulation:" ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [], "source": [ "N = 1e6\n", "p0 = 1./N\n", "p1 = 1-p1\n", "s = 0.01" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "100 2.73199429462e-06\n", "200 7.46377043595e-06\n", "300 2.03907499683e-05\n", "400 5.57055417224e-05\n", "500 0.000152172803824\n", "600 0.000415626408093\n", "700 0.0011346741306\n", "800 0.00309384842062\n", "900 0.00840733978591\n", "1000 0.0226391835409\n", "1100 0.0595165224224\n", "1200 0.147404271787\n", "1300 0.320805405939\n", "1400 0.563397103729\n", "1500 0.779025113177\n", "1600 0.905938943823\n", "1700 0.963387337251\n", "1800 0.986280125989\n" ] } ], "source": [ "f = [array([p0, p1])]\n", "tick = 0\n", "S = array([1,1-s])\n", "while f[-1][0] < 0.99:\n", " tick += 1\n", " f.append(f[-1] * S)\n", " f[-1] /= f[-1].sum()\n", " if tick%100 == 0: print tick,f[-1][0]\n", "f = array(f)" ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [], "source": [ "t = arange(0,tick+1)\n", "y = np.vectorize(p)(t)" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 56, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEECAYAAADAoTRlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcVHXiPvDnMMMwwzCAM4gEWBpa2fTNG1qbdxC0i6WW\nbFn6M2stNcvK2tTMyrW+rdr2bdN0N8xKayXz1kXzWpq0KQVe0BRKSyVAhtsIDMzMOb8/sMlxRhAZ\n5szleb9evpoDH2cePuHD4cw5nyNIkiSBiIiCSojcAYiIyPtY/kREQYjlT0QUhFj+RERBiOVPRBSE\nWP5EREFI2dyAJUuWIDc3F5GRkVi0aJHbMcuXL0deXh7CwsIwZcoUdO7c2eNBiYjIc5rd8x8yZAhm\nzZp10c//8MMPKCkpwZtvvolJkybhnXfeuaQXzs/Pv/SU1CzOp+dwLj2L8+lZnprPZsu/W7du0Gq1\nF/18Tk4OBg0aBADo2rUrampqUFlZ2ewL8xvCszifnsO59CzOp2d5rfybU15eDoPB4Ng2GAwoLy9v\n7dMSEVEb8sgbvlwhgojIvzT7hm9z9Ho9TCaTY9tkMkGv17uMy8/Pd/p1ZejQoa19aTpPRkaG3BEC\nBufSszifnpWamoqsrCzHttFohNFobPHztLr8k5OT8eWXX6Jfv344duwYtFotoqOjXca5C1hUVNTa\nl6dzdDodzGaz3DECAufSs/x5PoWaGkS89RZCDx+G8sgRKE+fBgA0hEXgmUmncfxEKF55pQp6vei1\nTPHx8R75gdps+b/xxhs4cuQIqqurMXnyZIwZMwZ2ux0AkJaWhl69eiE3NxfTpk2DWq3G5MmTWx2K\niEhuFRUCvv+vDmMXL4XS3uD0OVX9WeitJbhmWAzCwvzzsLcg55LO3PP3HH/eu/I1nEvP8tn5lCSE\nHjwI9ebNCPvmGxQuy8LXe6OxZ08Y9u1ToahIgV69GrDyyM3oeGZ/419RqWDr0gW2Ll1Q/fzzsCck\neD12fHy8R56n1Yd9iIj8ieKXX6BZuxaadesQ+tNPjo/P6X8Y1f2HYsCAeowbV4tu3axQKgH1xkdR\nYbPBev31sCUlAaGhMqb3HJY/EQUV3aSpCD+U6/Lx5fesQc1rvVw+brnzTm/E8jqWPxEFhQMHQpGZ\nqYX+pyewDBMAAKJWC8vw4bAMG4b6gQPlDehlLH8iCjx2OxTHj8OW1AXffKPC66/rcPq0AhMm1OK+\nmYNQO28kLOnpqE9Ph6TRyJ1WFix/IgoooTk5iJo9G+LpUgzv8iN+LdPhqafMuPPOOiiVAKBC5eLF\ncseUHcufiAKCUFGBqHnzEL56teNjL13xFjqveehc6dP5OCVE5PfCdu1C9BPToSgtcXxMDFOjR087\nathybnFaiMjvVR0tg+G84q8bPhzVL74Ie8eOMqbybbyTFxH5tc8+U6Pvm4/i4HV3wR4TA9O776Ii\nM5PF3wzu+RORX7Lbgb//XYf16zVYuaoCMR3n44zdDjEmRu5ofoHlT0R+RaithVnUYvLkdqirE/DF\nF2UwGERIaAf/XGVHHjzsQ0R+Q7N6NWJu6Y+n7yzHFVfY8dFHJhgM3ltRM5Bwz5+I/EL4e+8h+tz9\nxFfU346GJ9dCCo2TOZX/4p4/Efk8zerVjuIHgMhELQQF66s1uOdPRD5NvWEDop9+2rHd0LMnTB9+\nCCkyUsZU/o8/OonIp4n7DkA4d9sRq9EI08qVLH4PYPkTkc9qaADuOvoGVv1pEazdusH0n/9AcnOb\nWGo5HvYhIp8kScCsWVEID5cwcPm9OCPeHTA3UvEFLH8i8knLlmlx4IAK69aVQaEAoGDxexLLn4h8\niyhib44ab78dgS++KINWy0u32gKP+RORz1D9979ol5KO/5t0GosWVSIhwS53pIDFPX8i8glCZSWi\np06FsrgYm5W34KywFPVIlTtWwOKePxH5hKi5c6EsLgYAKHVhsF53ncyJAhvLn4hkF7ZlC8LXrHFs\nVy1YADEhQcZEgY/lT0SyEqqqEP3ss47t2lGjYLn1VhkTBQeWPxHJSoqMxOa+z6E2RAt7bCyq5s2T\nO1JQ4Bu+RCSr34oVGJf9BL78qD86h56G1K6d3JGCAsufiGT18stRGD++Fgn9E9GARLnjBA0e9iEi\n2ezerUJubiimTTsrd5Sgw/InIlk0NADPPx+Fl1+ugkbDq3i9jYd9iMjrtG+/jf8WXI2E+HuRllYv\nd5ygxD1/IvIqxa+/IvK11zB69URsKLkFIRXlckcKSix/IvIq3WuvQbBaAQAanZJn98iE5U9EXhN6\n4ADC1693bFc9/zwgCDImCl4sfyLyGt2CBY7HdcOHw9qnj4xpghvLn4i8QqiqgnDkGABAEgSY//pX\nmRMFt2bP9snLy8OKFSsgiiJSUlIwcuRIp89XV1fjn//8JyorKyGKIkaMGIHBgwe3VV4i8lNSVBSm\npe/HgBMfYVS3Q7Bdc43ckYJak+UviiIyMzMxZ84c6PV6zJw5E8nJyUhM/OMqvM2bN6Nz584YO3Ys\nqqurMX36dAwYMAAKhaLNwxOR/ygpCcHHG6Ix9au7Ud1+lNxxgl6Th30KCwsRFxeH2NhYKJVK9OvX\nDzk5OU5j2rVrh9raWgBAXV0ddDodi5+IXLz9dgTuvrsW7duLckchNFP+5eXlMBgMjm29Xo/ycudz\nclNTU3Hq1Ck88sgjeOaZZzBhwoQ2CUpE/qusLAQffxyOKVO4jIOvaPUVvuvWrUOnTp3w4osvori4\nGH/729+wYMECaDQap3H5+fnIz893bGdkZECn07X25ekclUrF+fQQzqVnqVQqrFqlwujRNnTtqpU7\nTkDIyspyPDYajTAajS1+jibLX6/Xw2QyObZNJhP0er3TmGPHjmHUqMbjd78fIioqKkJSUpLTOHcB\nzWZziwOTezqdjvPpIZxLzwkpK4N67AOwn3gUE9fdCrM5TO5Ifk+n0yEjI6PVz9PkYZ+kpCQUFxej\ntLQUNpsN2dnZSE5OdhoTHx+PgwcPAgAqKytRVFSEDh06tDoYEfm/8PffR2j+QSysmYre8x+UOw6d\np8k9f4VCgYkTJ2L+/PmOUz0TExOxdetWAEBaWhpGjRqFJUuW4JlnnoEoinjggQcQERHhlfBE5MMs\nFmhXrHBs1v75z/JlIReCJEmyraVaVFQk10sHHB6q8BzOpWeEf/QRomfMAADY4uNRmp0NhIbKnMr/\nxcfHe+R5eIUvEXmeJEH77387NmseeojF72NY/kTkcSGlpbDWNK7cKWkjUDt2rMyJ6EK8mQsReZzY\noQOmpeWhV9XXeGjICUiRkXJHoguw/InI42prBaxdp8WjW/rAet0ggO+h+Bwe9iEij9uwQYM+fRqQ\nkGCXOwpdBMufiDzugw/CMW5cjdwxqAksfyLyqP37Q2EyhWDwYN6Y3ZfxmD8ReUzU7Nk4+0NHPHrX\nA1AoeLGnL2P5E5FHhJSVIXzlSoyx2SAe/l+UPvQdRC714rN42IeIPEKzbh0Emw0AYOvRncXv41j+\nROQR4ectM1zrgVUnqW2x/Imo1ZSHDiH08GEAgBimRt2IETInouaw/Imo1dRff+14bLl1OK/o9QMs\nfyJqNfOUqbjtin04ftcjqL3/frnj0CXg2T5E1Gp796rwk647VItfQIMgdxq6FNzzJ6JWy8rSICOj\nFgKL329wz5+IWsViATZt0mDHjlK5o1ALcM+fiFplxw41jEYr4uJEuaNQC7D8ieiyabKykLPqJO66\nq07uKNRCPOxDRJclpKwM0U8/jaWiiLqS/0HFPesBtVruWHSJuOdPRJdF/dlnEMTGQz0hOg2L38+w\n/Inosmg2bHA8rrvrLhmT0OVg+RNRiylOn0bY3r0AAEmhgOWOO2RORC3F8ieiFlN/+qnjcX3//hBj\nYmRMQ5eD5U9ELWZJT8e7iTNRFXs1D/n4KZ7tQ0QtVqTtgifNf8Og76dCreL5/f6I5U9ELbZpkxqp\nqRaoNQIAhdxx6DLwsA8RtdimTRrcdptF7hjUCix/ImqR8nIB+/eHYvDgermjUCuw/Inokgk1Ndi6\nVY0BA+qh0Uhyx6FW4DF/Irpk+nHjMOpQNfoPuB0hZeN4iqcf454/EV2SkDNnoNq7F51qjiD5y0WA\nxD1/f8byJ6JLov7ySwjnCr/hppsgtm8vcyJqDZY/EV0S9RdfOB5bbr1VxiTkCSx/ImqWUFmJsD17\nHNt1LH+/1+wbvnl5eVixYgVEUURKSgpGjhzpMiY/Px/vvfce7HY7dDodXnzxxbbISkQyUZw+jcrY\nJLQrOoqGHj0gJiTIHYlaqcnyF0URmZmZmDNnDvR6PWbOnInk5GQkJiY6xtTU1CAzMxOzZ8+GwWBA\ndXV1m4cmIu+yGY2YfPMPSL/yEEYPKZY7DnlAk4d9CgsLERcXh9jYWCiVSvTr1w85OTlOY7755hvc\ndNNNMBgMAIDIyMi2S0tEsmhoaLxXb59xV8KanCx3HPKAJvf8y8vLHaUOAHq9HoWFhU5jfvvtN9jt\ndrz00kuoq6vDbbfdhoEDB7ZNWiKSxbffhiEpycabtAeQVl/kZbfbcfz4cbzwwguor6/H888/j65d\nu+KKK67wRD4i8gFbtqiRns61fAJJk+Wv1+thMpkc2yaTCXq93mmMwWCATqeDSqWCSqVCt27d8Msv\nv7iUf35+PvLz8x3bGRkZ0Ol0nvgaCIBKpeJ8egjn0pkkAdu3a/DJJ3WXNS+cT8/LyspyPDYajTAa\njS1+jibLPykpCcXFxSgtLYVer0d2djaeeOIJpzF9+vTB8uXLIYoirFYrCgoKcIebW7q5C2g2m1sc\nmNzT6XScTw/hXP5BeegQzB9sxf9YxyD+igSYzUKLn4Pz6Vk6nQ4ZGRmtfp4my1+hUGDixImYP3++\n41TPxMREbN26FQCQlpaGhIQEdO/eHTNmzIAgCEhNTXU6G4iI/JdmwwbErlyCz7EQZ+dNQvXcuXJH\nIg8RJEm+BTqKiorkeumAw70rz+Fc/qH94MEILSgAAJjefRf16ektfg7Op2fFx8d75Hl4hS8RuaU4\nftxR/KJajYYBA2RORJ7E8icit9TnDu8CQMOAAZA0GhnTkKex/InILfWWLY7Hlss43EO+jeVPRG6V\nvPAq5qpeQW2P3rCkpsodhzyM5U9Ebu0ovgHbkp9C5ecbIXboIHcc8jCWPxG5xat6AxvLn4hciCKw\nbZsaaWks/0DF8iciF3l5oWjXTkSnTna5o1AbYfkTkZOQsjJs3RLGQz4BrtWrehJRAJEkxIwYgRm/\nKdAwbCiEyumQoqPlTkVtgHv+ROSgPHoUyl9/RUfrcXT+6iNI4eFyR6I2wvInIofzr+qtHzwYUKnk\nC0NtiuVPRA68qjd4sPyJCAAQcuYMQnNzAQBSSAgsQ4bInIjaEsufiAAAISUlONPxRgBAQ9++kC64\nax8FFp7tQ0QAANsNN2BacjaGjD2BjJTTcsehNsY9fyICANhswI4datx8tx62y7gnLPkXlj8RAQD2\n7VOhY0cb4uNFuaOQF7D8iQgAsHWrGmlp9XLHIC9h+RMRgMZVPLmQW/Bg+RMFOeWPP8L+9HzcWPUN\n/qdbndxxyEtY/kRBTvP55+j4nyVYXz4E0bNnyR2HvITlTxTkws5f0qF/fxmTkDex/ImCWEhREVQH\nDwIApNDQxvV8KCiw/ImCmHrbNsfjhptvhhQZKWMa8iaWP1EQO38VT0tamoxJyNtY/kRBrOylV/Cs\n+v9gvnkgyz/IcG0foiD2zakkbOvWE9M/uUfuKORl3PMnCmLbtvHCrmDF8icKUpLEq3qDGcufKEgd\nPaqEKALdutnkjkIyYPkTBaGQkhJs26JCeroFgiB3GpID3/AlCkKG++/HUwUVqOqfipCSpyF26CB3\nJPIy7vkTBRnFyZMIPXIEMbZidM5eAykiQu5IJAOWP1GQOf/Crob+/SBptTKmIbmw/ImCzPkLuVmG\nDpUxCcmp2fLPy8vD9OnT8fjjj2P9+vUXHVdYWIh7770X3333nUcDEpHnCGYzwr791rHN8g9eTZa/\nKIrIzMzErFmz8Prrr2PPnj04deqU23GrVq1Cjx49IElSm4UlotYJOXMGpVcnw44QNNxwA8SEBLkj\nkUyaPNunsLAQcXFxiI2NBQD069cPOTk5SExMdBq3adMm3Hzzzfjpp5/aLikRtZr96qvxbPJWXH97\nKR6+lf9eg1mTe/7l5eUwGAyObb1ej/LycpcxOTk5SE9PBwAIPGmYyGeJYuNVvYNGh8N2/fVyxyEZ\ntfoN3xUrVmDs2LEQBAGSJPGwD5EP+/77UBgMIjp3tssdhWTW5GEfvV4Pk8nk2DaZTNDr9U5jfv75\nZ7zxxhsAALPZjLy8PCiVSiQnJzuNy8/PR35+vmM7IyMDOp2u1V8ANVKpVJxPDwnkudyxIwwjRohe\n/foCeT7lkpWV5XhsNBphNBpb/BxNln9SUhKKi4tRWloKvV6P7OxsPPHEE05j3nrrLcfjJUuWoHfv\n3i7Ff7GAZrO5xYHJPZ1Ox/n0kECdS0kCNm7UYOnSCpjNVq+9bqDOp1x0Oh0yMjJa/TxNlr9CocDE\niRMxf/58iKKIlJQUJCYmYuu584TTePMHIr8Qun8/GhZ/hP5n/4wbrjECUMkdiWQmSDIepC8qKpLr\npQMO9648JxDnMnLePEQsXQoAqBk/HlWvvuq11w7E+ZRTfHy8R56HV/gSBTpJgnrzZsemJTVVxjDk\nK1j+RAFOefQolCdOAABErRb1/fvLG4h8AsufKMCpN21yPK5PSQHUahnTkK9g+RMFOPWWLY7HluHD\nZUxCvoTlTxTgCt9chanqf6NmaDosKSlyxyEfwTt5EQW4z/cm4mTa/ahaepvcUciHcM+fKMB9+qkG\nI0bUyR2DfAzLnyiAmUwh2L8/FCkp9XJHIR/D8icKYF98ocaQIRZoNFxwkZyx/IkClCo7G19uDMGI\nERa5o5AP4hu+RAFI8csviBkzBusQBSHxdpwd/neA99qg83DPnygAaT77DAAQhSqoTCUsfnLB8icK\nQOpPP3U8rhsxQsYk5KtY/kQBRnHiBFQHDwIAJJUKlnO3WCU6H8ufKMBoztvrrx80CFJUlIxpyFex\n/IkCjLV7d+yKvgOiQslDPnRRPNuHKMAUdBqC0YoxyN13BIpIjdxxyEex/IkCzNq1GowYYYGig17u\nKOTDeNiHKIBIUmP5jx5dK3cU8nEsf6IAsn9/KOx2Ab16WeWOQj6O5U8UKKxWrF2rwd131/KaLmoW\ny58oAAi1teiQnIzbPngYEztsBERR7kjk41j+RAFAvWULFGVlGN2wGl3/9RKXc6BmsfyJAoAmK8vx\nuG70aJY/NYvlT+TnFCdPImzXLgCAJAiou+cemRORP2D5E/m58NWrIUiNN2upHzwY9oQEmRORP2D5\nE/k7mw21QjgAoPa++2QOQ/6C5U/k53YNm4ObOp5ExYKFsKSlyR2H/ASXdyDycx99FI47xtpRN5Z7\n/XTpWP5EfqymRsBnn2mwY0ep3FHIz/CwD5Ef27hRg759GxAXx4u6qGVY/kT+SJIgScCKFeEYP75G\n7jTkh1j+RH5IvXkzwobdgwEl6zC4P8ufWo7lT+SHtO++i5j8b7HkzL2IXPKW3HHID7H8ifyM8uhR\nhO3ZAwCQFArUjhkjcyLyRyx/Ij+jffddx2PLsGEQeUUvXYZLOtUzLy8PK1asgCiKSElJwciRI50+\nv3v3bmzcuBGSJEGj0eDhhx/GVVdd1SaBiYKZUFkJzSefOLZrHnxQxjTkz5rd8xdFEZmZmZg1axZe\nf/117NmzB6dOnXIa06FDB7z00ktYuHAh7r77bvzrX/9qs8BEwUx54gTOqtoBAKzXXYeGP/1J5kTk\nr5rd8y8sLERcXBxiY2MBAP369UNOTg4SExMdY6655hrH4y5dusBkMrVBVCKqv7EHBsUcw/L73kXX\nmyK4dDNdtmbLv7y8HAaDwbGt1+tRWFh40fE7duxAz549PZOOiJxs2xYGpUaJjrNHoZ69T63g0eUd\nDh06hJ07d2LevHkun8vPz0d+fr5jOyMjAzqdzpMvH9RUKhXn00N8eS6XLdPgqaesiIz0zXzu+PJ8\n+qus827eYzQaYTQaW/wczZa/Xq93OoxjMpmg1+tdxv3yyy9YtmwZZs+ejYiICJfPuwtoNptbHJjc\n0+l0nE8P8dW53LdPhd9+0yAlpQI+GO+ifHU+/ZVOp0NGRkarn6fZN3yTkpJQXFyM0tJS2Gw2ZGdn\nIzk52WlMWVkZFi5ciGnTpiEuLq7VoYjoPKIIWCx4880IPPLIWSi5HCN5QLPfRgqFAhMnTsT8+fMd\np3omJiZi69atAIC0tDSsWbMGNTU1eOeddxx/59VXX23b5ERBQv3ZZ9A8/zL62mbh3iV3AFDJHYkC\ngCBJ5+7/JoOioiK5Xjrg8Fdrz/GpubTb0T41FaEFBQCA6mefxdknnpA5VMv41HwGgPj4eI88D6/w\nJfJhmg0bHMUvRkSgZtw4mRNRoGD5E/kqmw261193bNb85S+Q3JxsQXQ5WP5EPip85Uoojx8HAIiR\nUTj7l7/InIgCCcufyEfVdzPioLo3AODs1CmQoqJkTkSBhCeNEfmo9wsHYV33dHz2/1bAMmyY3HEo\nwLD8iXyQ2Sxg4UId3nuvHJYb75I7DgUgHvYh8kFvvRWBgQPrceONVrmjUIDinj+Rjzl2TIkPPwzH\nli1n5I5CAYx7/kQ+QjCbETn7ebz2dD2eftqMK64Q5Y5EAYzlT+QjdAsWIGLFu/jP/u54JPJDueNQ\ngGP5E/kA1bffQrt8OQAgyl6BEP7LpDbGbzEimQlmM6KnT4dwbpkty5AhqLuLZ/hQ22L5E8ksau5c\nKM/dF9seFY3KhQt5e0Zqcyx/IjlJEsrU8bCf+6dY9eorEHlPDPICnupJJKM6SwhG7H0Nzz5yK+5W\nboCFh3vIS1j+RDKRJGDmzChcd50VaXOMMAstvw8r0eVi+RPJZNWqcOzfH4rPPy/jIX7yOpY/kZeF\nFBXh6586YcECHT75pAzh4bLdTI+CGN/wJfIi7fLlaD9gENb/JRtvv12BLl3sckeiIMXyJ/ISzZo1\niHzhBSgstVh1dhSG1H8pdyQKYix/Ii/QrF+P6CefdFzIZevRHfU33yxzKgpmLH+iNqb+9FNEP/44\nBLFxoTZrt24wvf8+oNHInIyCGcufqI0Vq69CnRgGALBeey1Mq1fzRuwkO5Y/URs6elSJtFlD8cm9\n78F6/fUwrV4N0WCQOxYRT/Ukaiu7d6vw2GPtMHduNVJH34Iz9s2AQiF3LCIALH8ij1IcPw7bVZ2w\neIkOmZlavP12BW65peHcJ1n85Dt42IfIExoaoPv73xE7aBDWDfsYmzer8fnnZ/4ofiIfwz1/olYK\n/eEHRD/3HELz8wEAj/74LEZnXQ3EJ8ucjOjiWP5El6uhAdF//SvCs7KcPmzv2xuKhFjw2l3yZSx/\nostUL4bi5PcVuPbctqgKg3nmc6h5+GHwPozk6/gdStRC9fXA+++Ho/+ADvjf2IUQlaGoGzYMZ3bu\nQM2kSSx+8gvc8ydqiihCtWcPVDk5OD3xKXz8cTiWLdPi2mttWLq0Ar17x+DM8Z2wd+4sd1KiFmH5\nE11IkhC6fz8069dDs3EjFMXFAICUpRORmBqFZcsq0KuX1TGcxU/+iOVPdAH1uHHQbdzo8vEvRiyC\nbeHLMiQi8jyWPwUfSYLi1ClAFGG/6ioAQF0d8N13Ydi5MwzGPb0xHX+Uv71du8Z7604YJ1diIo9r\ntvzz8vKwYsUKiKKIlJQUjBw50mXM8uXLkZeXh7CwMEyZMgWd+Wsw+ZCQ06eh3rYNoUePQnnsGJQ/\n/ghFRQV+7P8A3jC+je+/V+HwYSVuuMGKwYPrMfCVIRBnLIRl2DDUjRyJ+v79gdBQub8MIo9qsvxF\nUURmZibmzJkDvV6PmTNnIjk5GYmJiY4xP/zwA0pKSvDmm2+ioKAA77zzDubPn9/mwSnI1dVBYTIh\npKys8Y/JBEmtbtxDP0eSAJMpBGc/P45bXprl8hTq3O8R+ScRM2ZUo2dPKyIiGtfa12lvRPGwgyx8\nCmhNln9hYSHi4uIQGxsLAOjXrx9ycnKcyj8nJweDBg0CAHTt2hU1NTWorKxEdHR0G8YmnyZJQEND\n4/r1ViuEhgagvh6C3Q77lVe6jq+rg2bTJgj19RDq6iCYzRDOngVCQ2F+9lmX4cqCAsQOHuzy8VLD\nNZj53/E4dUqBkycb/2g0Enpf0RdbLxgrRkQgrncspj9WCSgv+GcQEsLip4DXZPmXl5fDcN7ys3q9\nHoWFhU2OMRgMKC8vv6Ty148f31gU50hqNSr+/W+XcUJdHdo9/LDLxyW1GhWZmW7H6x988Nwg5+cv\nf+899+PHj3cdr9Gg/IMP3I+//363eco//NDteMO997p9ftPq1a7ja2thyMhw+/ymNWvcjg+/4w6E\n2Z2vKZU0GpjWrnU7PsbN4TtJo0HZhg2u46ur0X7oUAh2OyCKgN0OwW6HqNWidO9e5+eQALG8Ch1v\nNLo8jzVKj12fHEFDg4D6egEWC9DQIABlNRj/9DSX8dXqGEw7+QrM5hCcPSugurrxvyhT4oTLaCC6\n8ldc07UBKSkiOna0IzHR3rg3L4Wi5tmxsHfqBOu118J27bWwJyTwfHwKah55w1c6r9BaQr19u9N2\nrSIC48Y53+RCkgCNzYzPd3/l8vdrFREYO1bvNBZoHL85e7fb8X/+s8Fp7O/jt36X7TK+JiQC99zj\nvPa6JAHhdjO27/vO7fhRowxOY4HG8V/9kON2/J13xuDC6dPYzPjmQK7b8bffHuM2T/ahXFy4ZmRN\nSASGD49xyvJ7nu+O5Lt9/rS09o0F3tjxsNsFaKxq5J8+7TL+bKUdXbvGQRQF2GyN4yVJQKQQgSqX\n0UB9dQOmTm2HsDDp3B9ApZIQrRQx3s34cGs1Bg+uh04nQqeToNOJ0Gol6KNtEJPVkKKjYY+JgRgT\nA9FggBgTgwfHVgBqtfMTCQKqFixw8wpEwavJ8tfr9TCZTI5tk8kE/QV3ILqUMQCQn5+P/Pw/Cic1\nNRWGC1rmK41hAAAFeUlEQVQvHMB2uKMG4PoDJhzAVy0c7/ojAQDC3I7XAnD9kdD0eNcfCU2Pd/2R\nAACGi453/ZHQ9PgDLRzv+iMBAPRux0cCqHE7PtLt+AgAP7od7/7/lxKA6+8D59TVAYDLD7yIi41v\nIZ1O56FnIoDz6Ukmkwnbz9txNhqNMBpdf9NuTpO/9yYlJaG4uBilpaWw2WzIzs5GcrLzSoXJycnY\ntWsXAODYsWPQarVuD/kYjUZkZGQ4/mzf7r7m6fJkXbC4GF0+zqVncT49a/v27U5dejnFDzSz569Q\nKDBx4kTMnz/fcapnYmIitm5tfPssLS0NvXr1Qm5uLqZNmwa1Wo3JkydfVhAiIvKeZo/59+zZEz17\n9nT6WFpamtP2Qw895NlURETUpmQ73eFyf1Uh9zifnsO59CzOp2d5aj4F6XJP1SEiIr/FE52JiIIQ\ny5+IKAh5fVXPS1kojlxNnToVGo0GISEhUCgUePXVV3H27Fn84x//QFlZGdq3b48nn3wSWq0WALBu\n3Trs3LkTISEhePDBB9G9e3eZvwJ5LVmyBLm5uYiMjMSiRYsA4LLm7+eff8bixYthtVrRs2dPPPj7\nleRBxN1cZmVlYceOHYiMjAQA3HfffY4TRTiXTSsrK8PixYtRVVUFQRCQmpqK2267re2/PyUvstvt\n0mOPPSaVlJRIVqtVmjFjhnTy5ElvRvBbU6ZMkcxms9PHPvjgA2n9+vWSJEnSunXrpJUrV0qSJEkn\nT56UZsyYIVmtVqmkpER67LHHJLvd7vXMvuTw4cPSzz//LD311FOOj7Vk/kRRlCRJkp577jmpoKBA\nkiRJeuWVV6Tc3FwvfyXyczeXWVlZ0qeffuoylnPZvIqKCun48eOSJElSXV2d9Pjjj0snT55s8+9P\nrx72OX+hOKVS6Vgoji6NdMF78+cvqjd48GDs27cPALBv3z7069cPSqUSsbGxiIuLc1mTKdh069bN\nsdf0u5bMX0FBASoqKmCxWNClSxcAwMCBA7H3grWNgoG7uQTcL/PCuWxedHQ0OnXqBABQq9VISEhA\neXl5m39/evWwz6UsFEfuCYKAefPmISQkBEOHDsXQoUNRVVXluJo6KioKVVWNK+pUVFSga9eujr/7\n+2J75Kyl86dUKp2WLtHr9ZzX82zevBm7du3C1VdfjfHjx0Or1XIuW6i0tBQnTpxA165d2/z7k3fy\n8hPz5s1Du3btUF1djXnz5iEhIcHp84IgNPn3m/t8sOP8tE56ejruueceAMDq1avx/vvv82r/FrJY\nLFi0aBEmTJgAjUbj9Lm2+P706mGfS10Ejly1a9cOABAZGYm+ffuisLAQUVFRqKysBNC4NxAVFQWA\n83ypWjJ/BoPBZU+K8/qHqKgoCIIAQRCQkpLi+I2ec3lpbDYbFi1ahIEDB6Jv374A2v7706vlfykL\nxZGr+vp61J1bxdJiseDAgQO48sorkZycjK+++goA8PXXX6NPnz4AGhfb27NnD2w2G0pLS1FcXOw4\nDkh/aOn8RUdHQ6PRoKCgAJIkYffu3Y5/qMGuoqLC8Xjv3r248txNeziXzZMkCUuXLkVCQgJuv/12\nx8fb+vvT61f45ubmOp3qOWrUKG++vF8qLS3FgnPr0YuiiP79+2PUqFFNngq2du1a7Ny5EwqFAhMm\nTECPHj3k/BJk98Ybb+DIkSOorq5GdHQ0MjIy0KdPnxbP3++n0jU0NKBnz56YOHGinF+WLC6cyzFj\nxuDw4cM4ceIEBEFA+/btMWnSJMfxas5l03788UfMnTsXV155pePwztixY9GlS5c2/f7k8g5EREGI\nV/gSEQUhlj8RURBi+RMRBSGWPxFREGL5ExEFIZY/EVEQYvkTEQUhlj8RURD6/8uhi85eensqAAAA\nAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot(t, y, 'b-', lw=1)\n", "plot(t, f[:,0], 'r--', lw=3)" ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [], "source": [ "Tfix = lambda N,s: -2*log(N)/log(1-s)" ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [], "source": [ "N = logspace(3,8,50)\n", "T = Tfix(N,s)" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEJCAYAAACQZoDoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8VPWd//HX9yQGQgyJiaHhIl4YFAwCYwMVo5GLSwWj\n0i6mSFvbFGSbdtdtHrZe2K2XIsvP7XJdU/VBd7t1tyrRcgfxAgECqSZYojaARil1A42RmRCIEHKZ\n7++PaMZZgSTkcibJ+/l4+HhkTs6ZeZ+PIe+cOTNnjLXWIiIi8hnH7QAiIhJeVAwiIhJCxSAiIiFU\nDCIiEkLFICIiIVQMIiISIrI1KwUCAR588EESEhJ48MEHqampYenSpRw9epSkpCRycnKIiYkBYM2a\nNeTn5+M4DllZWYwZMwaAgwcPkpubS319PV6vl6ysrM7bKxEROW+tOmLYvHkzQ4YMwRgDwNq1axk9\nejTLly9n1KhRrF27FoDy8nIKCwtZsmQJ8+fP59e//jWfv01i5cqVZGdns2LFCioqKigpKWlVwNLS\n0vPZrx5HcwjSLII0iyDNIqi9s2ixGHw+H3v37mXy5MnNv+T37NnDTTfdBMDEiRMpLi4GoLi4mLS0\nNCIjIxkwYADJycmUlZVRVVVFbW0tHo8HgPT0dIqKiloVUP+zm2gOQZpFkGYRpFkEdXox/Pa3v+U7\n3/kOjhNctbq6mvj4eADi4uKorq4GoKqqisTExOb1EhMT8fv9VFVVkZCQ0Lw8ISEBv9/fruDnqy0D\na826Z1untcvPdbuzf9A1i7M/dnvX1SxaXudMy1uzTLM48+2OnMU5i+Gtt96if//+XH755Zztyhmf\nP73UXegXwNkfu73rahYtr6NZnHu5iqFtyzprFuZc10p67rnnKCgowHEc6uvrOXXqFOPHj+fDDz/k\n0UcfJT4+nqqqKh577DGWLVvWfK5hxowZACxcuJDMzEySkpJ47LHHWLp0KQC7du1i3759zJs370uP\nWVpaGrKDN998c8jRhoiInJvP52Pr1q3Nt1NSUkhJSWn19ud8VdLs2bOZPXs2APv27WP9+vX8wz/8\nA//zP//D9u3bmTFjBjt27GDcuHEApKamsnz5cjIyMvD7/VRUVODxeDDGEB0dTVlZGR6Ph4KCAqZN\nm3bGxzzTDhw5cqTVO9RTxcbGcuLECbdjhAXNIkizCNIsggYNGkRmZuZ5b9+ql6t+7vOnjWbMmMHS\npUvJz89vfrkqwJAhQ5gwYQI5OTlEREQwZ86c5m3mzp1Lbm4udXV1eL1exo4de96hRUSk85zzqaRw\noSMG/TX0RZpFkGYRpFkEDRo0qF3b653PIiISQsUgIiIhVAwiIhJCxSAiIiFUDCIiEkLFICIiIVQM\nIiISQsUgIiIhVAwiIhJCxSAiIiFUDCIiEkLFICIiIVQMIiISQsUgIiIhVAwiIhJCxSAiIiFUDCIi\nEkLFICIiIVQMIiISQsUgIiIhVAwiIj2ErTlOYPWz7b6fyA7IIiIiLrIna7CvrcNu34y59vp235+K\nQUSkm7KnTmK3rsdu3YgZMx5n/mJMUnK771fFICLSzdjaU9htG7Gvr8ekXIvz0L9iBgzqsPtXMYiI\ndBP29Gns9s3YV1ZjRozG+dkizMAhHf44KgYRkTBn6+uwO7Zgt6wGzwic+x7HDL600x5PxSAiEqZs\nfT1216vYzS/BZR6cex/GDL2i0x9XxSAiEmZsQz22cCt204sw+FKcH8/HXDa8yx7/nMVQV1fHo48+\nSn19PQ0NDYwbN47Zs2eTl5fHtm3b6N+/PwB33XUXXq8XgDVr1pCfn4/jOGRlZTFmzBgADh48SG5u\nLvX19Xi9XrKysjp510REuhfb2Ih9Ix+7cRUMGIQz72eYYSO6PMc5iyEqKopHHnmEPn360NjYyMMP\nP8yBAwcwxpCRkUFGRkbI+uXl5RQWFrJkyRL8fj8LFixgxYoVGGNYuXIl2dnZeDweFi1aRElJCWPH\nju3UnRMR6Q5soBH75k7sxhfgootxsn6CuTLFtTwtPpXUp08fABoaGggEAsTExABgrf3SusXFxaSl\npREZGcmAAQNITk6mrKyMpKQkamtr8Xg8AKSnp1NUVKRiEJFezQYascW7mgrhwjic7/4YM2K027Fa\nLoZAIMADDzzAxx9/zNSpU7nkkkt444032LJlCzt37uSKK67g7rvvJiYmhqqqKoYPDz4PlpiYiN/v\nJzIykoSEhOblCQkJ+P3+ztkjEZEwZwMB2PsHAuufhz59cWbNg6vHYoxxOxrQimJwHIdf/vKXnDx5\nkoULF1JaWsrUqVOZOXMmAKtWreLZZ58lOzu708OKiHRn1lp4+00C656HiAicmVkw6tqwKYTPtfpV\nSf369cPr9fLhhx+SkhJ87mvy5Mk88cQTQNORgM/na/6ez+cjMTHxS0cIPp8v5Ajii0pLSyktLW2+\nnZmZSWxsbOv3qIeKiorSHD6jWQRpFkHhPAtrLQ0lb1L74m+wjY3EzJpD5Fev79RCyMvLa/46JSUl\n5Pd2S85ZDMePHyciIoKYmBjq6up49913mTlzJseOHSM+Ph6AoqIihg4dCkBqairLly8nIyMDv99P\nRUUFHo8HYwzR0dGUlZXh8XgoKChg2rRpZ3zMM+3AiRMnWr1DPVVsbKzm8BnNIkizCArHWVhrYX8J\ngXXPQe0pnNtng/c6ah0Hamo67XFjY2PJzMw87+3PWQzHjh0jNzeXQCCAtZb09HSuueYannzySQ4d\nOoQxhqSkJObNmwfAkCFDmDBhAjk5OURERDBnzpzmRpw7dy65ubnU1dXh9Xp14llEejR74J2mQqip\nxtx2Fyb1BozTPT7pwNgzvbwozBw5csTtCK4Lx7+G3KJZBGkWQeEyC/t+KYH1z0HVUcxtszDj0zFO\nRJdmGDSofRfU0zufRUQ6gP3wQFMhVP4VkzELc91ETETXFkJHUTGIiLSDPVTW9LLTw3/B3Hon5vqb\nMZHd+1dr904vIuIS+9HBpiOEv3yImX4nJvshzAUXuB2rQ6gYRETawB7+S9MRwof7Mbf8Lebv7sdc\nEOV2rA6lYhARaQX713LshuexB97BfP2bmB/kYD67ZFBPo2IQETkHW3kEu2EVtvSPmL+5A+fuv8f0\njXY7VqdSMYiInIH9pAK7aRX27SLMlNtwZj+Die7ndqwuoWIQEfkC6/sEuzkP+8dCzMRbcRY+g+l3\noduxupSKQUQEsFU+7OYXscUFmPSpOAuewlzY3+1YrlAxiEivZqursC+/hP1DPuaGv8FZ8CtMbJzb\nsVylYhCRXsmeqMZuWY3d9RpmwiScX+Ri4i5yO1ZYUDGISK9ia45jX12L3fkKZtyNOI/+O+aiRLdj\nhRUVg4j0CvZkDfa1ddj8zZivXo/z82WYxCS3Y4UlFYOI9Gj21Ens1vXYrRsxY8bj/NNiTFKy27HC\nmopBRHokW3sKu20j9vX1mJRrcR78V8xX2nc56t5CxSAiPYo9fRq7fRP2lTWYEaNxfrYIM3CI27G6\nFRWDiPQItu40gdfXYbesBs9InPsexwy+1O1Y3ZKKQUS6NVtfj931Kse3/B47dBjOvY9ghl7hdqxu\nTcUgIt2SbajHFm7FbsqDQZdy4X2Pc2qAziF0BBWDiHQrtrER+0Y+duMqGDAIZ979mGEjiIyNhTD4\nzOeeQMUgIt2CDTRi39yJ3fgCXHQxTtZPMFemuB2rR1IxiEhYs4EAds8u7Ibn4cI4nO/+GDNitNux\nejQVg4iEJRsIwN4/NH2MZp++OLPmwdVjMca4Ha3HUzGISFix1sLbbxJY9zxERODMzIJR16oQupCK\nQUTCgrUW/vQWgXXPQWMjzh2zYcx4FYILVAwi4iprLewvaSqE2lM4t88G73UYx3E7Wq+lYhAR19gD\n7zQVQk015ra7MKk3qBDCgIpBRLqcLdtHYN3vwP9JUyF8LR3jRLgdSz6jYhCRLmM/PEBg/XPw8RFM\nxrcwEyZjIlQI4eacxVBXV8ejjz5KfX09DQ0NjBs3jtmzZ1NTU8PSpUs5evQoSUlJ5OTkEBMTA8Ca\nNWvIz8/HcRyysrIYM2YMAAcPHiQ3N5f6+nq8Xi9ZWVmdv3ciEhbsobKml50ePoSZnolJm4KJvMDt\nWHIW5yyGqKgoHnnkEfr06UNjYyMPP/wwBw4cYM+ePYwePZo77riDtWvXsnbtWr797W9TXl5OYWEh\nS5Yswe/3s2DBAlasWIExhpUrV5KdnY3H42HRokWUlJQwduzYrtpPEXGB/egggQ3Pw6EyzPQ7MdkP\nYS5QIYS7Fs/y9OnTB4CGhgYCgQAxMTHs2bOHm266CYCJEydSXFwMQHFxMWlpaURGRjJgwACSk5Mp\nKyujqqqK2tpaPB4PAOnp6RQVFXXWPomIy+zhj2h86v8RWPEY5qpROAufwZl0q0qhm2jxHEMgEOCB\nBx7g448/ZurUqVxyySVUV1cTHx8PQFxcHNXV1QBUVVUxfPjw5m0TExPx+/1ERkaSkJDQvDwhIQG/\n39/R+yIiLrN/LcdueB574B3M17+J+cFPMH36uh1L2qjFYnAch1/+8pecPHmShQsX8qc//Snk+x39\n5pPS0lJKS0ubb2dmZhIbG9uhj9EdRUVFaQ6f0SyCwmUWjRWHqf39b2koKaLPrXfS58cPYfpGd2mG\ncJlFuMjLy2v+OiUlhZSU1l9wsNWvSurXrx9er5eDBw8SFxfHsWPHiI+Pp6qqiri4OKDpSMDn8zVv\n4/P5SExM/NIRgs/nCzmC+KIz7cAJXUqX2NhYzeEzmkWQ27Own1RgN63Cvl2EmXwbZuEz1Ef3o76+\nAeq7NpfbswgnsbGxZGZmnvf25zzHcPz4cT799FOg6RVK7777Lpdffjmpqals374dgB07djBu3DgA\nUlNT2b17Nw0NDVRWVlJRUYHH4yE+Pp7o6GjKysqw1lJQUMD48ePPO7SIuMv6PiHw37kEFt7XdAns\nx5/BuW0WJrqf29GkA5zziOHYsWPk5uYSCASw1pKens4111zD5ZdfztKlS8nPz29+uSrAkCFDmDBh\nAjk5OURERDBnzpzmp5rmzp1Lbm4udXV1eL1evSJJpBuyVT7s5hexRTsxN30dZ8FTmNj+bseSDmas\ntdbtEC05cuSI2xFcp8PkIM0iqKtmYaursC+/hP1DPuaGm5tOLPeP7/THbQv9XAQNGtS+jzjVO59F\n5KzsiWrsltXYXa9hJkzCeexJTPyZzw9Kz6FiEJEvsTXHsa+uxe58BTPuBpxHVmASLnY7lnQRFYOI\nNLMna7Cvrcdu34S59nqcny/DJCa5HUu6mIpBRLCnTmK3rsdu3YgZMw5n/mJMUrLbscQlKgaRXszW\nnsLmb8K+tg6T4sV58F8xX2nfiUvp/lQMIr2QPX0au2Mz9pU1mKuuwfnZv2AGXuJ2LAkTKgaRXsTW\n12F3bMFuWQ3DRuDk/AIz5DK3Y0mYUTGI9AK2vh676zXs5hfhMg/OvQ9jhl7hdiwJUyoGkR7MNjRg\nC1/HbsqDQZfi/Hg+5rLhLW4nvZuKQaQHso2N2DfysRtXwYCBOPPuxwwb4XYs6SZUDCI9iA00Yot2\nYje80HRxu6yfYK5s/eWWRUDFINIj2EAA+9Zu7Prn4cJYnO/+GDNitNuxpJtSMYh0Y02FUNj0ucpR\nfXC+NRdSvB3+AVrSu6gYRLohay28XUTNplUELDh/+z0Y9VUVgnQIFYNIN2KthT/9kcC630FjAzGz\n5nDqytEqBOlQKgaRbsBaC/tLCKx7Dk6dxLljNngncEFcHLX6DALpYCoGkTBn33u36QjhRDXmtrsw\nqWkYJ8LtWNKDqRhEwpT9YF/TEYKvsqkQxqdjIlQI0vlUDCJhxh58r6kQPj6MyfgW5rpJmEj9U5Wu\no582kTBh//JBUyEcPoSZnolJm4KJvMDtWNILqRhEXGY/Otj0PoRDZZjpd2KyH8JcoEIQ96gYRFxi\nD39EYP1z8OF+zC3fxNzzU0xUH7djiagYRLqa/Ws5dsPz2APvYL7+DcwPfoLp09ftWCLNVAwiXcRW\nHsFuWIUt/SPm5ttx7v57TN9ot2OJfImKQaST2U8qsJvysG+/iZl8G87sZzDR/dyOJXJWKgaRTmL9\nn2A3vYh9azdm0nScx5/BxFzodiyRFqkYRDqYPebDbn4RW1SASZ+K8/hTmAv7ux1LpNVUDCIdxFZX\nYbf8Hlu4DXPDzTi/yMX0j3c7lkibtVgMR48eJTc3l+rqaowxTJkyhenTp5OXl8e2bdvo37/pL6G7\n7roLr9cLwJo1a8jPz8dxHLKyshgzZgwABw8eJDc3l/r6erxeL1lZWZ24ayJdw56oxm5Zjd31GmbC\nJJzHnsTEJ7gdS+S8tVgMkZGRfO973+Oyyy6jtraWBx54gNGjmy7zm5GRQUZGRsj65eXlFBYWsmTJ\nEvx+PwsWLGDFihUYY1i5ciXZ2dl4PB4WLVpESUkJY8eO7bSdE+lMtuY49tW12J2vYMbdiPPov2Mu\nSnQ7lki7tVgM8fHxxMc3HQ737duXwYMH4/f7gc8uBfx/FBcXk5aWRmRkJAMGDCA5OZmysjKSkpKo\nra3F4/EAkJ6eTlFRkYpBuh17sgb72nrs9k2Ya6/H+fkyTGKS27FEOkybzjFUVlZy6NAhrrzySt57\n7z22bNnCzp07ueKKK7j77ruJiYmhqqqK4cOHN2+TmJiI3+8nMjKShITg4XVCQkJzwYh0B/bUSezW\nDditGzBjxuHMX4xJSnY7lkiHa3Ux1NbWsmTJEr7//e/Tt29fpk6dysyZMwFYtWoVzz77LNnZ2e0O\nVFpaSmlpafPtzMxMYmNj232/3V1UVJTm8JmunoWtPcXpV9ZyelMekaPH0XdBLhEDh3TZ45+Lfi6C\nNItQeXl5zV+npKSQkpLS6m1bVQwNDQ0sXryYG2+8kfHjxwMQFxfX/P3JkyfzxBNPAE1HAj6fr/l7\nPp+PxMTELx0h+Hy+kCOIc+3ACX1CFbGxsZrDZ7pqFvb0aeyOzdhX1mCuugbz04UEBl7CSYAw+X+h\nn4sgzSIoNjaWzMzM897eaWkFay1PP/00gwcP5tZbb21eXlVV1fx1UVERQ4cOBSA1NZXdu3fT0NBA\nZWUlFRUVeDwe4uPjiY6OpqysDGstBQUFzSUjEk5sfR2B19cT+Ke/w354ACfnFzjzfoYZeInb0US6\nRItHDO+99x4FBQUMHTqU+++/H2h6aeru3bs5dOgQxhiSkpKYN28eAEOGDGHChAnk5OQQERHBnDlz\nmj+ofO7cueTm5lJXV4fX69WJZwkrtr4eu+s17OYX4dJhOPf+HDN0mNuxRLqcsWd6aVGYOXLkiNsR\nXKfD5KCOnoVtaMAWbsVuyoNBl+Dc/m3M5cNb3jAM6OciSLMIGjRoULu21zufpdeyjY3YN7ZjN74A\nSck49/wU4xnpdiwR16kYpNexgUZs0U7shhcgPhEn6x8xV45yO5ZI2FAxSK9hAwHsW7ux65+HmAtx\nvvMjGDG6+RyYiDRRMUiPZwMB2PtG0+cqR/XBmXUPXD1WhSByFioG6bGstfB2UdPnKhsH55t3wzWp\nKgSRFqgYpMex1sKf3iKw7jlobMC5fTaM/ZoKQaSVVAzSY1hrYX9JUyGcOolz+11w7fUYp8X3cYrI\nF6gYpEew771LYN3v4Hg15rZZmHE3YJwIt2OJdEsqBunW7Af7CKz9Hfg/wdx2F2Z8OiZChSDSHioG\n6ZbswfeanjL6+DAm41uY6yZhIvXjLNIR9C9JuhX7lw+o2ZRH4C8fYG7NxFw/BRN5gduxRHoUFYN0\nC/Z//0xg/fNw6H36fOM7BO75GeYCFYJIZ1AxSFizhz8isOE5+GA/5pZvYu65jz6JF1Oni6WJdBoV\ng4QlW1GO3fACdv/bmKkzMFk/wfTp63YskV5BxSBhxVYewW5chX33LczNt+N890eYvv3cjiXSq6gY\nJCzYTyqwm/Kwb7+JmXwbzsJnMP1i3I4l0iupGMRV1v8JdtOL2Ld2YyZOw3n8GUzMhW7HEunVVAzi\nCnvMh938ErZoJyZ9Ks7jT2Eu7O92LBFBxSBdzB6vwr78e2zhNswNN+P8IhfTP97tWCLyBSoG6RL2\nxHHsK6uxu17DXDcR57EnMfEJbscSkTNQMUinsp+ewL66FrtjCyY1Defh5ZiEi92OJSLnoGKQTmFP\n1mBfX4/N34TxTsD5+VJM4gC3Y4lIK6gYpEPZ2pPY1zdgt27AjB6HM38xJinZ7Vgi0gYqBukQtvYU\nNn8z9rW1mKvH4jzwBCZ5sNuxROQ8qBikXezp09gdL2NfWY256hqcn/0LZuAlbscSkXZQMch5sfV1\n2J2vYF/+PQy7CifnF5ghl7kdS0Q6gIpB2sTW12N3v4bd/BIMvQLn3p9jhg5zO5aIdCAVg7SKbWjA\n/mEbdlMeDLwEJ/tBzOVXuh1LRDpBi8Vw9OhRcnNzqa6uxhjDlClTmD59OjU1NSxdupSjR4+SlJRE\nTk4OMTFNFz1bs2YN+fn5OI5DVlYWY8aMAeDgwYPk5uZSX1+P1+slKyurc/dO2s02NmLf2I7d+AIk\nJePMvQ/jGel2LBHpRC0WQ2RkJN/73ve47LLLqK2t5YEHHmD06NFs376d0aNHc8cdd7B27VrWrl3L\nt7/9bcrLyyksLGTJkiX4/X4WLFjAihUrMMawcuVKsrOz8Xg8LFq0iJKSEsaOHdsV+yltZAON2KKd\n2A0vQHwiTtY/Yq4c5XYsEekCTksrxMfHc9lllwHQt29fBg8ejN/vZ8+ePdx0000ATJw4keLiYgCK\ni4tJS0sjMjKSAQMGkJycTFlZGVVVVdTW1uLxeABIT0+nqKiok3ZLzpcNBAgUFxB45B+w21/G+c6P\ncH66UKUg0ou06RxDZWUlhw4dYvjw4VRXVxMf33Txs7i4OKqrqwGoqqpi+PDhzdskJibi9/uJjIwk\nISF4bZyEhAT8fn9H7IN0ABsIQMkbTZ+rHNUHZ9Y9cPVYjDFuRxORLtbqYqitrWXx4sV8//vfJzo6\nOuR7+uXRfVlr4Z1iAut+B8bB+ebdcE2q/p+K9GKtKoaGhgYWL15Meno648ePB5qOEo4dO0Z8fDxV\nVVXExcUBTUcCPp+veVufz0diYuKXjhB8Pl/IEcTnSktLKS0tbb6dmZlJbGzs+e1dDxIVFdWhc7DW\n0vB2MbUv/gZbX0e/zB9wQWpatyiEjp5Fd6ZZBGkWofLy8pq/TklJISUlpdXbtlgM1lqefvppBg8e\nzK233tq8PDU1le3btzNjxgx27NjBuHHjmpcvX76cjIwM/H4/FRUVeDwejDFER0dTVlaGx+OhoKCA\nadOmfenxzrQDJ06caPUO9VSxsbEdMgdrLex/m8D65+Dkpzi33wXXXs9px+F0TU0HJO18HTWLnkCz\nCNIsgmJjY8nMzDzv7Vsshvfee4+CggKGDh3K/fffD8Ds2bOZMWMGS5cuJT8/v/nlqgBDhgxhwoQJ\n5OTkEBERwZw5c5r/Cp07dy65ubnU1dXh9Xr1iqQuZt/7E4H1v4PqY5jbZmHG3YBxItyOJSJhxlhr\nrdshWnLkyBG3I7iuPX8N2Q/2EVj3HPgqMRmzMF+7CRPRfQtBfxkGaRZBmkXQoEGD2rW93vncg9mD\n7zU9ZVRxGJPxLcx1kzCR+l8uIuem3xI9kP3Lh02vMjp8CDM9E5M2BRN5gduxRKSbUDH0ILb8zwTW\nPQ+H3sdMm4nJfghzgQpBRNpGxdAD2MMfYTc8j/1gH+br38Tccx8mqo/bsUSkm1IxdGO2ohy74QXs\n/rcxU2c0Xc+oT1+3Y4lIN6di6IZs5RHsxlXYd9/C3Hw7znd/hOnbz+1YItJDqBi6EXv0Y+ymPGzJ\nG5hJGTgLn8H0i3E7loj0MCqGbsD6P+HkqpUE/rAdc9MtOI8/jYnRW/9FpHOoGMKYPebDbn4J++YO\nzJQMnAVPYWL7ux1LRHo4FUMYssersC+vxhZuxaRNwVmQS/TgoTToXZ0i0gVUDGHEnjiOfWU1tuBV\nzHUTcR57EhP/5SvQioh0JhVDGLCfnsC+uha7Ywtm3A04jyzHJCS5HUtEeikVg4vsyRrs6+ux+Zsw\n3gk4/7wEc/FX3I4lIr2cisEF9tRJ7NYN2K0bMNek4jz0b5gBA92OJSICqBi6lK09hc3fhH1tHWbk\nWJwHnsAkD3Y7lohICBVDF7CnT2N3bMa+sgZz5Sicny7EDBrqdiwRkTNSMXQiW1+H3fkK9uXfw7Cr\ncHIewwy53O1YIiLnpGLoBLa+Hrv7Nezml2DoFTj3/hwzdJjbsUREWkXF0IFsQwP2D9uwm/Jg4CU4\n2Q9hLh/udiwRkTZRMXQA29iIfXM7duMquPgrOHPvw3hGuh1LROS8qBjawQYascW7sBtegLh4nO/d\ni7lqlNuxRETaRcVwHmwggH2rELvheegXg/PtH8KI0Rhj3I4mItJuKoY2sNbC3jcIrH8OLojCyZwD\nKV4Vgoj0KCqGVrDWwjvFTYUAON+4G0anqhBEpEdSMZyDtRZK/0hg3XPQUI9z+2wY+zUVgoj0aCqG\nM7DWwoF3CKz7HZz8FHPbXZivXo9xHLejiYh0OhXD/2Hf/1NTIVQfw9w2CzPuBowT4XYsEZEuo2L4\njP1gf9M5hKMfYzK+hfnaREyECkFEep8Wi+FXv/oVe/fupX///ixevBiAvLw8tm3bRv/+TZ8/fNdd\nd+H1egFYs2YN+fn5OI5DVlYWY8aMAeDgwYPk5uZSX1+P1+slKyurs/apTeyf328qhL+WY27NxEyY\njIlUX4pI79Xib8BJkyYxbdo0nnzyyeZlxhgyMjLIyMgIWbe8vJzCwkKWLFmC3+9nwYIFrFixAmMM\nK1euJDs7G4/Hw6JFiygpKWHs2LEdv0etZD/6sOmk8v/+GTP9TsyP/wkTeYFreUREwkWLxTBy5Egq\nKyu/tNxa+6VlxcXFpKWlERkZyYABA0hOTqasrIykpCRqa2vxeDwApKenU1RU5Eox2PI/E1j/PBx8\nHzNtJuaHD2AuiOryHCIi4eq8nzPZsmULO3fu5IorruDuu+8mJiaGqqoqhg8PXjQuMTERv99PZGQk\nCQnBD7X/jZO4AAAJTklEQVRPSEjA7/e3L3kb2SMfYdc/jy0rxdzyt5i592Gi+nRpBhGR7uC8imHq\n1KnMnDkTgFWrVvHss8+SnZ3docE6iq0ox25Yhd1fgpk6AyfrHzF9+rodS0QkbJ1XMcTFxTV/PXny\nZJ544gmg6UjA5/M1f8/n85GYmPilIwSfzxdyBPFFpaWllJaWNt/OzMwkNja2zRkbKw5Tu/q/adj7\nBn2mz6RP9v2Y6H5tvp9wERUVdV5z6Ik0iyDNIkizCJWXl9f8dUpKCikpKa3e9ryKoaqqiosuugiA\noqIihg5t+pjK1NRUli9fTkZGBn6/n4qKCjweD8YYoqOjKSsrw+PxUFBQwLRp085432fagRMnTrQ6\nm/VVYjeuwpa8gZmUgXn8aer7xVDf0AhtuJ9wExsb26Y59GSaRZBmEaRZBMXGxpKZmXne27dYDMuW\nLWP//v0cP36c7Oxs7rzzTvbt28ehQ4cwxpCUlMS8efMAGDJkCBMmTCAnJ4eIiAjmzJnTfPmIuXPn\nkpubS11dHV6vt8NPPFv/J9jNL2L37MbcdAvO409jYvTXg4hIWxl7ppcXhZkjR46c9Xv2mA+7+SXs\nmzswN07FfP2bmNj+XZiua+ivoSDNIkizCNIsggYNGtSu7bvtO7ns8Srsy6uxhVsxaVNwFuRi+l/k\ndiwRkW6v2xWDPXEc+8pq7K7XMNdNxHnsSUz8mU9ki4hI23WbYrCfnsC+ug6742XMuBtwHl6OSbjY\n7VgiIj1OtyiGwPrnsPmbMN4JOP+8BHPxV9yOJCLSY3WLYuBoJc5D/4YZMNDtJCIiPV63KAbnBz9x\nO4KISK+hjyQTEZEQKgYREQmhYhARkRAqBhERCaFiEBGRECoGEREJoWIQEZEQKgYREQmhYhARkRAq\nBhERCaFiEBGRECoGEREJoWIQEZEQKgYREQmhYhARkRAqBhERCaFiEBGRECoGEREJoWIQEZEQKgYR\nEQmhYhARkRAqBhERCRHZ0gq/+tWv2Lt3L/3792fx4sUA1NTUsHTpUo4ePUpSUhI5OTnExMQAsGbN\nGvLz83Ech6ysLMaMGQPAwYMHyc3Npb6+Hq/XS1ZWVifuloiInK8WjxgmTZrE/PnzQ5atXbuW0aNH\ns3z5ckaNGsXatWsBKC8vp7CwkCVLljB//nx+/etfY60FYOXKlWRnZ7NixQoqKiooKSnphN0REZH2\narEYRo4c2Xw08Lk9e/Zw0003ATBx4kSKi4sBKC4uJi0tjcjISAYMGEBycjJlZWVUVVVRW1uLx+MB\nID09naKioo7eFxER6QDndY6hurqa+Ph4AOLi4qiurgagqqqKxMTE5vUSExPx+/1UVVWRkJDQvDwh\nIQG/39+e3CIi0knaffLZGNMROUREJEy0ePL5TOLi4jh27Bjx8fFUVVURFxcHNB0J+Hy+5vV8Ph+J\niYlfOkLw+XwhRxBfVFpaSmlpafPtKVOmMGjQoPOJ2ePExsa6HSFsaBZBmkWQZtHE5/OxdevW5tsp\nKSmkpKS0evvzOmJITU1l+/btAOzYsYNx48Y1L9+9ezcNDQ1UVlZSUVGBx+MhPj6e6OhoysrKsNZS\nUFDA+PHjz3jfKSkpZGZmNv/3xZ3rCHl5eR267tnWae3yc90+29cdRbM4+2O3d13NouV1zrS8Ncs0\nizPf/uLXW7duDfk92pZSAIh49NFHHz3XCsuWLSMvL4+jR4+ydetWYmJimDJlCuvWrWP16tV8+umn\nZGVlERUVRf/+/ampqeHpp59m9+7d/OAHP2DgwIEAXH755Tz11FNs3LgRj8fDtGnTWhWwtLS0zTvV\nkgEDBnToumdbp7XLz3X78687Yw5neuz2rqtZtLyOZnHu5a1Zplmc+XZHzcLYz19PGqby8vLIzMx0\nO4brNIcgzSJIswjSLILaO4uwf+dzZ/wF0B1pDkGaRZBmEaRZBLV3FmF/xCAiIl0r7I8YRESka6kY\nREQkhIpBRERCqBhERCTEeb3z2S2HDx9m8+bNnDhxgrFjxzJ58mS3I7mqtraWxx57jDvvvJNrr73W\n7TiuKS0tZdWqVVxyySWkpaVx9dVXux3JNdZaXnjhBU6dOsWwYcOaL3bZGx04cICCggIaGxs5fPgw\nCxYscDuSa3w+H7/5zW+IiYlh4MCBzJgx45zrd6tiGDx4MPfccw+BQIBly5b1+mJYv349EyZMcDuG\n64wxREdHU19ff9ZLrfQWxcXF+P1+YmNjQy5o2RuNGDGCESNGUFxc3Hxl597qo48+4mtf+xo33ngj\ny5Yta3F914vhTB8EBFBSUsJ//dd/EQgEmDx5cnPD7dmzh1dffZUpU6a4FbnTtGUW77zzDkOGDKGu\nrs7FxJ2nLbMYOXIkV199NdXV1fz2t7/l3nvvdTF5x2vLLI4cOcJVV13FzTffzJIlSxg1apSLyTte\nW39fAOzatYvs7Gw34naqtsziqquu4oknniA/P5/09PQW79v1cwxn+iCgQCDAf/zHfzB//nyWLFnC\n7t27KS8vB5quxzR//nx27NjhRtxO1ZZZ7Nu3j/fff5/du3fz+uuv09PejtKWWXx+hd+YmBgaGhrc\niNup2jKLxMTE5s9P6YlXPm7r74ujR4/Sr18/+vbt60bcTtWWWeTn5zNr1iwefvhh/vjHP7Z4364f\nMYwcOZLKysqQZR988AHJycnN1/1IS0tjz549HD9+nDfffJP6+voe+S7Htsxi1qxZAGzfvp3+/fv3\nuF8CbZnFkSNHKCkp4eTJk9xyyy1uxO1UbZnF9OnT+c///E/279/f6/+NDBkyhG3btjFp0iQ3ona6\ntszC6/Xy0ksvsWvXrlZdz8n1YjgTv98f8vxoQkICH3zwAVdffXWvO7F4tll8buLEiS6kcsfZZjFj\nxoyzXq23pzrbLKKiovjhD3/oYrKud65/I73t2klnm8Wll17Kfffd1+r7cf2pJBERCS9hWQxn+sCf\n3vpqE80iSLMI0iyCNIugjppFWBbDsGHDqKiooLKykoaGBgoLC0lNTXU7lis0iyDNIkizCNIsgjpq\nFq5fXXXZsmXs37+fEydOEBcXR2ZmJpMmTWLv3r0hL7n6xje+4WbMLqFZBGkWQZpFkGYR1JmzcL0Y\nREQkvITlU0kiIuIeFYOIiIRQMYiISAgVg4iIhFAxiIhICBWDiIiEUDGIiEgIFYOIiIRQMYiISIj/\nD6kM1Wde2958AAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot(N,T)\n", "xscale('log')\n", "grid(False,'minor')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python [default]", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.2" } }, "nbformat": 4, "nbformat_minor": 1 }