{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# Large amplitude pendulum" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Generate and write out pendulum data for large angle oscillations. (See [this paper](https://arxiv.org/pdf/physics/0510206.pdf) for the approximate formula used to calculate the period!). Here, $\\theta_0$ is the maximum displacement of the pendulum, assumed to occur at time $t = 0$.\n", "\n", "$$\n", "T = - 2\\pi \\sqrt {\\frac{L}{g}} \\frac{{\\log \\left[ {\\cos \\frac{{{\\theta _0}}}{2}} \\right]}}{{1 - \\cos \\frac{{{\\theta _0}}}{2}}}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This formula can be inverted." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs4AAACaCAYAAAC5bM7kAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae2dS87dxLqGf6LoNFF2kNDunjADEkawQ58GlxEAMwDRoxfBDIARcGnsPuwRcJlBON0tpB0idJpHynlf//4cLy/fr2X7KcnLdt1c9VT59edy2euVFy9e3OAgAIFjE/j888/vqYZ/HruW/WsnHq/0j01MCEAAAuMIoL2X3I6gvXcvq8QeBCBQR0An+6/y/1Trn+rCd+D3vsr4nsr/ww7KuloRxeOBDvajlofafr7agTkQBCCwCgGd11trN9pb09J71t47NfXBCwIQKBHQCf6Vdn/Xeq9Gs2vztpY9l991mN2pTX9Xpr6Z+Gb2zMkQAhDYlEAi2o321vSCPWsvI841DYoXBIKATu53te0Rg/8Ov52u76kuSY6oqlyfiKkvLh75fU2L97/U8h8tb2nxTcunWi/inLeWP7V8pOXrRQ5CphCAwKoEdC6not1ob0PLq412qb13VfA3VSc/yqhzPyj8vboA/G5uxGbrR0A0w4IE1L6eF+yRSJ/cSRqdfaqvsns6gkdWU3WvqYw2nG+0th69q3VmKGvtNvCNy9LOOvejjveTlkmslJ7pH0u31sT81UZo90SGKSdX+yah3bkWTNKThTmjvTWA1W5P5W0dv3IKe6U8VcMjPB9XFj+ixtUQELwjPL6vqRleJQI2mj3aufdRyMeqx/eleqW2WdYZl7WYUiL2vmH5bukC6zg+pqdslMsy6rDKi+kfo8itk0jtg3avg3rLo6Si3WhvRy9IVHu/ULGr9nDxflB5qsZXueB3VJNgcUrlERCNsRABtXE28qnss5HQhQ6zVrZvqz7JPjmq6I55XxivCl9rtN+j3E91PI94FyI5ppGUfpePIMfUdU9p3LYq7xGmXu0J+6plVRunpN1ob7/WT0p71YeuBsvk55pYP27KI87ex3UQELwkHgF1FJPg6QQ8YvGb2rsY/Zye5WY5uM/uxV2MOK9ZaLW1R4rd3h5tmMP5ZsUDErWP/OY4AHn0J4B292e185gpaTfa26Mz7U17MZx7NGolSiqPgCrFYncuAjqJfVfpUYvFXkibq6xd+agursdvXfFSCM/L6qkxa40w11Xbbf5AZfALipOc8pht+sekgpA4CKDdQeKga51zyWi3yoL2Dutnu9FeDOcBDZufCD4xd29QDaj2GaN+pkrv/fNz0W4ewf0xdlJd69yyoZqN9ObbmxRVx/ZNhg1e94E5nLXisfLNHvHNkSF5DCcg/vH4Hu0ejm9PKVLSbrR3QM/Zk/beHVAvot5+YeEoj+9pzxoC+QX2EKPNefU8x84v/i7mlP/cfz/6hfKca7pEUW/l2fffAj3HOjN2lWbqXOe4AXN9JuVVVISNMQRSenw/pvyk6SCgc9W6nZJ2o715mx1NezGcO07GCFbDxyOgI7wsFtVifU0gRhqvXg64jno8H/XzB6pV26jc94pzMe97gCjuApjq489wuqzuC3MYu+b5q/L8RMuiNzEuNO6SgJij3ZdIjrq3a+1WP0V7d6K9GM79JcQnZYwe9U9FzN0QkHD5RQ5fZP0t3y3n2c7CTHXw6MugaRpK4xfk/BmeszsbzP66xptaJs0Rd3otMf0Dw3n9noV2r8981SPq/EpKu1UetHd8D0heezGcezRufhKk9AioR6mJMoJA/NHGVyPSppjksQp1MTqcYiH7lCk/B5v+qKkzC6XvO00j8nIf8E3UB1omGc55hs5vlukfeX6sehBAu3tAOkaU1LQb7c371RG1F8O5n2h4xMLulI/vb6t+it/41vEhjE212OJz7FbsFZ4j/FAifGXEys9hfrnQo7nf1sWR/yCnPPzUwWk+0tI2dcVxOp3ymnv6R+cxiZARQLvP0RFS0260d2S/24P23hlZt9MkUyMm9QjoNOBXrmjezh4l2PpzaHPWfPfTTQxDbeOnPV5fGc32l/N0FLeb/3SkKU4WceCPb6DuKc/s+APT1kX3I0hP/Zgrv7pj4JcTEGe0+wS9IW/n1LQb7Z3W95LWXgzn7sZN7RFQd4mJMYZAtPMcL4ONOf6saXQx8YXk51kz3S4zjxq2fWXDhugSTwlifrina8zhPF3Dbq78bnPjt4lAnNPBvSke/vsmEO2chHajvbN0pqS1F8O5u41TewTUXWJijCEQX0uJE3ZMHimlcX2SuJBMgaKL0AOl9x+StBnGruvo+c8t5Ytj+iZksivVwdM/cMsTQLuXZ5zCEVLTbrR3eq9IWnsxnFsaWBc6P+pL7RFQS4kJmkAgM45Kxs2ErJZPqnLamGwz6DwlwF/I2Lvz/OKuOcaPFCeEdrb6il9M+zBLa8EczuWcc/rHHGU6XB55e6Hdh2vZ2gqtqt3qW2jvy2Y4pfbefVl/tmoIJPUIqKZ8m3pJQL5XAfzlgUlO+Qz94sGk41UT5xdZG0Z7MjQ9Mm4j7MpgzOuz+zl2eT0ead34eTyFeUTahuhF22l/rm8m23j2VJC5LhBuN1/oPV0jDHNt4mYmgHa3ANX5gXa38OkIQnsFSH3otNqL4dx+hqT2CKi9tCuG6qSxoemRuE2N3pmqnI1YKK8L42umvGfPRsz9qN8Gc5S7egz7H2HKiec1P6lWrrJ/Nb9ZfD5RnKsbikq6vru/KKKPcXWcvhlU4kW5mtquEp3dkQTQ7gZwaHcDmB7eaO8FpCtNPIv23r3AwE6VQHZxU2eIi101/Mz7ZrP7ObR5A8ZFdi8jgN+p3DbyP/JFUEt1dNn1aXuZLq/2fCuVwf3Bx/QohMsW618U1jhirHi1zvVSgL973JhWYX7akX1uTNs2lu1cdz9KneuPRmLutPOdnKfK5T9DUVbZTWdd2zkMN50A2t3MEO1uZtMVgvaKkDTs1Np7p6uXnDVcHcMXbi+7GIXcoJ38qPnbDY67xCH9GN5uF1+hUN98riVu5qLstzW4/bXhuHq/1TEf6vBfl9fabjR8ywWu2W79kobyfaE0ftwco8E22r3YKHD9X+RxtDvJecTZro7zbcjw37hBmzPP4aU4aAq1O9rd3rZodzufxlD1LbQX7b2529hDCMhGLIRhdQNkJ+g9TSMMgJ0UubGYHh2121tbm79HQsOIvsmNhtXroeP6D0NssISzQVvrFO9dBZQvQHXx/HfXb9QF2E9ha00RCpblujUVq6//3NM/+h73LPHQ7vaWRrvb+fQJRXv7UJoWJ1ntZcS5uWFtkNgdxTi8rc0MvzJafGEqjLUZstw6izCK4kTdujx9j+82CCMh0nh/q/nNHiWOY3vUN7hG2WzwekrFN1o8Olzr8jhJfHtXZXkehdR2481AxOm5Lk//6JmEaAMIoN0NsNSH0e4GNgO90d6BwIZGT1l77w6tzInix2PUXTy+X7ld/Pi90bDJxdlxPLLoedA/y2/y/FDlM7tTuWK0+UbbhZE0+4EaMtQx4y+dLcQ+vsvji5vn0nn/Xr7vf8Wrzim3kRpze7WZORsNXZ9vy6POvnI7xw3VVRkUZsPTf1v/mhZ/9aJp5OtjhTWONs9e6u4MfUPldvEyx430EtM/umtxnhhod3Nbo905G2kM2nvdT9DeayZXPhjOV0gKjzCo9jYKWVRgwQ0bPGEgXR3GYVpsJPnlrvgTgqt4iXjEKOJW7WxD96E4ZUa71hZzT1Mo5gbnfll4hVlmgCncnKM9PNJbF7eSdP5dHbcw7MvbcST5ZUan1k/kZ4PfI9QX/UNhrn/jTZnCrpzS+Fw1r6dabHDbgC/Kov2pLniGJkzNL/ravakZkb6WQLRTcK6NdFJPtPtlw6O9L1ncoL0lGB2bdzrCzxwcFzXEt9QLdHJ5NLTPqJtFKYy5Ug7JbcZFNoyj1Qoolj72j1qXj13H7ZniXfXDPJ39neYmz+8qnsNScnm5PfLsG4TgH0X0SLXDhjgz9Ii808V5OyR9V9xg6pHyyS6vf5aPtuPGbXK+ZFAQiD4Q7VYEnHlDfQ3tzjuAWKC9aO9o7WXEuUZJ85MqC9F22aipiX06L48Q9vmahkW6GDVNmFIYQzZOV3fqX1Uj0dw8Ilt2/oRZkxHgmxOnsfM65hhnHnU/ysuGhUd7PULrbb+oVy2HvBd1nuMcj0qzfqIyZFN7tO59zuVpyjdyviCW9+eoRPSNqpE/JW+3p/NborxTyrXrtOoPRRsN6Ue7rnT/wqPdJVbqH1XNQ3vR3lIPad682xx06pC4E2kyVhaHk18AYp6o5+z54u1RtbmNgqF18bSAVoNY4cGv14iz4sdnxJzORpOnIHwh/yy91v7smEXNRp7bxAw+lH9mYGltg8t5+KJpvxh91Ganiwttlldn7BkjqNwX/Uv7LovreMGtGq9SBBvKNkDtPPL8YbbV8KO8nL9fTntP21lf0tp8qxeRhhzm8XadtLieH2nt9jJ/G/P/GHGEjKPycN2W+BSfy2YXfeV2b9rvEnlOK9ExUof2XJxba1ZN/dD9BO0WBLFIUrtVrov+kbeZ9QPtHXaynFJ77wxjdJrYcYGMi9uqFc9P4q+09kR9L/4+ro2cX7VtA3ITp2P7otTHcHcZPYp5IU5NhVY8X2TiQuPvAL+tpRAwbXukJAw787DRV7SNtj2f1cb8D9r+m5aIK69OZ7G0e3a72vQ3a1uVvw/jKGjGSWncNve0LrhEhMraX7TwCLYXG5m+4egcpa7kMdeuj21n49k3P/6zlK7yZwniR/Hd9q533DwNYRfZ9F1HX+kbvy1enBvxxKMtLmH9CaDdNax0fqDdNVxKXmgv2lvqDu2bGM71fOJitpUxZYPiYlRXwmfD0kaF7+C3ch/owN/2OHjdPN0b1cEGTpPhH4bEWOPExvbFi2Y9yuko93vGWyNaLbe2A6vO7hNm59HaPkaj+T9TOo9S+2L6RNtDbjSUZB6n49rod5lddi9hSGuzn1MeftnJN1M2oD31ZImbgP/kpZmzr4S2hKHXr8LE6iKAdtcTQrvruYQv2hskeq7PrL13ezI6W7S4mNko2cLZuHmqjunR03IZbGhkL1TJPwzNNcvnUb0nPQ7o8l8Y/nkaG0dN6cOQaDJOnKddtM3tnn7FwmFjDaYw1Muci7xX3nA9mvi0FcX9woawxb/Lub6eBrNF/6krm+vrm0FP3RhTpm+ULpi9pu0lPnsYfSP6Sl09hvpFnlf9eWhGxL8gEDyD70XgCjs+h9HuS9B70G6097LN+uydVnvv9KFzwjhxgQxjbm0ENoRsSDSJf5RvtXKpLL4gtZUpK4vieRTTznUonPxtdDd+Kq1U17jwldP602Vf5R51db+Y2lEk7LfRZKj3Sz1TrJyv63bBrWf22U2D8uiT1nGijbLslW7wSG/PcnVG07E9UuxR55iq05mmHEHp/Sk/T9HxMiqPcn4bbNf15w2KcZhDBk+0O29SnRdod0v3zvmgvS2M6oLE7bTae7cOCH7bPr5Xh2yacpAZPAq3oZE5bdvo8WhjXDBuAy5/PW/440uvwXs2fFuniegYNnD9SNDOj8+9drkeaXHZu0ZEn+fxtbp1yiMTfe3FaGQ8is0iKNx1D6P6NtG433gcPy71yFQqv7mYWRiz5uYRqyEjpzaG+xjNLqX7lkeczTV4x4itw1d3KsvD1Q867IBLGGHR3+4PKwqxOwhsylN9Ge3OGyjXGOt2ktqt8qG9aG+HnNQHYzjXcwkj1IZFEi4/yW3sFKNq8ouRQou1y/yWlm+12FD1fmZMKV5haMtvrLNx1/rFAx3Hxt4Qg69aFhso1QufDUl/dcF1t4u2uZGftz1vOoQ5izDy5/nIdJOSqexum0ntozxc9q6bkqycedypN1GT6kzijED0N/dh3HwEgmfwnS/nkTnpnLOBhnYnpt1o78gOvf9koQ2hFYNrhOFcj6xqvNXHWtfXo71+HJ0ZprkYF/+Qpv1sOoPWnk7xWHH9732TDLKonvLJOpjW0eEiaO618w8D+UbH8yi3bwTsbFTbFeHa/kxxihuJLHT4z+iTZ/ihSLFTAkv3+51iSbLYaHepWaSPmb5pvXQfdv6FNqPdpUZgcwqBpfvtqLJhOLdji8ep7bEWDpUIeSqCDeLiMaC2qyOVnucbo70eeZ7zSwnvK78wYJesbWYcqx7xWTXPifYc2Butn2vxZnZh1LZvDsa+EOh8qi4M86r/1b6O7YvRv7R43dcV303um4B4yRAo+obbXsscYl7kmUwtj1UQtPu2PdHuY/Xrs9Wm0MmUtBfDub0bznGBbD9CR6g6i+fw3te68VG8wvwosFxWf6KrvF97FMXxVI9i1Lo20q2nDfbWx/vK60VL+qsgxX/lyvPlXDjX1/WuGv+u0708nW8Upo42l4vwf+Wdtm0d1+WYPC9X+Qxi1lYmwoYREPu6/teUyatNAfgnS6BT/5YuufoY2v0SchLaHcVBe4PE+usjaC+Gc32/CeOsPnQlX3UwT1V4Q+vySHP2OEx+5Xm9NmrLo6/FI7OmouZ5e3qHR3SzUd2muPLv/Dc25TfEEGk6VFzsPJpc9wUP332GUT33C22rnwszMWtiif98BP6aL6sip+jrhQcbsxBAuy8xot2XPLI9tLcGSppeSWqvjYUYTYl1mvi2KdV/bXPYbGqCR5Hf0gleHVW1MV2MxCrcFwqPbPxNS+Hk/1hL7ZcW5O+87X7TYiO10Smuj9dlWDemHxgQj1c9raFuhN3GRtw4YHgMhLt2dLWh28p/SjLKKf0cN2Ojjt2USGWaq98VTziU51zTP5qKfUZ/tBvtPmO/z+qM9rY2/VjtfT1yXX2ULQ68k/X/blHOvNP7ZUC/4Ff91JoN4pjL7OLZaK4bnbXhWWs4y9/zhf2ioQ0bf5rMUztsRNc5f01j7tHduuPYL4ySpmkhzxTHdS1uHJoyGuDvY97Tcn9AGqL2I+Cbrinf2O53lOVj0TeWZzz3EdDu289cot1z96x95If2LthONpz/yvOP9YKH213Wm4ivKHnahY1aG8VVVzVwbSBXR6VbR4hleMY0jzCsbRxX843jthnVEWeutcv1Zal81XxdRs/Lxu2DQPmF1X2UeL1SorfLska79W34lgGRuemj3XMTnZYf2tvMb6z2/hFZMuIcJBJaS+ze6Fscxb2a0iC/Yk50Wz4WVS0eca2drqEw+4dx3ZZVbZjS2/D3aK5dNrIrv/Jo+W1I/qswH6vxeAqv3iBcpB+541HsKOPILEhWJaC2elN+TTdjNwr3DZDn2Ls/fKv9xrgKx0FgFwTUj9HumpYSF7S7hssSXmKN9i4BtpQnhnMJxkk3bai+q5Otbp6lDfDWfwtsYqb8bBT5xcZs2oXWHkW38ZyaizJhPM/bMn6K0fYJQ/cH97slLqjz1uTljVX0lbnzJz8IjCGAdt9SQ7svew/ae8lj9r07s+dIhnsjEMZN3ahz4wuGbZWUMZTNnVacslF0X/vlL3+0ZbFmmEec7S7+yvvWi98JBLoeE3tUpPHpwoTjLpHUfdcu+srtHr8Q2JYA2n3LH+2+7Ido7yWP2fcwnGdHursMw3jxXWrhZPy2Pu4pItZv2AiPaSARI1VDKUYRGbWIlpq47tl3PMXo14mHWit59I3oK2sdl+NAoI0A2n1LJ87PNlanCEN712lmpmqswznZo+hE8xc2PL+0OuLc9binrU6e5/dLRMhPZh8nXkqMoBTWMYoYo4oplGlyGcTaF5PP8oz8BMDuQ/mvYfx5eo7nMLe5Rwps+npKW7otwmJEK/rKFmXgmBC4IOBzWQvandgXkdDei246dSdJ7WXEeWqzHiO9H/l5jrNHhcN5/mnr1zkiYs36Z/mVDVEbcDE6UhN9U68wJI82auHPDH6aL56rbqNvrRHe1j9dUJlsyLu/XdxIad/z4lN00TcuyptiQSnT6Qig3S/fQUil8dHe+VoiSe3FcJ6vgfecUxi12XSN3LAZbSQovQ3uZzaEcmPIhpKN6RRd/EFH2dBPsZxDy/SR2JefIngE2AZt+eZoaJ6d8fNjds1ldxmiz2V55v3kwq/zYOtFcP+1i75yu8cvBLYnEOcM2r19W0QJ0N4gMX2dpPYyVWN6w+4+BxktMR85DC0/Qq/+8cqgeirP4jG8tlMecY4bhDhBB9Uz4cjmX0yXqZZTbeL62sBtqre/iBJsqsnb9j26/WlTBOXpD/NnU0i0HSPMnu9so77xU4VN+a3kHzdVY3isVEQOc0YCOmfQ7mYN26pLoL3zkU9SezGc52vgvefkkYvss3Ra24B+MqZCEnKPJn6vdfY9U61tHPnf/jwXL0UXxlA8EkqxjIPLJN7Vf1e0mGftoDDX1cbtQy0WJn895Wst2SNG7U9xNoCfN2WgMD+NqE4BStVgjmrEzUX0lfBnDYEUCKDdKbRCXgZraaU4aG8FyIDdJLXXhvOreSViPaBOh4/698PX8GUFPVcuRgM9zaLR+HmZpHbL6Wyg+c9PbDzbSPuHliSdyumyZmXTunVubpIV6FEo1cs3M25bG8qZk5/F/EZrt5P/ft3tNMkpDx+j8bvfCn/R5wCK90qfeCvGCTZzGs7o7bINiHYP54t2D2fWmkJahva2EuoMTEl7X4/S2nD+K9+JdYSxvrn594kgxFw5jxBnRtWYuksobFxc/ZvhmLxWTOMy+87Wy5zG0YpVqD+U2sN18vzmh9r2hfEm1t6Wy/6aVX4W+BCpLKDuR/H8Ul+WT02451l+WOOfeSldagZxU1ELf5U5Rjz8BYOmehfxB2ygtwNgjYiKdg+Epv6Ndg9k1hY91w60tw1SS1iC2vtHFPdObLA+N4HcKIjpFN+djEbUO4ykQ1Q/Fx5/WcPGsQ0/j6gXdSxvq8L386Wx7orvOdHZ/OSmSIozp3HZdJg1/YPXoW6o1gTIsZYlkJ9zoWFo97K4e+Weayva24tWY6RktdcjzjgIBAFP17g5oPET9Wta+4sfF1MZmiLuxT8Xbr/gafH2aLKdnySUX9zzdvkLGI8dqc4pDz+J8MXZnMp5ZNEVbv+s/2Qex/kJdvFE5jg1oyZHIoB2J9KaaO9sDZGs9t6ZrYpkdAQCfqmh8VH7ESrYUIcYrXnUEL5Hb3+z2Yaw17H4M0nlEWHXN3uRRf42DD23/Z7WV07+foHPhvjFqHUp4geKU33prxS828238pKn+jnF3YKl4LMSQLtnxTkpM7R3Er4icbLay4hz0UZsyPCxURVG5JmAxGfb4g5393VXW/6tqxKKU7ws6Ljaz76E0pRO4X6R0lMWPLpsQzpz8qs1tiO8bq00fgznEXB/G9nH/Vl+KRre0ScYcVYj4dIkoHMH7U6kadQWaO88bZGs9jLiPE8Dk8uOCeQXnWwOq7ZjXtWOa7Ro0W3c+iXAsntfO0Onafwo1p5G4pGywYZ3+eBLbatsLpf7w9wvBi5VZPKFwKkI6Bz1DcNZtBvtTaR3Yzgn0hAUY3MCMdrZOM938xKmUQAbyG/qglW+wXhP+8Gvs5SK6xHr8pMN51Xe78xjpQjRFxhtXgk4h4HACAKhPXG+jshiF0nQ3kSaCcM5kYagGJsTiJfk9vYpvVXByei1gRvTNW6071HZbMRnYEFilCgb1VU+Y/IYeMjB0WNkfeho+uADkQACEBhN4BTajfaO7h+zJ8Rwnh0pGe6RgETJo4p+7OfRUFw7gfIjQ/+BSuOfntRlI9ZO7+9Bm/UXWlIcbXbRsxGsvLzex0EAAokROJl2o70J9D8M5wQagSIkQyD7BqqEOF5KSKZgiRUkHhl6tNjfiB40lcF8tXysxRcBvxwYI0bJVNNlVGFcP5cRBwEIpE3gLNqN9ibQDzGcE2gEipAMgRg5jUf0yRQspYLIqPQIsUfnPVo8ZorFN8rjXS9K/5rW2SfxtJ2Si/mSvlDhIACBtAmcQrvR3jQ64d00ikEpILA9AYnST1psEHr6wdWffGxfwqRK4BEec7r4pF2fEopxOU2qI7q+efLXNFItXx/UxIHAKQicTLvR3o17tUecPWLk76k+27gsKR3eI2D/1JLq3MuUWB2tLE9UIc+/jRHHo9Vvrvp4hMeG5eHOEdXJ0zS8LDkSbt21xvxTx/PNGm4eAmj3PBz3mMtZtBvtndY7x2qvpyTaVr6x4fxAi/8R7L4W3C0Bj6S9o8UXT9yJCMiI+TKvbnaCnKjqg6oqThaR8sjxoPSJR/4sL58vxEs566415h2xvLfUQU6YL9p9wkZ3lc+i3Wjv5A4+Vns9mGZbOTOcJ5eCDCBwMAI2nj0HF4OmpWHFZ8z85pYctw/K29xzr3/Q9vPtS0QJIACBAQROod1o74AesUBUjzjjIACBSwIx0hgjj5eh7B2ZgEcs7ZjjfsuBXwjsiQDavafWuizrbrQXw/my4diDwE0+0mjD6RNtM+p8rj7hm6Wvjziic65mpLZnJIB277rVd6O9vb+qoQ7peSExCvNI236Z8FP5H+7loF13PQo/CwH16y+1eJ6zT+bo97PkTSZpElB7f6KS+UbpUO2teqHdaXY5SrUAAbR7AagLZ7k37e014pwL71da+08LvPilIBvMv2qbrw8s3KnIfjMCMepswwN3YALSMRvM/i619e0wc5tVF/ddtPvAfZeq1RJAu2uxpOe5R+3tZTgLdXZBKSNXZd0xfYGJD4+Xg9mGwO4JqI/7G77+vNVXu68MFegi8I0i+DveS36CrqsMS4Sj3UtQJc+kCaDdSTdPtXC7096+hrNHlZ+qM3pUpuz8SSp/85YRuTIVtg9DQH3b0zUeaB0vLhymblTkloDa1vrm5b0DMkG7D9ioVKmbANrdzWjrGHvV3r6Gsw3k31XJpkeYVYN66/bg+BCYk4ANqi/U/+nnc1JNJy8/UXivRd/SKenwkqDdw5mR4jgE0O6023KX2nu3D1NdUJpGYrI/CFF44wuCCvNIXRgcr+XHe1K9SGnfcfxY8Xkex6uf5V/85W0pztM8zhtafy9/Xxwyp22XyaMszieO+4H8j/pnDVm9+VmOgPrOb1p8DvhPgsr9c7mDkvMqBNSu1gjPay40ZJUDr3SQvN/WHQ3trrAWBq8AAA5vSURBVKOC36EIqP+j3Ym26J61t5fhXMddlbbweopG4xvoiuM/ErAx/FDb2Z8laO07jH/ZT0vm5Od8ftXyobYzQ1lrG782ij3S7c7v43k+9dv20zpz2v7RYVr8FQRfBD0y+HYeHHH8uB0HgdEE1KcOaViNBnKQhGpX3widqm1VZ7T7IP2XanQTQLu7GW0RY8/a66kar+bQYt2XoY1Y/7tW/EXxkHQ2dG0sh3NeNpCL0WXt2wi284XNLo5XGM233pnhbmPZFwN/Ju+RtiNtHmX0y11/jwxYQwACEJiRwFC9nfHQhZai3XNSJS8IQGAPBMZq7+tROY84/5XvxDrCGtcyTD1qbEO3aQpHllbhNoRjBNnGrA3b6miw/W30Xoh4Ja3DbWj/rOXCKZ5Ho+3n6Rj+rvQzbf+ptUeRftQyxrhXssz9Oza61jqe6+GR9KrR3pbU8yp/K0fQ/ovyPtsQgMA+COjcfWVASXvr7YA8O6OqjGh3DSVxmazfaHcNWLwgsAKBlbT3j6jK4KkaKqDnLN/X+sIAjgyra8WzIHm6xn0t32qpGr8x8vwfhTW5iPO8KYL8bVzbeQrIZ1pimohHo/1PYItO11D+Llsx/UTbo5zyGXLxHXUMEkEAAucjIG1BuxuafQ79Rrsb4OINgYMRuDOkPhIGG6NvaF2MNGvbn+oKw/Yiuzz+/8jTc5U9uurR5+pUi9j3i35NLuLYCG9yHgHPyqG1R55dThuhNpg/0nYY1k3p8YcABCBwSALSP7T7kC1LpSAAgbUJ9Dacc8PzLa2rLwNakD09os55XvJ3SlN++cYjz5nL8/S+pys8uvW9/FWcd7U43CO6V6PcCvNLhHY+lo3ji+/tKtx/aGCDvTZ/+eMgAAEIJEtAGtY2YNBZbqW3LqLdnaSIAAEIQKCbQC/DWcLrkVwbpv6zE/99a7HIb+hf1FrE7XwxyEaItfYI9tWfTOg4nuIRc4D/oe335RfptZs5x/EXNcI4/0zb1QuN9yP8NhW/EIAABBInkGuZ39nwV4cGO6VDuwdTIwEEIACBZgJ95zj7JTsL8MVobp5tGLZ1R8nmG0u8P1Ggp1s4D48Av6bFhrgN8Hh58L+17/nIfnklvtPs8GyahtZ+CdD5eRrGc63tnJ/3wyi2v/P31AytMudjOd+Y7pF7s4IABCCQNgFrXa5dHrTwEtrXt+Bod19SxIMABCDQg0Avw1li3Tb/uPEwSmejupgPXYro6R4XUz7yC0LrC3yKY+O3MY7CbUCHEV06HJsQWI+A+qGfcMQIoftsGDvxtKTs55s/x/dN4sWXZeSHg8CN+oXf1/CUuMHOaQcnUgKlQ7vHgCNN0gTUr9HmpFtoH4XrZTjvoyqUEgLJEPCTGRse/kOfMJptjPwpv2daX8zV176fvhTxtI2DQJWA5yhnT+eqAexDAAK9CaDNvVERsYkAhnMTGfwhMJ6AjZyLJy3aj5Hl72qy9WcaMZxrwOCVjf56lKztc51gggAE+hFAm/txIlYLgTstYQRBAAIDCeQGsueVVt3j3KMuzEb1L9UE7EMgJ+B3NpjGQ3eAwAQCaPMEeCS9IIDhfIGDHQhMJuCvw/gF1aqL6Rl1c/D94isjzlVi7GcEMJrpCBCYhQDaPAtGMsFwpg9AYEYCMnLqDGMfwSPO/pOeKwO5Jc2MJSMrCEAAAucl0KKzaPN5u8WomnuOc3yp4tmoHI6ZyCOGr2vxC144CEwiIMH2VAzPU62b3zwpbxLvloB1958uvfrH1c3Ubmu1fcHR7u3bYDclQJt301RzFnSs9npQLPuqm0ecfVH3t5Pva8HdEvCbt+9oic+HwQUCUwi0zW+eki9p90vAumuNeUcX73v7rUZyJUe7k2uSpAuENifdPIsUbqz2uq/YVr6x4YyDAASWJdA2v3nZI5M7BCAAAQg0EUCbm8jg30gAw7kRDQEQmI2A71Rr5zfPdgQyggAEIACBoQTQ5qHEiM+IM30AAksSKM2ha3ppcMnDkzcEIAABCNQQQJtroODViwAjzr0wEQkCowkwh240OhJCAAIQWIwA2rwY2mNnjOF87PaldtsTyN7CVTEYcd6+LSgBBCAAgSCANgcJ1oMI3B0Um8gQgEAnAT0C/EKRHmjx1xLiyyzfy/+59n/Wmn+BEwgcBCAAgTUJoM1r0j7usTCcj9u21GwjAhLnTzc6NIeFAAQgAIEGAmhzAxi8BxFgqsYgXESGAAQgAAEIQAACEDgrAQzns7Y89YYABCAAAQhAAAIQGEQAw3kQLiJDAAIQgAAEIAABCJyVgA3nV/PKx/qsLOrq/fc6T/wgAAEITCSA3k4E2JEc7e4ARDAETkpgrPa+HrxsOP+V78Q6wljf3PwbCBCAAAQWIIDeLgC1lCXaXYLBJgQgUBAYq71/RA5M1QgSrCEAAQhAAAIQgAAEINBCAMO5BQ5BEIAABCAAAQhAAAIQCAIYzkGCNQQgAAEIQAACEIAABFoIYDi3wCEIAhCAAAQgAAEIQAACQQDDOUiwhgAEIAABCEAAAhCAQAsBDOcWOARBAAIQgAAEIAABCEAgCGA4BwnWEIAABCAAAQhAAAIQaCGA4dwChyAIQAACEIAABCAAAQgEAQznIMEaAhCAAAQgAAEIQAACLQQwnFvgEAQBCEAAAhCAAAQgAIEggOEcJFhDAAIQgAAEIAABCECghcArL168aAk+Z9Dnn3/+QDW/p/Vv5yRArSEAgaUJSF/e9DHQmflIiyXaPR9OcoLAIQlM1d47zkDLi1xwDglpRKWeKs2vYvJ4RFqSQAACEGglYN1VhF+9aPtea2QChxBAu4fQIi4ETkZgrPYq3UdaspFmpmq0dxouaO18CIUABCCQIgG0O8VWoUwQOAABDOcDNCJVgAAE9ktAoxjP91t6Sg4BCEBgnwTGai+Gc317x4Xsfn0wvhCAAAQmEUBbJuFrTIx2N6IhAAIQEIHJ2ovhTD+CAAQgAAEIQAACEIBADwIYzvWQntV74wsBCEAAAgkTQLsTbhyKBoEjEMBwrm/FeNx3rz4YXwhAAAKTCIS2hNZMyozEBYHgGXyLADYgAAEIiEBoQ2jFYCgYzvXIYtTitfpgfCEAAQhMIhDz7EJrJmVG4oJA8ES7CyRsQAACJQKTtRfDuUSztBl3InFnUgpiEwIQgMBkAqEtoTWTMySDjEDwDL5ggQAEIFAmENoQWlEO67WN4VyPKUYt4s6kPha+EIAABMYRiBHR0JpxuZCqSiB4ot1VMuxDAAImMFl7MZzrO1LcicSdSX0sfCEAAQiMIxDa8vu45KRqIIB2N4DBGwIQyAhM1l4M5/qe5L9ttWPU4pYDvxCAwLwEHuTZhdbMm/t5cwueaPd5+wA1h0Abgcnai+FcjzdGgQJwfSx8IQABCIwjEIZdaM24XEhVJRA80e4qGfYhAAETmKy9GM71HSnEN4b062PhCwEIQGAcgTDsQmvG5UKqKoHgiXZXybAPAQiYwGTtxXCu6Uj6//IQ3xttB+SamHhBAAIQGEUgDLtCa0blQqILAmj3BQ52IACBawKTtRfD+Rpq+MQFDcM5iLCGAAQmEyjdjD/X9vPJGZJBlQDaXSXCPgQgcDOX9mI4N3em3/IgDOdmRoRAAALDCYSmhIE3PAdStBFAu9voEAaB8xKYRXsxnJs70M950MPmKIRAAAIQGEzgzTzFT4NTkqAPAbS7DyXiQOB8BGbRXgzn5o4ToxaPmqMQAgEIQGAwgbfyFGHgDc6ABK0E0O5WPARC4LQEZtFeDOfm/vNLHhR3KM0xCYEABCDQn0BoCiPO/ZkNiYl2D6FFXAich8As2ovh3NBh8pd2sjmIpQnlDbHxhgAEINBNQFriN7ofaOHFwG5co2Kg3aOwkQgChyYwp/ZiOLd3lR/y4Mft0QiFAAQg0ItAaAmjzb1wjY6Edo9GR0IIHJLAbNqL4dzeP37Mg99uj0YoBCAAgV4EPshjfdsrNpHGEkC7x5IjHQSOSWA27cVwbukgGtr3qJC/s/puSzSCIAABCPQlkI16SFtiRLRvOuINIIB2D4BFVAicg8Bs2ovh3N1hvnMUCXFMKu9OQQwIQAACFQK5hniOM0Zzhc1Cu2j3QmDJFgJ7IjC39mI4d7f+93mUGObvTkEMCEAAAtcEYo4d0zSu2Szhg3YvQZU8IbA/ArNqL4ZzRwfQnUpM1/ioIyrBEIAABNoI+ObbX9NgxLmN0kxhaPdMIMkGAvsnMKv2Yjj36xBPFO2ehDjuWvqlIhYEIAABEZB2eKqXl68BsioBtHtV3BwMAmkRWEJ7MZx7tLHAf5lH+7hHdKJAAAIQqBL4LPewIYdbiQDavRJoDgOBdAnMrr0Yzv0b28bzuxJiv9yDgwAEINCLQK4Z/jLPD9r2V3pw6xJAu9flzdEgkASBpbQXw7l/88ZIUdy99E9JTAhA4MwE4v2IT88MYcO6o90bwufQENiQwCLai+Hcs0XzkSJf+D7J72J6piQaBCBwcgK+2f5auvH7yTlsUn20exPsHBQCKRBYRHsxnAc0rQTYj/x88WPUeQA3okLgrASkGZ+o7p7exWjzhp0A7d4QPoeGwAYEltReDOfhDRqjzg+GJyUFBCBwFgISbhvMX2j5WNvMbd6+4dHu7duAEkBgcQJLay+G88AmVIP4G6z+pNRXA5MSHQIQOBeBb1Tdn6QZfIIugXZHuxNoBIoAgXUILKq9GM4jGlEC7M/SPdA6Jp6PyIUkEIDAUQlIG/zNdy/vHbWOe6wX2r3HVqPMEOhPYA3tvVsqzlMdsLSbbfrzSQh/lcrtvrn8S3y+08Jj2HpG+ELgrAT8ROo9tCHJ5ke7k2wWCgWBWQhM1l7p9lOVpHE6rg1nv+zmEdQ6x1vgdVTkJ7C/abEA39eC4dzACW8InI2AdMFzmz2v+aez1X0P9UW799BKlBECwwnMqL1+N6XR/T+q2GDJSsK3EgAAAABJRU5ErkJggg==\n", "text/latex": [ "$\\displaystyle \\left[ 2 \\operatorname{acos}{\\left(- \\frac{2 \\pi \\sqrt{\\frac{L}{g}} W\\left(- \\frac{T \\sqrt{e^{- \\frac{T}{\\pi \\sqrt{\\frac{L}{g}}}}}}{2 \\pi \\sqrt{\\frac{L}{g}}}\\right)}{T} \\right)}, \\ 2 \\operatorname{acos}{\\left(- \\frac{2 \\pi \\sqrt{\\frac{L}{g}} W\\left(\\frac{T \\sqrt{e^{- \\frac{T}{\\pi \\sqrt{\\frac{L}{g}}}}}}{2 \\pi \\sqrt{\\frac{L}{g}}}\\right)}{T} \\right)}\\right]$" ], "text/plain": [ "⎡ ⎛ ⎛ ____________ ⎞ ⎞ ⎛ ⎛ \n", "⎢ ⎜ ⎜ ╱ -T ⎟ ⎟ ⎜ ⎜ \n", "⎢ ⎜ ⎜ ╱ ───────── ⎟ ⎟ ⎜ ⎜ \n", "⎢ ⎜ ⎜ ╱ ___ ⎟ ⎟ ⎜ ⎜ \n", "⎢ ⎜ ⎜ ╱ ╱ L ⎟ ⎟ ⎜ ⎜ \n", "⎢ ⎜ ⎜ ╱ π⋅ ╱ ─ ⎟ ⎟ ⎜ ⎜ \n", "⎢ ⎜ ___ ⎜ ╱ ╲╱ g ⎟ ⎟ ⎜ ___ ⎜ \n", "⎢ ⎜ ╱ L ⎜-T⋅╲╱ ℯ ⎟ ⎟ ⎜ ╱ L ⎜T⋅╲╱\n", "⎢ ⎜-2⋅π⋅ ╱ ─ ⋅W⎜────────────────────────⎟ ⎟ ⎜-2⋅π⋅ ╱ ─ ⋅W⎜────\n", "⎢ ⎜ ╲╱ g ⎜ ___ ⎟ ⎟ ⎜ ╲╱ g ⎜ \n", "⎢ ⎜ ⎜ ╱ L ⎟ ⎟ ⎜ ⎜ \n", "⎢ ⎜ ⎜ 2⋅π⋅ ╱ ─ ⎟ ⎟ ⎜ ⎜ \n", "⎢ ⎜ ⎝ ╲╱ g ⎠ ⎟ ⎜ ⎝ \n", "⎢2⋅acos⎜─────────────────────────────────────────⎟, 2⋅acos⎜───────────────────\n", "⎣ ⎝ T ⎠ ⎝ \n", "\n", " ____________⎞ ⎞⎤\n", " ╱ -T ⎟ ⎟⎥\n", " ╱ ───────── ⎟ ⎟⎥\n", " ╱ ___ ⎟ ⎟⎥\n", " ╱ ╱ L ⎟ ⎟⎥\n", " ╱ π⋅ ╱ ─ ⎟ ⎟⎥\n", "╱ ╲╱ g ⎟ ⎟⎥\n", " ℯ ⎟ ⎟⎥\n", "──────────────────⎟ ⎟⎥\n", " ___ ⎟ ⎟⎥\n", " ╱ L ⎟ ⎟⎥\n", " 2⋅π⋅ ╱ ─ ⎟ ⎟⎥\n", " ╲╱ g ⎠ ⎟⎥\n", "────────────────────⎟⎥\n", "T ⎠⎦" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import sympy as sym\n", "sym.init_printing()\n", "#\n", "T = sym.Symbol('T')\n", "L = sym.Symbol('L')\n", "g = sym.Symbol('g')\n", "theta = sym.Symbol('theta')\n", "#\n", "sym.solve(T + 2*sym.pi*sym.sqrt(L/g)*sym.log(sym.cos(theta/2))/(1 - sym.cos(theta/2)), theta)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The $W$ referred to above is the Lambert $W$ function:\n", "\n", "$$\n", "{\\theta _0} = 2\\arccos \\left( - {\\frac{{2\\pi }}{T}\\sqrt {\\frac{L}{g}} {\\rm{LambertW}}\\left( - {\\frac{{T\\exp \\left[ -{\\frac{T}{{2\\pi \\sqrt {\\frac{L}{g}} }}} \\right]}}{{2\\pi \\sqrt {\\frac{L}{g}} }}} \\right)} \\right).\n", "$$" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAGECAYAAADayDLFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABGCklEQVR4nO3dd3iUVfr/8fedEAhVepPehFAlNBtiY7E3XAVFsWHXn6vftayKumvZdd3V1V3LKoIrggVFRcWyioBILwIJTXoNoSckpJ3fH8/AxjBJJmEmU/J5XddcmXnqfWaSuXOe85xzzDmHiIiIxJa4cAcgIiIiwacELyIiEoOU4EVERGKQEryIiEgMUoIXERGJQUrwIiIiMUgJXioNMxtkZs7MGh7jcZaZ2eNBCivQczY1s6/NLNPMwt631czWm9n94Y4DwMweN7Mdvs92ZLjjKaroexXK987M7jez9aE4tkQfJXiJOGY21vdl7cws18zWmtlfzazmMR56FtAM2BWEMCva/UBzoBdeGSqEL3ku87OqL/CvioqjOGbWDRgN3Ir3vrwX3ogC8qv3zvd7PjSM8UiMqhLuAESK8S0wAkgATgPeAGoCt5XnYGaW4JzLAbYHLcKK1QFY4JxbHe5AAJxzO8Mdg08H38/JLkpG7Yqg905inGrwEqkOOee2O+c2OefeBcYDlwCY5/dm9ouZZZnZUjO75vCOZtbGVysaZmbfmVkWcIu/S/Rmdplv/0NmtsnM/mBmVmh9YzP7xHeeDWZ2Q2mBm1l73z7bfZfUF5rZBUW2uczMfvYdd7eZ/WBmTYo53nrgYuBaX/xjfcuPqvn5uRzszGyUmX3gi2Vt4ffKt01zMxtvZrvM7KCZLTazM3yXu0cDXQtdURlZzHlamdnHZnbA9/jIzFoUWv+4r2njKt/ndsDMJpfWXGJm3c3s20Lv01gzO+7wMYGPfZsWlNR0YWbPmtlK33HWm9lfzCzRT3zX+dZnmNlbZlbVzG73/W7sMrO/mVlcof3W+/Z9x7fPdivl8nvh987+dzn9A9/7u75wPEX2G2lmGUWW/d53zgwzexuo5ed815tZipllm9kqM7u3cBkkdulDlmiRhVebB/gTcCNwB5AEPAO8ZmbnF9nnGbxLoUnA5KIHNLNk4APgI6A78CDwEHBnoc3G4tUSz8b7B+NaoE0psdYCvgTOAXoCk4CPzKyz77xNgYnAOKALMBD4TwnH64t3ReN9vMvQ95Ry/qIeAz7xxfIeMMbMWvtiqQn84CvTpXjvw5O+/d4DngdW+s7r9xK47x+iyUAT4EzgDLzmhMmF/1nyneNK33kGAycCTxUXtJnVAKYCGUA/334nA2N8m/wVuNn3/HB8xckEbsB7v28HrgL+UGSbNnj/SF0AXA5cgfe+9fXFexNwly+Own4HpAK98f4hetrMLishlsL6+n7e7Iu/bwnb/oqZ/Rbvb2G079wrfbEU3uZm4Gm834EuwH3AA3jvgcQ655weekTUAy+pTin0uh+QjpdcauIl+9OK7PMC8IXveRvAAfcV2WaQb3lD3+vxwHdFtnkc2Ox73sm3/SmF1rcG8oHHy1im2cAjvue9fcdtXYb9pwBjiyxzwNAiy9YD9xfZ5plCr6sAB4FrfK9vBg4cfk/8nPdxYJmf5UfOg/ePTD7QptD6dkABcHah42QDxxXa5g/AmhLKfDOwD6jt5zPs4Hs91PsaK/Pv2K2Fz+2LL6tIfB8CO4GqhZZNA14u8j58U+TYbwAzS/hM/H1GRT/Ho953YCSQUej1LODfRbb5Flhf6PVGYESRbf4fkFLW90yP6HuoBi+RaojvsmM28BMwHa/2lAQkAlN96zN8ly1vA9oXOcb8Us7RBfixyLKZwPFmVse3vgCYe3ilc24DsLWkg5pZTd8l4BQz2+OLrw/QyrfJErwv4mVmNsnMbjOzRqXEeix+PvzEOZeHl7Qa+xadCPzsnEs/huN3AbY659YXOs9avPcpqdB2G5xz+wq93loojuKO+7Nz7kChZbPwPpMk/7v4Z2ZDzWzm4cvZwN/53+dx2MYi8e0AVjnv3o3Cy4rG/JOf12WKr5y6FHNuAHy/Uy3xrm4V/lt5lqP/ViQG6SY7iVTTgVFALl7yyAUws7a+9Rfi1U4Kyy3yOrOUcxhe7ckf51tfHn8FhuDd+b4ar8b8NlAVwDmXb2aDgQF4l35vBJ4xs9Odc0vKcB5/MSb42a7o++L4X/NcectYWGnvYyBxHMtxS2RmA/CaRJ4A7gX2AhfhfU6F+YvP37L4QM99DAoI7LMtyeH39la8f4ykklENXiLVQefcGufchsPJ3ScFOIR3eXtNkceGMp4jBTi1yLJT8S7RH8BrV42jULuombXCa18uyanA2865Sc65n4HNFKkxOc9PzrknfMffitc+XRY7KdTubN5NemXtQrcQ6FHCzW45lJ7QUvCuerQpFEs7vPcppYzxFD1uTzOrXWjZyXifSWoZjnMKsMU590fn3Dzn9URofQxxFTXAz+uyxJfL0e/xTqBJkXsYehXZJrWYcwPgnNsBbAHa+/lbWVOG+CRKqQYvUcU5d8DM/gr81fflNx3vprYBQIFz7vUyHO55YJ7vbux38RLtfcDDvnOtNLOpeJc4R+G10f7N97Mkq4BLzewTvC/v0XjNCsCRGuXZwFd4l3xPxLuUWtZk+B1wh5nNwmsDfxqvnbss3sW7uXCymT2E989Id+CAc+57vPbi1mbWG++KyQHn3KEix/gWr9lhvJndjVfzfAnvn4fvyhhPYePxat1vm9ljQD3gNeCjMiaoVXj/gFyNdwn7N8CwY4irqAG+9+5DvHsErgWuLsP+64GzzOwHvN4je/Da+usDD5vZRN9xi/aVfxHvvZnn234o0B/YXWibx4GXzGwv8AXeVYDewPHOuWfKEKNEIdXgJRo9ivfFdT+wHPgG767ndWU5iHNuId6d0pcDy/DaJp8FXi602Ujfcb8DPsNLiOtLOfTvgDRgBt7d9LN9zw/bh1ernIJ3Cf954I/OuXfKEj/ePyNr8b7cP8S7uSutLAdwzmUCp+PV9D7Dez+f4H+XwCfhJYb/4tUqj0qMzjmH18Ngpy+W7/HGG7jEt65cnHMH8ZJxHbz7ID7BS9CldlUscpzPgOfwbsT8Ge+mwMfKG5cffwN6AIvw7mp/zDn3YRn2vw+v58Em3zFwzqXi3VcyqlDMTxfeyTn3Ht7fwVO+/br7Yim8zRt479cIvH/CZviOWaa/FYlOdgx/fyIilZqv3/rLzrmi7fkiYacavIiISAxSghcREYlBukQvIiISg1SDFxERiUFK8CIiIjEopvrBN2zY0LVp0yZox1u9ezUd63cM2vHCLZbKE0tlAZUnksVSWUDliWTlKcuCBQvSnXP+h7oO92D4wXwkJye7YEp+LbjHC7dYKk8slcU5lSeSxVJZnFN5Ill5ygLMd5psRkREpPJQghcREYlBSvAiIiIxKKZusvMnNzeXzZs3k51d1jk44C+9/kJqalkmhYpssVSeomVJTEykRYsWJCSUdUZNEZHYFPMJfvPmzdSuXZs2bdrw65kXS+d2Oro06hKiyCpeLJWncFmcc+zatYvNmzfTtm3bUvYUEakcYv4SfXZ2Ng0aNChzcpfoYWY0aNCgXFdpRERiVcwneEDJvRLQZywi8muVIsGHW61atY5a9uqrr/L2228H/VwzZsyga9eu9OrVi6ysrKAf35/169fTrVu3CjmXiIgEJubb4CPVrbfeGpLjjh8/nvvvv5/rr78+oO3z8/OJj48PSSwiIhI+qsGHyeOPP85f//pXAAYNGsS9997LwIED6dKlC/PmzeOyyy6jY8eOPPLII0f2ueSSS0hOTqZr1668/vrrRx3zjTfe4P333+fJJ5/k6quvxjnH//3f/9GtWze6d+/Ol5O/BGDatGmcccYZDB8+nO7duzNt2jROP/10fvvb39KpUycefPBBxo8fT79+/ejevTu//PILACNHjuTDDz88cj5/VybWr1/PaaedRu/evenduzezZs0K6vsmIiKBqVQ1+Cc+W07K1v0Bb38w9yA1EvaVuE1S8zqMvrDrsYZG1apVmT59Oi+++CIXX3wxCxYsoH79+rRv3557772XBg0aMGbMGOrXr09WVhZ9+/bl8ssvp0GDBkeOcdNNNzFz5kwuuOAChg4dyqRJk1i8eDFLliwhPT2dXsm9GH7BcADmzp3LsmXLaNu2LdOmTWPJkiWkpqZSv3592rVrx0033cTcuXN58cUXeemll3jhhRcCKkfjxo355ptvSExMZPXq1QwbNoz58+cf8/sjIiJloxp8hLjooosA6N69O127dqVZs2ZUq1aNdu3asWnTJgD+8Y9/0LNnTwYMGMCmTZtYvXp1icecOXMmw4YNIz4+niZNmtD3pL7MmzcPgH79+v2qS1nfvn2PnLN9+/YMHjz4SDzr168PuBy5ubncfPPNdO/enSuuuIKUlJSyvA0iIjFr6rJtHMzJq7DzVaoafFlr2ik7U0hqlBSiaH6tWrVqAMTFxR15fvh1Xl4e06ZN49tvv+Wnn36iRo0aDBo0qNRuYd48BP7VrFnT7/mLxnD4/ABVqlShoKDgyLFzcnKOOu7f//53mjRpwpIlSygoKCAxMbHEGEVEKoP563dz6zsLuffsTtxzdsXMfqcafJTYt28f9erVo0aNGqxYsYLZs2eXus/AgQN57733yM/PZ+fOncyfPZ9+/fqVO4Y2bdqwYMECAD755BNyc3P9xtmsWTPi4uL4z3/+Q35+frnPJyISC3LyCnjoo6UcX7c6N51WcYNxKcFXgIMHD9KiRYsjj7/97W9lPsaQIUPIy8ujR48ePProowwYMKDUfS699FJ69OhBz549OfPMM7nvsfto2rRpeYoAwM0338wPP/xAv379mDNnzlFXAQBuv/12xo0bx4ABA1i1apXfbUREKpN/z1jL6rQMnry4KzWrVeCF8+LmkY3Gh7/54FNSUso8v+5hy9OWl3vfSBRL5fFXlmP5rMMtlua0di62yhNLZXFO5alo63ZmuE5/+MLd9s78UrfVfPAiIiJRwDnHo58sIyE+Lii9rcpKCV5ERCQEPlm8lRmr0/n9kBNoUqfibzhWghcREQmyvQdz+OOUFHq1rMvV/VuHJYZK1U1ORESkIjz75Qr2ZuXyn0u7Ex8XnsmwVIMXEREJornrdjNx3iZuOrUtSc3rhC0OJXgREZEgOZSXz8Mfe33eK2pAm+IowVeAp556iq5du9KjRw969erFnDlzAG/s+GAN5epv4pdwKDohjT9jx47lzjvvrKCIREQqzus/rGVNWgZ/uqQbNaqGtxVcbfAh9tNPPzFlyhQWLlxItWrVSE9PPzLE6xtvvBHm6EREJFjWpWfy0vdrOL9HM87o3Djc4agGH2rbtm2jYcOGR8Z2b9iwIc2bNwe8aWIPz7RWq1YtHnjgAZKTkzn77LOZO3cugwYNol27dnz66aeAV/O9+OKLGTJkCCeccAJPPPGE33M+99xz9O3blx49ejB69Gi/29SqVYv77ruP3r17c9ZZZ7Fz504AfvnlF4YMGUJycjKnnXYaK1asALya+d13383JJ59Mu3btjtTSnXPceeedJCUlcf7555OWlnbkHG3atCE9PR2A+fPnM2jQoKPiKG4K2tKmsN24bmMA776ISMVwzvHI5KVUi49j9AUVM4dJaSpXDf7LB2H70oA3b52bCQmlDLXatDuc+2yxqwcPHsyTTz5Jp06dOPvss7nyyis5/fTTj9ouMzOTQYMG8ec//5lLL72URx55hG+++YaUlBSuu+66I7PNHZ7mtUaNGvTt25fzzz+fPn36HDnO119/zerVq5k7dy7OOS666CKmT5/OwIEDjzpf7969ef7553nyySd54oknePnllxk1ahSvvvoqHTt2ZM6cOdx+++189913gPfPysyZM1mxYgUXXXQRQ4cO5eOPP2blypUsXbqUHTt2kJSUxA033BDoW1yikqawHf/meIb0GxKU84iIHKvJi7fw45pd/PGSbjQOQ593fypXgg+DWrVqsWDBAmbMmMH333/PlVdeybPPPsvIkSN/tV3VqlUZMsRLWN27d6datWokJCQcNV3rOeecc2QO+Msuu4yZM2celeC//vprTjzxRAAyMjJYvXr1UQk+Li6OK6+8EoBrrrmGyy67jIyMDGbNmsUVV1xxZLtDhw4deX7JJZcQFxdHUlISO3bsAGD69OlHpqRt3rw5Z5555jG+Y/9zeApb4KgpbD+d+mnQziMiciz2ZObwxympXp/3fq3CHc4RlSvBl1DT9mdDkKaLjY+PZ9CgQQwaNIju3bszbty4oxJ8QkICZl5fyeKmawWObFPca+ccDz30ELfcckuZYjQzCgoKqFu3LosXL/a7TeEpZV2hqWiLxnBY4elli5vatqQpaEuawlaz1IlIpHjmy1T2Z+XyzGXdiQtTn3d/1AYfYitXrmT16tVHXi9evJjWrcs/qtE333zD7t27ycrKYvLkyZxyyim/Wv+b3/yGMWPGkJGRAcCWLVt+1S5+WEFBwZG273fffZdTTz2VOnXq0LZtWz744APAS7hLliwpMZ6BAwcyceJE8vPz2bZtG99///2RdYWnl500aZLf/QOZglZEJFLNXruL9+dv5qbT2tGlWfj6vPtTuWrwYZCRkcFdd93F3r17qVKlCh06dOD1118v9/FOPfVURowYwZo1axg+fPivLs+D1+afmprKSSedBHhNBO+88w6NG//6js6aNWuyfPlykpOTOe6443jvvfcAGD9+PLfddht/+tOfyM3N5aqrrqJnz57FxnPppZfy3Xff0b17dzp16vSr+wtGjx7NjTfeyNNPP03//v397n/zzTdz8cUX069fP8466yxNLysiUeNQXj5/+HgpLepV556zwtvn3a/ippmLxkesTxf71ltvuTvuuKPc+xcuT82aNYMRUthoutjIFkvliaWyOKfyBNML36xyrR+Y4r5bsSMox9N0sSIiImG2dmcG//x+DRf0aMYZJ4S/z7s/ukQfRUaOHHnUzXnldbiNXkREysY5xx8+Xka1hDgeuzAy+rz7oxq8iIhIGXy0cAs/rd3Fg+d2pnHtyOjz7k+lSPCuUJcuiU36jEWkIuzOzOFPn6fQu1VdhvWNnD7v/sR8gk9MTGTXrl1KADHMOceuXbtITIzc/6RFJDY8/UUqB7LzeDrC+rz7E/Nt8C1atGDz5s1Hxlovi+0HtmPpkf0BlkUsladoWRITE2nRokUYIxKRWPfTL7v4cMFmbhvUns5NI6vPuz8xn+ATEhJo27ZtufYd8foI5o+aH+SIwieWyhNLZRGRyHe4z3ur+jW4+8wI7PPuR8wneBERkWP1yrRfWJueybgb+lG9any4wwlIzLfBi4iIHIs1aRn86/tfuKhnc07v1Cjc4QRMCV5ERKQYXp/3pSQmxPFohMzzHigleBERkWJ8uGAzc9bt5qHzutCodrXSd4ggIUvwZjbGzNLMbFkJ2wwys8VmttzMfii0fL2ZLfWt051UIiJS4XZn5vD0F6n0aV2PK/u0DHc4ZRbKm+zGAi8Db/tbaWZ1gX8BQ5xzG82s6GC+Zzjn0kMYn4iISLGe+jx6+rz7E7IavHNuOrC7hE2GAx855zb6tj960nIREZEwmLUmnUkLN3PL6e3o1KR2uMMpFwvlCG9m1gaY4pzr5mfdC0AC0BWoDbzonHvbt24dsAdwwGvOuWInUDezUcAogMQGicldn+4atPhT01Pp0rBL0I4XbrFUnlgqC6g8kSyWygIqTyBcQRX2bLkNnFGvxb+wuLygHr845SnLglsWLHDO9fG7srh5ZIPxANoAy4pZ9zIwG6gJNARWA51865r7fjYGlgADAzmfv/ngj4XmTY5csVQW51SeSBZLZXFO5QnE81+vdK0fmOKmr0oL+rFLEkvzwW8GpjrnMp3X1j4d6AngnNvq+5kGfAz0C1uUIiJSaaxJy+CVaWu4pFdzTusYPX3e/Qlngv8EOM3MqphZDaA/kGpmNc2sNoCZ1QQGA8XeiS8iIhIMBQWOhz9eSo2qVXgkyvq8+xOyu+jNbAIwCGhoZpuB0Xht7jjnXnXOpZrZVOBnoAB4wzm3zMzaAR+b2eH43nXOTQ1VnCIiIuD1eZ+7bjfPXtadhrWiq8+7PyFL8M65YQFs8xzwXJFla/FdqhcREakI6RmHeOqLVPq1qc9vo7DPuz8ayU5ERCq9pz9P5WBOHk9d2i0q+7z7owQvIiKV2o9r0vlo0RZuPb09HaO0z7s/SvAiIlJpZed687y3aVCDO87oEO5wgkrzwYuISKX1z+/XsH7XQcbf1J/EhOiY5z1QqsGLiEiltHrHAV794RcuO/F4TunQMNzhBJ0SvIiIVDqH+7zXrFaFP5wfO0P3FqYELyIilc778zcxb/0eHj63Cw1ioM+7P0rwIiJSqew8cIinv0ilX9v6XNGnRbjDCRkleBERqVSe+jyFrNx8nr60O75RU2OSEryIiFQaM1bvZPLirdw2qAMdGtcKdzghpQQvIiKVQnZuPo9MXkbbhjW5fVD7cIcTcuoHLyIilcJL361mw66DvBuDfd79UQ1eRERi3qodB3jth7Vc1vt4To7BPu/+KMGLiEhMKyhwPPzRUmonVuGR86N/nvdAKcGLiEhMe2/+JuZv2MPD53Whfs2q4Q6nwijBi4hIzEo7kM0zX6TSv219hibHbp93f5TgRUQkZv1pSirZuQU8FeN93v1RghcRkZj0w6qdfLpkK7cNah/zfd79UYIXEZGYk5WTzyOTl9KuYU1uPyP2+7z7o37wIiISc/7x3Wo27c5iws0DqFYl9vu8+6MavIiIxJQV2/fz7+lrGZrcgpPaNwh3OGGjBC8iIjGjcJ/3h8+LzXneA6UELyIiMWPCvI0s3LiXR85PqlR93v1RghcRkZiQdiCbZ79cwUntGnBZ7+PDHU7YKcGLiEhMePKzFA7lFvDUpd0qXZ93f5TgRUQk6n2/Mo0pP2/jjjM60K5R5evz7o8SvIiIRLX0jEM88OHPtG9Uk1sHtQt3OBFD/eBFRCRqOWfc+95i9mXlMu6GfpW2z7s/SvAiIhK1svaexow96TxzWXe6NKsT7nAiii7Ri4hIVJq9dhcH95zBxb2ac1XfluEOJ+IowYuISNRJzzjE3RMWEZ+wu1LOFBcIJXgREYkqBQWOe99bzN6sXGo3fp9a1dTa7I8SvIiIRJVXfviFGavTefzCrlSptiPc4UQsJXgREYkac9bu4vmvV3JRz+YM66d295IowYuISFRIzzjE3RMX0bpBTZ6+TO3upVGCFxGRiHe43X3PwVxeHn6i2t0DoAQvIiIR73C7++gLk+ja/LhwhxMVlOBFRCSizV23m+e/XsmFPZszvF+rcIcTNZTgRUQkYu3KOMRdExbSqn4NntYscWWiBC8iIhGpoMBx7/tLfO3uvamdmBDukKKKEryIiESkV6f/wvRVO3nsgiS6Ha9297JSghcRkYjjtbuv4oIezbi6v9rdy0MJXkREIsou3zjzLetV5xn1dy83dSQUEZGIUVDg+N37S9idmcNHt5+sdvdjoBq8iIhEjNemr+WHVTt59EK1ux8rJXgREYkI89bv5q9fr+T8Hs24Ru3ux0wJXkREwm53Zg53vbuIFvWq86za3YNCbfAiIhJWXrv7YrW7B1nIavBmNsbM0sxsWQnbDDKzxWa23Mx+KLR8iJmtNLM1ZvZgqGIUEZHwe33GWqat3MmjF3RRu3sQhfIS/VhgSHErzawu8C/gIudcV+AK3/J44J/AuUASMMzMkkIYp4iIhMn89bt57quVnN+9GdcMaB3ucGJKyBK8c246sLuETYYDHznnNvq2T/Mt7wescc6tdc7lABOBi0MVp4iIhMfuzBzumuC1uz9zudrdg82cc6E7uFkbYIpzrpufdS8ACUBXoDbwonPubTMbCgxxzt3k224E0N85d2cx5xgFjAJIbJCY3PXprkGLPzU9lS4NuwTteOEWS+WJpbKAyhPJYqksEDnlcc7Yv304uVntqHv8G1Sptq1cx4mU8gRDecqy4JYFC5xzffyudM6F7AG0AZYVs+5lYDZQE2gIrAY64V2qf6PQdiOAlwI5X3Jysgum5NeCe7xwi6XyxFJZnFN5IlkslcW5yCnPq9PWuNYPTHHjZq07puNESnmCoTxlAea7YnJiOO+i3wykO+cygUwzmw709C1vWWi7FsDWMMQnIiIhsGDDbv7y1UrO696UEWp3D5lw9oP/BDjNzKqYWQ2gP5AKzAM6mllbM6sKXAV8GsY4RUQkSPZk5nDnu4s4vm51nr28h9rdQyhkNXgzmwAMAhqa2WZgNF6bO865V51zqWY2FfgZKMC7LL/Mt++dwFdAPDDGObc8VHGKiEjFKChw3PfBEnZl5DDptpOpo/7uIRWyBO+cGxbANs8Bz/lZ/gXwRSjiEhGR8Pj3jLV8tyKNJy7qSvcW6u8eahqqVkREQu5wu/u53Zpy7Ulqd68ISvAiIhJSe3zjzB9ftzp/Hqp294qisehFRCRkDre7p6vdvcKpBi8iIiHzxkyv3f0P53dRu3sFU4IXEZGQWLBhD3+ZupIhXdXuHg5K8CIiEnR7D+Zw17sLaVY3Ue3uYaI2eBERCSrnHPe9v4SdGYeYdNvJHFdd7e7hoBq8iIgE1Rsz1vHfFWn84bwu9GhRN9zhVFpK8CIiEjQLN+7hz1NXMKRrU647uU24w6nUlOBFRCQovHb3RTQ9Tu3ukUBt8CIicsycc9z/wRLSDmTz4a1qd48EqsGLiMgxe3PmOr5NTePh87rQs2XdcIcjKMGLiMgxWrRxD89+uYLfdG3CSLW7RwwleBERKbe9B7353Zsel8hfLu+pdvcIojZ4EREpF6/d/WfSDmTzwa0nc1wNtbtHEtXgRUSkXLx29x08dG4XeqndPeIowYuISJkt3rSXP09dweCkJlx/SptwhyN+KMGLiEiZ7DuYyx3jF9KkTiLPDVW7e6RSG7yIiATMOcf9Hy5Ru3sUUA1eREQCNubH9XyTsoMH1e4e8ZTgRUQkIIs37eXZL1M5J6kJN6jdPeIpwYuISKn2HczlzncX0rh2In9Vu3tUUBu8iIiUyDnH/324hO37svng1pPU7h4lVIMXEZES/f2bVXydsoMHz+3Mia3qhTscCZASvIiIFOs/szfwj+/WcGWfltx4attwhyNlUKYEb2Y1zSw+VMGIiEjkmLpsG499soyzuzTmqUu7qd09ypSY4M0szsyGm9nnZpYGrAC2mdlyM3vOzDpWTJgiIlKR5qzdxd0TF3Niy7q8NKw3VeJ1wTfalPaJfQ+0Bx4CmjrnWjrnGgOnAbOBZ83smhDHKCIiFWjF9v3c9PZ8WtarzpvX9aV6VV24jUal3UV/tnMut+hC59xuYBIwycx0O6WISIzYsjeLkWPmUaNqPG/f2J96NauGOyQppxJr8IeTu5kNMLPah5ebWW0z6194GxERiW57MnO49s05ZObkMe6Gfhxft3q4Q5JjEGijyitARqHXmb5lIiISA7Jy8rlx3Dw27cni39f2oXPTOuEOSY5RoAnenHPu8AvnXAEaJEdEJCbk5Rdw14RFLNq0lxev7MWAdg3CHZIEQaAJfq2Z3W1mCb7HPcDaUAYmIiKh55zjkcnL+DZ1B09e1JVzuzcLd0gSJIEm+FuBk4EtwGagPzAqVEGJiEjF+Pu3q5k4bxN3ndmBESe1CXc4EkQBXWZ3zqUBV4U4FhERqUDvzN7AP/67mt/2acHvzukU7nAkyAKqwZtZJzP7r5kt873uYWaPhDY0EREJlanLtvPYJ8s4q3Njnr60u0api0GBXqL/N95gN7kAzrmfUY1eRCQqzV23m7snLqJny7q8PFyj1MWqQD/VGs65uUWW5QU7GBERCa28nMbcNG4eLepVZ4xGqYtpgSb4dDNrDzgAMxsKbAtZVCIiEnRb9maxf9s1VK8az9s39NModTEu0L7sdwCvA53NbAuwDtAY9CIiUWLvwRyuGzMX56oy7oZ+tKhXI9whSYgFehf9WuBsM6sJxDnnDoQ2LBERCZbs3HxuHDefjbsOUqfJBDo3HRrukKQCBHoX/T1mVgc4CPzdzBaa2eDQhiYiIscqL7+AO99dxMKNe3jhql4kVN8Q7pCkggTaBn+Dc24/MBhoDFwPPBuyqERE5Jg553j0k/+NUneeRqmrVAIei9738zzgLefckkLLREQkAr3w7WomzN3EnWdolLrKKNAEv8DMvsZL8F/5po4tCF1YIiJyLN6ZvYEXfaPU3TdYo9RVRoHeRX8j0AtY65w7aGYN8C7Ti4hIhDk8St2ZGqWuUiuxBm9mbcCbHtY5t9A5t9f3epdz7mfztAh9mCIiEoh56/83St0/NUpdpVZaDf45M4sDPgEWADuBRKADcAZwFjAab4Y5EREJo5XbD3DjWG+Uujc1Sl2lV2KCd85dYWZJwNXADUAzvK5yqcAXwFPOuWx/+5rZGOACIM05183P+kF4/zis8y36yDn3pG/deuAAkA/kOef6lLVgIiKVyda9WVw3Zi6JCd4odfU1Sl2lV2obvHMuBfhDOY49FngZeLuEbWY45y4oZt0Zzrn0cpxXRKRS2Xswh2vHzCXzUB7v33qSRqkTIPC76MvMOTcd2B2q44uIiDdK3U2+Uepev7YPXZrVCXdIEiHMORe6g3s36U0p4RL9JLz2+63A/c655b5164A9eJPbvOace72Ec4wCRgEkNkhM7vp016DFn5qeSpeGXYJ2vHCLpfLEUllA5YlkkVwW5+I4sONKcg52onbjD6hWK6XUfSK5POURS+UpT1kW3LJgQbHN2M65kD2ANsCyYtbVAWr5np8HrC60rrnvZ2NgCTAwkPMlJye7YEp+LbjHC7dYKk8slcU5lSeSRWpZCgoK3IOTfnatH5jixv64LuD9IrU85RVL5SlPWYD5rpicWGIbvJn1Lmm9c25h4P9nHLXv/kLPvzCzf5lZQ+dcunNuq295mpl9DPQDppf3XCIiscYbpW4jd5zRnutObhPucCQClXaT3fO+n4lAH7zatAE9gDnAqeU9sZk1BXY455yZ9cO7H2BX4RnrfM8HA0+W9zwiIrFm/BxvlLorkltw/+ATwh2ORKjSusmdAWBmE4FRzrmlvtfdgPtL2tfMJgCDgIZmthmvv3yC77ivAkOB28wsD8gCrvIl+ybAx76Rl6oA7zrnppa7hCIiMeSr5dt5dLI3St0zl2mUOileoEPVdj6c3AGcc8vMrFdJOzjnhpWy/mW8bnRFl68FegYYl4hIpTFv/W7unrCIHi3q8vLwEzVKnZQo0ASfamZvAO/g3dl+Dd5gNyIiUgFW7fBGqTu+bnXGjOxLjaqBfn1LZRXob8j1wG3APb7X04FXQhKRiIj8SuFR6sZplDoJUEAJ3jmXbWb/BL7Fq8GvdM7lhjQyERFh38Fcrhszl4zsPN675SRa1tcodRKYgBK8b1CaccB6vLvoW5rZdc4brU5EREIgOzefm96ex4ZdBxl7Q1+SmmuUOglcoJfonwcGO+dWAphZJ2ACkByqwEREKrO8/ALumrCI+Rv28PKw3pzcvmG4Q5IoE+gtmAmHkzuAc24Vvi5vIiISXM45Hv1kOd+k7GD0BUmc36NZuEOSKBRoDX6+mb0J/Mf3+mq8+eFFRCTIXvyvN0rd7YPaM/KUtuEOR6JUoAn+NuAO4G68NvjpwL9CFZSISGX17pyNvPDtaoYmt+D/fqNR6qT8Ar2L/hDwN99DRERC4Ovl23lk8lLOOKGRRqmTY1baZDPvO+d+a2ZL8brH/YpzrkfIIhMRqUTmr9/NXRMW0b1FXf55dW8SNEqdHKPSavCHB7a5INSBiIhUVqt2HOAG3yh1b2mUOgmS0iab2WZm8cCbzrmzKygmEZFKI3Xbfq4dM5dqGqVOgqzUa0DOuXzgoJkdVwHxiIhUGvPW7+a3r/1EvBnv3tRfo9RJUAV6HSgbWGpm3wCZhxc65+4OSVQiIjHuv6k7uH38Qo6vW523b+xHi3pK7hJcgSb4z30PERE5RpMWbOb3k36ma/M6vDWyLw1qVQt3SBKDAu0mN87MqgOtCo9oJyIiZfPGjLX86fNUTunQgNdG9KFWNd1QJ6ERUD8MM7sQWAxM9b3uZWafhjAuEZGY4pzjz1NX8KfPUzmve1PGjOyr5C4hFWhHy8eBfsBeAOfcYkDjJ4qIBCAvv4CHPlrKK9N+YXj/Vrw0rDfVqsSHOyyJcYH++5jnnNtXZFSlowa+ERGRX8vOzeeeiYv4avkO7j6zA/ee00kj1EmFCDTBLzOz4UC8mXXEG5N+VujCEhGJfgeycxn19gJ+WruL0Rcmcb0mjpEKFOgl+ruArsAh4F1gH/D/QhSTiEjUS884xLB/z2be+t28eFUvJXepcKWNRZ8I3Ap0AJYCJznn8ioiMBGRaLVp90FGvDmH7fuz+fd1fTjjhMbhDkkqodIu0Y8DcoEZwLlAF1RzFxEp1srtBxjx5hwO5RUw/qYBJLeuF+6QpJIqLcEnOee6A5jZm8Dc0IckIhKdFmzYzfVvzaN61Xg+uPUkOjWpHe6QpBIrLcHnHn7inMvTnZ8iIv59vyKN28YvoNlx1Xn7hn4aV17CrrQE39PM9vueG1Dd99oA55yrE9LoRESiwORFW7j/gyV0blabsdf3o6GGnpUIUNp0sRqJQUSkBGNmruPJKSmc3L4Br41IpnZiQrhDEgEC7wcvIiKFOOd4/utVvPz9GoZ0bcoLV/UiMUF1IokcSvAiImWUX+B4ZPIyJszdyLB+LfnTJd2Jj9M9ShJZlOBFRMrgUF4+/2/iYr5ctp07zmjP/YNP0NCzEpGU4EVEApRxKI9Rb89n1i+7ePSCJG48VaPTSeRSghcRCUB6xiGuf2seKdv287ff9uSy3i3CHZJIiZTgRURKsXnPQa59cy5b92Xx72uTObNzk3CHJFIqJXgRkRLk5TTi8ldmkZWTzzs39qdPm/rhDkkkIIHOJiciUuks2LCHfVtvwDl4/9aTlNwlqqgGLyLix7SVadz2zkIsLotJt52voWcl6qgGLyJSxCeLt3DTuPm0a1STus3fVHKXqKQELyJSyNgf13HPxMUkt67HhFEDiKuSGe6QRMpFl+hFRPCGnv37N6v4x3drGJzUhH8MO1FDz0pUU4IXkUovv8Dx2CfLGD9nI1f2aclTl3ajSrwucEp0U4IXkUrtUF4+v3tvCZ8v3cZtg9rz+99o6FmJDUrwIlJpZRzK49b/LGDmmnQeOb8LN53WLtwhiQSNEryIVEq7M3O4/q25LNu6n+ev6MnlyRp6VmKLEryIVDpb9mYx4s05bNmTxWvXJHN2koaeldijBC8ilcqatAOMeHMuGYfyeOem/vTV6HQSo5TgRaTSWLRxD9ePnUdCfBzv33ISXZrVCXdIIiGjBC8ilcL0VTu59Z0FNKpdjf/c0J9WDTQ6ncQ2JXgRiXmfLdnK795fTIfGtRl3Q18a104Md0giIReykRzMbIyZpZnZsmLWDzKzfWa22Pd4rNC6IWa20szWmNmDoYpRRGJbQYHj5e9Wc/fERZzYqh7v3TJAyV0qjVDW4McCLwNvl7DNDOfcBYUXmFk88E/gHGAzMM/MPnXOpYQqUBGJPTsPHOJ37y9mxup0LurZnL8M7aGhZ6VSCVmCd85NN7M25di1H7DGObcWwMwmAhcDSvAiEpBZa9K5573F7M/K5ZnLunNV35YanU4qHXPOhe7gXoKf4pzr5mfdIGASXi19K3C/c265mQ0FhjjnbvJtNwLo75y7s5hzjAJGASQ2SEzu+nTXoMWfmp5Kl4Zdgna8cIul8sRSWUDlCRbnjIN7Tidr7+nEJ+yidpP3qVI17ZiOqc8mssVSecpTlgW3LFjgnOvjd6VzLmQPoA2wrJh1dYBavufnAat9z68A3ii03QjgpUDOl5yc7IIp+bXgHi/cYqk8sVQW51SeYNixL8td9dpPrvUDU9y97y1yGdm5QTmuPpvIFkvlKU9ZgPmumJwYtrvonXP7Cz3/wsz+ZWYN8Wr0LQtt2gKvhi8i4teM1Tu5973FZB7K57mhPbiiT8vSdxKJcWFL8GbWFNjhnHNm1g/vjv5dwF6go5m1BbYAVwHDwxWniESuvPwCXvh2Nf+ctoaOjWsx4ebedGxSO9xhiUSEkCV4M5sADAIamtlmYDSQAOCcexUYCtxmZnlAFnCV73JDnpndCXwFxANjnHPLQxWniESnbfuyuGfCYuau381v+7TgiYu6Ub2q7pIXOSyUd9EPK2X9y3jd6Pyt+wL4IhRxiUj0+35lGr97bzGH8gr4+5U9ufREzQQnUpRGshORqJGbX8Bfv17Jaz+spXPT2vzz6t60b1Qr3GGJRCQleBGJClv2ZnH3hEUs2LCH4f1b8dgFSRq4RqQESvAiEvG+TdnBfR8sIb/A8dKwE7mwZ/NwhyQS8ZTgRSRi5eQV8JepK3hj5jq6Nq/DP4f3pk3DmuEOSyQqKMGLSETatPsgd05YxJJNe7nupNY8dF4XXZIXKQMleBGJOFOXbef/PlwCwCtX9+bc7s3CHJFI9FGCF5GIcSgvn2e+WMHYWevp2eI4XhrWm1YNaoQ7LJGopAQvIhFhw65M7nx3EUu37OOGU9ry4LmdqVolLtxhiUQtJXgRCbvPf97Gg5N+xgxeH5HM4K5Nwx2SSNRTgheRsMnOzedPn6fwzuyNnNiqLi8NO5EW9XRJXiQYlOBFJCzW7szgjncXkbptP7cMbMf9vzmBhHhdkhcJFiV4EalwnyzewsMfLaVqlTjGjOzDmZ2bhDskkZijBC8iFSYrJ58nPlvOxHmb6NumHv8YdiLNjqse7rBEYpISvIhUiDVpB7hj/CJW7jjA7YPa87tzOlFFl+RFQkYJXkRC7sMFm3l08jJqVI1n3A39OL1To3CHJBLzlOBFJGQO5uTx6OTlTFq4mQHt6vPiVSfSpE5iuMMSqRSU4EUkJFZuP8Ad7y7kl50Z3H1WR+45qyPxcRbusEQqDSV4EQkq5xzZ+3tz0cszqZ2YwDs39ueUDg3DHZZIpaMELyJBk3Eoj0c+XkpG+sWc0qEef7+yF41r65K8SDgowYtIUKRs3c+d7y5k/a5MatT7L2/f8LwuyYuEkfqoiMgxcc7xzuwNXPKvH8nMyePdmwdQo950JXeRMFMNXkTKLW1/No9/tpwvlm5nYKdG/P23PWlQq1q4wxIRlOBFpBxy8gp468d1/OO/q8nNd/x+yAncOrA9caq1i0QMJXgRKZMfVu3kic+Ws3ZnJmd1bsyjFyTRpmHNcIclIkUowYtIQDbtPsiTU1L4JmUHbRrU4K2RfTmjc+NwhyUixVCCF5ESZeXk88q0Nbw6fS1V4ozfDzmBG09tS7Uq8eEOTURKoAQvIn4555i6bDt/+jyVLXuzuLBncx4+r7NmfxOJEkrwInKU1TsO8Phny/lxzS46N63NxFEDGNCuQbjDEpEyUIIXkSP2Z+fy4rerGTdrPTWqxvPERV25un8rTesqEoWU4EWEggLHR4u28OyXK9iVeYir+rbk/sEnqE+7SBRTghep5JZu3sdjny5j0ca99GpZlzEj+9CjRd1whyUix0gJXqSS2p2Zw3NfrWDivE00qFmV54b24PLeLTRYjUiMUIIXqWTy8gsYP2cjz3+9ksycfG44pS33nN2ROokJ4Q5NRIJICV6kEpmzdhejP13Oiu0HOKVDAx6/sCsdm9QOd1giEgJK8CKVwPZ92Tz9RSqfLtnK8XWr88rVvRnSrSlmuhwvEquU4EVi2KG8fN6cuY6Xv1tDXoHj7rM6ctvp7aleVaPQicQ6JXiRGPX9ijSenJLCuvRMzklqwqPnJ9GqQY1whyUiFUQJXiTGbNiVyZOfpfDfFWm0a1iTcTf04/ROjcIdlohUMCV4kRhxMCePf36/hn9PX0dCvPHQuZ25/pS2VK2iUehEKiMleJEo55zj86XbeOrzVLbty+bSE4/nwXM706ROYrhDE5EwUoIXiWIrtx/g8U+X89PaXSQ1q8M/hp1I3zb1wx2WiEQAJXiRKLQvK5cXvl3F2z9toHZiFf54STeG92tFvEahExEfJXiRKFJQ4PhwwWb+PHUFuw/mMLxfK+4ffAL1alYNd2giEmGU4EWixOJNexn96XKWbNpLcut6jLuoH92OPy7cYYlIhFKCF4lw6RmH+MvUFbw/fzONalfjb7/tyaUnHq9R6ESkRErwIhFqT2YO4+ds4LXpa8nKyWfUwHbcdWYHamtSGBEJgBK8SIRZl57JmzPX8uGCzWTnFnBm58Y8fF4XOjSuFe7QRCSKKMGLRADnvJne/j1jHf9dsYOEuDguObE5N57ajhOaarY3ESm7kCV4MxsDXACkOee6lbBdX2A2cKVz7kPfsvXAASAfyHPO9QlVnCLhlJtfwBdLt7Fv6yiufH029WokcNcZHRhxUhsa1a4W7vBEJIqFsgY/FngZeLu4DcwsHvgz8JWf1Wc459JDE5pIeO3PzmXi3I2M/XE9W/dlE59Qjacu7cZlJ7bQTG8iEhQhS/DOuelm1qaUze4CJgF9QxWHSCTZtPsgb/24nvfmbSQzJ58B7erzx0u68fvpF3B1/+vCHZ6IxBBzzoXu4F6Cn+LvEr2ZHQ+8C5wJvOnb7vAl+nXAHsABrznnXi/hHKOAUQCJDRKTuz7dNWjxp6an0qVhl6AdL9xiqTzRVpbc7OPJ2ncyOZlJgKNarWVUP+4nqlTbBkRfeUoTS+WJpbKAyhPJylOWBbcsWFBsM7ZzLmQPoA2wrJh1HwADfM/HAkMLrWvu+9kYWAIMDOR8ycnJLpiSXwvu8cItlsoTDWXJyy9wXy7d6i7/14+u9QNTXLfRU93TX6S4rXsPHrVtNJSnLGKpPLFUFudUnkhWnrIA810xOTGcd9H3ASb6ButoCJxnZnnOucnOua0Azrk0M/sY6AdMD1+oIoHLPJTHB/M3MebH9WzcfZCW9asz+sIkrujTklrV1HFFRCpG2L5tnHNtDz83s7F4l+gnm1lNIM45d8D3fDDwZJjCFAnY9n3ZjPtpPeNnb2B/dh69W9XloXM7M7hrU00CIyIVLpTd5CYAg4CGZrYZGA0kADjnXi1h1ybAx76afRXgXefc1FDFKXKslm/dx5sz1vHpkq0UOMeQbk258dR2JLeuF+7QRKQSC+Vd9MPKsO3IQs/XAj1DEZNIsBQUOKatSuONGeuY9csualSNZ8RJrbn+5La0alAj3OGJiGgkO5GyyM7N5+NFW3hjxlp+2ZlJ0zqJPHRuZ67q14rjqmuMeBGJHErwIgFIzzjEf37awDuzN7ArM4euzevwwpW9OL9HMxLi48IdnojIUZTgRUqwescB3py5jo8WbSEnr4CzuzTmxlPbMaBdfU3XKiIRTQlepAjnHLN+2cW/Z6xl2sqdVKsSx9DkFtx4alvaN9KMbiISHZTgRXxy8gr4bMlW3pi5jtRt+2lYqyq/O6cT1wxoTf2aVcMdnohImSjBS6W392AO4+dsZNys9aQdOESnJrX4y+U9uKhXcxITNPGLiEQnJXiplPILHPPX7+azn7cyacEWsnLzOa1jQ567oicDOzZU+7qIRD0leKk0snPzmbk6na+Wb+e/K9LYnZlD1SpxXNSzOTed1pbOTeuEO0QRkaBRgpeYtu9gLt+t3MFXy3bww6qdZOXmUzuxCmd2bszgpKacfkIjjQ8vIjFJ32wSc7bty+KblB18tXw7c9buJq/A0bh2NS5PPp7BSU0Z0K4BVauo77qIxDYleIl6zjnWpGXwtS+p/7x5HwDtG9Xk5oHtGJzUhJ4t6hKnCV9EpBJRgpeo5JyxYMMevk7ZztfLd7AuPROAXi3r8vshJzA4qSkdGqvPuohUXkrwEjUO5eXz0y+7+DplB7s33sflr8yiSpxxUvsG3HBqWwYnNaFJncRwhykiEhGU4CWiHcjOZdrKnXy1fDvTVu4k41AeNavGk5C4gecuvIZBJzTWJC8iIn4owUvESTuQzbcpaXy1fDuzfkknN9/RoGZVLujRjMFdm3By+4ac+tYjXNzrgXCHKiISsZTgJSKsS8/kq+Xb+Xr5dhZt2otz0Kp+DUae3IbBXZvSu1U94nWTnIhIwJTgJSyccyzdss+X1HewOi0DgG7H1+HeszsxuGsTTmhSWyPKiYiUkxK8VJjc/ALmrtvNV8u3803KDrbtyyY+zujXpj7D+7finKQmtKhXI9xhiojEBCV4CamDOXn8sHInX6fs4L+pO9ifnUdiQhwDOzbivsEncFbnxtTTTG0iIkGnBC9Bk5NXwC87M0jZup+UbftJ2bqfhRv3cCivgLo1EjgnqSmDuzZhYMdGVK+qWdpEREJJCV7KZd/BXC+J+xJ56rb9rE47QG6+A6BalTg6N6vD8P6tGJzUlL5t6lElXsPDiohUFCV4KVFBgWPTnoNHknjKtv2kbjvAlr1ZR7ZpVLsaSc3qMLBTI5Ka1yGpWW3aNKiphC4iEkZK8HJEdm4+K7cf8CVxr2a+YvsBMg7lARBn0L5RLfq0qceIZq1JalaHLs3q0Kh2tTBHLiIiRSnBV1JpB7JJ3XbgVzXztTszKPCusFOrWhW6NKvNZb2PJ6lZHZKa16FTk9okJqjtXEQkGijBx7i8/ALWpWeSsm0/mbvO4doxc0nZup/0jENHtjm+bnW6NKvDed2a+i6xH0eLetU1+5qISBRTgo8hB7JzWbH9wJHL66nbvEvsh/IKfFsMIL3GIQad0Iguzer4LrHXpm4NdVMTEYk1SvBRwjlHdm4B+7NzOZCdy76sPNIzDnlt5lv3k7p9Pxt2HTyyfd0aCSQ1q8OIAa29ZN68Dtd8eg5f3DInjKUQEZGKogRfQZxzZObksz8r15ek844835/lPT9wqMiyItsd7oJWVNuGNenavA5XJLc4ksyb1kk8aphXs/yKKKqIiEQAJfgAFRS4UhPw/5b973nhbQr85+cjEhPiqJOYQJ3qCdRJrEK9GlVp3aAmdRKrUKd6ArUTq/xqfd0aVenQuBa1quljFBGRX1NmKMa3f7uel/cVsPCJk8kvcOS7X2dnA47zPQ6LNyM+zntUOfwzPo74GkaVWr9e7j2P+9WyOH8Tq2T7HnuPvUyvpW+Gt84/9gNFgFgqC6g8kSyWygIqT9g17Q7nPlshp1KCL0ZCvGGWR53qCf9L1oUTc7wdtdzQXeciIhIhnHMx80hOTnbBlPxacI8XbrFUnlgqi3MqTySLpbI4p/JEsvKUBZjvismJGktUREQkBinBi4iIxCAleBERkRikBC8iIhKDlOBFRERikBK8iIhIDFKCFxERiUFK8CIiIjFICV5ERCQGKcGLiIjEICV4ERGRGKQELyIiEoOU4EVERGKQuSLznEczM9sJbAjiIRsC6UE8XrjFUnliqSyg8kSyWCoLqDyRrDxlae2ca+RvRUwl+GAzs/nOuT7hjiNYYqk8sVQWUHkiWSyVBVSeSBbssugSvYiISAxSghcREYlBSvAlez3cAQRZLJUnlsoCKk8ki6WygMoTyYJaFrXBi4iIxCDV4EVERGJQpU/wZjbEzFaa2Roze9DPejOzf/jW/2xmvcMRZ6ACKM/VvnL8bGazzKxnOOIMVGnlKbRdXzPLN7OhFRlfWQRSFjMbZGaLzWy5mf1Q0TGWRQC/a8eZ2WdmtsRXnuvDEWcgzGyMmaWZ2bJi1kfb90Bp5Ym274ESy1Nou2j4Hii1LEH7HnDOVdoHEA/8ArQDqgJLgKQi25wHfAkYMACYE+64j7E8JwP1fM/PjfbyFNruO+ALYGi44z6Gz6YukAK08r1uHO64j7E8DwN/9j1vBOwGqoY79mLKMxDoDSwrZn3UfA8EWJ6o+R4IpDy+bSL+eyDAzyZo3wOVvQbfD1jjnFvrnMsBJgIXF9nmYuBt55kN1DWzZhUdaIBKLY9zbpZzbo/v5WygRQXHWBaBfD4AdwGTgLSKDK6MAinLcOAj59xGAOdctJfHAbXNzIBaeAk+r2LDDIxzbjpefMWJpu+BUssTZd8DgXw+EB3fA4GUJWjfA5U9wR8PbCr0erNvWVm3iRRljfVGvFpJpCq1PGZ2PHAp8GoFxlUegXw2nYB6ZjbNzBaY2bUVFl3ZBVKel4EuwFZgKXCPc66gYsILumj6HiirSP8eKFUUfQ8EImjfA1WCGFQ0Mj/LinYrCGSbSBFwrGZ2Bt4f9qkhjejYBFKeF4AHnHP5XkUxYgVSlipAMnAWUB34ycxmO+dWhTq4cgikPL8BFgNnAu2Bb8xshnNuf4hjC4Vo+h4IWJR8DwTiBaLjeyAQQfseqOwJfjPQstDrFni1jbJuEykCitXMegBvAOc653ZVUGzlEUh5+gATfX/UDYHzzCzPOTe5QiIMXKC/a+nOuUwg08ymAz2BSEzwgZTneuBZ5zUkrjGzdUBnYG7FhBhU0fQ9EJAo+h4IRLR8DwQiaN8Dlf0S/Tygo5m1NbOqwFXAp0W2+RS41ncX7QBgn3NuW0UHGqBSy2NmrYCPgBERWjMsrNTyOOfaOufaOOfaAB8Ct0foH3Ugv2ufAKeZWRUzqwH0B1IrOM5ABVKejXi1EMysCXACsLZCowyeaPoeKFWUfQ+UKoq+BwIRtO+BSl2Dd87lmdmdwFd4d2COcc4tN7Nbfetfxbsj8zxgDXAQr1YSkQIsz2NAA+Bfvv9281yETtQQYHmiQiBlcc6lmtlU4GegAHjDOVdit6BwCfCz+SMw1syW4l3ifsA5F5GzfpnZBGAQ0NDMNgOjgQSIvu8BCKg8UfM9AAGVJ2qUVpZgfg9oJDsREZEYVNkv0YuIiMQkJXgREZEYpAQvIiISg5TgRUREYpASvIiISAxSgheJImZ2qZk5M+scxGMOMrMpvucXmW9mODO7xMySynG8aWZWpi5XZvaCmQ0s67mKHCPD97ORr5uRSKWmBC8SXYYBM/EGlgk659ynzrlnfS8vAcqc4MvKzOoDA3yTcBRdF1/W4znndgLbzOyUYMQnEq2U4EWihJnVAk7BGzv8qkLLB5nZD2b2vpmtMrNnzZvve66ZLTWz9r7txprZq2Y2w7fdBX7OMdLMXjazk4GLgOd881K3L1wzN7OGZrbe97y6mU00b27x9/DGzz58vMFm9pOZLTSzD3xlKGooMLXQPuvN7DEzmwlcYWY3m9k88+aVn+Qb3QvfKHo/+db9scgxJwNXl/lNFokhSvAi0eMSYKpvaNHdZta70LqewD1Ad2AE0Mk51w9vrPG7Cm3XBjgdOB941cwS/Z3IOTcLb3jW/3PO9XLO/VJCXLcBB51zPYCn8CbKwMwaAo8AZzvnegPzgd/52f8UYEGRZdnOuVOdcxPxps7s65zriTdk542+bV4EXnHO9QW2F9l/PnBaCTGLxDwleJHoMQxv3nV8P4cVWjfPObfNOXcI+AX42rd8KV5SP+x951yBc2413rjwwWjLHwi8A+Cc+xlviE2AAXiX+H80s8XAdUBrP/s3A3YWWfZeoefdfFcdluLVyrv6lp8CTPA9/0+R/dOA5mUuiUgMqdRj0YtECzNrgDftajczc3jjvzsz+71vk0OFNi8o9LqAX/+dFx2buixjVefxv0pB0Zq/v+MY8I1zbpifdYVl+TleZqHnY4FLnHNLzGwk3jjeJZ33cHxZpZxXJKapBi8SHYYCbzvnWvtmzWoJrKPs83hfYWZxvnb5dsDKErY9ANQu9Ho9vsvvvngOm46vvdvMugE9fMtnA6eYWQffuhpm1snPeVKBDiXEURvvprkEft2u/iP/uxehaHt7JyAiJ+oRqShK8CLRYRjwcZFlk4DhZTzOSuAH4EvgVudcdgnbTgT+z8wW+f4h+Ctwm5nNwptz+7BXgFpm9jPwe3zzvfvuZh8JTPCtm43/JoHP+XWtvKhHgTnAN8CKQsvvAe4ws3nAcUX2OcN3XJFKS7PJiVQSZjYWmOKc+zDcsRTlu2P+Aufc3iAdbzpwsXNuTzCOJxKNVIMXkUhwH9AqGAcys0bA35TcpbJTDV5ERCQGqQYvIiISg5TgRUREYpASvIiISAxSghcREYlBSvAiIiIxSAleREQkBv1/fOgjmn8CAeYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "#\n", "def periodLima(L, thetaZero):\n", " if not hasattr(periodLima, \"period.g\"):\n", " periodLima.g = 9.81 # m/s^2\n", " periodLima.small = 1E-9\n", " #\n", " T = -2*np.pi*np.sqrt(L/periodLima.g)*np.log(np.cos(thetaZero/2))/(1 - np.cos(thetaZero/2))\n", " #\n", " return T\n", "#\n", "L = 0.5 # m\n", "thetaBot = 0.0001\n", "thetaTop = np.pi/2\n", "nTheta = 10\n", "thetaArr = np.linspace(thetaBot, thetaTop, nTheta + 1)\n", "periodArr = periodLima(L, thetaArr)\n", "#\n", "plt.figure(figsize = (8, 6))\n", "plt.title('Period as function of amplitude', fontsize = 14)\n", "plt.xlabel('Amplitude (rad)')\n", "plt.ylabel('Period (sec)')\n", "plt.plot(thetaArr, periodArr, label = \"Lima formula\")\n", "plt.plot(thetaArr, 2*np.pi*np.sqrt(L/periodLima.g)*np.ones(nTheta + 1), \n", " label = \"Simple pendulum\")\n", "plt.grid(color = 'g')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Create Monte Carlo data of a series of measurements of a pendulum of a fixed length but with a range of amplitudes, including large values where the small angle approximation is not valid.\n", "Create both small and large datasets, so the influence of statistical precision in making inferences can be investigated." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "Debug = False\n", "thetaBot = 0.1\n", "thetaTop = np.pi/2\n", "nTheta = 16\n", "L = 0.5\n", "thetaArr = np.linspace(thetaBot, thetaTop, nTheta) # rad\n", "sigma = 0.5\n", "exact = np.zeros(nTheta)\n", "#\n", "nMeas = 10\n", "periodTabSm = np.zeros((nTheta, nMeas))\n", "#\n", "for n in range(0, nTheta):\n", " exact[n] = periodLima(L, thetaArr[n])\n", " periodTabSm[n, 0:nMeas] = np.random.normal(exact[n], scale = sigma, size = nMeas) # seconds\n", "if Debug: print(periodTab)\n", "#\n", "np.savetxt(\"thetaTable.csv\", thetaArr, delimiter = ',', )\n", "np.savetxt(\"periodTableSmall.csv\", periodTabSm, delimiter = ',')\n", "#\n", "nMeas = 1000\n", "periodTabLg = np.zeros((nTheta, nMeas))\n", "#\n", "for n in range(0, nTheta):\n", " periodTabLg[n, 0:nMeas] = np.random.normal(exact[n], sigma, nMeas) # seconds\n", "if Debug: print(periodTab)\n", "#\n", "np.savetxt(\"periodTableLarge.csv\", periodTabLg, delimiter = ',')" ] } ], "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": 4 }