Eden_model.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3006434b-2456-4572-a0ac-ae50dc1bda6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "###################################################################\n",
    "def FFGN():\n",
    "    B = 2\n",
    "    a = 6\n",
    "    θ = 0.4\n",
    "    T = 1\n",
    "    n = 1\n",
    "    \n",
    "    #decay rates\n",
    "    def u():   \n",
    "        return a * pow(B, -n)\n",
    "\n",
    "    #uncorrelated uniformly deviates\n",
    "    def r():\n",
    "        return pow(e, -u())\n",
    "\n",
    "    #uncorrelated uniformly distributed noise   ？：取值范围如何界定，为什么\n",
    "    uu_noise = np.random.uniform(0,1,T)\n",
    "\n",
    "    #weight function\n",
    "    def W():       \n",
    "        return np.sqrt((12 * (1 - r()**2) * (pow(B,0.5-θ) - pow(B,θ-0.5)) * pow(a * pow(B,-n), 1-2*θ)) / math.gamma(2-2*θ))\n",
    "\n",
    "    #auxiliary functions 辅助方程，不知道什么意思\n",
    "    def X():\n",
    "                                    # The current time\n",
    "        if T == 1:\n",
    "            x1 = (uu_noise - 0.5) / np.sqrt(1 - r()**2)               \n",
    "        elif T > 1:\n",
    "            x1 = r() * X() + (uu_noise - 0.5)\n",
    "        return x1\n",
    "\n",
    "    N = 20\n",
    "    #desired correlation noise\n",
    "    def η():\n",
    "        for n in range(1,N):\n",
    "            x = W() * X()\n",
    "        return x\n",
    "    def η1():\n",
    "        if η() >= 0:\n",
    "            return 1\n",
    "        elif η() < 0:\n",
    "            return 0\n",
    "    return η1()\n",
    "#########################################################################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "03b446a5-370b-4354-ad3d-b0c1f6e3570f",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'x' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-5-b6da2dd453f2>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      9\u001b[0m         \u001b[0mi\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mmath\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msqrt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen_a\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 11\u001b[0;31m \u001b[0mLen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m128\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2000\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     12\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     13\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-5-b6da2dd453f2>\u001b[0m in \u001b[0;36mLen\u001b[0;34m(L, t)\u001b[0m\n\u001b[1;32m      4\u001b[0m     \u001b[0ml\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m         \u001b[0mlen_1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmath\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msqrt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mx_first\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0my_first\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      7\u001b[0m         \u001b[0mlen_a\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mlen_1\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0maverage\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen_1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m         \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0maverage\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen_1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'x' is not defined"
     ]
    }
   ],
   "source": [
    "import math\n",
    "def Len(L,t):\n",
    "    i=1\n",
    "    l=0\n",
    "    for i in range(0,n):\n",
    "        len_1 = math.sqrt((x[i]-x_first)**2+(y[i]-y_first)**2)\n",
    "        len_a += (len_1-np.average(len_1))**2\n",
    "        print(np.average(len_1))\n",
    "        i += 1\n",
    "    return (math.sqrt(len_a/n))\n",
    "Len(128,2000)\n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2502a3bb-6705-4236-a054-aaa76a247370",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "first point\n",
      "[96, 74]\n",
      "-----------\n",
      "iterations\n",
      "2000\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fb241377880>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQEAAAEICAYAAABf40E1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUMElEQVR4nO3df7BcdX3G8fdDIiBgNGkuMSSBxDajoi3i3EHQWikRjUgNnSljmFIvmpo6g4o/ZjRRp7TTapn6Y3TaaicCEpXCRMQSGRXSqHWYVvRG0BJiSCRAApfkUkEUKhL59I9zLpxs9t7du+ec3b37fV4zO7t7zu6ez+69+5zP9+zZPYoIzCxdR/S6ADPrLYeAWeIcAmaJcwiYJc4hYJY4h4BZ4hwCZolzCPSIpKWSQtLsXtfSjunUK+kiSbdUvPyjJH1B0qOS7pJ0VpWPnzKHQIUk3SPp/yT9qnD65x7V8RtJ8xum356/kZd2u6YKrAeWAM8H3gd8VdKcKheQB82VedA8KOl9VT5+v3IIVO9PIuK4wumdPapjD3DBxBVJvw88u0e1lCLpCGAN8LGIeDwibgR+DLy54kX9DbAcOAn4Y+ADklZWvIy+4xDoEkmzJH1C0kOS7gbe2DD/uZKukDQm6X5Jfy9pVj7vIkm35Pd/WNIeSW9oscgvAW8pXB8BvthkmV+UNC7pXkkfyd9wpeqdxmvyGUl78zXvNkmvnuSmJwCLgP8qTPsBcNp0lteGtwB/FxEPR8QO4PPARRUvo+84BLrn7cC5wKnAMPBnDfM3AgeB38tv8zrgLwvzXwHsBOYD/whcIUlTLO/7wBxJL87fnG8Gvtxwm38Cngu8AHgN2ZvgrRXV244fAi8D5gH/BnxF0tFNbvd8YDwifl2Ydl8+/TCS1kl6ZLLTJPeZSxY2Py5M/jHwkmk+p5knInyq6ATcA/wKeKRwens+79vAOwq3fR0QwGxgAfAE8OzC/AuA7+SXLwJ2F+Ydk9/3+VPU8VrgI8A/ACuBLfmyAlgKzMqXeXLhfn8FfLeiem/p4PV7GDilyfRh4KmG1/Vx4OsV/u2W5M/v6MK0s4F7ev1/VfdpRmyZnmHOi4j/aDL9BGBv4fq9hcsnAc8Cxgor9yMabv/gxIWIeDy/3XEtavkS8D1gGQ1DAbKO4siGOu4la7urqLclSe8n6x5OIHsDzsnrajQGPBIR8wr3vYRq19K/ys/nAL8uXP5lhcvoSw6B7hkjW9tMOLFweS/ZmnV+RBysaoERca+kPcA5ZBvWih4CniR7Q99ZqOn+btSbj/8/CKwAtkfEU5IeBpoNcR4AfilpfkQ8lE97Idl2gWaP/SHgQ5MtOyIOC8+IeFjSGHAKWddEfnl7m09pxvI2ge7ZBLxb0uJ8/LluYkZEjAE3A5+UNEfSEZJ+V9JrKljuGuCsiHisODEifpvX9FFJz5F0EtlHbxPbDequ9zlk2xTGgdmS/ppszQuApDMlRb68AK4EPibpGEmvINtQuanZA0fEx+LQT2gOOU1R0xeBj0iaK+lFZNtFrprGc5qRHALV+3rDfgJfy6d/HriJbGPTj4DrG+73FrL2/E6ysfF1wMKyxUTEzyJidJLZ7wIeA+4GbiHbOHdll+q9CfgmcBfZUOPXHDqcWAL8d+H6Zfny9pMF1dsi4tFpLK8dlwI/y+v5T+DjEfGtipfRd5RvADHrK5IuB74SETf1upZB5xAwS1xtwwFJKyXtlLRb0rrW97BBI+nVDUOjp0+9rs2eUUsnkO+cchfZ56z7yHYKuSAi7pzyjmbWdXV9RHga2c4tdwNIuhZYxTMfRR1i/vz5sXTp0ppKMTOAbdu2PRQRQ43T6wqBRRy6pXcf2W6vT5O0FlgLcOKJJzI6OtkGbDOrgqR7m02va5tAsx0+Dhl3RMSGiBiOiOGhocPCycy6pK4Q2Mehe5stJtvry8z6TF0h8ENguaRlko4EVgOba1qWmZVQyzaBiDgo6Z1ke4XNAq6MiIHfB9tsJqrtC0QR8Q3gG3U9vplVw98dMEucQ8AscQ4Bs8Q5BMwS5xAwS5xDwCxxDgGzxDkEzBLnEDBLnEPALHEOAbPEOQTMEucQMEucQ8AscQ4Bs8Q5BMwS5xAwS5xDwCxxDgGzxDkEzBLnEDBLnEPALHEOAbPEOQTMEucQMEtcxyEgaYmk70jaIWm7pEvy6fMkbZG0Kz+fW125Zla1Mp3AQeD9EfFi4HTgYkknA+uArRGxHNiaXzezPtVxCETEWET8KL/8S2AHsAhYBWzMb7YROK9kjWZWo0q2CUhaCpwK3AosiIgxyIICOH6S+6yVNCppdHx8vIoyzKwDpUNA0nHAV4H3RMSj7d4vIjZExHBEDA8NDZUtw8w6VCoEJD2LLACujojr88n7JS3M5y8EDpQr0czqVObTAQFXADsi4lOFWZuBkfzyCHBD5+WZWd1ml7jvq4C/AP5H0u35tA8BlwGbJK0B7gPOL1WhmdWq4xCIiFsATTJ7RaePa2bd5T0GzRLnEDBLnEPALHEOAbPEOQTMEucQMEucQ8AscQ4Bs8Q5BMwS5xAwS5xDwCxxDgGzxDkEzBLnEDBLnEPALHEOAbPEOQTMEucQMEucQ8AscQ4Bs8Q5BMwS5xAwS5xDwCxxDgGzxDkEzBJXxVGJZ0m6TdKN+fV5krZI2pWfzy1fppnVpYpO4BJgR+H6OmBrRCwHtubXzaxPlT00+WLgjcDlhcmrgI355Y3AeWWWYWb1KtsJfBr4APBUYdqCiBgDyM+Pb3ZHSWsljUoaHR8fL1mGmXWq4xCQdC5wICK2dXL/iNgQEcMRMTw0NNRpGWZWUseHJgdeBbxJ0jnA0cAcSV8G9ktaGBFjkhYCB6oo1Mzq0XEnEBHrI2JxRCwFVgPfjogLgc3ASH6zEeCG0lWaWW3q2E/gMuBsSbuAs/PrZtanygwHnhYR3wW+m1/+X2BFFY9rZvXzHoNmiXMImCXOIWCWOIeAWeIcAmaJcwiYJc4hYJY4h4BZ4hwCZolzCJglziFgljiHgFniHAJmiXMImCXOIWCWOIeAWeIcAmaJcwiYJc4hYJY4h4BZ4hwCZolzCJglziFgljiHgFniHAJmiSsVApKeJ+k6ST+VtEPSGZLmSdoiaVd+PreqYs2semU7gc8A34qIFwGnADuAdcDWiFgObM2vm1mf6jgEJM0B/gi4AiAifhMRjwCrgI35zTYC55Ur0czqVKYTeAEwDnxB0m2SLpd0LLAgIsYA8vPjm91Z0lpJo5JGx8fHS5RhZmWUCYHZwMuBz0XEqcBjTKP1j4gNETEcEcNDQ0MlyjCzMsqEwD5gX0Tcml+/jiwU9ktaCJCfHyhXopnVqeMQiIgHgb2SXphPWgHcCWwGRvJpI8ANpSo0s1rNLnn/dwFXSzoSuBt4K1mwbJK0BrgPOL/kMsysRqVCICJuB4abzFpR5nHNrHu8x6BZ4hwCZolzCJglziFgljiHgFniHAJmiXMImCXOIWCWOIeAWeIcAmaJcwhYV0hCUq/LsCYcAmaJcwhYz7TbHbRzO3cZnXMImCWu7O8JWAIm1rIR0XSNGxHTeqzG2zeb1rhsq49DwNpWfEO2+8ZvfBMXg2TiMaYzLJju8q01DwfMEudOwFoqrrEnNFsr90vr3thp2NTcCZglzp1Awoob5KZae3b68VzxsSZuM1k3MV3N6m18PHcE7XEIJG6yNw6Uf/PUNTxo9SnFVMt1MBzOwwGzxLkTSNhkH9dNmGpeL002pJisvnaHM6l2B+4EzBLnTiBx7az9+qUDmK5iB9PYzczU51SHUp2ApPdK2i7pDknXSDpa0jxJWyTtys/nVlWsmVWv4xCQtAh4NzAcES8FZgGryQ5PvjUilgNbmcbhyq37JnbZnTgN0ri4+Jz8ewaTK7tNYDbwbEmzgWOAB4BVwMZ8/kbgvJLLsBpFxJSftQ+CVs+p8TVITZlDk98PfILsyMNjwC8i4mZgQUSM5bcZA45vdn9JayWNShodHx/vtAwzK6nMcGAu2Vp/GXACcKykC9u9f0RsiIjhiBgeGhrqtAybhqla4pTXhqkPFcoMB14L7ImI8Yh4ErgeeCWwX9JCgPz8QPkyzawuZULgPuB0Sccoi9EVwA5gMzCS32YEuKFciVa1xo2BxVOKUu6CoMR+AhFxq6TrgB8BB4HbgA3AccAmSWvIguL8Kgq1cibb136qN36/7SlYl2a/nJRSKJTaWSgiLgUubZj8BFlXYGYzgPcYTJh/wfdQrb6HMKjdgb87YJY4dwIDKrW1eBWm2kYyqF0AuBMwe1qqwekQMEucQ8CsDYPcJTgEzBLnDYMDbtC/IWjluRMwS5w7gQHT6Y9tWmuDukuxQ2DA+Lf06jNob/4JHg6YJc6dwIDxmr8+gzoccCdgljh3AgPEXUC9Bq0DmOBOYMCk/is5dRnk19QhYJY4DwcGSKufC7PODfKPjLgTMEucQ8CS1GwtPtWafbJ5g/ArzR4ODJhUfiG4rGavz1SvWavbN86fSUMFdwJmiXMnMKDcEXTXTFrzN3InYJY4dwIDqkwH4C5i+mby9wpadgKSrpR0QNIdhWnzJG2RtCs/n1uYt17Sbkk7Jb2+rsLNrBrtDAeuAlY2TFsHbI2I5cDW/DqSTgZWAy/J7/NZSbMqq9amVNWBRQfhY686TbVr9kw8uGvLEIiI7wE/b5i8CtiYX94InFeYfm1EPBERe4DdwGnVlGqtTPxzzsSWdCaZ6k0+E/8GnW4YXBARYwD5+fH59EXA3sLt9uXTDiNpraRRSaPj4+MdlmFmZVX96UCzeGwaiRGxISKGI2J4aGio4jJspq2NBsVMGwpA5yGwX9JCgPz8QD59H7CkcLvFwAOdl2dmdes0BDYDI/nlEeCGwvTVko6StAxYDvygXIlWNXcJnSmO9wfp9Wu5n4Cka4AzgfmS9gGXApcBmyStAe4DzgeIiO2SNgF3AgeBiyPitzXVbh0qfqY901rXXhmkN32jliEQERdMMmvFJLf/KPDRMkWZWfd4j8EB1c4a3l1A+4qv1VRfQ252vId+7yL83QGzxLkTGHD+HkBr032NplrL9/tavxl3AmaJcycwYBrXZu4AWkv9NXInMKAG7bNsq49DwCxxHg4MGG8IrN+gdVjuBMwS5xBIxKCtvaw6Hg4MqHYOpNHukGGmDzGq/o7ETP49wWbcCZglzp1AQjpdc83UDmBCFfUPylq/GXcCZolzJ5CwZmu3mb7Wb6b4PFs9v5m+/aMTDgEbOO18rXeyN3k70wdtaODhgFni3AnYIVoNEfr548XJ1tCT/QhI40d97Q4VBo07AbPEuROwlopryqnWlo1r1OKas+4OYqqf/JrOfaYyaDsJTXAnYJY4dwLW0lTbBIpj6yrXlM06gOKy69jWkOpPsDsEbFqmenO2Yzqf2bezbCvPwwGzxLkTsJbq2hA23fa+k41/ZWqqczn9pGUnIOlKSQck3VGY9nFJP5X0E0lfk/S8wrz1knZL2inp9TXVbWYVaWc4cBWwsmHaFuClEfEHwF3AegBJJwOrgZfk9/mspFmVVWsDqfGjxGa/5T/Ia+JeaxkCEfE94OcN026OiIP51e+THYIcYBVwbUQ8ERF7gN3AaRXWa4np1Zs/peCpYsPg24Bv5pcXAXsL8/bl0w4jaa2kUUmj4+PjFZRhZp0oFQKSPkx2CPKrJyY1uVnTOI2IDRExHBHDQ0NDZcqwAdQva+J+qaNOHX86IGkEOBdYEc+8SvuAJYWbLQYe6Lw8M6tbR52ApJXAB4E3RcTjhVmbgdWSjpK0DFgO/KB8mWZWl5adgKRrgDOB+ZL2AZeSfRpwFLAl/zz1+xHxjojYLmkTcCfZMOHiiPhtXcWbWXnqh/HO8PBwjI6O9roMs4EmaVtEDDdO927DZolzCJglziFgljiHgFniHAJmiXMImCXOIWCWuL7YT0DSOPAY8FCvawHm4zqKXMehZnIdJ0XEYV/U6YsQAJA02mxHBtfhOlxHvXV4OGCWOIeAWeL6KQQ29LqAnOs4lOs41MDV0TfbBMysN/qpEzCzHnAImCWuL0JA0sr8OAW7Ja3r4nKXSPqOpB2Stku6JJ8+T9IWSbvy87ldqGWWpNsk3djDGp4n6br8mBI7JJ3Rozrem/897pB0jaSju1XHJMfZmHTZdR1no5vH++h5COTHJfgX4A3AycAF+fELuuEg8P6IeDFwOnBxvux1wNaIWA5sza/X7RJgR+F6L2r4DPCtiHgRcEpeT1frkLQIeDcwHBEvBWaRHcuiW3VcxeHH2Wi67JqPs9GsjnqO9zHxa6q9OgFnADcVrq8H1veolhuAs4GdwMJ82kJgZ83LXUz2z3UWcGM+rds1zAH2kG8sLkzvdh0TP1s/j+zn724EXtfNOoClwB2tXoPG/1XgJuCMuupomPenwNVV1NHzToBpHKugTpKWAqcCtwILImIMID8/vubFfxr4APBUYVq3a3gBMA58IR+WXC7p2G7XERH3A58A7gPGgF9ExM3drqPBZMvu5f9uR8f7aKYfQqDtYxXUVoB0HPBV4D0R8WiXl30ucCAitnVzuU3MBl4OfC4iTiX7LkfXts9MyMfbq4BlwAnAsZIu7HYdberJ/26Z43000w8h0NNjFUh6FlkAXB0R1+eT90tamM9fCByosYRXAW+SdA9wLXCWpC93uQbI/g77IuLW/Pp1ZKHQ7TpeC+yJiPGIeBK4HnhlD+oommzZXf/fLRzv488j7/3L1tEPIfBDYLmkZZKOJNvAsbkbC1b2e+lXADsi4lOFWZuBkfzyCNm2glpExPqIWBwRS8me+7cj4sJu1pDX8SCwV9IL80kryH46vqt1kA0DTpd0TP73WUG2gbLbdRRNtuyuHmejtuN91LmRZxobQM4h29r5M+DDXVzuH5K1TT8Bbs9P5wC/Q7ahbld+Pq9L9ZzJMxsGu14D8DJgNH89/h2Y26M6/hb4KXAH8CWyY1x0pQ7gGrJtEU+SrWHXTLVs4MP5/+1O4A0117GbbOw/8b/6r1XU4d2GzRLXD8MBM+shh4BZ4hwCZolzCJglziFgljiHgFniHAJmift/rJPgTlAFaPgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import copy\n",
    "import random\n",
    "from math import e\n",
    "#When will we consider that the simulation has ended?\n",
    "#When all the positions in the colony have been occupied?\n",
    "\n",
    "L = 128\n",
    "axes = plt.gca()\n",
    "\n",
    "\n",
    "colony = np.ones((L,L), dtype = int)\n",
    "col_coord = np.array([[(i, j) for j in range(L)] for i in range(L)])\n",
    "\n",
    "#define a set of active sites, avoiding multiple instances of the same active site. This set will initialy be empty:\n",
    "active_coord = set()\n",
    "\n",
    "#Add the first active position:\n",
    "x_first = np.random.randint(0,L)\n",
    "y_first = np.random.randint(0,L)\n",
    "active_coord.add((x_first,y_first))\n",
    "colony[x_first,y_first] = 0 \n",
    "\n",
    "def func_veins(coord, LL):    #计算相邻位置\n",
    "    veins_coord =[]\n",
    "    if coord[1] + 1 <= LL - 1:\n",
    "        veins_coord.append((coord[0],coord[1]+1))\n",
    "    if coord[1] - 1 >= 0:\n",
    "        veins_coord.append((coord[0],coord[1]-1))\n",
    "    if coord[0] - 1 >= 0:\n",
    "        veins_coord.append((coord[0]-1,coord[1]))\n",
    "    if coord[0] + 1 <= LL-1:\n",
    "        veins_coord.append((coord[0]+1,coord[1]))\n",
    "    return veins_coord\n",
    "\n",
    "\n",
    "n = 0 \n",
    "while n < 2000:                                           #Main loop of Eden model\n",
    "    #We must save the unoccupied positions of the neighbors of the active positions in a list.\n",
    "    #To do this we first simply calculate the neighboring positions of the active positions\n",
    "    veins = []                                           #initialize\n",
    "    for i in range(len(active_coord)):\n",
    "        one_active_coord = list(active_coord)[i]\n",
    "        veins_function = func_veins(one_active_coord,L)\n",
    "        veins = veins + veins_function\n",
    "    veins = list(set(tuple(x) for x in veins))\n",
    "    inact_veins = [veins[j] for j in range(len(veins)) if colony[veins[j]] == 1]#delete occupied data grid points\n",
    "    inact_veins = set(inact_veins)\n",
    "    new_active = random.sample(inact_veins,1)[0]         #select a random coordinate of these neighbors \n",
    "    colony[new_active] = 0                              #occupy it\n",
    "    active_coord.add(new_active)\n",
    "    n = n + 1\n",
    "\n",
    "print('first point')\n",
    "print([x_first,y_first])\n",
    "print('-----------')\n",
    "print('iterations')\n",
    "print(n)\n",
    "\n",
    "plt.title('Eden Model_a,θ = 0')\n",
    "plt.imshow(colony, interpolation=\"none\", cmap=\"gray\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0eea07fe-1ce0-4c29-a0f0-06c32d2ffa36",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import copy\n",
    "import random\n",
    "#When will we consider that the simulation has ended?\n",
    "#When all the positions in the colony have been occupied?\n",
    "\n",
    "L = 40\n",
    "axes = plt.gca()\n",
    "axes.set_xlim(0,L)\n",
    "axes.set_ylim(0,L)\n",
    "\n",
    "colony = np.ones((L,L), dtype = int)\n",
    "col_coord = np.array([[(i, j) for j in range(L)] for i in range(L)])\n",
    "\n",
    "#define a set of active sites, avoiding multiple instances of the same active site. This set will initialy be empty:\n",
    "active_coord = set()\n",
    "\n",
    "#Add the first active position:\n",
    "x_first = np.random.randint(0,L)\n",
    "y_first = np.random.randint(0,L)\n",
    "active_coord.add((x_first,y_first))\n",
    "colony[x_first,y_first] = 0 \n",
    "\n",
    "def func_veins(coord, LL):    #计算相邻位置\n",
    "    veins_coord =[]\n",
    "    if coord[1] + 1 <= LL - 1:\n",
    "        veins_coord.append((coord[0],coord[1]+1))\n",
    "    if coord[1] - 1 >= 0:\n",
    "        veins_coord.append((coord[0],coord[1]-1))\n",
    "    if coord[0] - 1 >= 0:\n",
    "        veins_coord.append((coord[0]-1,coord[1]))\n",
    "    if coord[0] + 1 <= LL-1:\n",
    "        veins_coord.append((coord[0]+1,coord[1]))\n",
    "    return veins_coord\n",
    "\n",
    "\n",
    "n = 1 \n",
    "while n < 500:                                           #Main loop of Eden model\n",
    "    #We must save the unoccupied positions of the neighbors of the active positions in a list.\n",
    "    #To do this we first simply calculate the neighboring positions of the active positions\n",
    "    veins = []                                           #initialize\n",
    "    for i in range(len(active_coord)):\n",
    "        one_active_coord = list(active_coord)[i]\n",
    "        veins_function = func_veins(one_active_coord,L)\n",
    "        veins = veins + veins_function\n",
    "    veins = list(set(tuple(x) for x in veins))\n",
    "    inact_veins = [veins[j] for j in range(len(veins)) if colony[veins[j]] == 1]#delete occupied data grid points\n",
    "    inact_veins = set(inact_veins)\n",
    "    new_active = random.sample(inact_veins,1)[0]         #select a random coordinate of these neighbors \n",
    "    colony[new_active] = 0                               #occupy it\n",
    "    active_coord.add(new_active)\n",
    "    n = n + 1\n",
    "\n",
    "print('first point')\n",
    "print([x_first,y_first])\n",
    "print('-----------')\n",
    "print('iterations')\n",
    "print(n)\n",
    "\n",
    "plt.title('Eden Model_a')\n",
    "plt.imshow(colony, interpolation=\"none\", cmap=\"gray\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c55a3b10-8f82-42c1-a3bc-9668123d9514",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAZpUlEQVR4nO3dX5Bk9VUH8O+hJxiXFtBaHKrSzCyUohIUzSxIgn9odk0hUmxZpqxYk1SipVMOCUGDZUG2fMgDpSFTUaxIWRTgy65OEoSgqSQk6JAqH0BpAskCSSS4OzTECjxEbbaUYnN8uN07vT23d7r7nl/fc/p+P1W3dvr29Jlv/7ln79y+fUZUFUREFNcZZQcgIqJi2MiJiIJjIyciCo6NnIgoODZyIqLg5sr4obt379Y9e/ZMdNvXXnsNZ511lm0gY94zes8H+M/oPR/AjBa85Wu1Wq+q6nnbrlDVqS9LS0s6qY2NjYlvOy3eM3rPp+o/o/d8qsxowVs+AE9oTk/loRUiouDYyImIgmMjJyIKjo2ciCg4NnIiouDYyCvk8GFgzx6g1cr+PXy47ETp9e7zGWdU7z5X6XmuulLOI6fpO3wYWFkBjh/PLh87ll0GgOXl8nKlxPtcjftM3COvjIMHtzbunuPHs/Wzivc5M+v3mdjIK2Nzc7z1s4D3eef1NBvYyCtiYWG89bOA93nn9TQb2Mgr4vbbgV27Tl23a1e2flbxPmdm/T4TG3llLC8Dd98NLC5mlxcXs8uz/AZY/30Wqd59Bqpxn4lnrVTK8nK2PPoocPRo2Wmmo3efq6SKz3PVcY+ciCg4NnIiouDYyImIgmMjJyIKjo2ciCg4NnIiouDMGrmI1ETkqyLyOauaVWc+uS/BWLwbbwTm5rLztOfmsssW9Votm3rWD6J5Ptg/htbPs3k+spf3hzwnWQB8GMDfAvjcTt/LP768s0OHVHftUgW2ll27svVFC26srRkUVF1dPTVfb1ldLV5vbW2jcD3rB9E8n9o/htbPs3m+Ad63Z2/5MOSPL1s18QaAfwJwDRu5TcbFxfwNaHGxeMGTG3ihgqq1Wn7GWq14vf5GOWk96wfRPJ/aP4bWz7N5vgHet2dv+YY1csmuK0ZE7gfwpwB+CMAfqer1Od+zAmAFAObn55fW19cn+lmdTgf1er1A2vQsMrZaw69bWipWsNNooN5uFyxon7G/XqPRQbu99RgWvc/bTFDQPB/8P8/m+QZ435695Ws2my1V3bvtirzuPs4C4HoAd3W/vhrcI+ceOffIJ6rJPXJ/vOXDkD1yizc7rwJwg4gcBbAO4BoROWRQt9LMp9glGIvX+8szo66fdj3r+2yeL0XNAPeZEsjr7pMu4B65qtplPHQo25ESyf4t8L7kKQU31taMCmZvevX22mq14m+C9eqtrW2Y1LN+EM3zqf1jaP08m+fr43179pYPKd/sPFmMjVxV/Wf0nk/Vf0bv+VSZ0YK3fMMauekYW1V9FMCjljWJiOj0+MlOIqLg2MiJiIJjIyciCo6NnIgoODZyIqLg2Mgd8z4V75SaViMarTN6zxcho/kYzjRTJCst75zE1AvPI9+Z96l4gzUtpguaZ/SeL0JG8zGcaaZIpuKt32AaHwgadWEj35n3GRyDNU0Gwlhn9J4vQkbzoT9pZtak4q3fDGvkPLTi1IkT463f0ebmeOvLqFm1eilqeq+HBK9tYiP3qlYbb/2OFhbGW19GzarVS1HTez0keG0TG7lX3qfiJalZtXopanqvB05UTCLveEvqhcfIR+N9Kl5/TbMRjdYZveeLkNF8DGeaKZIpeOs34Jud0+M9o/d8qv4zes+nyowWvOUb1sh5aIWIKDg2ciKi4NjIiYiCYyMnIgqOjZyIKDg2ciKi4Ao3chG5QEQ2RORZEXlGRG62CEZERKOx2CN/A8AtqnoJgCsBfEBELjGoG475hFPv400jZPSeL0LGBGNsrTOaj3yOJu/k8iILgIcA/MrpvmcWPxDUP+2zN9Gt0LRP7+NNI2T0ni9CxgRjbK0zmo987uOt32Aan+wEsAfAJoCzT/d9s9jI+6d99o/mnHjap/fxphEyes8XIWOCMbbWGc1HPvfx1m+GNXLJritOROoAvgLgdlV9IOf6FQArADA/P7+0vr4+0c/pdDqo1+tFoibRam193Wh00G5vZVxaKlhw0CQF++p1Gg3U2+1i9QZqbuMho/d8AzW38ZDROt9ATYuMKSL2eOs3zWazpap7t12R193HXQC8CcDDAD48yvdzj3zMgh731CJk9J4vQkbukRcvYgipZq2IiAC4F8BzqvqJovWiMp/2GWAcqfuM3uulqOm9XoKaHIuL4nvkAH4BgAL4GoCnust1p7vNLO6Rq25N+1xb27CZ9ul9vGmEjN7zRciYYIytdUbzkc9d3voNOMZ2erxn9J5P1X9G7/lUmdGCt3zDGjk/2UlEFBwbORFRcGzkRETBsZETEQXHRk5EFBwbORFRcGzkhjj90GFG7/kiZAww/dA6Y4qnOam8cxJTL7N4HjmnHzrM6D1fhIwBph9aZzTflg2BHwhKi7NWHGb0ni9CxgCzVqwzmm/LhoY1ch5aMbK5Od76qRc0D5igZtXqpajpvV6Kmsb1Utzl1NjIjSwsjLd+6gXNAyaoWbV6KWp6r5eipnG9FHc5NTZyI5x+aFCzavVS1PReL0VN43op7nJyecdbUi+zeIxcldMPXWb0ni9CxgDTD60zmm/LRsA3O6fHe0bv+VT9Z/SeT5UZLXjLN6yR89AKEVFwbORERMGxkRMRBcdGTkQUHBs5EVFwbORERMGZNHIRuVZEvikiz4vIrRY1I+L0Q4cZveeLkJHTD2d/+iGAGoBvA7gIwJkAngZwyeluM4vnkXP6ocOM3vNFyMjph9WYfgjg7QAe7rt8G4DbTnebWWzknH7oMKP3fBEycvphiOmHkl03ORF5F4BrVfV3u5ffC+DnVfWDA9+3AmAFAObn55fW19cn+nmdTgf1er1Q5hRara2vG40O2u2tjEtLBQsOmqRgX71Oo4F6u12s3kDNbTxk9J5voOY2HjJa5xuo6TGj+bZsqNlstlR177Yr8rr7OAuAdwG4p+/yewF88nS34R75mAU97qlFyOg9X4SM3CMPsUdu8WbnSwAu6Lvc6K6rFE4/NKhZtXopanqvl6Impx+a7JHPAXgBwIXYerPzrae7zSzukaty+qHLjN7zRcjI6Ycu3uhUHb5HXriRZ7VxHYBvITt75eBO3z+rjbzHe0bv+VT9Z/SeT5UZLXjLN6yRzxnt1X8ewOctahER0Xj4yU4iouDYyImIgmMjJyIKjo2ciCg4NnIiouDYyA1x+qHDjN7zRciYYPrhjTcCc3NZxLm57LKziLHknZOYepnF88g5/dBhRu/5ImRMMP1wdTX/I/Crq24inuSt3yDlB4LGXWaxkXPWisOM3vNFyJhg1kqtlr+t1GpuIp7krd8Ma+Q8tGJkc3O89VMvaB4wQc2q1UtR03s9ACdOjLd+JymelmjYyI0sLIy3fuoFzQMmqFm1eilqeq8HoFYbb/1OUjwt0bCRG+H0Q4OaVauXoqb3egBWVsZbv5OQ0wqt5R1vSb3M4jFyVU4/dJnRe74IGRNMP1xdzY6Jr61taK02+RudCSOqqr9+A77ZOT3eM3rPp+o/o/d8qsxowVu+YY2ch1aIiIJjIyciCo6NnIgoODZyIqLg2MiJiIJjIyciCq5QIxeRj4vIN0TkayLyoIica5QruRTT0qwnurmfihcho/d8ETIGmH5YeXnnJI66AHgngLnu1x8D8LFRblf2eeQppqVZT3RzPxUvQkbv+SJkDDD9MKUo55GbfcgHwK8DODzK95bdyFNMS7Oe6OZ+Kl6EjN7zRcgYYPphSlEauWTXFSci/wjgU6p6aMj1KwBWAGB+fn5pfX19op/T6XRQr9cnzglkv84Ns7RUvGaj0UG7vZVxoprWIfvqdRoN1NvtYvUGam7jIaP3fAM1t/GQMcHGYr6tJGTRbyw1m82Wqu7ddkVed+9fADwC4EjOcqDvew4CeBDI/mPYaeEeeQkhuTfpL1+EjNwjLzvCKTDpR/RVdb+qXpqzPAQAIvJ+ANcDWO7+IPdSTEuznugWYYqd+4ze66Wo6b0eEmwrVOwYOYBrATwL4Lxxblf2Hrlqmmlp1hPd3E/Fi5DRe74IGQNMP0wlyh550Ub+PIAXATzVXf56lNt5aOQpec/oPZ+q/4ze86kyowVv+YY18rmCe/M/VuT2RERUHD/ZSUQUHBs5EVFwbORERMGxkRMRBcdGTkQUXJhGbj50LsX4Q+9T5zi5z1++CBkDbCspIoaSd05i6mXc88j7B7D1PtJbaABbivGH3qfOcXKfv3wRMgbYVlJE7IlyHnmIRt4/7qF/NsPE4x5SDFvxPuOCc0L85YuQMcC2kiJiT5RGHuLQyubmeOunXzBBTe/1UtSsWr0UNb3XS1AzRcRoQjTyhYXx1k+/YIKa3uulqFm1eilqeq+XoGaKiNGEaOTmA9g4xa54vRQ1q1YvRU3v9RLUTBExnLzjLamXSYZm9Qawra1t2AxgSzH+0PvUOU7u85cvQsYA20qKiKpxjpGHaeQ93h7YPN4zes+n6j+j93yqzGjBW75hjTzEoRUiIhqOjZyIKDg2ciKi4NjIiYiCYyMnIgqOjZyIKDiTRi4it4iIishui3p5OP3QsB4n9/nJFyFjgG2F0w8LnhMO4AIADwM4BmD3KLfh9MMSMnJyn798ETIG2FY4/dDgA0EA7gdwGYCjqRo5px8aZOTkPn/5ImQMsK1w+qFCsusmIyIHAFyjqjeLyFEAe1X11SHfuwJgBQDm5+eX1tfXR/45rdbW141GB+12/eTlpaUJgvcXHDRRwVNrdhoN1NvtYjWtM1rnG6i5jYeM3vMN1NzGQ8YA20qKiD2dTgf1en3nb5ySZrPZUtW9267I6+79C4BHABzJWQ4AeBzAOd3vOwrukfvdE+LepL98ETIG2Fa4Rz7CR/RVdb+qXjq4AHgBwIUAnu7ujTcAPCki5xf9X2cQpx86rJeiZtXqpajpvV6Cmpx+iJ33yEddkHCPXJXTD11OxYuQ0Xu+CBkDbCucfhikkfd4e2DzeM/oPZ+q/4ze86kyowVv+YY18jnDPfs9VrWIiGh0/GQnEVFwbORERMGxkRMRBcdGTkQUHBs5EVFwlW3kKaal3XgjMDeXfWR4bi677Cnk/v2ASJZPJLtcWNUm91Vw+mHvddNbLF431tuK84cwvbxzElMvZZ9HnmJa2urqVq3+MQKrqz5C7tuXn2/fvgnzJcjofnJfBacf9r9u+pcirxvrbSXhQ2gzbdUQUn8gaJyl7EaeYjZDrZb/4qzVfITsL9GfD5gwX4KM7ueEVHDWSl60oq8b620l4UNoM9vJ0LBGXslDK5ub460fxYkT463fUYqQ1qwzVq1eipoBXjfW20oFH8JtKtnIFxbGWz+KWm289TtKEdKadcaq1UtRM8DrxnpbqeBDuE0lG3mKaWkrK+Ot35FxyH37xls/Eu+T9rzXS1EzwOvGeltx/hBOR97xltRL2cfIVdNMS1tdzY7zra1taK1W4I3ORCF7b1z1jvsVeqMzUUb3k/sqOP1w8A1Pi9eN9baS6CG0m7ZqBHyzc3q8Z/SeT9V/Ru/5VJnRgrd8wxp5JQ+tEBHNEjZyIqLg2MiJiIJjIyciCo6NnIgoODZyIqLgCjdyEblJRL4hIs+IyB0WoYiIaHSFGrmINAEcAHCZqr4VwJpJqqDMR18az+as4ATWk+NSRYxHCxs+iOZjYr2PAo6QMcHznOJhPCnv5PJRFwCfBrB/3NvN4geCzEdfGs/mTDGa0/v40P5xqf2LxWhhqzG25mNivY8CjpAxwfNsFREpPtkJ4CkAHwXwOICvALh8lNvNYiM3H31pPJszxWhO7+ND+8el9i8Wo4Wtxtiaj4n1Pgo4QsYEz7NVxGGNXLLrhhORRwCcn3PVQQC3A9gA8CEAlwP4FICLNKeoiKwAWAGA+fn5pfX19fF/fQDQ6XRQr9cnum1KrdbW141GB+32VsalpYIFB01Q0DwfzCOaZ7TO11+w02ig3m4XLBggo3lA+M/o+HluNpstVd277Yq87j7qAuCLAJp9l78N4Lydbsc98jELco98Itwjd7i3GyFjwD3yometfBZAEwBE5GIAZwJ4tWDNkMxHXxoXrOAEVvejhYEEY2K9Pykpanqvl6bkqfK6+6gLssZ9CMARAE8CuGaU283iHrlqgtGXxrM5U4zm9D4+tDcutbcnbjVa2HKMrfmYWO+jgCNkTPA8W0QEx9hOj/eM3vOp+s/oPZ8qM1rwlm9YI+cnO4mIgmMjJyIKjo2ciCg4NnIiouDYyImIgmMjJyIKjo3csaTT0mgiESZIWoc0nyBpHzHC8MO08s5JTL3wPPKdpRg6Z5kvNY8ZI0yQtJ7cZz5BUu0fR+9TOC2B55HHcvAgcPz4qeuOH8/WUzlSPCfmNY0L3n33eOtHYX2fvdebBjZypzY3x1tP6aV4TsxrGhc8cWK89aOwvs/e600DG7lTCwvjraf0Ujwn5jWNC9Zq460fhfV99l5vGtjInUo+LY3GFmEQoHVB8wmS8D+sMOK2x0bu1PJydhxycTE7W2BxMbu8vFx2surqf04Am+fE/Hk2DnnXXcDq6tYeeK2WXb7rrgnz2Uc0fwxTPM/J5b0DmnrhWSvl8p5P1X9G7/lUmdGCt3zgWStERLOJjZyIKDg2ciKi4NjIiYiCYyMnIgqOjZyIKLhCjVxEflZEHhORp0TkCRG5wioYUVV4n37IKZz+zRW8/R0APqqqXxCR67qXry6ciqgiDh/OPiXZG9J07NjWpyYn+gCKcUHzfJRE0UMrCuDs7tfnAHi5YD2iSvE+/TDiJMAqkuzDQhPeWOSnADwMQJD9p/AOVT025HtXAKwAwPz8/NL6+vpEP7PT6aBer08WeEq8Z/SeD/Cf0SpfqzX8uqWlYgU7jQbq7Xahgub5BlTlebbSbDZbqrp32xV5H/fsXwA8AuBIznIAwF8C+I3u9/0mgEd2qqf8iH7pvOdT9Z/RKt/iYv4fblhcLF7w5B+WKFDQPN+AqjzPVjDpR/RVdb+qXpqzPATgfQAe6H7rZwDwzU6iMXiffhhxEmAVFT1G/jKAX+5+fQ2Afy9Yj6hSvE8/5BTOGIqetfJ7AO4UkTkA/4vuMXAiGt3ysnFj7BV89FHg6FGzcuRXoUauqv8CwOAtDyIimhQ/2UlEFBwbORFRcGzkRETBsZETEQXHRk5EFBwbOdGMMR5+SAEUPY+ciBzhtMJq4h450QzhtMJqYiMnmiGbm+Otp9nARk40QxYWxltPs4GNnGiGcFphNbGRE80Q4+GHFATPWiGaMcbDDykA7pETEQXHRk5EFBwbORFRcGzkRETBsZETEQUnqjr9HyryCoBjE958N4BXDeOk4D2j93yA/4ze8wHMaMFbvkVVPW9wZSmNvAgReUJV95ad43S8Z/SeD/Cf0Xs+gBkteM/Xw0MrRETBsZETEQUXsZHfXXaAEXjP6D0f4D+j93wAM1rwng9AwGPkRER0qoh75ERE1IeNnIgouLCNXERuEpFviMgzInJH2XmGEZFbRERFZHfZWfqJyMe7j9/XRORBETm37EwAICLXisg3ReR5Ebm17DyDROQCEdkQkWe7r72by86UR0RqIvJVEflc2VnyiMi5InJ/9zX4nIi8vexMg0TkD7vP8RER+TsReXPZmYYJ2chFpAngAIDLVPWtANZKjpRLRC4A8E4AHv/Q1pcBXKqqPwPgWwBuKzkPRKQG4K8A/CqASwD8lohcUm6qbd4AcIuqXgLgSgAfcJgRAG4G8FzZIU7jTgBfVNWfBHAZnGUVkbcA+BCAvap6KYAagHeXm2q4kI0cwCqAP1PV/wMAVf1uyXmG+XMAfwzA3TvKqvolVX2je/ExAI0y83RdAeB5VX1BVV8HsI7sP2w3VPU7qvpk9+v/QdaA3lJuqlOJSAPArwG4p+wseUTkHAC/BOBeAFDV11X1e6WGyjcH4AdFZA7ALgAvl5xnqKiN/GIAvygij4vIV0Tk8rIDDRKRAwBeUtWny84ygt8B8IWyQyBriC/2XW7DWZPsJyJ7APwcgMdLjjLoL5DtQHy/5BzDXAjgFQB/0z38c4+InFV2qH6q+hKy3/Q3AXwHwH+p6pfKTTWc278QJCKPADg/56qDyHL/CLJfbS8H8GkRuUinfC7lDhk/guywSmlOl09VH+p+z0FkhwsOTzNbdCJSB/D3AP5AVf+77Dw9InI9gO+qaktEri45zjBzAN4G4CZVfVxE7gRwK4A/KTfWFhH5YWS/DV4I4HsAPiMi71HVQ6UGG8JtI1fV/cOuE5FVAA90G/e/isj3kQ23eWVa+YDhGUXkp5G9AJ4WESA7bPGkiFyhqv9Zdr4eEXk/gOsB7Jv2f4JDvATggr7Lje46V0TkTcia+GFVfaDsPAOuAnCDiFwH4M0AzhaRQ6r6npJz9WsDaKtq7zeZ+5E1ck/2A/gPVX0FAETkAQDvAOCykUc9tPJZAE0AEJGLAZwJRxPKVPXrqvqjqrpHVfcge+G+bZpNfCcici2yX79vUNXjZefp+jcAPy4iF4rImcjeXPqHkjOdQrL/me8F8JyqfqLsPINU9TZVbXRfd+8G8M/Omji628GLIvIT3VX7ADxbYqQ8mwCuFJFd3ed8H5y9IdvP7R75Du4DcJ+IHAHwOoD3OdmjjOSTAH4AwJe7vzU8pqq/X2YgVX1DRD4I4GFkZwncp6rPlJkpx1UA3gvg6yLyVHfdR1T18+VFCukmAIe7/2G/AOC3S85ziu4hn/sBPIns0ONX4fjj+vyIPhFRcFEPrRARURcbORFRcGzkRETBsZETEQXHRk5EFBwbORFRcGzkRETB/T8QpaY/2kVa/gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAYmUlEQVR4nO3dfbBc9V3H8fcntw01pK15uNII3HvpGMemTrV1pbbT1oeEKeBIcLQKXjpB6dySiIPT0Wmc6zAOTkasOj6MKL2DnUnLWtriA7GmVkirVSvIpVIgrZA0JiE0hUtoS2mUlObrH+esnCy7e3f3nH08n9fMnT0Pv93z42T53N/9nidFBGZmNv5WDLoDZmbWHw58M7OScOCbmZWEA9/MrCQc+GZmJeHANzMriUICX9LFkh6RdFDSzgbrr5X0kKQHJP2rpE1FbNfMzNqnvOfhS5oAHgUuAo4B9wFXRsQXMm1eERHPpNOXATsi4uJcGzYzs44UMcK/EDgYEYci4hRwO7A126AW9qmzAV/tZWbWZy8p4DPOBR7LzB8D3ljfSNIvA+8BVgI/sdyHrl+/PmZmZgronplZedx///1PRcRko3VFBH5bIuJm4GZJvwD8JrCtvo2kOWAOYGpqisXFxX51z8xsLEg60mxdESWdx4HzM/PnpcuauR24vNGKiFiIiEpEVCYnG/6CMjOzLhUR+PcBGyVdIGklcAWwJ9tA0sbM7E8CBwrYrpmZdSB3SScinpd0HfBJYAL4QETsl3QjsBgRe4DrJG0BvgV8lQblHDMz661CavgRsRfYW7fshsz09UVsx8zMuucrbc3MSsKBb2Y2LKpVmJmBFSuS12q10I/v22mZZmbWwo4dcMstULv7wZEjMDeXTM/OFrIJj/DNzAatWoU///MXwr7m5EmYny9sMw58M7NBu/ba5uuOHi1sMw58M7NB2rEDnn22+fqpqcI25cA3MxuUWimnlV27CtucA9/MbFBalXIAzj67sAO24MA3MxuMLVtal3IA3v/+QjfpwDcz67ctW2DfvtZttm8vdHQPDnwzs/7asWP5sD/7bPizPyt8077wysysX9oZ2UPhpZwaj/DNzHqtWoWzzmov7HtQyqnxCN/MrJfaHdUDbN7ck1JOjUf4Zma90MmoHpKwv/vunnbJgW9mVqRqFVavhquuglOn2nvP9u09D3twScfMrDidlG9qelzGyfII38wsr07LNzV9KONkOfDNzLrVTfkG4GUvg9tu62vYQ0GBL+liSY9IOihpZ4P175H0BUkPStonabqI7ZqZDcyOHUnQf/Obnb1v82b4n//p2amXreQOfEkTwM3AJcAm4EpJm+qa/SdQiYjXAXcA78u7XTOzgaiN6pe7y2W91asHMqrPKmKEfyFwMCIORcQp4HZga7ZBRHw6Ik6ms/cA5xWwXTOz/upmVF8r33zjGwMZ1WcVcZbOucBjmfljwBtbtL8G+EQB2zUz658dOzof1ff5oOxy+nrQVtJVQAX4vSbr5yQtSlpcWlrqZ9fMzBrrpoQzBOWbRooI/MeB8zPz56XLziBpCzAPXBYRzzX6oIhYiIhKRFQmJycL6JqZWZeyZ+C0W8KpBf0QlG8aKSLw7wM2SrpA0krgCmBPtoGk1wPvJwn7JwvYpplZsapVWL8epOSn01r99u1DG/Q1uWv4EfG8pOuATwITwAciYr+kG4HFiNhDUsJZDXxMEsDRiLgs77bNzHKrVuHd7+789MoaCT70oaEO+ppCbq0QEXuBvXXLbshMbyliO2ZmhermQGzWxATs3j0SYQ++0tbMyipv2K9ePVJhD755mpmVTd4SzurVcMstIxX0NR7hm1k51A7KdnM7BIB164b6DJx2eIRvZuNvx45kVB7R2fvWrYM//uORDfh6DnwzG2+d1upHuGSzHJd0zGz8VKswM5OcMtlJ2I/AufR5eIRvZuOl2/LN9u19e/LUoHiEb2bjo1a+6STsawdjxzzswSN8MxsH1Spcfz2cONH+e0owoq/nEb6ZjaZsnf6qqxz2bfAI38xGT7d1egmuvbaUYQ8e4ZvZKKldPNVpnR5gejq5yVlJwx48wjezUdHNqL7kI/p6HuGb2XDrdlS/bl3pR/T1PMI3s+FVrcLcHJw82f57PKpvyiN8MxtO1Sps29ZZ2LtO35JH+GY2fDqp14/ZDc56yYFvZsOjkwuoXLrpmAPfzIZDJ/V6j+q7UkgNX9LFkh6RdFDSzgbr3ybpc5Kel/SzRWzTzMZIu/X6iYnkvjdPPeWw70LuwJc0AdwMXAJsAq6UtKmu2VHgauAv827PzMZMbWT/7W+3bieN3DNkh00RJZ0LgYMRcQhA0u3AVuALtQYRcThdd7qA7ZnZOLn++uVH9rV6vcM+lyJKOucCj2Xmj6XLzMxaq1aXP0DrC6gKM1QHbSXNAXMAU1NTA+6NmfVUrW7fzMSESzgFK2KE/zhwfmb+vHRZxyJiISIqEVGZnJwsoGtmNnRqt0q46qrWdXuHfeGKCPz7gI2SLpC0ErgC2FPA55rZuKkdoG2njOOwL1zuwI+I54HrgE8CXwQ+GhH7Jd0o6TIAST8s6RjwDuD9kvbn3a6ZjZh2T71ctSo5x94KV0gNPyL2Anvrlt2Qmb6PpNRjZmXTydWzExOwsODRfY8M1UFbMxsznVw9u2qVw77HfLdMM+uNTu52uW6dw74PPMI3s+K1e/WsT73sK4/wzaw41SrMzCSnXLZzcNZh31ce4ZtZfp0cmAXf7XJAHPhmlk8nB2Zdwhkol3TMLJ/5+fbPwnHYD5QD38y6V63CkSPLt5ue9lk4Q8AlHTPrTq2U04rPrR8qHuGbWWfaPRPH59YPHY/wzax97R6gve02B/0Q8gjfzNrXzgHa6WmH/ZBy4JtZe9o5QLtqFeza1Z/+WMcc+Ga2vHYO0PpMnKHnGr6ZNVetJmWcViN7n4kzMhz4ZtZYuwdoHfYjw4FvZi9Wu7Xxcne79AHakeIavpmdqd1bG/sA7chx4JtZopNbG/sA7UgqJPAlXSzpEUkHJe1ssP4sSR9J198raaaI7ZpZQWqj+nZOu7ztNjh82GE/gnIHvqQJ4GbgEmATcKWkTXXNrgG+GhHfA/wh8Lt5t2tmBehkVO8HjI+8Ikb4FwIHI+JQRJwCbge21rXZCuxOp+8ANktSAds2s07VQl6Cd76zvbtd+tbGY6GIwD8XeCwzfyxd1rBNRDwPfB1YV8C2zawT9aWbiOXf43r92Biqg7aS5iQtSlpcWloadHfMxkcnpZsa1+vHThGB/zhwfmb+vHRZwzaSXgK8EnjRwy8jYiEiKhFRmZycLKBrZtb2Adksj+rHUhGBfx+wUdIFklYCVwB76trsAbal0z8LfCqinb8lzSy3dh9BCB7Vj7ncgZ/W5K8DPgl8EfhoROyXdKOky9JmfwGsk3QQeA/wolM3zaxAtRLOihXLj+xr5094VD/2Crm1QkTsBfbWLbshM/2/wDuK2JaZLaPde+BAEvK7djnkS8L30jEbF+3c2bLGd7gsJQe+2Thod1QvwdSUR/Ul5cA3G1W1Ef3Ro0mtvp07Wx4+3Jeu2XBy4JuNovoRve9saW0YqguvzKyF7Jk327a1f6qlz76xlEf4ZqOg0xE9+MCsvYhH+GajoN2LpyYmkgOzHtVbAw58s2HVycVT8MIdLU+f9pWy1pBLOmbDqN3TLCcmkoD3qZbWBge+2TBqp4TjGr11yCUds2F09Gjzda7RW5cc+GbDIFuvn5mBtWsbt5uedo3euuaSjtmg1dfrjxyBl74UVq6EU6deaOeLpywnj/DNBq1Rvf5b34KXvzwZ0buEYwXxCN9s0JrV659+Gp56qr99sbHmEb7ZoE1NdbbcrEsOfLN+qz9Ae+mlSX0+y/V66wEHvlk/ZR8oHpG87t6d3AzN9XrrMdfwzfqp0QHakydh717fq956LtcIX9JaSXdJOpC+rmnS7h8kfU3Sx/Nsz2zkNTtA2+pCK7OC5C3p7AT2RcRGYF8638jvAe/MuS2z0dPuBVU+QGt9kDfwtwK70+ndwOWNGkXEPuAbObdlNloa1eufeSa5oCrLB2itT/IG/jkRcTyd/gpwTs7PMxsfvqDKhsyyB20l3Q28qsGq+exMRISkyNMZSXPAHMCU/8S1UZJ9oHjtVsW+oMqGzLKBHxFbmq2T9ISkDRFxXNIG4Mk8nYmIBWABoFKp5PrlYdY3je6FMzeX1OtPnHhxew9mbEDylnT2ANvS6W3AnTk/z2z0NDvVEnxBlQ2VvIF/E3CRpAPAlnQeSRVJt9YaSfoX4GPAZknHJL0953bNhker0s3Cguv1NjRyXXgVESeAzQ2WLwLvysy/Nc92zIZGo1r91FTjZ85OTSXh7oC3IeFbK5i1q9FplnNzvheOjQwHvlm7Wt0WwaUbGwGKGM6TYSqVSiwuLg66G2YvWLEiGdnXk5LHDpoNAUn3R0Sl0TqP8M3a5fvW24hz4Ju1a9cu1+ptpDnwzdo1O+tavY00B75Zvfo7XFarL6ybnU3uW3/6dPLqsLcR4gegmGU1u00CONxt5HmEb5bV7NTL+fnG7c1GiAPfLMtPpLIx5sA3y/KplzbGHPhmWT710saYA98sy6de2hhz4Fv5tDrtEnzqpY0tn5Zp5eLTLq3EPMK3cvFpl1ZiDnwrF592aSXmwLdy8WmXVmIOfCsXn3ZpJZYr8CWtlXSXpAPp65oGbX5Q0r9L2i/pQUk/n2ebZrn4tEsrsVxPvJL0PuDpiLhJ0k5gTUS8t67N9wIREQckfTdwP/CaiPhaq8/2E6/MzDrXyydebQV2p9O7gcvrG0TEoxFxIJ3+MvAkMJlzu2Zm1qG8gX9ORBxPp78CnNOqsaQLgZXAl5qsn5O0KGlxaWkpZ9fMzCxr2cCXdLekhxv8bM22i6Q21LQ+JGkD8CHgFyOi4ROfI2IhIioRUZmc9B8B1oXlrqI1K7Flr7SNiC3N1kl6QtKGiDieBvqTTdq9Avh7YD4i7um6t2at+Cpas5bylnT2ANvS6W3AnfUNJK0E/gb4YETckXN7Zs35KlqzlvIG/k3ARZIOAFvSeSRVJN2atvk54G3A1ZIeSH9+MOd2zV7MV9GatZTr5mkRcQLY3GD5IvCudPo24LY82zFry9RUUsZptNzMfKWtjRFfRWvWkgPfxoevojVryffDt/EyO+uAN2vCI3wzs5Jw4JuZlYQD38ysJBz4ZmYl4cA3MysJB76ZWUk48M3MSsKBb2ZWEg58M7OScOCbmZWEA99Gg59kZZab76Vjw89PsjIrhEf4Nvz8JCuzQjjwbfj5SVZmhXDg2/Br9sQqP8nKrCMOfBt+fpKVWSFyBb6ktZLuknQgfV3ToM20pM+lDy/fL+naPNu0EvKTrMwKoYjo/s3S+4CnI+ImSTuBNRHx3ro2K9PtPCdpNfAw8OaI+HKrz65UKrG4uNh138zMykjS/RFRabQub0lnK7A7nd4NXF7fICJORcRz6exZBWzTzMy6kDd8z4mI4+n0V4BzGjWSdL6kB4HHgN9tNrqXNCdpUdLi0tJSzq6ZmVnWsoEv6W5JDzf42ZptF0ltqGF9KCIei4jXAd8DbJPU8BdDRCxERCUiKpOTk13859iyfMWqWWktG/gRsSUivr/Bz53AE5I2AKSvTy7zWV8mqeG/tYjOj41+hXDtitUjRyDihStWHfpmpZC3pLMH2JZObwPurG8g6TxJ35FOrwHeAjySc7vjo58hPExXrPovDbO+yxv4NwEXSToAbEnnkVSRdGva5jXAvZI+D/wz8PsR8VDO7Y6PfoZw0Vesdhva/kvDbCBynZbZS6U5LXPFiiT06klw+nSx25qZScK13vQ0HD7c2WfV39AMkouh2jk/vsh+mNkZenlapuXVz9sGFHnFap6/THxvHLOBcOAPWj9vG1DkFat5Qtv3xjEbCAf+oPX7tgGzs0nZ5PTp5LXb7eQJbd8bx2wgHPjDoKgQ7qc8oe1745gNhJ94Zd2phfP8fFLGmZpKwr7d0J6ddcCb9ZkD37rn0DYbKS7pmJmVhAPfzKwkHPhmZiXhwDczKwkHvplZSTjwzcxKwoFvZlYSDnwzs5Jw4JuZlYQD38ysJBz4ZmYl4cA3MysJB76ZWUnkCnxJayXdJelA+rqmRdtXSDom6U/zbNPMzLqTd4S/E9gXERuBfel8M78NfCbn9szMrEt5A38rsDud3g1c3qiRpB8CzgH+Mef2zMysS3kD/5yIOJ5Of4Uk1M8gaQXwB8CvLfdhkuYkLUpaXFpaytk1MzPLWvaJV5LuBl7VYNV8diYiQlI0aLcD2BsRxyS13FZELAALAJVKpdFnmZlZl5YN/IjY0mydpCckbYiI45I2AE82aPYm4K2SdgCrgZWSno2IVvV+MzMrWN5n2u4BtgE3pa931jeIiP9/6Kmkq4GKw97MrP/y1vBvAi6SdADYks4jqSLp1rydMzOz4ihiOEvllUolFhcXB90NM7ORIun+iKg0Wucrbc3MSsKBb2ZWEg58M7OSGL/Ar1ZhZgZWrEheq9VB98jMbCjkPS1zuFSrMDcHJ08m80eOJPMAs7PN32dmVgLjNcKfn38h7GtOnkyWm5mV3HgF/tGjnS03MyuR8Qr8qanOlpuZlch4Bf6uXbBq1ZnLVq1KlpuZldx4Bf7sLCwswPQ0SMnrwoIP2JqZMW5n6UAS7g54M7MXGa8RvpmZNeXANzMrCQe+mVlJOPDNzErCgW9mVhJD+wAUSUvAkRwfsR54qqDuFMn96oz71Rn3qzPj2K/piJhstGJoAz8vSYvNnvoySO5XZ9yvzrhfnSlbv1zSMTMrCQe+mVlJjHPgLwy6A024X51xvzrjfnWmVP0a2xq+mZmdaZxH+GZmljFygS/pYkmPSDooaWeD9WdJ+ki6/l5JM5l1v5Euf0TS2/vcr/dI+oKkByXtkzSdWfdtSQ+kP3v63K+rJS1ltv+uzLptkg6kP9v63K8/zPTpUUlfy6zr5f76gKQnJT3cZL0k/Una7wclvSGzrpf7a7l+zab9eUjSZyX9QGbd4XT5A5IW+9yvH5P09cy/1w2ZdS2/Az3u169n+vRw+p1am67r5f46X9Kn0yzYL+n6Bm169x2LiJH5ASaALwGvBlYCnwc21bXZAdySTl8BfCSd3pS2Pwu4IP2ciT7268eBVen09lq/0vlnB7i/rgb+tMF71wKH0tc16fSafvWrrv2vAB/o9f5KP/ttwBuAh5usvxT4BCDgR4B7e72/2uzXm2vbAy6p9SudPwysH9D++jHg43m/A0X3q67tTwGf6tP+2gC8IZ1+OfBog/8ne/YdG7UR/oXAwYg4FBGngNuBrXVttgK70+k7gM2SlC6/PSKei4j/Bg6mn9eXfkXEpyOi9sDde4DzCtp2rn618Hbgroh4OiK+CtwFXDygfl0JfLigbbcUEZ8Bnm7RZCvwwUjcA3ynpA30dn8t26+I+Gy6Xejf96ud/dVMnu9m0f3q5/freER8Lp3+BvBF4Ny6Zj37jo1a4J8LPJaZP8aLd9b/t4mI54GvA+vafG8v+5V1Dclv8JqXSVqUdI+kywvqUyf9+pn0T8c7JJ3f4Xt72S/S0tcFwKcyi3u1v9rRrO+93F+dqv9+BfCPku6XNDeA/rxJ0uclfULSa9NlQ7G/JK0iCc2/yizuy/5SUm5+PXBv3aqefcfG7wEoQ07SVUAF+NHM4umIeFzSq4FPSXooIr7Upy79HfDhiHhO0rtJ/jr6iT5tux1XAHdExLczywa5v4aapB8nCfy3ZBa/Jd1f3wXcJem/0hFwP3yO5N/rWUmXAn8LbOzTttvxU8C/RUT2r4Ge7y9Jq0l+yfxqRDxT5Ge3Mmoj/MeB8zPz56XLGraR9BLglcCJNt/by34haQswD1wWEc/VlkfE4+nrIeCfSH7r96VfEXEi05dbgR9q97297FfGFdT9ud3D/dWOZn3v5f5qi6TXkfwbbo2IE7Xlmf31JPA3FFfKXFZEPBMRz6bTe4GXSlrPEOyvVKvvV0/2l6SXkoR9NSL+ukGT3n3HenFgolc/JH+RHCL5E792oOe1dW1+mTMP2n40nX4tZx60PURxB23b6dfrSQ5SbaxbvgY4K51eDxygoINXbfZrQ2b6p4F74oUDRP+d9m9NOr22X/1K230fyQE09WN/ZbYxQ/ODkD/JmQfU/qPX+6vNfk2RHJd6c93ys4GXZ6Y/C1zcx369qvbvRxKcR9N919Z3oFf9Ste/kqTOf3a/9lf63/5B4I9atOnZd6ywnduvH5Ij2I+ShOd8uuxGklEzwMuAj6Vf/v8AXp1573z6vkeAS/rcr7uBJ4AH0p896fI3Aw+lX/iHgGv63K/fAfan2/808H2Z9/5Suh8PAr/Yz36l878F3FT3vl7vrw8Dx4FvkdRIrwGuBa5N1wu4Oe33Q0ClT/truX7dCnw18/1aTJe/Ot1Xn0//nef73K/rMt+ve8j8Qmr0HehXv9I2V5OcyJF9X6/311tIjhE8mPm3urRf3zFfaWtmVhKjVsM3M7MuOfDNzErCgW9mVhIOfDOzknDgm5mVhAPfzKwkHPhmZiXhwDczK4n/A9xriPhUsWizAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import random as rd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "\n",
    "def d(x,y):\n",
    "    r = pow(x**2 + y**2, 0.5)\n",
    "    return r\n",
    "\n",
    "plt.figure()\n",
    "\n",
    "t=100     # 设定时间步长\n",
    "\n",
    "point = np.array([[0,0]])\n",
    "step = np.array([[0,1],[0,-1],[1,0],[-1,0]])\n",
    "edge = step\n",
    "\n",
    "dis=[]\n",
    "dis = np.array(dis)\n",
    "wt = []\n",
    "for i in range(0,t):\n",
    "    growsite = rd.randint(0,len(edge)-1)\n",
    "    point = np.append(point,np.array([edge[growsite]]),axis=0)\n",
    "    edge = np.append(edge,np.array([edge[growsite]]) + step,axis=0)\n",
    "    j=0\n",
    "    for m in edge:\n",
    "        for n in point:\n",
    "            if m[0]==n[0] and m[1]==n[1]:\n",
    "                edge = np.delete(edge,j,axis=0)\n",
    "                j=j-1\n",
    "        j=j+1\n",
    "    # print(i)\n",
    "    for j in edge:\n",
    "        dis = np.append(dis,d(j[0], j[1]))\n",
    "    w = np.std(dis)\n",
    "    wt.append(w)\n",
    "wt = np.array(wt)\n",
    "\n",
    "plt.scatter(edge[:,0],edge[:,1],c='b')\n",
    "plt.scatter(point[:,0],point[:,1],c='r')\n",
    "plt.grid(True)\n",
    "plt.figure()\n",
    "plt.scatter(np.log10(np.arange(1,len(wt)+1)),np.log10(np.array(wt)),c='r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bb6e568c-1739-4631-b6d3-a077bb61b1d2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "β = 0.3756537272820365\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAN4UlEQVR4nO3cX4jdd5nH8ffHpKWKxYodRfIHsxLHzUUL/c+iu2PLrklvglDYtmLZooSyVvayZS/0ojcrsuCKrWEoofTGXKxF4xItwnKsULvWhTZtWhJmU7aZTaHUipIKW9I8e3HOco7jpPPLzG9m0vN9v2BgzjnfmTx5SN45+c3MSVUhSZp+79vsASRJG8PgS1IjDL4kNcLgS1IjDL4kNcLgS1IjVgx+kkNJXk/y4gUeT5LvJFlIcizJdf2PKUlaqy7P8B8D9r7L4/uA3aO3A8D31j6WJKlvKwa/qp4C3nyXI/uBx2voGeCqJB/va0BJUj+29vA5tgGnJ24vju57benBJAcY/i+AK6644vqdO3f28Mu/950/f573vc8vp4C7mOQuxtzF2MmTJ9+oqpnVfGwfwc8y9y37eg1VNQ/MA8zOztaJEyd6+OXf+waDAXNzc5s9xiXBXYy5izF3MZbkv1f7sX38k7kI7Ji4vR0408PnlST1qI/gHwHuGX23zi3A76rqTy7nSJI214qXdJJ8H5gDrk6yCHwDuAygqg4CR4HbgQXgD8C96zWsJGn1Vgx+Vd21wuMFfLW3iSRJ68Ive0tSIwy+JDXC4EtSIwy+JDXC4EtSIwy+JDXC4EtSIwy+JDXC4EtSIwy+JDXC4EtSIwy+JDXC4EtSIwy+JDXC4EtSIwy+JDXC4EtSIwy+JDXC4EtSIwy+JDXC4EtSIwy+JDXC4EtSIwy+JDXC4EtSIwy+JDXC4EtSIwy+JDXC4EtSIwy+JDXC4EtSIwy+JDXC4EtSIwy+JDXC4EtSIzoFP8neJCeSLCR5cJnHP5Tkx0meT3I8yb39jypJWosVg59kC/AwsA/YA9yVZM+SY18FXqqqa4E54J+TXN7zrJKkNejyDP8mYKGqTlXV28BhYP+SMwVcmSTAB4E3gXO9TipJWpOtHc5sA05P3F4Ebl5y5rvAEeAMcCXwt1V1fuknSnIAOAAwMzPDYDBYxcjT5+zZs+5ixF2MuYsxd9GPLsHPMvfVktufB54DbgU+CfwsyS+q6vd/9EFV88A8wOzsbM3NzV3svFNpMBjgLobcxZi7GHMX/ehySWcR2DFxezvDZ/KT7gWeqKEF4BXg0/2MKEnqQ5fgPwvsTrJr9IXYOxlevpn0KnAbQJKPAbPAqT4HlSStzYqXdKrqXJL7gSeBLcChqjqe5L7R4weBh4DHkrzA8BLQA1X1xjrOLUm6SF2u4VNVR4GjS+47OPH+GeBv+h1NktQnf9JWkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEQZfkhph8CWpEZ2Cn2RvkhNJFpI8eIEzc0meS3I8yc/7HVOStFZbVzqQZAvwMPDXwCLwbJIjVfXSxJmrgEeAvVX1apKPrtO8kqRV6vIM/yZgoapOVdXbwGFg/5IzdwNPVNWrAFX1er9jSpLWasVn+MA24PTE7UXg5iVnPgVclmQAXAn8S1U9vvQTJTkAHACYmZlhMBisYuTpc/bsWXcx4i7G3MWYu+hHl+Bnmftqmc9zPXAb8H7gl0meqaqTf/RBVfPAPMDs7GzNzc1d9MDTaDAY4C6G3MWYuxhzF/3oEvxFYMfE7e3AmWXOvFFVbwFvJXkKuBY4iSTpktDlGv6zwO4ku5JcDtwJHFly5kfAZ5NsTfIBhpd8Xu53VEnSWqz4DL+qziW5H3gS2AIcqqrjSe4bPX6wql5O8lPgGHAeeLSqXlzPwSVJF6fLJR2q6ihwdMl9B5fc/hbwrf5GkyT1yZ+0laRGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGdAp+kr1JTiRZSPLgu5y7Mck7Se7ob0RJUh9WDH6SLcDDwD5gD3BXkj0XOPdN4Mm+h5QkrV2XZ/g3AQtVdaqq3gYOA/uXOfc14AfA6z3OJ0nqydYOZ7YBpyduLwI3Tx5Isg34AnArcOOFPlGSA8ABgJmZGQaDwUWOO53Onj3rLkbcxZi7GHMX/egS/CxzXy25/W3ggap6J1nu+OiDquaBeYDZ2dmam5vrNuWUGwwGuIshdzHmLsbcRT+6BH8R2DFxeztwZsmZG4DDo9hfDdye5FxV/bCPISVJa9cl+M8Cu5PsAv4HuBO4e/JAVe36//eTPAb8m7GXpEvLisGvqnNJ7mf43TdbgENVdTzJfaPHD67zjJKkHnR5hk9VHQWOLrlv2dBX1d+tfSxJUt/8SVtJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGGHxJaoTBl6RGdAp+kr1JTiRZSPLgMo9/Mcmx0dvTSa7tf1RJ0lqsGPwkW4CHgX3AHuCuJHuWHHsF+KuqugZ4CJjve1BJ0tp0eYZ/E7BQVaeq6m3gMLB/8kBVPV1Vvx3dfAbY3u+YkqS12trhzDbg9MTtReDmdzn/ZeAnyz2Q5ABwAGBmZobBYNBtyil39uxZdzHiLsbcxZi76EeX4GeZ+2rZg8nnGAb/M8s9XlXzjC73zM7O1tzcXLcpp9xgMMBdDLmLMXcx5i760SX4i8COidvbgTNLDyW5BngU2FdVv+lnPElSX7pcw38W2J1kV5LLgTuBI5MHkuwEngC+VFUn+x9TkrRWKz7Dr6pzSe4HngS2AIeq6niS+0aPHwS+DnwEeCQJwLmqumH9xpYkXawul3SoqqPA0SX3HZx4/yvAV/odTZLUJ3/SVpIaYfAlqREGX5IaYfAlqREGX5IaYfAlqREGX5IaYfAlqREGX5IaYfAlqREGX5IaYfAlqREGX5IaYfAlqREGX5IaYfAlqREGX5IaYfAlqREGX5IaYfAlqREGX5IaYfAlqREGX5IaYfAlqREGX5IaYfAlqREGX5IaYfAlqREGX5IaYfAlqREGX5IaYfAlqREGX5IaYfAlqREGX5Ia0Sn4SfYmOZFkIcmDyzyeJN8ZPX4syXX9jypJWosVg59kC/AwsA/YA9yVZM+SY/uA3aO3A8D3ep5TkrRGXZ7h3wQsVNWpqnobOAzsX3JmP/B4DT0DXJXk4z3PKklag60dzmwDTk/cXgRu7nBmG/Da5KEkBxj+DwDgf5O8eFHTTq+rgTc2e4hLhLsYcxdj7mJsdrUf2CX4Wea+WsUZqmoemAdI8uuquqHDrz/13MWYuxhzF2PuYizJr1f7sV0u6SwCOyZubwfOrOKMJGkTdQn+s8DuJLuSXA7cCRxZcuYIcM/ou3VuAX5XVa8t/USSpM2z4iWdqjqX5H7gSWALcKiqjie5b/T4QeAocDuwAPwBuLfDrz2/6qmnj7sYcxdj7mLMXYytehep+pNL7ZKkKeRP2kpSIwy+JDVi3YPvyzKMddjFF0c7OJbk6STXbsacG2GlXUycuzHJO0nu2Mj5NlKXXSSZS/JckuNJfr7RM26UDn9HPpTkx0meH+2iy9cL33OSHEry+oV+VmnV3ayqdXtj+EXe/wL+DLgceB7Ys+TM7cBPGH4v/y3Af6znTJv11nEXfwF8ePT+vpZ3MXHu3xl+U8Admz33Jv65uAp4Cdg5uv3RzZ57E3fxj8A3R+/PAG8Cl2/27Ouwi78ErgNevMDjq+rmej/D92UZxlbcRVU9XVW/Hd18huHPM0yjLn8uAL4G/AB4fSOH22BddnE38ERVvQpQVdO6jy67KODKJAE+yDD45zZ2zPVXVU8x/L1dyKq6ud7Bv9BLLlzsmWlwsb/PLzP8F3warbiLJNuALwAHN3CuzdDlz8WngA8nGST5zyT3bNh0G6vLLr4L/DnDH+x8AfiHqjq/MeNdUlbVzS4vrbAWvb0swxTo/PtM8jmGwf/Muk60ebrs4tvAA1X1zvDJ3NTqsoutwPXAbcD7gV8meaaqTq73cBusyy4+DzwH3Ap8EvhZkl9U1e/XebZLzaq6ud7B92UZxjr9PpNcAzwK7Kuq32zQbButyy5uAA6PYn81cHuSc1X1ww2ZcON0/TvyRlW9BbyV5CngWmDagt9lF/cC/1TDC9kLSV4BPg38amNGvGSsqpvrfUnHl2UYW3EXSXYCTwBfmsJnb5NW3EVV7aqqT1TVJ4B/Bf5+CmMP3f6O/Aj4bJKtST7A8NVqX97gOTdCl128yvB/OiT5GMNXjjy1oVNeGlbVzXV9hl/r97IM7zkdd/F14CPAI6NntudqCl8hsOMumtBlF1X1cpKfAseA88CjVTV1Ly3e8c/FQ8BjSV5geFnjgaqaupdNTvJ9YA64Oski8A3gMlhbN31pBUlqhD9pK0mNMPiS1AiDL0mNMPiS1AiDL0mNMPiS1AiDL0mN+D9l2QhVNNNvxAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAqoElEQVR4nO3de5zWc/7/8cerURRLUUtKje/+WE05D/JtnTassKzFspvDIqXEYhc5rMWuw7IqSwcj58b6opZKTuWQRJrSUQ5JEUkhUekwvX5/vK+pabpOc811mGuu5/12u24z1+fzvj7X++OjXr0Pr/fb3B0REZFYGuW6AiIiUr8pUIiISFwKFCIiEpcChYiIxKVAISIicSlQiIhIXAoUIkkws2IzczPbKtd1Eck2BQopWGa2wMxWm9kP1V735qAefc2swszWmNnDKV6jiZk9HbknN7Mj01pJKWgKFFLofu3u21V79c1BHb4A/gE8mKhgJBAUxzg9ETgL+DJ9VRNRoBCJysyKzOxfZrbMzOYDJ9Q4v4OZPWBmi83sczP7h5kVRc790cwmRj7/rZl9YmbdYn2Xu49092eAr1Otr7uvdfeB7j4RqEz1OiLRKFCIRHchcCKwP1AKnFbj/CPAeuD/RcocC/Sodv4Q4AOgJXAH8ICZWYbrLJIRChRS6J4xs+XVXhdGjv8OGOjun7n7N8BtVR8ws52BbsBl7r7S3b8CBgBnVrvuQne/390rCUGlNbBzVu5IJM00g0MK3W/cfVyU47sCn1V7v7Da7+2BxsDiao2ERjXKbxwncPdVkXLb1bZyZtYOmFnt0PbATDPbEHnfx90fr+11RWpDgUIkusXAbtXet6v2+2fAGqClu6/PZCXc/VOgedV7M1sAHOnuCzL5vSLVqetJJLongUvNrK2ZtQD6VZ1w98XAS8BdZra9mTUys5+Z2RGpfJGZbWVm2wBFQJGZbZNKvoaZbR25DkCTyHU0LiJ1pkAhhW50jTyK/0aO3w+8CMwApgEja3zuHKAJ8B7wLfA0YRwiFdcDqwnB6KzI79encJ0PIp9tQ6j7akI3mUidmDYuEhGReNSiEBGRuBQoREQkLgUKERGJS4FCRETiapB5FC1btvTi4uJcV0NEJG9MnTp1mbu3inauQQaK4uJiKioqcl0NEZG8YWYLY51T15OIiMSlQCEiInEpUIiISFwKFCIiEpcChYiIxKVAISKS55o0AbNNryZN0nt9BQoRkTxmBuvWbX5s3br0BgsFChGRPNSiRQgSsdQMHnWR00BhZseZ2QdmNs/M+sUoc6SZTTezOWb2erbrKCJS35jB8uXZ+76cBQozKwIGETapLwF+b2YlNco0BwYDJ7l7R+D0bNdTRKS+qBqDiOYURnIHV2bke3PZojgYmOfu8919LfAEcHKNMn8ARkb2Dcbdv8pyHUVEci5egNiFxTzNqYzkVLoynmasTPv35zJQtCFsUl9lUeRYdXsCLczsNTObambnxLqYmfU0swozq1i6dGkGqisikl3xAgQ4FzCMuXTgeMZyNbdzCJNZxbYANG2avnrkclHAaLdfc1/WrYADga5AU+AtM3vb3T/c4oPuZUAZQGlpqfZ3FZG8VVQEGzbEPv8z5lFGT37Jq7zKkfSkjHnssfF848awalX66pPLFsUiYLdq79sCX0Qp84K7r3T3ZcAEYN8s1U9EJKuqZjLFChJFrOdK7mAWe3MA07iQMroyfrMgMXw4rF2b3nrlskUxBdjDzHYHPgfOJIxJVPcscK+ZbQU0AQ4BBmS1liIiWRBvqivAfrzLA1zAAbzLSE6hL/eymF03K+MZ6kvJWYvC3dcDfYEXgbnAk+4+x8wuMrOLImXmAi8AM4F3gGHuPjtXdRYRSafy8kTjELANq7mNfkzhIFqzmFMjQ9fVg0Tv3pkLEgDmmbx6jpSWlro2LhKR+ixRCwLgCF7jfi5kD+YxjAu4kjtZTouN5xs3Tl83k5lNdffSaOeUmS0ikkVHH504SOzAcu6jJ69xFI3YwC8Zz4UM2xgkGjUKLYh0j0XE0iC3QhURqY+SaUX8hv8yiIvZmSXcwZXcyI2sphkQprymczZTshQoRESyIFGQ2JkvuZe+nMYIprMvv2Y00zgQSG8XUyrU9SQikmHxg4RzPg8wlw6cyBj6cRsHMWVjkMhmF1MsChQiIhkUL0j8Dx8zjqN5gB7MZB/2YSb/pB/raczw4ZmdyVQb6noSEcmQWEGiiPVczgBu4m+sozE9uY9h9MBpVG+CQ3VqUYiIZECsILEv03mbztzJVbzEsZTwHvfTk4t6188gAWpRiIikXbQgsQ2ruYGbuZI7WUZLTuMpRnAqYPU2QFRRi0JEJE1i5UgczuvMYF+u4XYe4VxKeI8RnEY+BAlQi0JEJC2aNNly+9Ht+Y47uIpelPEx/0NXxvEKXYGQNFdZmYOKpkAtChGROjLbMkicxLO8Rwk9GMad/IW9mbUxSDRunD9BAhQoRETqpGZX0858yZOczrP8hmW05BAmcxV3bsyuznXyXCoUKEREUrR5kHD+yEO8Rwm/ZjTXcCulVDCVTevsNW2af0ECFChERFJSPUjsznxe5hge4nxm04l9mcHtXMN6Gm8sU1KSm3Wa0kGBQkSklqqCRBHruYK7mE0nDmIKvRjKkbzGh/x8s/LDh8OcOTmoaJpo1pOISC1UBYl9mMEwenAQFTzLSfRhMF/QZovy+TD9NRG1KEREkmQGW/Mj/+A6KiilHZ9yOk/yG55psEECFChERBKq2rL0MCYwg325jlsZzll0YC5Pczqw+dSnpk0bTpAAdT2JiMTVogVsWP4dg+lHb4Yyn905mpcZz9FRyzdvDt9+m906ZppaFCIiMZjBYctHMYeO9KSMu7iCvZkVM0jsumvDCxKgFoWISFQ72xKe4FLO4ElmsA+n8F8qOChm+V13hc8/z2IFs0iBQkSkOnf+2OgR5nIF27KS6/gHd3DVZjkRNTVt2nCDBKjrSURkk/nzebnRsTzMecyhI/syg1u5LmGQyNdEumQpUIiIVFZC//6s/NneHMJkejOYI3idD9gr7se6dm34QQLU9SQihW7mTN4/rAd7rZjCeH5NHwbzOW0TfqwhTX9NRC0KESlMP/4I11/Pun0PZMcVCziDJziZZxMGiebNCytIQI4DhZkdZ2YfmNk8M+sXp9xBZlZpZqdls34i0kC98Qbstx/ccgvldKcDc3mSM6iZOFdT164Nc/prIjkLFGZWBAwCugElwO/NrCRGuX8CL2a3hiLS4KxYAX36wOGH88kHaziWFzmPh/mGnRJ+tHdvGDcuC3Wsh3LZojgYmOfu8919LfAEcHKUcpcAI4Cvslk5EWlgRo+GkhIqh9xHfy6nE7N5mWOT+mjv3jB4cIbrV4/lMlC0AT6r9n5R5NhGZtYGOAUYmuhiZtbTzCrMrGLp0qVpraiI5LElS+CMM+Ckk5j1eQsO5S3+TH9WsW1SHx8+vLCDBOR21lO0zsCaQ0QDgavdvdJq7jdY84PuZUAZQGlpaYENNYnIFtzh0Ufh8stZ8+1K/s7fuYOrWEeTWl1CchsoFgG7VXvfFviiRplS4IlIkGgJHG9m6939mazUUETy0yefQK9e8PLLTKQLF3I/79Mh6Y/n477WmZTLrqcpwB5mtruZNQHOBEZVL+Duu7t7sbsXA08DfRQkRCSmykoYMAA6deLH196iD4M4nAm1ChL5uq91JuWsReHu682sL2E2UxHwoLvPMbOLIucTjkuIiGw0cyb06AFTpjBu6xM4b90QFm3WaZFYQ17Yry5ympnt7mOBsTWORQ0Q7v7HbNRJRPLMjz/CLbfA7bdDixac0/g/PLYmcU5ETcOHQ/fumalivtMSHiKSvyZOhAsvhPffh3POodVj/VnmiXMiquvatXDzI5KlJTxEJP+sWAEXXwyHHQarV3PPCS9gjz5SqyBRtV2pgkRiChQikl/GjIGOHWHIELjsMlotmc2lz/2qVpfo3bswVn1NF3U9iUh++Oor+NOf4IknoFMnXujxNN1uPKTWl1FuRO2pRSEi9VtV4lyHDjByJNx8M/tVTq11kKjqapLaU4tCROqvBQtC4txLL0GXLnD//TQ7sAOrV9fuMpr2WjdqUYhI/VNZCQMHhrGISZPg3nu5eO8JWEntg0Tv3goSdaUWhYjUL7NmhcS5d96BE06AIUNosc9uLF9eu8s0ahTijdSdWhQiUj+sWQM33AAHHADz58Pjj9Nnt9FYu9oHia5dFSTSSS0KEcm9N98MrYj334ezz4b+/WmxR0u1IuoJtShEJHe+/x769t2YOMfzz1P+q0exVrUPEs2bK0hkigKFiOTGc89BSUnYFejSS2H2bPqMOo6zzqr9pQp1L+tsUdeTiGTX0qUhce4//wmzmp56Cjp3pmNHeO+92l3KDB57TIv5ZZpaFCKSHe7hb/UOHeDpp+Gmm2DaNMo/7oxZ7YNE166wYYOCRDaoRSEimbdwYUice/FFOPRQGDYMSkpSakU0bap1mrJNLQoRyZzKSrj77tDF9OabcM89MHEi5e+WpNSK0GJ+uaEWhYhkxpw5cMEFMHkydOsGQ4dCu3ZqReQhtShEJL3WrIG//Q323x8+/hjKy+G55yh/o11KrYjmzRUkck0tChFJn0mTQuLc3Llw1lkwYAC0bMnRR8P48bW/nHafqx8UKESk7r7/Hq69FgYNgt12g7FjQ3cT0KYNfPFF7S6nrqb6RV1PIlI3Y8eGwepBg+CSS2D2bOjWjfLykOdQ2yDRtauCRH2jFoWIpGbpUrjsMnj88ZBh/eabYeorpNTV1Ly5sqvrK7UoRKR23GH48JA499RTcOONMG0aHHoo5eVhYb7aBonevRUk6jO1KEQkeQsXhr/Vn38eOncOiXMdOwJqRTRkalGISGKVlfDvf4egMGFC+H3iROjYkT59wlhEbYOEFvLLHzkNFGZ2nJl9YGbzzKxflPPdzWxm5DXJzPbNRT1FCtqcOfCLX4SF/A47DObMoXzHS9huhyLMYMiQ2l3OLPRcadpr/shZ15OZFQGDgGOARcAUMxvl7tXTcT4BjnD3b82sG1AGHJL92ooUoDVr4Lbb4NZbYfvtYfhwyv0PnL+nsXZtapfcdVftX52PcjlGcTAwz93nA5jZE8DJwMZA4e6TqpV/G2ib1RqKFKq33gqJc++9B92783SXAXQ/v1XKAQLC0MbgwemromRPLrue2gCfVXu/KHIslguA5zNaI5FC98MPYROhLl1CEt1zz3H0l8M5vU/qQWLXXcNEKQWJ/JXLQGFRjnnUgmZHEQLF1TEvZtbTzCrMrGLp0qVpqqJIAXnhhTBYfe+90Lcvlx87Bzvh+JSW3qjSu7e6mhqCXHY9LQJ2q/a+LbBFDqeZ7QMMA7q5+9exLubuZYQxDEpLS6MGHBGJYtmykDhXXs7iFh04zScy6Z7/rdMltUZTw5LLFsUUYA8z293MmgBnAqOqFzCzdsBI4Gx3/zAHdRRpuNxDVnWHDqz/z5PcyN8o/vZdJpF6kOjaNVxWQaJhyVmLwt3Xm1lf4EWgCHjQ3eeY2UWR80OBG4CdgMFmBrDe3UtzVWeRBuPTT/n8pN60mTGWtzmEHgxjDp1SvpwGqhs2c294vTSlpaVeUVGR62qI1D8bNjDlvMHs9eg1GM613MogLmYDRbW+1FZbwcMPa8/qhsLMpsb6h7gys0UKQJ8+0NHeY1LRLzjo0Ut4ky50Yjb3cGlKQaJ3b1i3TkGiUChQiDRA5eXQsmXIgm5ia2k55GamsT978iFn8RjdeJ6FFNf6ur17a6prIYo5RmFmA4E3gUnurgluIvVceXlYZePranMDD+FthtGDTszhcX7PZQxkKT+t1XUbNYJevRQcClm8FsU84BTgTTNbYGaPm9nFZra/maklIlIPVG85nHXWpiCxLT8wgMuYxP+yPSs4gTF05/FaBYmq1kNlpYJEoYvZonD3e4F7AcysNdAF+F/gcuCnwPbZqKCIRG8txHIsL3IfvWjHpwymD9dwGz/wk6S+xwwuukiBQTYXd3qshTmpexMCRBeghNDSeCzzVROR8vLQ7bNyZeKyO7GM/lzBOTzGXPbiMN5gEl2S+p7ttoOhQzU4LdHFG6N4mdBqmE5YkO9Wd5+bpXqJFKTatBw2cc7kCe7mT7TgW27mr9zCdaxl67if0tiDJCtei2I+sC+wB/A1sMzMlrr7sqzUTKRApBYcgrZ8xhB6cyLPMZmD6cp4ZrN3zPJqOUgq4o1R9AIws+2BzoTup4vNrBUw293PzU4VRRqOugSF6owN9GYIt9OPRmzgMgZwD5dskROhVoOkQzJLeKwBVgGrI7+3BZpkslIiDUW6AkN1ezGXYfSgC5N4kWPpxX2b5UQoOEi6xZzmamYDzGwysBi4GfgJcB/wc3eP3bYVKTDl5VBcHGYMbbVV+Nmo0ZZTVuuqMWu5nr8znf3Yi/c5h0c4jhdYSDE77RS2F9V0VsmEeC2KT4By4F13r8xSfUTySp8+oc+/asm0ysiflHQvoXYwkxlGD/ZmNiObnMmGAXfzaJ+f8mh6v0YkqnhjFP/OZkVE8kF5OVx3HSxcGFoMmV5TsxkruWub6+m15m6sTRsYMprfnnhiZr9UpAZlWIskoSoD+qyzQpCAzAaJnXaCcVe/xMriTlz040Csd2+YMwcUJCQHFChE4qgeINI5IF1do8ifwvbtI+MMy75m2Qnn0vWfv4Ktt4Y33oBBg2B7LYYguZEwUJjZzmZ2QGSNp52zUSmRXMtkgKg+8Fw1+OwOCz5xuhc9AR06hJ3nrr8epk+HX/wivRUQqaV4mdn7AUOBHYCq1WPbmtlyoI+7T8t47URyoLwcevaEVatSv0b18YuddoK7706Q5PbZZ2FkfMwYOOigsJfoPvukXgGRNIo36+lhoJe7T65+0Mw6Aw8RsrZFGpTycjj33E2zl2ojqYBQ04YNYdpUv37hS/v3h0svhaLabyYkkinxup62rRkkANz9bWDbzFVJJLuq50GcfXbyQaLm2MKyZbUMEu+/D4cfDhdfDJ07w+zZcPnlChJS78RrUTxvZs8BjwKfRY7tBpwDvJDpiolkQ81upmRmMqXUcqhu7Vq44w74+99h223DxtPnnBMilUg9FC+P4lIz6wacDLQBDFgEDHL3sVmqn0hGXXdd8mMRdQ4QAO+8Az16wKxZcMYZ4YI7a46I1G9x13py9+eB57NUF5GsKi/flBMRT1ERPPJIHQPEypXw17+GwNC6NTz7LJx0Uh0uKJI9KeVRmFlZuisikk1VXU6JNGuWhiDx0kvQqRMMGBBW65szR0FC8kq86bE7xjoFHJ+Z6ohkVvUlOGKpmtravj3ccksdgsTXX8Of/xwizZ57woQJcNhhKV5MJHfidT0tBRYSAkMVj7xPfod2kXoi2fyIxx6rYwvCHZ58Mkxz/eabEJmuvx622aYOFxXJnUQ73HV1909rnjCzz6KUF6nXkhm4bt++jkFi0aKQODd6NJSWhm6nfZVyJPkt3hjFQKBFjHN3pOPLzew4M/vAzOaZWb8o583M/h05P9PMDkjH90ph+nSLf/Jsrlmz0NWUkg0bYMgQKCkJWdV33QVvvaUgIQ1CzEDh7oPcfUaMc/fU9YvNrAgYBHQDSoDfm1lJjWLdCHt27wH0BIbU9Xul8FQl1MXLkWjfHsrKUmxNvP8+HHFEaEkcckhInLviirCLkUgDkPD/ZDP7bZTD3wGz3P2rOnz3wcA8d58f+Z4nCDkb71UrczLwqLs78LaZNTez1u6+uA7fKwUg2X0jmjWrQ4BYty4kzt18c0ice+ihsP6HEuekgUnmnzwXAIcCr0beHwm8DexpZje7+2MpfncbNmV8Q0jmOySJMm0I27Nuxsx6ElodtGvXLsUqSUOQbLZ1nWY1TZkCF1wQEud+9zv497+VOCcNVjKBYgPQwd2XQFh2nNAFdAgwAUg1UET7Z1fNP9LJlAkH3cuAMoDS0tIM7zsm9Vkyg9ZmsGBBChdfuRJuuAEGDoRddlHinBSEZBLuiquCRMRXwJ7u/g2wrg7fvYiwdlSVtsAXKZQRATaNRSSTbZ1So3PcONh777DCa8+e8N57ChJSEJIJFG+Y2RgzO9fMzgVGARPMbFtgeR2+ewqwh5ntbmZNgDMj165uFHBOZPZTZ+A7jU9INFXdTckEiVrPbvrmGzjvPDjmGGjcGF5/Pcxw2mGHlOsrkk+S6Xq6GPgt8AtCV9AjwIjIAPNRqX6xu683s77Ai0AR8KC7zzGziyLnhwJjCVng84BVwHmpfp80bIm6m1LKtnaHp56CSy4JweLaa8N6TUqckwKTMFC4u5vZRGAtYXzgnUiQqLPIKrRjaxwbWv27CYFKJK54ORIpDVp//nmY7jpqFBx4oBLnpKAls2f274B3gNOA3wGTzey0TFdMJJGqMYlGjTZtIlRT+/Zh0DrpIFG141xJCbz8MvzrX/D22woSUtCS6Xq6DjioKmfCzFoB44CnM1kxkXhqToGNtitdrcciPvggXHTCBOjaFe67D372s7TUVySfJTOY3ahGYt3XSX5OJGNijUkUFYXxiFplWq9bB7feGloNM2fCgw+G1oSChAiQXIviBTN7EfhP5P0Z1BhXEMm2WGMSGzaEV9IqKsKOczNmwOmnh8S5XXZJSx1FGoqELQN3v5KQyLYPsC9Q5u5XZ7piIvHEyoNIOj9i1Sr4y1/C2kxLl8Izz4SlwRUkRLaQ1Kpl7j4CGJHhuogk7ZZbttxbIukxifHjw4fnzw87zv3zn8qJEIkjZovCzL43sxVRXt+b2YpsVlKkpu7dwxhE+/a1GJP49ls4/3w4+ugwmPHaa2GGk4KESFzxlhn/ibtvH+X1E3ffPpuVFIHNp8MWF4djCxaEMYm4U2Dd4emnoUMHePRRuOaaMCZxxBFZqbdIvtOC+ZIXak6HXbgwvIcErYjPP4eLLw6L9x1wALzwAuy3X6arK9KgaJqr5IVo02FXrQrHo9qwIfRFlZSErOo774TJkxUkRFKgQCH1VvWupliL/UWdJvvhh/DLX4aB6tLSsGfEX/6iHedEUqRAIfVG9cDQsmUYd164MP4WpptNh123Dm6/HfbZJ4xBPPBAWBpciXMidaJ/Ykm9UHMM4uuvE39ms+mwU6eGxLnp0+HUU+Gee6B160xVV6SgqEUhWVdz9lLV/taJdqWrstl02FNWwVVXwcEHw5IlMHJkmOGkICGSNmpRSFbFmr2UbJCoWg0WgFdegX16wscfw4UXwh13QPPmGai1SGFTi0KyKtbspaKixJ/d2NX07behm6lr19C8ePXV0LxQkBDJCAUKyaia3UyxZi9VVoZAUF3jxrDTTjW6mrYZEaa8PvwwXH11WO31yCMzexMiBU6BQjKm+j7W7uGnWfSyVYGg+pIcDz0Ey5ZFMq8nfUH3Eb+F006DXXeFKVPCDKemTbN7UyIFSGMUkjHRupncN+1fXaWqS6l79yhZ1hs2hGmuV14Ja9aEBfyuuEI5ESJZpBaFZEysPSPck1zM76OPwjhEz55h+Y1Zs8IMJwUJkazSnzjJmHbtoo9JbDZzKZp166B/f7jxRth6a7j/frjggtj9ViKSUWpRSELR8h6SOX/LLVsOUCfcM2LatLCZUL9+cPzxMHdumOGkICGSO+7e4F4HHnigS2LDh7u3b+9uFn4OHx69TLNm7qHDKLyaNdtUNpnzib7D3d1XrnS/6ir3oiL3XXZxHzEi7fcrIrEBFR7j71TzeAvp5KnS0lKvqKjIdTXqtZqJbxD+tV9zvCDWlNaq7qNE55Py6qshYe7jj0Pr4c47lRMhkmVmNtXdS6OdU9dTgUp22e5YA9JVxxOdj2v58hAgfvnL8P6VV8J4hIKESL2iQNFAJBpHqCnZv+A3W501yvFE52MaOTLsOPfQQ2Em06xZcNRRCT4kIrmQk0BhZjua2ctm9lHkZ4soZXYzs1fNbK6ZzTGzP+WirvkgWmJbz57xg0Wyf8EnGpCu9YD14sVhdddTT4VddoF33gm5EUqcE6m/Yg1eZPIF3AH0i/zeD/hnlDKtgQMiv/8E+BAoSeb6hTaY3b795oPJVa/27WN/JtEgdM2y8Qakkxqw3rDB/f773XfYwX2bbdxvv9197drUblhE0o76NphtZh8AR7r7YjNrDbzm7j9P8JlngXvd/eVE1y+0wexGjaJv7mMWEptjqVre+9NPQ0uiKjs67ebNC02cV18N6zKVlcEee2Tgi0QkVfVxMHtnd18MEPn503iFzawY2B+YHKdMTzOrMLOKpUuXprOu9V6q4wTdu4eZSRs2hJ9pDxLr14elv/feO+RHlJXB+PEKEiJ5JmOBwszGmdnsKK+Ta3md7YARwGXuviJWOXcvc/dSdy9t1apVXaufV1JKbMu0d98NmwldfTV06wbvvRdmODXS/AmRfJOxJTzc/ehY58xsiZm1rtb19FWMco0JQaLc3UdmqKp5r6olkJVupERWr4abboJ//QtatQq7zZ16ag4qIiLpkqu1nkYB5wK3R34+W7OAmRnwADDX3ftnt3r5J+rKq9n22muh1TBvXlib6c47ocUWE9pEJM/kqh/gduAYM/sIOCbyHjPb1czGRsp0Ac4Gfmlm0yOv43NTXYlr+fIwWH3UUWHAY/x4GDZMQUKkgchJi8Ldvwa6Rjn+BXB85PeJgFaCq+/++1+4+GJYsiTsGXHjjVsOmIhIXtMy45KaL7+Evn1hxAjYd18YPRoOPDDXtRKRDNAUFKkd97DjXIcOMGYM3HZb2JZUQUKkwVKLQpL38cdhLOKVV+Dww8MCfnvumetaiUiGqUUhia1fH2Yw7b03VFTAffeFLGsFCZGCoBaFxDd9epjqOm0anHwyDBoEbdrkulYikkVqUUh0q1fDNddAaSl8/jk89VSY4aQgIVJw1KKQLb3+ekic++gjOP/80O204465rpWI5IhaFLLJ8uXQq1dY4bWyEsaNCzOcFCRECpoChQTPPAMlJSGj+i9/CTvOdd0iJ1JECpACRaH78ks4/XQ45ZSwiN/kyaGrSdnVIhKhQFGo3OHBB0Pi3OjRcOutYepradR9S0SkgGkwuxBVT5w77LCQOPfzuBsMikgBU4six8rLobg47OdTXBzeZ8z69WGfiKrEuaFDw9LgChIiEodaFDlUXh7+Yb9qVXi/cGF4DxnYW2L6dOjRA6ZOhZNOgsGDlRMhIklRiyKHrrtuU5CosmpVOJ42P/4I114bxh4++wyefDLMcFKQEJEkqUWRQ59+WrvjtTZhQkic+/BDOO+80O2knAgRqSW1KHKoXbvaHU/ad9/BRRfBEUfAunXw8sthhpOChIikQIEih265Zct0hWbNwvGUPftsSJy7/374859D4tzRR9epniJS2BQocqh7dygrg/btwSz8LCtLcSB7yRL43e/gN7+Bli3h7bdDV9O226a72iJSYDRGkWPdu9dxhpM7PPxwaD2sWhWaI1deCY0bp6uKIlLgFCjy2fz5YRG/ceOUOCciGaOup3y0fj3cdRd06hTWZhoyRIlzIpIxalHkm5kzw45zFRXw61+HxLm2bXNdKxFpwNSiyBc//hgy8Q48MCRa/N//hRlOChIikmFqUeSDN94IiXMffAB//GOYzbTTTrmulYgUCLUo6rMVK6B3bzj8cFizBl56CR56SEFCRLIqJ4HCzHY0s5fN7KPIzxZxyhaZ2btmNiabdcy50aND4lxZGVxxBcyeDccck+taiUgBylWLoh8w3t33AMZH3sfyJ2BuVmpVHyxZAmecEVZ43XFHeOutMMNJiXMikiO5ChQnA49Efn8E+E20QmbWFjgBGJadauWQOzzySNhx7pln4B//CDObDj441zUTkQKXq0Cxs7svBoj8/GmMcgOBq4ANiS5oZj3NrMLMKpYuXZq2imbFJ5/Ar34VBqo7doQZM8IMpyZNcl0zEZHMBQozG2dms6O8Tk7y8ycCX7n71GTKu3uZu5e6e2mrVq3qVPesqayEAQNC4tzbb4eciNdfh732ynXNREQ2ytj0WHePuWSpmS0xs9buvtjMWgNfRSnWBTjJzI4HtgG2N7Ph7n5WhqqcXTNnhh3npkyBE08MQWK33XJdKxGRLeSq62kUcG7k93OBZ2sWcPdr3L2tuxcDZwKvNIgg8eOP8Ne/hsS5BQvgiSdg1CgFCRGpt3IVKG4HjjGzj4BjIu8xs13NbGwuKlReDsXF0KhR+FlenoEvmTgR9t8/DFT/4Q8wd26Y4WSWgS8TEUmPnGRmu/vXQNcox78Ajo9y/DXgtUzVp7wcevbctH/1woXhPdRxCfAqK1bANdeE7qXiYnjxRTj22DRcWEQk85SZTZhgVBUkqqxaFY7X2ZgxYSbTkCFw2WVhxzkFCRHJIwoUhDX2anM8KV99Bb//fVjhtXnzkDg3YABst10dLioikn0KFEC7drU7Hpc7PPpoSJwbORJuvhmmToVDDqlTHUVEckWBgrB7aLNmmx9r1iwcr5UFC+C44+Dcc0OgmD49zHBS4pyI5DEFCsKAdVkZtG8fJiC1bx/eJz2QXVkJAweGsYhJk+Dee2HChBAsRETynPajiOjePcUZTrNmhcS5d96BE04Ig9bKiRCRBkQtilStWQM33AAHHADz58Pjj4elwRUkRKSBUYsiFW++GVoR778PZ58N/ftDy5a5rpWISEaoRVEb338PffvCYYfB6tXw/PNhhpOChIg0YAoUyXruubDj3ODBcOmlYce5447Lda1ERDJOgSKRpUvDukwnngg77BBmNQ0cqMQ5ESkYChSxuMNjj4Uprk8/DTfdBNOmQefOua6ZiEhWaTA7moULoVevsHjfoYfCsGGh20lEpACpRVFdZSXcfXdInHvzTbjnnrA0uIKEiBQwtSiqfPstdOsGkyeHn0OHprjYk4hIw6IWRZXmzeFnPwubUzz3nIKEiEiEWhRVzDK0rZ2ISH5Ti0JEROJSoBARkbgUKEREJC4FChERiUuBQkRE4lKgEBGRuBQoREQkLgUKERGJy9w913VIOzNbCixM8eMtgWVprE4uNZR7aSj3AbqX+qih3AfU7V7au3uraCcaZKCoCzOrcPfSXNcjHRrKvTSU+wDdS33UUO4DMncv6noSEZG4FChERCQuBYotleW6AmnUUO6lodwH6F7qo4ZyH5Che9EYhYiIxKUWhYiIxKVAISIicRVkoDCz48zsAzObZ2b9opw3M/t35PxMMzsgF/VMRhL3cqSZfWdm0yOvG3JRz0TM7EEz+8rMZsc4n0/PJNG95Msz2c3MXjWzuWY2x8z+FKVMXjyXJO8lX57LNmb2jpnNiNzLTVHKpPe5uHtBvYAi4GPgf4AmwAygpEaZ44HnAQM6A5NzXe863MuRwJhc1zWJezkcOACYHeN8XjyTJO8lX55Ja+CAyO8/AT7M4z8rydxLvjwXA7aL/N4YmAx0zuRzKcQWxcHAPHef7+5rgSeAk2uUORl41IO3geZm1jrbFU1CMveSF9x9AvBNnCL58kySuZe84O6L3X1a5PfvgblAmxrF8uK5JHkveSHy3/qHyNvGkVfNWUlpfS6FGCjaAJ9Ve7+ILf+HSaZMfZBsPQ+NNFOfN7OO2ala2uXLM0lWXj0TMysG9if867W6vHsuce4F8uS5mFmRmU0HvgJedveMPpetUv1gHrMox2pG42TK1AfJ1HMaYQ2XH8zseOAZYI9MVywD8uWZJCOvnomZbQeMAC5z9xU1T0f5SL19LgnuJW+ei7tXAvuZWXPgv2bWyd2rj4ml9bkUYotiEbBbtfdtgS9SKFMfJKynu6+oaqa6+1igsZm1zF4V0yZfnklC+fRMzKwx4S/WcncfGaVI3jyXRPeST8+lirsvB14DjqtxKq3PpRADxRRgDzPb3cyaAGcCo2qUGQWcE5k50Bn4zt0XZ7uiSUh4L2a2i5lZ5PeDCc/866zXtO7y5ZkklC/PJFLHB4C57t4/RrG8eC7J3EsePZdWkZYEZtYUOBp4v0axtD6Xgut6cvf1ZtYXeJEwa+hBd59jZhdFzg8FxhJmDcwDVgHn5aq+8SR5L6cBvc1sPbAaONMj0yLqEzP7D2HWSUszWwT8jTBIl1fPBJK6l7x4JkAX4GxgVqQ/HOBaoB3k3XNJ5l7y5bm0Bh4xsyJCMHvS3cdk8u8wLeEhIiJxFWLXk4iI1IIChYiIxKVAISIicSlQiIhIXAoUIiISlwKFSJLM7IfEpWJ+tm9kJU+vnsSV7CqfZnZtqt8tUlcKFCLZ8SYhMWphjePdCMtE7AH0BIbE+LwCheSMAoVILUVaAXea2Wwzm2VmZ0SONzKzwZE9AsaY2VgzOw3A3d919wVRLpdwlU8zux1oamGPhPIM357IFgouM1skDX4L7AfsC7QEppjZBEL2bzGwN/BTwlLWDya4VqxVPjcut+Du/cysr7vvl57qi9SOWhQitfcL4D/uXunuS4DXgYMix59y9w3u/iXwahLXyqvVV6UwKVCI1F60v9zjHY8nb1ZflcKlQCFSexOAMyKbx7QibH36DjARODUyVrEzYWHARJJd5XNdZJlskaxToBCpvf8CMwl7lL8CXBXpahpBaCHMBu4j7KD2HYCZXRpZSbYtMNPMhkWuNRaYT1jl836gT4zvLIt8ToPZknVaPVYkjcxsu8gOaTsRWhldIkFEJG9p1pNIeo2JbCrTBPi7goQ0BGpRiIhIXBqjEBGRuBQoREQkLgUKERGJS4FCRETiUqAQEZG4/j88nqgDNsQihQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import random as rd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "from scipy import stats\n",
    "\n",
    "def d(x,y):\n",
    "    r = pow(x**2 + y**2, 0.5)\n",
    "    return r\n",
    "\n",
    "plt.figure()\n",
    "\n",
    "t=1000     # 设定时间步长\n",
    "\n",
    "point = np.array([[0,0]])\n",
    "step = np.array([[0,1],[0,-1],[1,0],[-1,0]])\n",
    "edge = step\n",
    "\n",
    "dis=[]\n",
    "dis = np.array(dis)\n",
    "wt = []\n",
    "for i in range(0,t):\n",
    "    growsite = rd.randint(0,len(edge)-1)\n",
    "    point = np.append(point,np.array([edge[growsite]]),axis=0)\n",
    "    edge = np.append(edge,np.array([edge[growsite]]) + step,axis=0)\n",
    "    j=0\n",
    "    for m in edge:\n",
    "        for n in point:\n",
    "            if m[0]==n[0] and m[1]==n[1]:\n",
    "                edge = np.delete(edge,j,axis=0)\n",
    "                j=j-1\n",
    "        j=j+1\n",
    "    # print(i)\n",
    "    for j in edge:\n",
    "        dis = np.append(dis,d(j[0], j[1]))\n",
    "    w = np.std(dis)\n",
    "    wt.append(w)\n",
    "wt = np.array(wt)\n",
    "\n",
    "plt.grid(True)\n",
    "plt.figure()\n",
    "plt.scatter(np.log10(np.arange(1,len(wt)+1)),np.log10(np.array(wt)),c='blue', label='OD')\n",
    "x=np.log10(np.arange(1,len(wt)+1))\n",
    "y=np.log10(np.array(wt))\n",
    "res = stats.linregress(x, y)\n",
    "plt.plot(x, res.intercept + res.slope*x, color='r', label='fitted line',)\n",
    "print('β =',res.slope) \n",
    "plt.title('Eden 1+1')  \n",
    "plt.ylabel(\"log10 W\")\n",
    "plt.xlabel(\"log10 t\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54eb5e6d-9d9b-44fc-a65b-524f48af0fc5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}