obsidian_backup/补课/多体动力学/18_Simulation_and_Visualization.ipynb

1375 lines
124 KiB
Plaintext

{
"cells": [
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"import sympy as sm\n",
"import sympy.physics.mechanics as me\n",
"me.init_vprinting(use_latex='mathjax')"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left( \\frac{l u_{1}}{2}\\mathbf{\\hat{a}_y}, \\ l u_{1}\\mathbf{\\hat{a}_y}, \\ u_{1}\\mathbf{\\hat{n}_z}, \\ u_{2}\\mathbf{\\hat{a}_x} + u_{1}\\mathbf{\\hat{n}_z}\\right)$"
],
"text/plain": [
"⎛l⋅u₁ ⎞\n",
"⎜──── a_y, l⋅u₁ a_y, u₁ n_z, u₂ a_x + u₁ n_z⎟\n",
"⎝ 2 ⎠"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m, g, kt, kl, l = sm.symbols('m, g, k_t, k_l, l')\n",
"q1, q2, q3 = me.dynamicsymbols('q1, q2, q3')\n",
"u1, u2, u3 = me.dynamicsymbols('u1, u2, u3')\n",
"\n",
"N = me.ReferenceFrame('N')\n",
"A = me.ReferenceFrame('A')\n",
"B = me.ReferenceFrame('B')\n",
"\n",
"A.orient_axis(N, q1, N.z)\n",
"B.orient_axis(A, q2, A.x)\n",
"\n",
"A.set_ang_vel(N, u1*N.z)\n",
"B.set_ang_vel(A, u2*A.x)\n",
"\n",
"O = me.Point('O')\n",
"Ao = me.Point('A_O')\n",
"Bo = me.Point('B_O')\n",
"\n",
"Ao.set_pos(O, l/2*A.x)\n",
"Bo.set_pos(O, l*A.x)\n",
"\n",
"O.set_vel(N, 0)\n",
"Ao.v2pt_theory(O, N, A)\n",
"Bo.v2pt_theory(O, N, A)\n",
"\n",
"Ao.vel(N), Bo.vel(N), A.ang_vel_in(N), B.ang_vel_in(N)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle - q_{3} u_{1} \\cos{\\left(q_{2} \\right)}\\mathbf{\\hat{b}_x} + u_{3}\\mathbf{\\hat{b}_y} + q_{3} u_{2}\\mathbf{\\hat{b}_z} + l u_{1}\\mathbf{\\hat{a}_y}$"
],
"text/plain": [
"-q₃⋅u₁⋅cos(q₂) b_x + u₃ b_y + q₃⋅u₂ b_z + l⋅u₁ a_y"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Q = me.Point('Q')\n",
"Q.set_pos(Bo, q3*B.y)\n",
"Q.set_vel(B, u3*B.y)\n",
"Q.v1pt_theory(Bo, N, B)\n",
"\n",
"Q.vel(N)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left( - \\frac{l u_{1}^{2}}{2}\\mathbf{\\hat{a}_x} + \\frac{l \\dot{u}_{1}}{2}\\mathbf{\\hat{a}_y}, \\ - l u_{1}^{2}\\mathbf{\\hat{a}_x} + l \\dot{u}_{1}\\mathbf{\\hat{a}_y}, \\ \\dot{u}_{1}\\mathbf{\\hat{n}_z}, \\ \\dot{u}_{2}\\mathbf{\\hat{a}_x} + u_{1} u_{2}\\mathbf{\\hat{a}_y} + \\dot{u}_{1}\\mathbf{\\hat{n}_z}\\right)$"
],
"text/plain": [
"⎛ 2 ↪\n",
"⎜-l⋅u₁ l⋅u₁̇ 2 ↪\n",
"⎜─────── a_x + ──── a_y, -l⋅u₁ a_x + l⋅u₁̇ a_y, u₁̇ n_z, u₂̇ a_x + u₁⋅u₂ a_y + u ↪\n",
"⎝ 2 2 ↪\n",
"\n",
"↪ ⎞\n",
"↪ ⎟\n",
"↪ ₁̇ n_z⎟\n",
"↪ ⎠"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Ao.acc(N), Bo.acc(N), A.ang_acc_in(N), B.ang_acc_in(N)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle (q_{3} u_{1} u_{2} \\sin{\\left(q_{2} \\right)} + q_{3} u_{1} \\sin{\\left(q_{2} \\right)} \\dot{q}_{2} - q_{3} \\cos{\\left(q_{2} \\right)} \\dot{u}_{1} - u_{1} u_{3} \\cos{\\left(q_{2} \\right)} - u_{1} \\cos{\\left(q_{2} \\right)} \\dot{q}_{3})\\mathbf{\\hat{b}_x} + (- q_{3} u_{1}^{2} \\cos^{2}{\\left(q_{2} \\right)} - q_{3} u_{2}^{2} + \\dot{u}_{3})\\mathbf{\\hat{b}_y} + (q_{3} u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)} + q_{3} \\dot{u}_{2} + u_{2} u_{3} + u_{2} \\dot{q}_{3})\\mathbf{\\hat{b}_z} - l u_{1}^{2}\\mathbf{\\hat{a}_x} + l \\dot{u}_{1}\\mathbf{\\hat{a}_y}$"
],
"text/plain": [
" ↪\n",
"(q₃⋅u₁⋅u₂⋅sin(q₂) + q₃⋅u₁⋅sin(q₂)⋅q₂̇ - q₃⋅cos(q₂)⋅u₁̇ - u₁⋅u₃⋅cos(q₂) - u₁⋅cos( ↪\n",
"\n",
"↪ ⎛ 2 2 2 ⎞ ⎛ 2 ↪\n",
"↪ q₂)⋅q₃̇) b_x + ⎝- q₃⋅u₁ ⋅cos (q₂) - q₃⋅u₂ + u₃̇⎠ b_y + ⎝q₃⋅u₁ ⋅sin(q₂)⋅cos(q₂ ↪\n",
"\n",
"↪ ⎞ 2 \n",
"↪ ) + q₃⋅u₂̇ + u₂⋅u₃ + u₂⋅q₃̇⎠ b_z + -l⋅u₁ a_x + l⋅u₁̇ a_y"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Q.acc(N)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle (2 q_{3} u_{1} u_{2} \\sin{\\left(q_{2} \\right)} - q_{3} \\cos{\\left(q_{2} \\right)} \\dot{u}_{1} - 2 u_{1} u_{3} \\cos{\\left(q_{2} \\right)})\\mathbf{\\hat{b}_x} + (- q_{3} u_{1}^{2} \\cos^{2}{\\left(q_{2} \\right)} - q_{3} u_{2}^{2} + \\dot{u}_{3})\\mathbf{\\hat{b}_y} + (q_{3} u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)} + q_{3} \\dot{u}_{2} + 2 u_{2} u_{3})\\mathbf{\\hat{b}_z} - l u_{1}^{2}\\mathbf{\\hat{a}_x} + l \\dot{u}_{1}\\mathbf{\\hat{a}_y}$"
],
"text/plain": [
" ⎛ 2 2 ↪\n",
"(2⋅q₃⋅u₁⋅u₂⋅sin(q₂) - q₃⋅cos(q₂)⋅u₁̇ - 2⋅u₁⋅u₃⋅cos(q₂)) b_x + ⎝- q₃⋅u₁ ⋅cos (q₂ ↪\n",
"\n",
"↪ 2 ⎞ ⎛ 2 ⎞ ↪\n",
"↪ ) - q₃⋅u₂ + u₃̇⎠ b_y + ⎝q₃⋅u₁ ⋅sin(q₂)⋅cos(q₂) + q₃⋅u₂̇ + 2⋅u₂⋅u₃⎠ b_z + -l⋅u ↪\n",
"\n",
"↪ 2 \n",
"↪ ₁ a_x + l⋅u₁̇ a_y"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = me.dynamicsymbols._t\n",
"\n",
"qdot_repl = {q1.diff(t): u1,\n",
" q2.diff(t): u2,\n",
" q3.diff(t): u3}\n",
"\n",
"Q.set_acc(N, Q.acc(N).xreplace(qdot_repl))\n",
"Q.acc(N)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"R_Ao = m*g*N.x\n",
"R_Bo = m*g*N.x + kl*q3*B.y\n",
"R_Q = m/4*g*N.x - kl*q3*B.y\n",
"T_A = -kt*q1*N.z + kt*q2*A.x\n",
"T_B = -kt*q2*A.x"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"I = m*l**2/12\n",
"I_A_Ao = I*me.outer(A.y, A.y) + I*me.outer(A.z, A.z)\n",
"I_B_Bo = I*me.outer(B.x, B.x) + I*me.outer(B.z, B.z)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"v_Ao_1 = Ao.vel(N).diff(u1, N)\n",
"v_Bo_1 = Bo.vel(N).diff(u1, N)\n",
"v_Q_1 = Q.vel(N).diff(u1, N)\n",
"\n",
"v_Ao_2 = Ao.vel(N).diff(u2, N)\n",
"v_Bo_2 = Bo.vel(N).diff(u2, N)\n",
"v_Q_2 = Q.vel(N).diff(u2, N)\n",
"\n",
"v_Ao_3 = Ao.vel(N).diff(u3, N)\n",
"v_Bo_3 = Bo.vel(N).diff(u3, N)\n",
"v_Q_3 = Q.vel(N).diff(u3, N)\n",
"\n",
"w_A_1 = A.ang_vel_in(N).diff(u1, N)\n",
"w_B_1 = B.ang_vel_in(N).diff(u1, N)\n",
"\n",
"w_A_2 = A.ang_vel_in(N).diff(u2, N)\n",
"w_B_2 = B.ang_vel_in(N).diff(u2, N)\n",
"\n",
"w_A_3 = A.ang_vel_in(N).diff(u3, N)\n",
"w_B_3 = B.ang_vel_in(N).diff(u3, N)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"F1 = v_Ao_1.dot(R_Ao) + v_Bo_1.dot(R_Bo) + v_Q_1.dot(R_Q) + w_A_1.dot(T_A) + w_B_1.dot(T_B)\n",
"F2 = v_Ao_2.dot(R_Ao) + v_Bo_2.dot(R_Bo) + v_Q_2.dot(R_Q) + w_A_2.dot(T_A) + w_B_2.dot(T_B)\n",
"F3 = v_Ao_3.dot(R_Ao) + v_Bo_3.dot(R_Bo) + v_Q_3.dot(R_Q) + w_A_3.dot(T_A) + w_B_3.dot(T_B)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}- \\frac{7 g l m \\sin{\\left(q_{1} \\right)}}{4} - \\frac{g m q_{3} \\cos{\\left(q_{1} \\right)} \\cos{\\left(q_{2} \\right)}}{4} - k_{t} q_{1}\\\\\\frac{g m q_{3} \\sin{\\left(q_{1} \\right)} \\sin{\\left(q_{2} \\right)}}{4} - k_{t} q_{2}\\\\- \\frac{g m \\sin{\\left(q_{1} \\right)} \\cos{\\left(q_{2} \\right)}}{4} - k_{l} q_{3}\\end{matrix}\\right]$"
],
"text/plain": [
"⎡ 7⋅g⋅l⋅m⋅sin(q₁) g⋅m⋅q₃⋅cos(q₁)⋅cos(q₂) ⎤\n",
"⎢- ─────────────── - ────────────────────── - kₜ⋅q₁⎥\n",
"⎢ 4 4 ⎥\n",
"⎢ ⎥\n",
"⎢ g⋅m⋅q₃⋅sin(q₁)⋅sin(q₂) ⎥\n",
"⎢ ────────────────────── - kₜ⋅q₂ ⎥\n",
"⎢ 4 ⎥\n",
"⎢ ⎥\n",
"⎢ g⋅m⋅sin(q₁)⋅cos(q₂) ⎥\n",
"⎢ - ─────────────────── - kₗ⋅q₃ ⎥\n",
"⎣ 4 ⎦"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Fr = sm.Matrix([F1, F2, F3])\n",
"Fr"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"TAs = -(A.ang_acc_in(N).dot(I_A_Ao) + me.cross(A.ang_vel_in(N), I_A_Ao).dot(A.ang_vel_in(N)))\n",
"TBs = -(B.ang_acc_in(N).dot(I_B_Bo) + me.cross(B.ang_vel_in(N), I_B_Bo).dot(B.ang_vel_in(N)))\n",
"\n",
"F1s = v_Ao_1.dot(-m*Ao.acc(N)) + v_Bo_1.dot(-m*Bo.acc(N)) + v_Q_1.dot(-m/4*Q.acc(N))\n",
"F1s += w_A_1.dot(TAs) + w_B_1.dot(TBs)\n",
"\n",
"F2s = v_Ao_2.dot(-m*Ao.acc(N)) + v_Bo_2.dot(-m*Bo.acc(N)) + v_Q_2.dot(-m/4*Q.acc(N))\n",
"F2s += w_A_2.dot(TAs) + w_B_2.dot(TBs)\n",
"\n",
"F3s = v_Ao_3.dot(-m*Ao.acc(N)) + v_Bo_3.dot(-m*Bo.acc(N)) + v_Q_3.dot(-m/4*Q.acc(N))\n",
"F3s += w_A_3.dot(TAs) + w_B_3.dot(TBs)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}- \\frac{19 l^{2} m \\dot{u}_{1}}{12} - \\frac{l m q_{3} u_{1}^{2} \\cos{\\left(q_{2} \\right)}}{4} + l \\left(- \\frac{m \\left(- q_{3} u_{1}^{2} \\cos^{2}{\\left(q_{2} \\right)} - q_{3} u_{2}^{2} + \\dot{u}_{3}\\right) \\cos{\\left(q_{2} \\right)}}{4} + \\frac{m \\left(q_{3} u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)} + q_{3} \\dot{u}_{2} + 2 u_{2} u_{3}\\right) \\sin{\\left(q_{2} \\right)}}{4}\\right) + \\frac{m \\left(2 q_{3} u_{1} u_{2} \\sin{\\left(q_{2} \\right)} - q_{3} \\cos{\\left(q_{2} \\right)} \\dot{u}_{1} - 2 u_{1} u_{3} \\cos{\\left(q_{2} \\right)}\\right) q_{3} \\cos{\\left(q_{2} \\right)}}{4} + \\left(- \\frac{l^{2} m \\left(- u_{1} u_{2} \\sin{\\left(q_{2} \\right)} + \\cos{\\left(q_{2} \\right)} \\dot{u}_{1}\\right)}{12} + \\frac{l^{2} m u_{1} u_{2} \\sin{\\left(q_{2} \\right)}}{12}\\right) \\cos{\\left(q_{2} \\right)}\\\\- \\frac{l^{2} m u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)}}{12} - \\frac{l^{2} m \\dot{u}_{2}}{12} + \\frac{l m q_{3} \\sin{\\left(q_{2} \\right)} \\dot{u}_{1}}{4} - \\frac{m \\left(q_{3} u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)} + q_{3} \\dot{u}_{2} + 2 u_{2} u_{3}\\right) q_{3}}{4}\\\\- \\frac{l m \\cos{\\left(q_{2} \\right)} \\dot{u}_{1}}{4} - \\frac{m \\left(- q_{3} u_{1}^{2} \\cos^{2}{\\left(q_{2} \\right)} - q_{3} u_{2}^{2} + \\dot{u}_{3}\\right)}{4}\\end{matrix}\\right]$"
],
"text/plain": [
"⎡ 2 2 ⎛ ⎛ 2 2 2 ⎞ ↪\n",
"⎢ 19⋅l ⋅m⋅u₁̇ l⋅m⋅q₃⋅u₁ ⋅cos(q₂) ⎜ m⋅⎝- q₃⋅u₁ ⋅cos (q₂) - q₃⋅u₂ + u₃̇⎠⋅ ↪\n",
"⎢- ────────── - ────────────────── + l⋅⎜- ──────────────────────────────────── ↪\n",
"⎢ 12 4 ⎝ 4 ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎣ ↪\n",
"\n",
"↪ ⎛ 2 ⎞ ⎞ ↪\n",
"↪ cos(q₂) m⋅⎝q₃⋅u₁ ⋅sin(q₂)⋅cos(q₂) + q₃⋅u₂̇ + 2⋅u₂⋅u₃⎠⋅sin(q₂)⎟ m⋅(2⋅q₃⋅u₁ ↪\n",
"↪ ─────── + ────────────────────────────────────────────────────⎟ + ────────── ↪\n",
"↪ 4 ⎠ ↪\n",
"↪ ↪\n",
"↪ 2 2 2 ⎛ ↪\n",
"↪ l ⋅m⋅u₁ ⋅sin(q₂)⋅cos(q₂) l ⋅m⋅u₂̇ l⋅m⋅q₃⋅sin(q₂)⋅u₁̇ m⋅⎝q₃⋅u ↪\n",
"↪ - ──────────────────────── - ─────── + ───────────────── - ─────── ↪\n",
"↪ 12 12 4 ↪\n",
"↪ ↪\n",
"↪ ⎛ 2 2 ↪\n",
"↪ l⋅m⋅cos(q₂)⋅u₁̇ m⋅⎝- q₃⋅u₁ ⋅cos (q₂) ↪\n",
"↪ - ────────────── - ───────────────────── ↪\n",
"↪ 4 4 ↪\n",
"\n",
"↪ ⎛ 2 ↪\n",
"↪ ⋅u₂⋅sin(q₂) - q₃⋅cos(q₂)⋅u₁̇ - 2⋅u₁⋅u₃⋅cos(q₂))⋅q₃⋅cos(q₂) ⎜ l ⋅m⋅(-u₁⋅u₂⋅ ↪\n",
"↪ ───────────────────────────────────────────────────────── + ⎜- ───────────── ↪\n",
"↪ 4 ⎝ ↪\n",
"↪ ↪\n",
"↪ 2 ⎞ ↪\n",
"↪ ₁ ⋅sin(q₂)⋅cos(q₂) + q₃⋅u₂̇ + 2⋅u₂⋅u₃⎠⋅q₃ ↪\n",
"↪ ──────────────────────────────────────── ↪\n",
"↪ 4 ↪\n",
"↪ ↪\n",
"↪ 2 ⎞ ↪\n",
"↪ - q₃⋅u₂ + u₃̇⎠ ↪\n",
"↪ ────────────── ↪\n",
"↪ ↪\n",
"\n",
"↪ 2 ⎞ ⎤\n",
"↪ sin(q₂) + cos(q₂)⋅u₁̇) l ⋅m⋅u₁⋅u₂⋅sin(q₂)⎟ ⎥\n",
"↪ ───────────────────── + ──────────────────⎟⋅cos(q₂)⎥\n",
"↪ 12 12 ⎠ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎦"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Frs = sm.Matrix([F1s, F2s, F3s])\n",
"Frs"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}\\dot{u}_{1}\\\\\\dot{u}_{2}\\\\\\dot{u}_{3}\\end{matrix}\\right]$"
],
"text/plain": [
"⎡u₁̇⎤\n",
"⎢ ⎥\n",
"⎢u₂̇⎥\n",
"⎢ ⎥\n",
"⎣u₃̇⎦"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"u = sm.Matrix([u1, u2, u3])\n",
"ud = u.diff(t)\n",
"ud"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}- \\frac{l^{2} m \\cos^{2}{\\left(q_{2} \\right)}}{12} - \\frac{19 l^{2} m}{12} - \\frac{m q_{3}^{2} \\cos^{2}{\\left(q_{2} \\right)}}{4} & \\frac{l m q_{3} \\sin{\\left(q_{2} \\right)}}{4} & - \\frac{l m \\cos{\\left(q_{2} \\right)}}{4}\\\\\\frac{l m q_{3} \\sin{\\left(q_{2} \\right)}}{4} & - \\frac{l^{2} m}{12} - \\frac{m q_{3}^{2}}{4} & 0\\\\- \\frac{l m \\cos{\\left(q_{2} \\right)}}{4} & 0 & - \\frac{m}{4}\\end{matrix}\\right]$"
],
"text/plain": [
"⎡ 2 2 2 2 2 ⎤\n",
"⎢ l ⋅m⋅cos (q₂) 19⋅l ⋅m m⋅q₃ ⋅cos (q₂) l⋅m⋅q₃⋅sin(q₂) -l⋅m⋅cos(q₂) ⎥\n",
"⎢- ───────────── - ─────── - ────────────── ────────────── ─────────────⎥\n",
"⎢ 12 12 4 4 4 ⎥\n",
"⎢ ⎥\n",
"⎢ 2 2 ⎥\n",
"⎢ l⋅m⋅q₃⋅sin(q₂) l ⋅m m⋅q₃ ⎥\n",
"⎢ ────────────── - ──── - ───── 0 ⎥\n",
"⎢ 4 12 4 ⎥\n",
"⎢ ⎥\n",
"⎢ -l⋅m⋅cos(q₂) -m ⎥\n",
"⎢ ───────────── 0 ─── ⎥\n",
"⎣ 4 4 ⎦"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Md = Frs.jacobian(ud)\n",
"Md"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left\\{ \\dot{u}_{1} : 0, \\ \\dot{u}_{2} : 0, \\ \\dot{u}_{3} : 0\\right\\}$"
],
"text/plain": [
"{u₁̇: 0, u₂̇: 0, u₃̇: 0}"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ud_zerod = {udr: 0 for udr in ud}\n",
"ud_zerod\n"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}- \\frac{7 g l m \\sin{\\left(q_{1} \\right)}}{4} - \\frac{g m q_{3} \\cos{\\left(q_{1} \\right)} \\cos{\\left(q_{2} \\right)}}{4} - k_{t} q_{1} + \\frac{l^{2} m u_{1} u_{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)}}{6} - \\frac{l m q_{3} u_{1}^{2} \\cos{\\left(q_{2} \\right)}}{4} + l \\left(- \\frac{m \\left(- q_{3} u_{1}^{2} \\cos^{2}{\\left(q_{2} \\right)} - q_{3} u_{2}^{2}\\right) \\cos{\\left(q_{2} \\right)}}{4} + \\frac{m \\left(q_{3} u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)} + 2 u_{2} u_{3}\\right) \\sin{\\left(q_{2} \\right)}}{4}\\right) + \\frac{m \\left(2 q_{3} u_{1} u_{2} \\sin{\\left(q_{2} \\right)} - 2 u_{1} u_{3} \\cos{\\left(q_{2} \\right)}\\right) q_{3} \\cos{\\left(q_{2} \\right)}}{4}\\\\\\frac{g m q_{3} \\sin{\\left(q_{1} \\right)} \\sin{\\left(q_{2} \\right)}}{4} - k_{t} q_{2} - \\frac{l^{2} m u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)}}{12} - \\frac{m \\left(q_{3} u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)} + 2 u_{2} u_{3}\\right) q_{3}}{4}\\\\- \\frac{g m \\sin{\\left(q_{1} \\right)} \\cos{\\left(q_{2} \\right)}}{4} - k_{l} q_{3} - \\frac{m \\left(- q_{3} u_{1}^{2} \\cos^{2}{\\left(q_{2} \\right)} - q_{3} u_{2}^{2}\\right)}{4}\\end{matrix}\\right]$"
],
"text/plain": [
"⎡ 2 ↪\n",
"⎢ 7⋅g⋅l⋅m⋅sin(q₁) g⋅m⋅q₃⋅cos(q₁)⋅cos(q₂) l ⋅m⋅u₁⋅u₂⋅sin(q₂)⋅cos(q ↪\n",
"⎢- ─────────────── - ────────────────────── - kₜ⋅q₁ + ──────────────────────── ↪\n",
"⎢ 4 4 6 ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ g⋅m ↪\n",
"⎢ ─── ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎣ ↪\n",
"\n",
"↪ 2 ⎛ ⎛ 2 2 2⎞ ⎛ ↪\n",
"↪ ₂) l⋅m⋅q₃⋅u₁ ⋅cos(q₂) ⎜ m⋅⎝- q₃⋅u₁ ⋅cos (q₂) - q₃⋅u₂ ⎠⋅cos(q₂) m⋅⎝q ↪\n",
"↪ ── - ────────────────── + l⋅⎜- ────────────────────────────────────── + ──── ↪\n",
"↪ 4 ⎝ 4 ↪\n",
"↪ ↪\n",
"↪ 2 2 ⎛ 2 ↪\n",
"↪ ⋅q₃⋅sin(q₁)⋅sin(q₂) l ⋅m⋅u₁ ⋅sin(q₂)⋅cos(q₂) m⋅⎝q₃⋅u₁ ⋅sin(q₂)⋅c ↪\n",
"↪ ─────────────────── - kₜ⋅q₂ - ──────────────────────── - ─────────────────── ↪\n",
"↪ 4 12 ↪\n",
"↪ ↪\n",
"↪ ⎛ 2 2 ↪\n",
"↪ g⋅m⋅sin(q₁)⋅cos(q₂) m⋅⎝- q₃⋅u₁ ⋅cos (q₂) - q₃⋅u₂ ↪\n",
"↪ - ─────────────────── - kₗ⋅q₃ - ──────────────────────────── ↪\n",
"↪ 4 4 ↪\n",
"\n",
"↪ 2 ⎞ ⎞ ↪\n",
"↪ ₃⋅u₁ ⋅sin(q₂)⋅cos(q₂) + 2⋅u₂⋅u₃⎠⋅sin(q₂)⎟ m⋅(2⋅q₃⋅u₁⋅u₂⋅sin(q₂) - 2⋅u₁⋅u₃⋅ ↪\n",
"↪ ────────────────────────────────────────⎟ + ──────────────────────────────── ↪\n",
"↪ 4 ⎠ 4 ↪\n",
"↪ ↪\n",
"↪ ⎞ ↪\n",
"↪ os(q₂) + 2⋅u₂⋅u₃⎠⋅q₃ ↪\n",
"↪ ──────────────────── ↪\n",
"↪ 4 ↪\n",
"↪ ↪\n",
"↪ 2⎞ ↪\n",
"↪ ⎠ ↪\n",
"↪ ── ↪\n",
"↪ ↪\n",
"\n",
"↪ ⎤\n",
"↪ cos(q₂))⋅q₃⋅cos(q₂)⎥\n",
"↪ ───────────────────⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎦"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gd = Frs.xreplace(ud_zerod) + Fr\n",
"gd"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left( \\left[\\begin{matrix}q_{1}\\\\q_{2}\\\\q_{3}\\end{matrix}\\right], \\ \\left[\\begin{matrix}u_{1}\\\\u_{2}\\\\u_{3}\\end{matrix}\\right], \\ \\left[\\begin{matrix}\\dot{q}_{1}\\\\\\dot{q}_{2}\\\\\\dot{q}_{3}\\end{matrix}\\right], \\ \\left[\\begin{matrix}\\dot{u}_{1}\\\\\\dot{u}_{2}\\\\\\dot{u}_{3}\\end{matrix}\\right]\\right)$"
],
"text/plain": [
"⎛⎡q₁⎤ ⎡u₁⎤ ⎡q₁̇⎤ ⎡u₁̇⎤⎞\n",
"⎜⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎟\n",
"⎜⎢q₂⎥, ⎢u₂⎥, ⎢q₂̇⎥, ⎢u₂̇⎥⎟\n",
"⎜⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎟\n",
"⎝⎣q₃⎦ ⎣u₃⎦ ⎣q₃̇⎦ ⎣u₃̇⎦⎠"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"q = sm.Matrix([q1, q2, q3])\n",
"qd = q.diff(t)\n",
"q, u, qd, ud"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left( \\left[\\begin{matrix}-1 & 0 & 0\\\\0 & -1 & 0\\\\0 & 0 & -1\\end{matrix}\\right], \\ \\left[\\begin{matrix}u_{1}\\\\u_{2}\\\\u_{3}\\end{matrix}\\right]\\right)$"
],
"text/plain": [
"⎛⎡-1 0 0 ⎤ ⎡u₁⎤⎞\n",
"⎜⎢ ⎥ ⎢ ⎥⎟\n",
"⎜⎢0 -1 0 ⎥, ⎢u₂⎥⎟\n",
"⎜⎢ ⎥ ⎢ ⎥⎟\n",
"⎝⎣0 0 -1⎦ ⎣u₃⎦⎠"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Mk = -sm.eye(3)\n",
"gk = u\n",
"Mk, gk"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left( \\left[\\begin{matrix}- \\frac{l^{2} m \\cos^{2}{\\left(q_{2} \\right)}}{12} - \\frac{19 l^{2} m}{12} - \\frac{m q_{3}^{2} \\cos^{2}{\\left(q_{2} \\right)}}{4} & \\frac{l m q_{3} \\sin{\\left(q_{2} \\right)}}{4} & - \\frac{l m \\cos{\\left(q_{2} \\right)}}{4}\\\\\\frac{l m q_{3} \\sin{\\left(q_{2} \\right)}}{4} & - \\frac{l^{2} m}{12} - \\frac{m q_{3}^{2}}{4} & 0\\\\- \\frac{l m \\cos{\\left(q_{2} \\right)}}{4} & 0 & - \\frac{m}{4}\\end{matrix}\\right], \\ \\left[\\begin{matrix}- \\frac{7 g l m \\sin{\\left(q_{1} \\right)}}{4} - \\frac{g m q_{3} \\cos{\\left(q_{1} \\right)} \\cos{\\left(q_{2} \\right)}}{4} - k_{t} q_{1} + \\frac{l^{2} m u_{1} u_{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)}}{6} - \\frac{l m q_{3} u_{1}^{2} \\cos{\\left(q_{2} \\right)}}{4} + l \\left(- \\frac{m \\left(- q_{3} u_{1}^{2} \\cos^{2}{\\left(q_{2} \\right)} - q_{3} u_{2}^{2}\\right) \\cos{\\left(q_{2} \\right)}}{4} + \\frac{m \\left(q_{3} u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)} + 2 u_{2} u_{3}\\right) \\sin{\\left(q_{2} \\right)}}{4}\\right) + \\frac{m \\left(2 q_{3} u_{1} u_{2} \\sin{\\left(q_{2} \\right)} - 2 u_{1} u_{3} \\cos{\\left(q_{2} \\right)}\\right) q_{3} \\cos{\\left(q_{2} \\right)}}{4}\\\\\\frac{g m q_{3} \\sin{\\left(q_{1} \\right)} \\sin{\\left(q_{2} \\right)}}{4} - k_{t} q_{2} - \\frac{l^{2} m u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)}}{12} - \\frac{m \\left(q_{3} u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)} + 2 u_{2} u_{3}\\right) q_{3}}{4}\\\\- \\frac{g m \\sin{\\left(q_{1} \\right)} \\cos{\\left(q_{2} \\right)}}{4} - k_{l} q_{3} - \\frac{m \\left(- q_{3} u_{1}^{2} \\cos^{2}{\\left(q_{2} \\right)} - q_{3} u_{2}^{2}\\right)}{4}\\end{matrix}\\right]\\right)$"
],
"text/plain": [
"⎛ ↪\n",
"⎜⎡ 2 2 2 2 2 ⎤ ↪\n",
"⎜⎢ l ⋅m⋅cos (q₂) 19⋅l ⋅m m⋅q₃ ⋅cos (q₂) l⋅m⋅q₃⋅sin(q₂) -l⋅m⋅cos(q₂) ⎥ ↪\n",
"⎜⎢- ───────────── - ─────── - ────────────── ────────────── ─────────────⎥ ↪\n",
"⎜⎢ 12 12 4 4 4 ⎥ ↪\n",
"⎜⎢ ⎥ ↪\n",
"⎜⎢ 2 2 ⎥ ↪\n",
"⎜⎢ l⋅m⋅q₃⋅sin(q₂) l ⋅m m⋅q₃ ⎥, ↪\n",
"⎜⎢ ────────────── - ──── - ───── 0 ⎥ ↪\n",
"⎜⎢ 4 12 4 ⎥ ↪\n",
"⎜⎢ ⎥ ↪\n",
"⎜⎢ -l⋅m⋅cos(q₂) -m ⎥ ↪\n",
"⎜⎢ ───────────── 0 ─── ⎥ ↪\n",
"⎝⎣ 4 4 ⎦ ↪\n",
"\n",
"↪ ⎡ 2 ↪\n",
"↪ ⎢ 7⋅g⋅l⋅m⋅sin(q₁) g⋅m⋅q₃⋅cos(q₁)⋅cos(q₂) l ⋅m⋅u₁⋅u₂⋅sin(q₂)⋅cos ↪\n",
"↪ ⎢- ─────────────── - ────────────────────── - kₜ⋅q₁ + ────────────────────── ↪\n",
"↪ ⎢ 4 4 6 ↪\n",
"↪ ⎢ ↪\n",
"↪ ⎢ ↪\n",
"↪ ⎢ g ↪\n",
"↪ ⎢ ─ ↪\n",
"↪ ⎢ ↪\n",
"↪ ⎢ ↪\n",
"↪ ⎢ ↪\n",
"↪ ⎢ ↪\n",
"↪ ⎢ ↪\n",
"↪ ⎣ ↪\n",
"\n",
"↪ 2 ⎛ ⎛ 2 2 2⎞ ↪\n",
"↪ (q₂) l⋅m⋅q₃⋅u₁ ⋅cos(q₂) ⎜ m⋅⎝- q₃⋅u₁ ⋅cos (q₂) - q₃⋅u₂ ⎠⋅cos(q₂) m⋅ ↪\n",
"↪ ──── - ────────────────── + l⋅⎜- ────────────────────────────────────── + ── ↪\n",
"↪ 4 ⎝ 4 ↪\n",
"↪ ↪\n",
"↪ 2 2 ⎛ 2 ↪\n",
"↪ ⋅m⋅q₃⋅sin(q₁)⋅sin(q₂) l ⋅m⋅u₁ ⋅sin(q₂)⋅cos(q₂) m⋅⎝q₃⋅u₁ ⋅sin(q₂) ↪\n",
"↪ ───────────────────── - kₜ⋅q₂ - ──────────────────────── - ───────────────── ↪\n",
"↪ 4 12 ↪\n",
"↪ ↪\n",
"↪ ⎛ 2 2 ↪\n",
"↪ g⋅m⋅sin(q₁)⋅cos(q₂) m⋅⎝- q₃⋅u₁ ⋅cos (q₂) - q₃⋅ ↪\n",
"↪ - ─────────────────── - kₗ⋅q₃ - ────────────────────────── ↪\n",
"↪ 4 4 ↪\n",
"\n",
"↪ ⎛ 2 ⎞ ⎞ ↪\n",
"↪ ⎝q₃⋅u₁ ⋅sin(q₂)⋅cos(q₂) + 2⋅u₂⋅u₃⎠⋅sin(q₂)⎟ m⋅(2⋅q₃⋅u₁⋅u₂⋅sin(q₂) - 2⋅u₁⋅u ↪\n",
"↪ ──────────────────────────────────────────⎟ + ────────────────────────────── ↪\n",
"↪ 4 ⎠ 4 ↪\n",
"↪ ↪\n",
"↪ ⎞ ↪\n",
"↪ ⋅cos(q₂) + 2⋅u₂⋅u₃⎠⋅q₃ ↪\n",
"↪ ────────────────────── ↪\n",
"↪ 4 ↪\n",
"↪ ↪\n",
"↪ 2⎞ ↪\n",
"↪ u₂ ⎠ ↪\n",
"↪ ──── ↪\n",
"↪ ↪\n",
"\n",
"↪ ⎤⎞\n",
"↪ ₃⋅cos(q₂))⋅q₃⋅cos(q₂)⎥⎟\n",
"↪ ─────────────────────⎥⎟\n",
"↪ ⎥⎟\n",
"↪ ⎥⎟\n",
"↪ ⎥⎟\n",
"↪ ⎥⎟\n",
"↪ ⎥⎟\n",
"↪ ⎥⎟\n",
"↪ ⎥⎟\n",
"↪ ⎥⎟\n",
"↪ ⎥⎟\n",
"↪ ⎥⎟\n",
"↪ ⎦⎠"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Md, gd"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left\\{g, k_{l}, k_{t}, l, m, t\\right\\}$"
],
"text/plain": [
"{g, kₗ, kₜ, l, m, t}"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Mk.free_symbols | gk.free_symbols | Md.free_symbols | gd.free_symbols"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}g\\\\k_{l}\\\\k_{t}\\\\l\\\\m\\end{matrix}\\right]$"
],
"text/plain": [
"⎡g ⎤\n",
"⎢ ⎥\n",
"⎢kₗ⎥\n",
"⎢ ⎥\n",
"⎢kₜ⎥\n",
"⎢ ⎥\n",
"⎢l ⎥\n",
"⎢ ⎥\n",
"⎣m ⎦"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p = sm.Matrix([g, kl, kt, l, m])\n",
"p"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"eval_eom = sm.lambdify((q, u, p), [Mk, gk, Md, gd])"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(array([0.43633231, 0.08726646, 0.1 ]), numpy.ndarray, (3,))"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"q_vals = np.array([\n",
" np.deg2rad(25.0), # q1, rad\n",
" np.deg2rad(5.0), # q2, rad\n",
" 0.1, # q3, m\n",
"])\n",
"q_vals, type(q_vals), q_vals.shape"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(array([0.1, 2.2, 0.3]), numpy.ndarray, (3,))"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"u_vals = np.array([\n",
" 0.1, # u1, rad/s\n",
" 2.2, # u2, rad/s\n",
" 0.3, # u3, m/s\n",
"])\n",
"u_vals, type(u_vals), u_vals.shape"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(array([9.81, 2. , 0.01, 0.6 , 1. ]), numpy.ndarray, (5,))"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p_vals = np.array([\n",
" 9.81, # g, m/s**2\n",
" 2.0, # kl, N/m\n",
" 0.01, # kt, Nm/rad\n",
" 0.6, # l, m\n",
" 1.0, # m, kg\n",
"])\n",
"p_vals, type(p_vals), p_vals.shape"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(array([[-1, 0, 0],\n",
" [ 0, -1, 0],\n",
" [ 0, 0, -1]]),\n",
" array([[0.1],\n",
" [2.2],\n",
" [0.3]]),\n",
" array([[-0.60225313, 0.00130734, -0.1494292 ],\n",
" [ 0.00130734, -0.0325 , 0. ],\n",
" [-0.1494292 , 0. , -0.25 ]]),\n",
" array([[-4.48963535],\n",
" [-0.02486744],\n",
" [-1.1112791 ]]))"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Mk_vals, gk_vals, Md_vals, gd_vals = eval_eom(q_vals, u_vals, p_vals)\n",
"Mk_vals, gk_vals, Md_vals, gd_vals"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.1, 2.2, 0.3])"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"qd_vals = np.linalg.solve(-Mk_vals, np.squeeze(gk_vals))\n",
"qd_vals"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([-7.46056427, -1.06525862, 0.01418834])"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ud_vals = np.linalg.solve(-Md_vals, np.squeeze(gd_vals))\n",
"ud_vals"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"def eval_rhs(t, x, p):\n",
" \"\"\"Return the right hand side of the explicit ordinary differential\n",
" equations which evaluates the time derivative of the state ``x`` at time\n",
" ``t``.\n",
"\n",
" Parameters\n",
" ==========\n",
" t : float\n",
" Time in seconds.\n",
" x : array_like, shape(6,)\n",
" State at time t: [q1, q2, q3, u1, u2, u3]\n",
" p : array_like, shape(5,)\n",
" Constant parameters: [g, kl, kt, l, m]\n",
"\n",
" Returns\n",
" =======\n",
" xd : ndarray, shape(6,)\n",
" Derivative of the state with respect to time at time ``t``.\n",
"\n",
" \"\"\"\n",
"\n",
" # unpack the q and u vectors from x\n",
" q = x[:3]\n",
" u = x[3:]\n",
"\n",
" # evaluate the equations of motion matrices with the values of q, u, p\n",
" Mk, gk, Md, gd = eval_eom(q, u, p)\n",
"\n",
" # solve for q' and u'\n",
" qd = np.linalg.solve(-Mk, np.squeeze(gk))\n",
" ud = np.linalg.solve(-Md, np.squeeze(gd))\n",
"\n",
" # pack dq/dt and du/dt into a new state time derivative vector dx/dt\n",
" xd = np.empty_like(x)\n",
" xd[:3] = qd\n",
" xd[3:] = ud\n",
"\n",
" return xd"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"def euler_integrate(rhs_func, tspan, x0_vals, p_vals, delt=0.03):\n",
" \"\"\"Returns state trajectory and corresponding values of time resulting\n",
" from integrating the ordinary differential equations with Euler's\n",
" Method.\n",
"\n",
" Parameters\n",
" ==========\n",
" rhs_func : function\n",
" Python function that evaluates the derivative of the state and takes\n",
" this form ``dxdt = f(t, x, p)``.\n",
" tspan : 2-tuple of floats\n",
" The initial time and final time values: (t0, tf).\n",
" x0_vals : array_like, shape(2*n,)\n",
" Values of the state x at t0.\n",
" p_vals : array_like, shape(o,)\n",
" Values of constant parameters.\n",
" delt : float\n",
" Integration time step in seconds/step.\n",
"\n",
" Returns\n",
" =======\n",
" ts : ndarray(m, )\n",
" Monotonically equally spaced increasing values of time spaced apart\n",
" by ``delt``.\n",
" xs : ndarray(m, 2*n)\n",
" State values at each time in ts.\n",
"\n",
" \"\"\"\n",
" # generate monotonically increasing values of time.\n",
" duration = tspan[1] - tspan[0]\n",
" num_samples = round(duration/delt) + 1\n",
" ts = np.arange(tspan[0], tspan[0] + delt*num_samples, delt)\n",
"\n",
" # create an empty array to hold the state values.\n",
" x = np.empty((len(ts), len(x0_vals)))\n",
"\n",
" # set the initial conditions to the first element.\n",
" x[0, :] = x0_vals\n",
"\n",
" # use a for loop to sequentially calculate each new x.\n",
" for i, ti in enumerate(ts[:-1]):\n",
" x[i + 1, :] = x[i, :] + delt*rhs_func(ti, x[i, :], p_vals)\n",
"\n",
" return ts, x"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
"x0 = np.empty(6)\n",
"x0[:3] = q_vals\n",
"x0[3:] = u_vals\n",
"\n",
"t0 = 0.1"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 0.1 , 2.2 , 0.3 , -7.46056427, -1.06525862,\n",
" 0.01418834])"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"eval_rhs(t0, x0, p_vals)"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [],
"source": [
"tf = 2.0\n",
"ts, xs = euler_integrate(eval_rhs, (t0, tf), x0, p_vals)"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.1 , 0.13, 0.16, 0.19, 0.22, 0.25, 0.28, 0.31, 0.34, 0.37, 0.4 ,\n",
" 0.43, 0.46, 0.49, 0.52, 0.55, 0.58, 0.61, 0.64, 0.67, 0.7 , 0.73,\n",
" 0.76, 0.79, 0.82, 0.85, 0.88, 0.91, 0.94, 0.97, 1. , 1.03, 1.06,\n",
" 1.09, 1.12, 1.15, 1.18, 1.21, 1.24, 1.27, 1.3 , 1.33, 1.36, 1.39,\n",
" 1.42, 1.45, 1.48, 1.51, 1.54, 1.57, 1.6 , 1.63, 1.66, 1.69, 1.72,\n",
" 1.75, 1.78, 1.81, 1.84, 1.87, 1.9 , 1.93, 1.96, 1.99])"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ts"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 4.36332313e-01, 8.72664626e-02, 1.00000000e-01,\n",
" 1.00000000e-01, 2.20000000e+00, 3.00000000e-01],\n",
" [ 4.39332313e-01, 1.53266463e-01, 1.09000000e-01,\n",
" -1.23816928e-01, 2.16804224e+00, 3.00425650e-01],\n",
" [ 4.35617805e-01, 2.18307730e-01, 1.18012770e-01,\n",
" -3.48867554e-01, 2.13303503e+00, 2.99423633e-01],\n",
" [ 4.25151779e-01, 2.82298781e-01, 1.26995479e-01,\n",
" -5.72475739e-01, 2.09463526e+00, 2.97361542e-01],\n",
" [ 4.07977506e-01, 3.45137839e-01, 1.35916325e-01,\n",
" -7.91929289e-01, 2.05197871e+00, 2.94628416e-01],\n",
" [ 3.84219628e-01, 4.06697200e-01, 1.44755177e-01,\n",
" -1.00444789e+00, 2.00382451e+00, 2.91554565e-01],\n",
" [ 3.54086191e-01, 4.66811935e-01, 1.53501814e-01,\n",
" -1.20715867e+00, 1.94877129e+00, 2.88335792e-01],\n",
" [ 3.17871431e-01, 5.25275074e-01, 1.62151888e-01,\n",
" -1.39708735e+00, 1.88551903e+00, 2.84974403e-01],\n",
" [ 2.75958810e-01, 5.81840645e-01, 1.70701120e-01,\n",
" -1.57117194e+00, 1.81314156e+00, 2.81247424e-01],\n",
" [ 2.28823652e-01, 6.36234892e-01, 1.79138543e-01,\n",
" -1.72630399e+00, 1.73132870e+00, 2.76708292e-01],\n",
" [ 1.77034532e-01, 6.88174753e-01, 1.87439792e-01,\n",
" -1.85939863e+00, 1.64055711e+00, 2.70722509e-01],\n",
" [ 1.21252574e-01, 7.37391466e-01, 1.95561467e-01,\n",
" -1.96749077e+00, 1.54215738e+00, 2.62531540e-01],\n",
" [ 6.22278505e-02, 7.83656188e-01, 2.03437413e-01,\n",
" -2.04785035e+00, 1.43826241e+00, 2.51334050e-01],\n",
" [ 7.92340069e-04, 8.26804060e-01, 2.10977435e-01,\n",
" -2.09810628e+00, 1.33164459e+00, 2.36370890e-01],\n",
" [-6.21508484e-02, 8.66753398e-01, 2.18068561e-01,\n",
" -2.11636717e+00, 1.22547077e+00, 2.17000797e-01],\n",
" [-1.25641864e-01, 9.03517521e-01, 2.24578585e-01,\n",
" -2.10132625e+00, 1.12301845e+00, 1.92757262e-01],\n",
" [-1.88681651e-01, 9.37208074e-01, 2.30361303e-01,\n",
" -2.05234005e+00, 1.02739894e+00, 1.63382255e-01],\n",
" [-2.50251853e-01, 9.68030042e-01, 2.35262771e-01,\n",
" -1.96947297e+00, 9.41324864e-01, 1.28837671e-01],\n",
" [-3.09336042e-01, 9.96269788e-01, 2.39127901e-01,\n",
" -1.85350427e+00, 8.66944026e-01, 8.92990821e-02],\n",
" [-3.64941170e-01, 1.02227811e+00, 2.41806873e-01,\n",
" -1.70589822e+00, 8.05745012e-01, 4.51378231e-02],\n",
" [-4.16118116e-01, 1.04645046e+00, 2.43161008e-01,\n",
" -1.52874283e+00, 7.58527111e-01, -3.10326579e-03],\n",
" [-4.61980401e-01, 1.06920627e+00, 2.43067910e-01,\n",
" -1.32466575e+00, 7.25419993e-01, -5.47370444e-02],\n",
" [-5.01720374e-01, 1.09096887e+00, 2.41425799e-01,\n",
" -1.09673819e+00, 7.05937072e-01, -1.08961098e-01],\n",
" [-5.34622520e-01, 1.11214698e+00, 2.38156966e-01,\n",
" -8.48377776e-01, 6.99048783e-01, -1.64889059e-01],\n",
" [-5.60073853e-01, 1.13311845e+00, 2.33210294e-01,\n",
" -5.83259685e-01, 7.03266085e-01, -2.21585204e-01],\n",
" [-5.77571643e-01, 1.15421643e+00, 2.26562738e-01,\n",
" -3.05242994e-01, 7.16728855e-01, -2.78102687e-01],\n",
" [-5.86728933e-01, 1.17571830e+00, 2.18219657e-01,\n",
" -1.83153063e-02, 7.37297421e-01, -3.33524242e-01],\n",
" [-5.87278392e-01, 1.19783722e+00, 2.08213930e-01,\n",
" 2.73444426e-01, 7.62647870e-01, -3.87002450e-01],\n",
" [-5.79075060e-01, 1.22071666e+00, 1.96603856e-01,\n",
" 5.65888142e-01, 7.90372680e-01, -4.37795193e-01],\n",
" [-5.62098415e-01, 1.24442784e+00, 1.83470001e-01,\n",
" 8.54811324e-01, 8.18087698e-01, -4.85291412e-01],\n",
" [-5.36454076e-01, 1.26897047e+00, 1.68911258e-01,\n",
" 1.13596271e+00, 8.43544507e-01, -5.29023016e-01],\n",
" [-5.02375194e-01, 1.29427680e+00, 1.53040568e-01,\n",
" 1.40505807e+00, 8.64744480e-01, -5.68660838e-01],\n",
" [-4.60223452e-01, 1.32021914e+00, 1.35980743e-01,\n",
" 1.65780471e+00, 8.80047808e-01, -6.03995480e-01],\n",
" [-4.10489311e-01, 1.34662057e+00, 1.17860878e-01,\n",
" 1.88994245e+00, 8.88268730e-01, -6.34907053e-01],\n",
" [-3.53791037e-01, 1.37326863e+00, 9.88136667e-02,\n",
" 2.09730514e+00, 8.88747420e-01, -6.61330213e-01],\n",
" [-2.90871883e-01, 1.39993105e+00, 7.89737603e-02,\n",
" 2.27590461e+00, 8.81389870e-01, -6.83221797e-01],\n",
" [-2.22594745e-01, 1.42637275e+00, 5.84771064e-02,\n",
" 2.42203585e+00, 8.66668771e-01, -7.00537408e-01],\n",
" [-1.49933669e-01, 1.45237281e+00, 3.74609842e-02,\n",
" 2.53239896e+00, 8.45580400e-01, -7.13220652e-01],\n",
" [-7.39617005e-02, 1.47774023e+00, 1.60643646e-02,\n",
" 2.60422940e+00, 8.19554682e-01, -7.21205288e-01],\n",
" [ 4.16518154e-03, 1.50232687e+00, -5.57179403e-03,\n",
" 2.63542555e+00, 7.90319326e-01, -7.24427319e-01],\n",
" [ 8.32279481e-02, 1.52603645e+00, -2.73046136e-02,\n",
" 2.62466048e+00, 7.59725678e-01, -7.22841890e-01],\n",
" [ 1.61967762e-01, 1.54882822e+00, -4.89898703e-02,\n",
" 2.57146487e+00, 7.29553925e-01, -7.16439450e-01],\n",
" [ 2.39111708e-01, 1.57071483e+00, -7.04830538e-02,\n",
" 2.47627011e+00, 7.01325433e-01, -7.05256667e-01],\n",
" [ 3.13399812e-01, 1.59175460e+00, -9.16407538e-02,\n",
" 2.34040397e+00, 6.76154998e-01, -6.89379799e-01],\n",
" [ 3.83611931e-01, 1.61203925e+00, -1.12322148e-01,\n",
" 2.16603688e+00, 6.54671168e-01, -6.68940608e-01],\n",
" [ 4.48593037e-01, 1.63167938e+00, -1.32390366e-01,\n",
" 1.95608296e+00, 6.37018712e-01, -6.44106831e-01],\n",
" [ 5.07275526e-01, 1.65078994e+00, -1.51713571e-01,\n",
" 1.71406620e+00, 6.22938813e-01, -6.15070314e-01],\n",
" [ 5.58697512e-01, 1.66947811e+00, -1.70165680e-01,\n",
" 1.44396681e+00, 6.11906749e-01, -5.82036029e-01],\n",
" [ 6.02016516e-01, 1.68783531e+00, -1.87626761e-01,\n",
" 1.15006549e+00, 6.03298605e-01, -5.45214730e-01],\n",
" [ 6.36518481e-01, 1.70593427e+00, -2.03983203e-01,\n",
" 8.36802201e-01, 5.96558484e-01, -5.04821104e-01],\n",
" [ 6.61622547e-01, 1.72383102e+00, -2.19127836e-01,\n",
" 5.08663225e-01, 5.91343485e-01, -4.61078217e-01],\n",
" [ 6.76882444e-01, 1.74157133e+00, -2.32960183e-01,\n",
" 1.70104973e-01, 5.87631757e-01, -4.14227908e-01],\n",
" [ 6.81985593e-01, 1.75920028e+00, -2.45387020e-01,\n",
" -1.74482332e-01, 5.85786780e-01, -3.64545604e-01],\n",
" [ 6.76751123e-01, 1.77677388e+00, -2.56323388e-01,\n",
" -5.20773373e-01, 5.86577261e-01, -3.12357072e-01],\n",
" [ 6.61127922e-01, 1.79437120e+00, -2.65694100e-01,\n",
" -8.64475606e-01, 5.91156371e-01, -2.58054046e-01],\n",
" [ 6.35193654e-01, 1.81210589e+00, -2.73435722e-01,\n",
" -1.20127673e+00, 6.01006553e-01, -2.02105739e-01],\n",
" [ 5.99155352e-01, 1.83013609e+00, -2.79498894e-01,\n",
" -1.52678495e+00, 6.17857094e-01, -1.45063860e-01],\n",
" [ 5.53351803e-01, 1.84867180e+00, -2.83850810e-01,\n",
" -1.83647709e+00, 6.43581213e-01, -8.75598764e-02],\n",
" [ 4.98257491e-01, 1.86797924e+00, -2.86477606e-01,\n",
" -2.12567032e+00, 6.80078039e-01, -3.02944186e-02],\n",
" [ 4.34487381e-01, 1.88838158e+00, -2.87386438e-01,\n",
" -2.38953235e+00, 7.29142906e-01, 2.59803906e-02],\n",
" [ 3.62801410e-01, 1.91025587e+00, -2.86607027e-01,\n",
" -2.62314224e+00, 7.92327527e-01, 8.04846656e-02],\n",
" [ 2.84107143e-01, 1.93402569e+00, -2.84192487e-01,\n",
" -2.82160885e+00, 8.70790639e-01, 1.32439931e-01],\n",
" [ 1.99458878e-01, 1.96014941e+00, -2.80219289e-01,\n",
" -2.98024641e+00, 9.65140541e-01, 1.81104004e-01],\n",
" [ 1.10051485e-01, 1.98910363e+00, -2.74786169e-01,\n",
" -3.09479777e+00, 1.07527459e+00, 2.25812638e-01]])"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xs"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x1e335aa5870>,\n",
" <matplotlib.lines.Line2D at 0x1e335aa5960>,\n",
" <matplotlib.lines.Line2D at 0x1e335aa5a50>,\n",
" <matplotlib.lines.Line2D at 0x1e335aa5b40>,\n",
" <matplotlib.lines.Line2D at 0x1e335aa5c30>,\n",
" <matplotlib.lines.Line2D at 0x1e335aa5d20>]"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAAsTAAALEwEAmpwYAABaA0lEQVR4nO3dd3hb1fnA8e/VsGzZ8t57ZC+yCCQkBEgCYY+wwoaWWVZpKVAo8CulQKEUKKVllx02LSEQMlgNZJG9vfdeki1rn98f13EcyLalK9nn8zz3kSxd6b6+vn51dKYihECSJEkKXTqtA5AkSZL6RiZySZKkECcTuSRJUoiTiVySJCnEyUQuSZIU4gxaHDQxMVHk5uZqcWhJkqSQ9eOPPzYJIZJ++rgmiTw3N5e1a9dqcWhJkqSQpShK+b4el1UrkiRJIU4mckmSpBAnE7kkSVKIk4lckiQpxMlELkmSFOJkIpckSQpxMpFLkiSFOE36kUuSdGDetja6Nm7EsWMniikMQ1wc+rg49LGxGFJSMKakaB2iFERkIpekICDcbqyff07nDyvp2rABV2npAfePOOooYs47j+jTTkVvsQQoSilYKVosLDF58mQhR3ZKEgiPh/b/fkrTc8/hrqpCHxtLxPjxPVv4mDHg9eBta8Pb2oqntRVXcTFtn3yCq6gYJTwcy8lziLvoIsyTJmn960h+pijKj0KIyT97XCZySQo84fNhXfQ5Tc8+i6usjPBRo0i89RaiZs5EUZSDv14IHJs30/bRR1g/W4TPZiP2gvNJvusu9FFRAfgNJC3IRC5JQcJdXU3V7b/GsXkzpqFDSbz1FiyzZx9SAt8Xn8NB0z+eo/nllzGmppL25z8Teewx/Ry1FAz2l8hlrxVJCqCOFSsonXc+rtJS0v/yGHn/+YToOXOOOIkD6MLDSf7NHeS89SaK0UjFVVdR9/Cf8XV19WPkUjCTiVySAkAIQdMLL1J57XUYkhLJ++B9Ys46C0XXf/+C5gkTyPv4I+Iuu4zWN96g7KKL8TQ19dv7S8FLJnJJ8jNvRyfVt95G45NPEj33FHIXLCDMT/Px68xmUu+7l6wXX8RVWUn55Vfgrq/3y7Gk4CETuST5kddmo+LKK7EtX07yXXeR/te/oouM9Ptxo2ZMJ/ulF/E0NKjJvKbG78eUtCMTuST5ia+zk8rrrsexaxeZz/6dhKuv6lNd+OEyT5pE9isv421tpfyyy3FVVgbs2FJgyUQuSX7gczio/NXNdG3cSMYTT2A58URN4og46iiyX30VX2cn5ZddjvMgA42k0CQTuST1M+FyUXXbbdhXrSL90UeIPuVkTeOJGDOa7NdfQ7jdVF53Pd62Nk3jkfqfTOSS1I+Ex0P1b++k85tvSX3wQWLOOkvrkAAIHz6crOf+gaeujuo7foPweLQOSepHMpFLUj+q//OfsX35JSn33E3cRRdqHc5eIsaPJ+X+P9D5/fc0/O1vWocj9SM5aZYk9ZO2Tz6h9e13iL/6auKvvFLrcPYp7oILcG7fTsvLrxA+ahQxp5+udUhSPwipEnn51mY2Lq+kodyK1+vTOhxJ6uHYto26Bx7EPGUKyb+5Q+twDijl7ruJmDSJ2nvvw7F9u9bhSP0gpOZa+frtnWz9thoAQ5iOlNxoUvNjSBsaS1pBDGHh8guGFHjetjZKz78A4XaT99GHGBIStA7poDxNTZTOOx9Fryf3ww8wxMVpHZJ0CAbMpFkdrQ5qi9upK26nrqSdxsoOhE+g6BSSsi1kDIslY1gcaUNkYpf8T3i9VN5wI50rV5L7xutEjB+vdUiHrGvzZsovvYzIGTPIfPbvAe3jLh2Z/SXykMt0UXHhDJ0cztDJ6gopbqeXupJ2qne1UrOrjY3LKln/ZQU6vUJqfgzZo+PJGhlPUpYFRScvVKl/Nf3jH3R+9x2pDz4QUkkcIGLsWJJuu5WGx5/AuvAzYs48Q+uQpCMUciXyg3E7vdQVt1O5o4XK7S00VXYAEB5lJGd0ArnjEskeHS9L61KfdXz3PyqvvZaYc88l7c8Ph2SJVni9lF9yKa6yMvIXfoohKUnrkKQD8FvViqIoWcDrQAoggBeEEE8f6DWBnI+8s91J1Y5WKrY2U761GWenB51BIXNYHLnjEimYmIw5OiwgsUgDh9dqpeTMs9BFRZH34QfowsO1DumIOUtKKD3nXCKPn0Hm32UVSzDzZyJPA9KEEOsURbEAPwLnCCG27e81Wi0s4fP6qCtpp3RjE6Ubm2hv7EJRIGN4HEMnp5A/IYnwSGPA45JCT83d99D+6afkLniHiLFjtQ6nz5pffoWGxx8n/fHHZRVLEAtYY6eiKP8BnhVCLNnfPsGwQpAQgpbaTorWNlC4pp72xi50eoXsUfGMPC6dnLEJ6PUh1TtTChDb8uVU3fQrEm64nuTbb9c6nH4hq1gCwN0FtZsgZRSYjmzB7IAkckVRcoFvgTFCCOtPnrsOuA4gOzt7Unl5eb8dt6+EEDRVdlC4pp5dq+vobHdhjg5jxNQ0Rh6XRmyyWesQpSDhaW2l5MyzMCQkkPf+eyhhA6darqeKRfZi6TuPExq2Q816qFkH1euhYRsIL1zyPgw7svl3/J7IFUWJAr4BHhZCfHSgfYOhRL4/Pq+P8q0tbPtfDeVbmhE+QeaIOMbPziZ7dLy8uAe56jt+g/XLL8l7/z3CR47UOpx+t7uKJeOpvxE9d67W4QQ/IaCjARq3Q90WqNusbk07wdc9n014LGRMhPSJkD4BcqaBOf6IDufXRK4oihFYCCwWQjx5sP2DOZH31tnmZPv3tWz5tprONidxaZGMn53F8Cmp6I2y2mWwsX6xmOrbbyfx1ltIuukmrcPxC+H1UjrvfLzWdgoWLQrpRtx+JQTYaqFpFzQVQuMOtcTdsB26WvbsF5UKqWO7tzFq4o7Lg34qAPqzsVMBXgNahBC3H8prQiWR7+b1+Cj6sYH1SyporuogIjqM8bOyGHtCJkaTXuvwpADwtLZSctrpGNPTyV3wDopx4DaKd65eTcUVVw7oD6z96mqDlmJoLum+LVK3pkJwdezZzxQNSSMgeaS6JY2AlNEQlezX8PyZyKcD3wGbgd0ToPxeCLFof68JtUS+mxCCqp2trP+ygsptLURYjEyam8vo49MxGGVCH8hq7ruP9k/+Q95HHxI+bJjW4fhd1e2/puPrryn4fBHGtDStw+k/Ph901EFrObSWQktp922Jer936RoFYrMgvgASh0Hi0O7bYWBJ7bdS9uEYMEP0g0VtcTur/ltC9c5WouJMTD4tlxHT0mRPlwGoa8MGyi6eT/w115Dyuzu1DicgXFXVlJx+OpY5c8h44nGtwzk8XW3QVq4m671uy9T7XueefRUdxGSq1R/x+RCfpybuhCEQlwvG4KpakoncT6p2tLDqvyXUlViJSzUz/YKhZI8O/kmTpEMjvF5KL7gAb1Mz+YsWoY/y/8LJwaLxmWdoeu6f5Lz9FuaJE7UOZw+HFdoq1OTcVqFurb3uO9v33j88BmJz1MQct/s2F2JzITYbDKHT80gmcj8SQlC2uZkVHxTS3tBF7rhEjjt/iOy2OAC0vPUW9Q/9iYy/PUn0qadqHU5A+ex2ik89DUNiIrnvv4eiC9C3TadtT1Lu2cr3JGxH2977G81qoo7NVhN1TNaehB2bAxGxgYk7AGQiDwCv28fG5ZWsXVSG1+tj/KxsJp+WKxtEQ5SnuZniuacSMXYMWS+/PCi7nrZ/upCaO+8k7eE/ETtvXv+8qavz50m6d4l6r3pqwBChJundW1xOr59z1a58g+RvIxN5AHW2O/nh42J2rqzDkhDOiZeNIGvkkfUblbRTc8/vaV+4kPz/fIIpP1/rcDQhhFBHfFZUULD4C/RRUQd/kccF7ZVqnfReSbr7vr1p7/31pn0k6Jw9t5GJgyZRH4xM5BqoKWzjqzd30FZvZ+S0NKbNGyLncgkR9nXrKL/kUhKuvTboV/zxt65Nmyi78KI93RGFgI76PY2HPQm7TN2sNajz53XTGdXeH70T9O5qj9hsiEyCQFXbhDiZyDXicXtZ81kZ67+sICLKyPHzh1Ewwb99TaW+EV4vpefNw2u1UvDZQnTmQdjW4XHuSdKtpVQ+sQD7rgaGXBGJvqsCPF17729J7y5R5+xdPx2XA5Y00Mnqxf4wYBaWCDUGo56p5xQwZGIyy9/YzhfPb2HYMSkcf/FwTBHy9Aej9o8/xrlzJxlP/W1gJ3GXvbv/9E+3UrBW07tUnZRjoXSzhebtZpIv+MWeRB2fpzYuBlk3vcFGZpIAScq2cP7dk/nx83LWLiqjprCNOVePIn2oXCsxmPg6O2l8+hkiJkzAcsopWofTdx6XWu3RVLhnlGJLCTQXg61m733NiWpizpm2p091XB7E5xEemUQ0v6Nl2TLiH7ldzo4YZGQiDyC9XseUM/LIHh3P0le28fGT65l4cjZTzsxHb5B1hMGg+dV/42lsJOOZp0Orl4q9pXsekF175gNpKlSrRoR3z37mBHXAS/7M7oEv+eptfJ7a3/oAkm7+FdbPP6fphRdJvff3/v19pMMiE7kGUvNiuPDeo1nxQRHrFldQsa2FU345htiUAfw1PgS4GxpofvllLHPnYp4wQetwfq73THuNO3ttO/buCaI3qSMTU8fA6HPV+4lDIaEAIo78G2BYbi4x555D24IFJFxz9cAauh/iZGOnxko2NLL8je14PYITLhnO8GNStQ5p0Kr9wx9o++Q/FHy2kLDsbG2D6WqF+m3qHNa7Z9lr3K4+vlt4jDpZU+IwSBoOicPVhB2b7bfGRXdNDcWnzCXmnHNIe+iPfjmGtH+ysTNI5Y9PIinbwpJXtrL01W1U72xlxkXD5CCiAHPs3EXbhx8Rf/nlgU3iHhc0F6pzWddvURN3/ba9669NMeoMe6POhqSRkDxCTeBRKQHvX21MTyf2ootofecdEn75C8JycgJ6fGnfZIk8SPi8PlYvLOXHL8qJS43klF+OJiHjEAZfSP2i4trr6Nq4kSFfLkYfG+ufg3Q2Q333wgO7E3fjTvC51ef1YWrJOnm0uhxY8mg1gUenB9WAGE9jI0VzTlYn1Hr8L1qHM6jIEnmQ0+l1HHt2ARlD41jy6lY+eHQtMy8Zzoipsh7S3zpWrKDzu+9Ivuuu/kniQqhTo+5eLaZ2k3rbu5RtSVPnrx4yW12EIGW0WpetD/4BY4akJOIvu5Tml18h8cYbMeXnaR3SoCdL5EGos93Jkpe3Ur2rjZHHpXH8RcMwhMmqFn8QPh+l887HZ7ORv+gzdIe7BqfHpTY21m2Guk1q0q7fAs7uJWsVvVrK3r1qTMoY9TYysf9/mQDyNDdTdNIsos84nfSHH9Y6nEFDlshDSGSMibNuG8/qT9WqloZyG3Ovlb1a/MH2xRc4t28n/fG/HDyJO6xqku4pZW9SGyF3V40YzWqiHnchpI5TE3byqAE5WMaQkEDsBRfQumABSb/6Fcb0dK1DGtRkiTzIlW1uYum/t+HzCk66fCRDJsnh/f1FeDyUnHEmitFA3iefoOi7v/XsXp9xdyl7dxVJS8meF5sTIW2cmrB338bnD6qh6O6aGopOPoW4+fNlv/IAkSXyEJU7NpGL7p3C4he3sPjFLdQWZzLtvCFyAFE/aP/Pf3GVlZHx0G9Rtry/J2HXbwF7854d4/LU0vVR8/ckbktaUDVAasGYnk7MmWfS9v77JN5wPYYEuaCKVmSJPER4PT5++KiYjcsrSc6N5pRrRxOdEKF1WKFDCLDVQcNWqN+Gr3oTJY//gN7oJHdOo5qT9Sa1l0jPKujdddrh0VpHH7ScJSWUnH4GCdddR/Kvb9c6nAFPzn44QBSva2D569tRdAqzrxpF7rjQbjTzC3uL2gC5eyBNQ/fAml6DaVoq0qn/HrJ+NZOoE2Z39xoZCnr5JfVwVd12O50rVjDkq+XoLRatwxnQZNXKAFEwMZmEzCgWv7iFz57bxPg52Rx79iCdq2V3wm7cAY271JGPDTvUVdJ3C4tSS9kjz1KTdfIofNEFNJ1zEebJuUTe/M9BX0XSVwnXXYtt8WJa336HxOuv0zqcQUkm8hAUm2xm3u8m8b/3i9iwpIKaXa2c/MsxxCQNwKoWnw+sVepEUI279p4YqrNxz37GSHV4esFJ3SMfu0dAxmT9LFG3vvQS3sYmkp56KrQmxgpSEaNHEzljBi2vvUb8FZejixiA12GQk1UrIa54XQNfvbkDn09wwqXDGXZ0iM7V4rDumWa1qVAdtt5UqE632nsRg/DY7nlFhqnD1JOGq1t05iGtMuO12SiePYfwcePIfvEF//0+g4x97VrKL7uclHvvJf7yy7QOZ8CSVSsDVMHEZJJyLCx5eRtLXt5G5fZWZlw4lLDwIPzTet3qtKq7E3ZzETQVqUm7o37PfopOXbQgcSjkzYTEId0TQg3r8/qNLf9+DW97O0m33db330fqYZ48mYhJk2h+9RXi5l+MYgjC628Ak2d7AIhOiODc30zomaulprCNWVeOJH1IbOCDEUKt8tg9J3bvUvbP5sZOVKdWHTJHTdYJQ9QGx/g8MJj6PTRvezstr72GZc5sIsaM7vf3H+wSfnENVTf9CtuXXxJ92mlahzOoyEQ+QOyeqyV7VDzLXtvOx39dx4TZ2Uw5Kw+D0Q+DVLweNTE3dc+J3VS4J3k72/fsZwjfe27sxKHqz/H5YI7v/7gOoOX1N/B1dJD4q18F9LiDRdQJJxCWk0Pzq//Gcuqpsv3hJ9rq7Xz91g5OumIk0Yn9247QL4lcUZRXgDOABiHEmP54T+nIpA+N46L7prDiwyLWL6mgfGszs68aRVL2EXYL8/nUCaB6d+Vr3KlWh3hde/azpKul6nEXqFUguxczOMS6a3/zWq20vP46ljmzCR8xQutwBiRFpyP+qiup+78/0rVuHeZJk7QOKWiUbmpi6Stb0el1dLY5+z2R90tjp6IoxwMdwOuHkshlY2dglG9pZvkb23HY3Bw1O4ujT8878DznDivUb+2eT2STOtVqw/a9Gxtjs9X5Q3oaG0eoCTvIB800PvsPmp59lryPPyJ85EitwxmwfHY7RSeehHnK0WT+/e9ah6M5n0+wZmEpaxeVkZRtYe71Y/o0kM+vjZ1CiG8VRcntj/eS+k/OmATm33+MWjr/soLCtfUcf9Ew8o5KUpN27Uao3QA169Wt91wiEXHqyMbJV6uJO3mU2jvEFHpzpHutVlpee42o2bNkEvczndlM7MUX0/zCC7jKywf1whOOTjdLXtlKxdYWRk5L4/j5w/xTzUkA68gVRbkOuA4gW+tltAaR8Egjsy4bxsgRTr75uI5F/9xMrmUbM8KfJlrfoO4UkwVpR8FRl6jziKSMCbrFDPqi5fU38NlsJN10k9ahDApxl15C8yuv0PL6G6T+4T6tw9FETWErS/+9nc42JydcOpxR09P92mbQb/3Iu0vkC2XVShDwOKF6HZSvgPLvoXI1uGx4hZ6NrgtZ034uKDrGTVSYcOYYwpNDtO/5IfBarRTNnoP56KPJ+sezWoczaNTcfQ/WxYsZ+tVy/624FITcLi8rPylm0/IqopMimHPNKFLzYvrt/WU/8oHM64bqH6HkGyj9FqrWgNepPpc0AsaeD1nHoM+awsT4fIa0OFj1nxLWraln69ZdTDzFydgTMzEOwMUrWt54A5/VSuJNNx72a4VP4HJ6cXV51M3hRfh8CAEIEIACGMP1hIUbem4NYbpB32Mj/uqraP/kE1rfe5/E667VOpyAqC1uZ9lr22hv6GLsCZlMPbcgYGvvyhJ5KBJC7ZtduARKvlZL3q4OQFGrRnJnQM40yDoWIvc/tWhTlY2V/ymhfHMz5pgwJs3NZeS0tAGz8LPXZqNo1mzMkyeT9dw/fva8o9NNW72dtno7thYHHa1OOlqddLap951dHjVbHyadXiEyxkRkrLpFxZqwJIQTkxxBbIqZ6IRwdHrte/L4W8U11+AsKmbI0iUoh7vyUghx2t2s+ayMjcsrscSFc9IVI8gc4Z+utX6d/VBRlHeAE4BEoB54QAjx8v72l4n8CLg61dJ24RIoWgJtFerj8QWQfwLkz1QT+BH0za4pbGPlJ8XUFrdjMhsYNT2dsSdkYokP7ZVtmv75Txqffob0t96jMyaL5qoOmqo7aK3ppLXejqPDvdf+EdFhRPVKvuFRRsIiDJgiDIRFGAgL16PoFRQARb0VQuB2enE5vLgdaqnd0emms91JZ5uTzjYXHW1OPM49A6F0eoWYpAji0yJJzIoiIdNCYmYUUXGmAVWS7/juOyqvvY70xx4l5uyztQ6n33m9PrZ+W8OahaU47G5GT09n2rwhfh1VLaexDUUdDbDzc9i5CIq/UqtLjGZ12PrQ2eqIyLj+6RUghKCuuJ2NyyspWd8IikLBhCRGH59B+tBYdLrQSDB2q4vGShsNRc2UvbOYjpgc7MqenjYms4H49EjiUszEpkQSm2omLsWMJT4cvdE/pWQhhFr6r7PT1qB+A2its9NS00l7456unaZIA0lZFlJyo0nOjSYlN5rI2P4f4RooQghKzjgTnclE7ocfDJgPKSEEZZub+f7DItrq7WQMj+W4eUOPfKzGYZCJPFS0lsG2/8D2hWpdNwJismHEaTBsrlpl4ofh671Zm7vY/HU12/5Xg6vLgzk6jIKJyQyZlExaQQxKkCT1znYnjeU2NXGX22iqtNHR6ux5PsJeT+rYLFLGZpKQEUViZhSRscFV6nU5PDRXd9JUqcbfUGGjuboT4VP/LyNjwkjJjyElL5rUvBiSciwh1ZbRumABdQ/+Hzlvv4154gStw+kT4ROUbmpiw5IKaovbiU0xM+28AnLHJQbsmpKJPJi1lKjJe+snar9uULsDjjgDhp+mzqOtQfJxO72Ub2mmaG09ZVua8bp9RMaEkT06gbQhMaQNiSUmKcLvF7HX66Otzk5TVQfN1R1qFUlVB3Zr98hSRZ3aNynbQnKOhYQUE9ZfXYplSDbZr7zi19j8wePy0lTVQX2ZlfpSK/VlVqzdJXdFp5CQEUlKXgwpudGk5EUTl2IOmg/Xn/J1dlJ4wolEzZhBxpN/1TqcI+Jxedmxso6Nyyppq7djiQ9nwsnZjJqRjj7AbR2y10qwsdbC1o9g8/vqYByAjEkw5yEYdRbE5WoaHoDRpGfIJLUk7nJ4KN/cTNG6Bko2NrL9+1pArVdOL4ghLj2SmKQIYpLMxCRFEGExHnKCFz5BV4cbu9VJZ7sLW1MXbfVdtDWq1RDWJkdPCVVnUIhPiyRrZDxJ2RaSsi0kZkXtVS/Z8uZbdDRUk/DEn/v/pASAIUxPan4Mqfl7uq112VzUl1qpK22nvtTKrtV1bP22GoCwCAPJORaSc6JJzrGQlGPBEh8eFN88dJGRxJ53Hi1vvUVyfQPGlNBYPFwIQVNlB7vW1LPjh1ocHW6Scyyc/MvRFExICrrGalkiDyRHu1ry3vw+lH4HCEgbD2Pmwehz1OHvIUD4BK11dmqK2qgtbqOuuB1bs4Pel5IhTNfdQGjAaNJjNOkxGHV4PT48bh9ejw+v24ery4Pd5u5J1L1fH5NsJjbZTGxyBPHpkSRkRBGbaj5gKUi4XBSdfArGjAxy3nwjKJKZP/h8grY6O/Vl7T2l9pbqTnzd5zE80khSdhQJGVEkZKq38amRfmsHOBBXRQXFp8wl8cYbSLr11gPu6/b6cLi9ONw+3F4fJoOOiDA94QZ9QNpp2urt7FpTT+Gaetrq7ej0CrljEzlqViZpQ2I1v55kiVwrPq/aULnxbdjxGXgcak+TmXep/bsTh2od4WFTdArx6ZHEp0cy5vgMQF0c2tbsoK3BjrWpC2uzA3eXB5fTi9vpxe3w4rB7MBjVBK836NAbdISF6zHHhBEZY8IcHYY5xoQl3nTEddlt//kPnro60h56SPN/On/S9fobjJyWDoDH7aW5upPGcqvaZlDVweavq/F6fID6d4vt7gK5e4tLjTzsb1CHKyw7m6iZM2lZ8C7Vp1/MzhYH1a1d1LY7qLM6qG130GB1YHd58fj2X7AMM+iwmAykxoSTFhNBemw4qTHh5MRHMio9mpx482En+64OF9U726ja0ULVjla18VmBjGGxjJ+dRcHEZMIjjX09BX4nS+T+0rgL1r8Bm95T15CMiIMx58NR8yFj4oAZ/h5MhMdD8amnoY+JIff99wZ0Ij9UPq+P9sYutX2hqoPWOjut9XbaG+z4vHv+9w1GHZaEcCwJEUQnhHf3ge/+gI0xERkThinSeMiJ0unxsqmqndWlLWyuakf34ypu+eLv/GXSfL7KmoROgWSLmohTo8NJiTYRaTIQbtQTYdQTbtRh0OtweXx0ub043F663F6sXR7q2tUPgZq2LqwOT88xzWF6hqdaGJkWzcTsOKYWJJARu2eCKkenu6d9pam6g8ZyG83VHYA6qCtjaCyZI+IpmJhMVFxw9haSJfJAcHWqDZbrXofKlaAzwNCT1eQ97BS/9zYZ7KyLFuGurCTl7rtkEu+m0+uIS40kLjWSoZNTeh73eX1Ymx09bRC2ZvVblK3ZQX1ZO85Oz8/fTIFws5HwKCMRFiPhkcaePvbGCANNTjclbXaKWjrZ0dyJ3efDhSA5Lpz8MZNo2zSEm62beeA3d5AZb8bQD/XMnU4PxfU2tpW2UVhlpaLOxuZVtWz7tpr/+hRSjUZSDQYiXOCz7/mdIixGErMsDJmUT+aIOJJzLEFX7304ZIm8P9RsgB9fhc0fgsumzsU98Qo1gUeFRuNOqBM+HyVnnIliMJD3yccoQTAHeijzuL3Y2110trvobHNitzrp6nDj6HDTZXPj6HTRZXNjtblwdXnQeQ4vj+j0Sk/1mt6goDeq0xro9AqKTkFRFJTuP6HwAQiEUNtnPO7udpbdW3fV0c8o4AxTaPJ5acVHq0GQlGXh2ImpzJ2cQUJU6BWsZIm8v7k6YcuHsPYVtdeJIUJdAWfiFZB9rKw6CTDbl0twlZSQ/tcnZBLvBwajnujEiJ8tgODx+vhfUROL1lWzpKWVLr2XqAQDM4clMntIElMyY4lQ9OpoV6cHt8OLx+XF4/bhsnVR//fn0BcMJ/LEWWqDd69N+NRE7fMJRPeGoqAoqN+wum8NYTr0Rh2G7s1o0qvfDiKNmMwGTGZjd5VQGDqdgtcn2FTVxuKt9Xy+pZZPv9jOHxZv59j8BM6flMlpY9MI99P0soEiS+SHq3EnrHkJNi4ApxWSRsLka2DchRARq3V0g5IQgtJ58xD2LvI/W4iiD+1/ymC0vdbKR+uq+GRDDY02J7FmI6eNTeOU0akcmx+PyXBo57zu4T/TumABQ5cvw5CU5Oeof04IwbZaK59vrmPhphrKmu1Ywg2cOyGDi4/OZlR6cC+QIkvkfeH1wM7PYPWLUPYd6MNg1Dlw9C8g6xhZ+tZY53ff4dy2nbSHH5ZJvB85PV4+21TL6z+Us6GyDaNe4cThyZw3MZMTRyQdcvLuLf7SS2h94w1a33uPJA3WTlUUhdHpMYxOj+E3Jw9jZUkLC9ZUsGBNJa//UM5RmTFcMz2P08em9UsdfqDIEvmBdDTCj/9Wq09sNepQ+clXq9UnkYlaRyehlrDKL70Md10tQxYvRjEGf1exYFfd1sVbK8t5d00lzZ0u8pMiueyYHM6ZkEF8ZN9nMaz45bU4CwsZsmwpiiE4ypJtdhcfr6/mjR/KKWnqJCM2gmum53HR0VlEmYIjRhggJfL31lby7a5GPF6B2+vD7RO4PT4UBYx6HWEGdTPpdUSFG4g1hxFnNhJnDiMuMoy0mHAyYiOIPNgfpmY9rHperQP3uqDgJDj9r2rPE50s8QUT+5o1dK1bR8p998kk3kdbqtt5/tsSPttUA8DskSlcMTWX44Yk9GsvoLhL5lN106+wffUV0XPm9Nv79kWsOYyrj8vjyqm5LNvRwAvfFvPQwm08vXQXl0/N4ZfT84nrhw8xfwmpRF7d2sW2GitGvQ6DXsGo12HUqxdYp9OD0+PD5fXh8viwOTy0d7n3+T5xZiOZcWYy4yIYkhzFsBQLw5IiKGhajmHN81C5CoyRMPFKmHIdJA0L5K8pHYbm519An5hI7PnztA4lJAkh+L64mX99U8x3hU1EmQz8ckY+V07L3asPdn+KmjkTQ3oabe+8EzSJfDedTmHOqBTmjEphXUUrL35bwnNfF/Pa9+VcMz2PX0zPIyYi+AoMA7pqxeP10d7lptXuptXuoqati6rWLqrbuqhu7aKyxU5zcwMX6ZZzpeFLMpRmanSprEk+H+eY+YzOz2J4iiWk6soGk67Nmym74EKSf/sbEn75S63DCSlCCJbvaODpZYVsqmonMcrENdNzufSYnIAkqqZ/PU/jU0+Rv2gRpvw8vx+vL3bW2fjbkl18sbWO6HAD1x2fz1XH5WlS5SJnP/yp5mJY9S/E+rdQ3J00Jk7hu/jz+dQxjg1VNlrtamk+wqjnqKwYphUkMq0ggaOyYjHKxB4Uqm65hc5VqxmyfBn6qKiDv0DqSeBPLS1kc3U72fFmbjyhgHMnZAS0C56nqYnCE08i/pL5pNxzT8CO2xdbqtt5aukulm5vIDEqjDvmDOfCyZkBLejJRA7qEmkVK+GHZ9V5T3QGdb6TY29Up43t2U1Q3mxnQ2UbGyrbWFPWwrZaK0Kow4Cn5MUzfUgis0emkJsYGfjfQ8JZWEjJmWeReNNNJN16i9bhBD0hBF/vauSpJbvYWNVOVnwEt5w4lHMnZmhWMKm+4zd0fPcdQ7/5Gp3ZrEkMR2J9RSuPLNrB6rIWhqdYuPf0kRw/LDBdKQdEIv/69RfZtOxL9QchEL0WVFQUXfdoMEW9r9Oh0+vVW50OndeBztGGztuFTqdDb0lEF5OO3hSBTm9AZzCgNxjQ6dVbvcGAzmDsue8WOmpsLspbXZS0Omjs9OBV9CREmxmdFc+47ASGp8cRFm7CEBaGwRim3ppMGE3hGE3h6IOkhX4gqL7zd9iWLWPIsqUY4uK0DieorS1r4S9f7GR1WUtQJPDd7GvXUn7Z5aT96SFizz9f01gOlxCCxVvr+POiHVS02Jk5LIn7Th/J0BT/rhI0IBL5rlUrqNm5DXXtcrpHfSkIIdTELgRC+NSVzn0Cn9uJr3EXomE7PqcdnzEKX1w+vqh0vD6Bz+vF5/Xg9aibz+PB63Hj83p7PebG63bj9Xqhj+dKpzcQFh5OmNmMKcKs3pojMUVGEWGJ7tnM0TFExMRgiU8kKj4evSH4Gle05KqspPiUucRfeSUpd/1O63CC1vZaK08s3smyHQ0kWUzcOmsoF03OIswQHFWDQghKzzobxWgM2aXgnB4vb/xQztPLCulyefnF9DxunTX04D3jjtCASOSHrKMRVr8Aa16ErlbInALH3aqutnOE3QeFEAifD6/HjdetJnxvd+K3dTpYU9zA97vqWF/ahNftIs6kY2JGJBPSzKRF6vA4nbidDlwOB64uO067HZe9E6fdjqOzA0eHFVdX1z6PHRGtJvXopCRiUtKISUklIiEOY0I0uhgzXT4Hdrcdu0fdutxddHm6cHgd6q3HgcPjwOl14vA6cHqcOH1OXF5Xz+b2uXF5XXiFF7fPjcfnwePz4BM+vMKLEKLndl90iq5n0yt69Do9Rp1R3fTGnvsmvQmTwUS4PhyT3kSEIQKz0YzZYO65jTRGYgmzEBUWhcWo3saExWAJs6DX6am9/wHaP/mEgqVLMCbLuWx+qrqti78u3snHG6qxmAzccEIBV03LxRwWfN8IW995h7r/+yO5771LxLhxWodzxJo7nDz2xQ7eW1tFWkw4fzhjFKeOSd3rw6mjtYWiNSsZMe14wo+wTWdwJPKWEvj+77DhbfA41cR93K3q3CcB0uXy8tXOBj7bVMuyHfU43D5yE8ycNzGTs8anEhflo8PdQYerA5vLRqe7E5vbRqerE1tXOx3WVuzt7TitNtxWG15rF3Q40Xd4COvwEdGpYPDtuTg8Oh9tUW7aLG5aLW7aLC6aol04TXsmEgrXh2MymDDp1QQablCTqFFnJEwfpm66sJ6Eq1f0GHQGDDoDekXfk5wVRUGn6FDYu+QkEHsle5/w9XwQuH3ung8Gl9e154Ok+0PF4XGoH0BuO27fvruL9pbliOSxZ9pZPyWBlZeOIy48jvjweOLD44kLjyMxIrFnizPFoR9E/f6tDjfPfVXMKytKUYCrj8vjxpkFxJiD9xudt6OTouOPx3LyyaQ/+ojW4fTZj+Ut3PfJVrbXWpkxNJG7piXSVbSRwtXfU7trBwBn3nEPw4457ojef2An8poNsOIpdfUdnQGOuhim3tKv/b+FEDi8Dmwu215bh7vjZ4/Z3Optu9NGna2NVkc7bp8dRe866HH0ip6osCiijFFEGiN7bndvZkMEZqeBcKtA1+6EZjuehjac9S242m097xOVmEhKwVAyho4gbchwUguGYQgL3gENAG6vG7vH3vNBZ3VZ1Q88tw2r00q7q53sV5aSv3Qn/75/CuWRdlqdrbQ6WnF6nT97P52iI84UR7I5mRRzCknmpJ77KZEppEamkmpOxWwMnYa2fXF5fLy1qpxnlhXS1uXm3PEZ/OaU4X7rB97f6v74R9o++JAh33w9INo7Wupqeffd/1L14w/EO5sASMotYNgx0xg6ZRoJmVlH/N4DL5ELoc578t2TUPIVmKLVyauOuQGi0/b5EqfXidVpxeayYXXtfbt7++nPu5Oy1WXF49vHHM29GHQGLEbLz6oFLGEW8IZT0eRje42bFpueSGMkM/IzOX1MHiNSkrGEWYg0RhKuP/K1FrtsVpoqyqgrKaKuuJC6ol1YG+sB0BuNpA0ZTuaoMWSNGkva0OEYTeFHdByteFpaKDppFtGnnEL6Y4/2PC6EoMvTRXNXM82OZpq6mvbaGuwNNNgbaOxqpMXR8rP3jQ6LJi0yjbSoNNIj00mP2rNlRmUSHRYdlPW3Qgi+3FbPI4u2U9Zs57ghCdxz6kjGZMQc/MVBxLFrF6VnnU3ynXeS8ItrtA7niHS2tbLzh+/YseIbagt3ApBYMIztpnw+syaQlZ3JY/PGMTazb3+bAZHIPR437RXF2NtXY/vxRaxNO7GZ47EOORFb+jisPjdWl3WfSdrqtOLyHbhEbNKbsIRZ9tqijdFEhUURHRa91+NRRjVB7348KizqkJKwzyf4oaSZt1dX8OXWOtxewZTceC6fmsPcMan93pPA3t5GbdFOKrdtoWrbFhpKixHCh95oJHPkGHKPmkje+EnEZ2QFZbLqreFvT9H8wgvkf7YQU37+Eb2Hy+uiwd5AXWcddfY69bZ7q+6oprazlk53516viTJGkWnJJCMqg8yoTLIsWWRZssi0ZJIWlYZRF/iqiy3V7Ty0cBurSlsYkhzFvaeN5IThSUH/N9yfsksvw9vURP7ni0JmGmKPy0Xxj6vY+s0yyjauQ/h8JOXkMeK4mQyfOoOYZHUhjy+21HH/f7bQ1OHkmuPyuOPkYUfcXjEgEvlnV87Esr2BO67V4wz7+QWrV/R7EnB3gu19G22K/llC7v2zSR/YieabOpx8+GMVb62qoKLFTrLFxCXHZHPJlGySo/1TWnbaO6neuY2KzRsp27iO5qoKACyJSeSNn8SQo6eSPWZc0PWU8VqtFJ00i8jjjiPz6af8dhwhBFaXlZqOGmo6aqjqqKLKVkV1RzXVHdVU2ar2KhDoFT3pUelkW7LJsmSRE51DdnQ2OdE5pEel93uSr2t38PjinXy0vop4cxi3zxnG/KOzQn70cfunn1Jz5+/IfvUVIqdO1Tqc/RJCUFu4k63fLGXnD9/h7OwkKiGRUTNOZNSME0nI3PcC6laHm8c+38Fbqyp4/vJJnDI69YiO79dErijKXOBpQA+8JIR49ED7H2kiX/WPK4j++xpaZg/F+dtbiA6P2ytJmw3mkCyR+HyCb3Y18toPZXy9sxGDTuHUsWlcc1wuE7L9W2dobWqgbOM6Stf/SPnmDbgdXYRFmMmfeDRDjj6WvAmTCQvXvq616V//ovGpp8n76EPCR43SLA6f8NFgb6DKVkWlrbJnK7eWU2mrpMPd0bOvQTGQYckg26Im9tzoXHJi1NtkczI65dCTr93l4flvSnjh2xK8PsE10/O46cQCosOD6wP3SPmcTopmnoB5yhQyn3la63B+xm5tZ/t3X7F5+Zc0V1VgCDMx9JhpjD5+FlljxqI7xEb1nXU2hqVEHXGe8lsiVxRFD+wC5gBVwBpgvhBi2/5ec8R15B4ntX98mLYPPyLv/fc0/Yf2l7KmTt5YWc57ayqxOT1MyI7lF9PzmDs61e+lLo/LRcWWjRSu/oHiH1fRZW3HEGYif+LRjJh2PLkTJmEMC/zyWD67naKTZhF+1Diyn38+4Mc/VEIIWp2tVFgrKLeW77VV2Cro8uzpXhphiCDbkk1uTC650bnkxeSRG5NLXnTeXo2vPp/go/XVPL54B/VWJ6ePS+PuuSPIig/tBtp9qf/L47S89hpDli/HmKJ9t1Lh81G+ZSObly2maM1KfF4PaUOGM+akOQyfejwmDUaj+jORTwUeFEKc0v3zPQBCiP32JepLY6e3vZ3i007HmJZG7rsLBuxCAh1ODx+sreTV78sob7aTHhPOldNyuXhKdkAmNfL5vNTs2M7Old+xa+UK7O1thEVEMGTysYycfgLZ48Yfcimkr5r//W8aHn2MnLffxjxxQkCO2d92l+TLreWUtZdRZu3e2suo6azBJ/Z0F002J5MXk4eZNDaWhFHVaGFEQgEPnjaVo3MH7jz4rrIyiueeStJtt5J4442axWFvb2PL10vZvGwxbfW1hEdZGHX8SYw9cQ6J2bmaxQX+TeTnA3OFEL/s/vly4BghxM0/2e864DqA7OzsSeXl5Ud8zPaFn1Hz29+Scu+9xF9+2ZEHHwK8PsFXOxp4+X+l/FDSTGSYnouOzubq43IDVirzeb1Ubt3Mju+/pXD1CrVeMC6eEdNPYPTMWSRm5fjv2E4nxXNOJiw3l5zXX/PbcbTk9DqptFZSai2lrL2MzQ2FrK7eSYe3GkW/p1tlhCGC3Og9Jfjdpfic6BwiDNpXf/WHimuuwVlaxpClSwJaSBNCULl1MxuXfk7R6h/weT1kjhzDuNlzGTplWtB03dU8kffW1+6HQggqr72OrvXryf9sIcbUI2s4CDVbqtt56bsSFm6qxScEp45N49oZ+YzPig1YDB63m5J1q9n27XJK16/F5/WSnFfAmBNmM3L6iUc8Ym1/do/8y37lZSKnTevX9w42LZ0unllWyJsrywk36rlhZj7nTLZQa6+ktL1U3bqTfU1HzV5zDaVHpv+smiY3OpcUc0pItRtZF39J9W23kfnP57CceKLfj+fo6GDrN8vYuPRzWmuqCI+MYtTMWYybNbdP/b39ZcBUrezmqqyk5IwziTr+eDL//kyf3ivU1LZ38e8VZby9ugKbw8OU3HiuOz6fk0Yko9MF7p/Wbm1nx4pv2Pr1MhrKitEbjQydMo2xJ51M1qixfe5GJlwuiubOxZiUTM6Cd0IqIR0Oh9vLKytK+efXxdhdXi4+OovbZw8jybL/9giHx6FW01jLepL87qoau8fes1+EIWJPQ2t0Tk+Cz47OJjos+BYaFm632h4yahRZz//Lb8epK9rFhiWL2Pn9d3hcTtKGDueoOacxbOp0TdqBDpU/E7kBtbFzFlCN2th5iRBi6/5e018jO5teeJHGJ58k87l/YDnppD6/X6jpcHp4d00lr/yvlOq2LvKTIrl2Rn7A55YGaCgrYfPyL9n+v69wdnYSk5zC2JNOYcyJc4iMPbKeN20ffkjtvfeR9fy/iJo5s58j1p7XJ/jgx0r+tqSQOquD2SOTuWvuiD7NoCeEoLGrUU3sB6mLjw+PV7tLdveqyYrOIseidp+MNGo3PXPjM8/Q9M9/MWTpEowZGf32vm6Hgx3ff8vGJYuoLynCaApn5IwTOGrOaSTnHtm4hEDzd/fD04CnULsfviKEePhA+/dXIhduN6XnzcPb0UHBwk/RRQ7OucE9Xh+fba7lxe9K2FJtJSEyjMun5nD5sTkkRAW2dOF2OSla/QObl39J5dZN6PR6CiYdw7hZp5AzbsIhl9KFx0Pxaaejt1jI/eD9AVUaF0KwbHsDf1m8g131HYzPiuX3p41kSl68X4/r8rqoslVRZi3r6U1TZi2jwlpBY1fjXvvGh8f39I3Piu4eABWVSaYlk4Tw/l3D86fcNTUUzZ5DwnXXknz77X1+v+aqCjYu/Zxt3yzHae8kITOb8SefzsgZJ2rS86QvBsSAoH2xr1tP+SWXEH/VVaTcfVe/vGeoEkIdNfrSd6Us39GAyaDjvImZ/HJGHgVJgV9Bp6Wmms3LF7P166V02axEJ6UwbtahldLb//tfan53F5nP/h3L7NkBitj/VhQ18fjinWyobCMvMZLfnTKcuT+ZJU8LdredSlslFTa162SVrYoKWwWVtkrqO+v3qo+PMESQHplOpiWT9Kh0MqIy9kxrEJlOrCm2z79P5Y030bV5M0OXL0M5goZGj9tN0erv2bj0c6q2bUGnNzDs2OM4as6pZIwYrfn5PlIDNpED1N7/AG0ffkjeB+8TPnJkv71vKCtqsPHSd6V8tL4al8fHSSOS+cX0PKYV+Lc0tS8et5vitSvZuOTznlL6kKOnMm72XLJHj/tZKV14vZSceRaKwUDeJx+HzJDtA/mxvIUnFu/ih5Jm0mPCuXXWUOZNytR8cYdD4fQ6qbZV7zXKtcpWRVVHFTUdNXsNggJ1ts3UyNSe+WtSzanqZGWRKT2Tlh1s/pqOb76h8vobyHjqKaLnnnLIsbbV1bJp+WK2fLWELms7MSmpjJs1lzEnzMYcE3ukpyBoDOhE3tO3PCOD3HfeHrB9y49Eo83JmyvLeXNlOc2dLkakWrjmuDzOGp8e8Hp0gJaaKjYt/YKt3yzD0WEjLi2dsbPmMnrmLMzR6oRC1s8/p/rXd5DxtyeJPvXUgMfYn34sb+GZZUV8s6uRxCgTN59YwMVTsjU59/6ye0qDaps6V83ura6zjtrOWpq7mvcq0YM6r1FiRCJJEUkkmZNIikgiISKBhPAEEiISiDfGEnbxbZhycsn9978PmPS9HjdFa1axadkXVGzegKLTkT9xCuPnnHpY1XmhYEAncoD2TxdSc+edpNz/B+IvuaRf33sgcLi9/HdjDa/8r5QddTYSIsO4eEoWlx6TQ7oG0516XC52rVrBpqWfU71jG3qDgaHHHMe4Wafguuc+cHvIX/hpSH4o767i+vuyIn4oaSY+MoxrZ+Rz5bScoFzcwd/cPjdN9ibq7fXU2euo76zvmZWy9+1PS/bnrvAx/1sfv7kxHGd6PHGmOGJNsXvmTOrQE7atBc+mSnydDsJio0mfNomCGTNISEwjwhhBhCGCcH24pvPS+4RPnX+/e4GX2PDYI+73P+ATuRCCyl/8gq5Nm8n/7LOgGOIbjIQQ/FDczKvfl7Fsez2KojBnZApXTMthan7gq10AmirK2LRsMdu+VRujIh0uxkydwcSbbu0ppYcCn0/w9a4GnvuqmLXlrSRZTFx/fD6XHJM9KBP44XJ4HLQ6Wml2NNPiaKGtupSCa/5C2dwxrDpvOK3OVqy2VkzF7SQWeUho1uNTBFXJXezMslGT5EDs5/I16U09W5g+rOf+7sVTehZS0enRoetZQEVRFBQUfPhAgA8fQgh18RThwevz4vF58AqvutqWr3vFLa8bp8+J06MuoNLbv2b/i+My5MIS++UqL6fkzLOImnUSmX/7W7+//0BT2WLnrVUVLFhTQZvdzZDkKOZPyea8CRnERQZ+JJvL0cWKS+dTiotWgw69wcCQKdMYN+uUfumX7i8Ot5eP1lXz8v9KKG7sJC0mnBtPKODCyVkDqgpFC1W33kbn6tWYX/oX21d8w84fvsPV1UVcWjpjTjyZgmnH4Ys00OHuoNPd2bMoSe8lD3dvTq+6MpXL6+q57xGentWsPD41Me9O1gI1YQtET3IHUBSlJ/nvXtbQoBgw6o2E6cJ6Vt0y6ozqNwJDeM/KXBGGCKalTyM1MghnPzxc/lyzs+mf/6Tx6WfIeuF5oo4/3i/HGGh2V7u8vaqCDZVthBl0nDYmlflTspmSFx+wUrpt+XKqbvoVaQ8/jOfoiWxavpjt336Fo7ODmJRUxsyczaiZs4hOTApIPAfTYHXw9uoK3vhBbX8YnR7NtTPyOX1cWkg0YgY7a2MD6998lW3fLMNuCsNoCmfoMeqAs1DuedIXgyaR+1wuSs85F+F0kr/wU3QRA2MOikDZVmNlwZoKPl5fjc3hITfBzLkTMjl3QgbZCf7rcyuEoHTePHwdnRQs+gzFoFZFuF1OilZ9z5avl1CxZRMoCjljxzPmxDkMmXxswOfA8PkE3xc389aqcpZsq8fjE8wakcwvZ+RzbH7gPvQGqq4OG4UrV7B9xddUbdsCQKLbR254FFNfeS0oplTW0qBJ5AD2NWsov/wKEn75C5J/+1u/HWcg63J5WbS5lg/XVfFDSTNCwNG5cZw7IZPTxqYSa+7fBGpbupSqm28h7ZFHiD33nH3u095Qx5avl7H166XYmhsxmSMZduxxjJxxIpkjRvu16qXe6uCT9dW8s7qCsmY7cWYjF0zOYv6UbPISB+dAtP7idjooWbeW7f/7unv+Hg9x6ZmMnD6TUTNOwvPfT2l44q/kL/wU05AhWoerqUGVyAFq7ruP9o8/URciGD7cr8ca6KrbuvhkfTUfr6+mqKEDg05hakECp41N4+RRKX0ePSp8PkrPm4fo6iL/s4U9pfH98fm8VG7ZzLbvllO46nvcTgfRScmMnH4Cw6cdT2JWTr+UjK0ON19sruOTDdU9H2ZTcuO55Jhs5o5JlfXffeB2Oijd8CM7f/gfJetW43E6iYyLZ8S04xk5/QSS8wp6/oae5mYKTziRuPkXk/r732scubYGXSL3trWpfcuzMsl9552gbSgLJUIItlRb+WxzLZ9vqaW82Y5OgWPyEpg1MpkThidTkBR52EnU+uWXVN96G+mPPUrM2Wcf1mvdDgdFa1ey7buvKN+4HiF8xGdkMXzqdIZPPf6wZ7Br7nCyfEcDS7fX89XORlweHzkJZs4en8HZ49M1GSE7UDjtnZSuX0vhmpWUrluD2+kgIjqGYcdMY9ix08kcNWa/c9xX3/EbOv73P4Z++w268NBaNLw/DbpEDnuGecu+5f1PCMG2Wiufb67j8y21FDeqCxZnxkUwc1gSM4clMSUv/qBVMMLnU9s0XC613/hBSuMHYm9vY9eq79n1w3dUbt8CQpCYlcOQKVMZMvnYvUp5u/l8gu11Vr7e2ciy7fWsr2xDCEiJNnHqmDTOmZDBUZkxsu77CNlamij5cTVFa1ZSsWUTPq+HiOgYhk6ZyvCpM8gcOQbdIYwV6Fy1moorrzxg1dtgMCgTuexbHjiVLXa+2dXI1zsb+b64CbvLC8CwlCgm58ZzdG4ck7LjyYqP2CspWr9YTPXtt5P++F+IOfPMfouno7WFXStXULhqBdU7tiGED0tiEvmTjkGXM4ZCXRJrKqysLm3B6vAAMCYjmlkjUpg9MoUxGQceQi7tm8/rpaZwB6Xr11K64Ucay0oAiE1NY8jR6gdq2rDhh726lBCCktNORx8dTe67C/wRekgYlIkcevUtP/FEv66+Lu3h8vhYV9HK2rIW1pS1sq68FZtTTZYWk4HhqRZGpkUzIiWSo/7vZoz4KPj0UwzG/hs04/H6qLc5KWvqZGtJDRUb1uIp3UxcawkG4cWlGGiNycGcP5qRR09h5qQRpMYM3q/sR0oIQVt9LRWbN1KxZSPlm9fj7OxE0enIGDGKvPGTyZ94NAmZ2X3+YGx57TXqH3mUvE8+JnzEiH76DULLoE3ksGcF9kCtOiLtzesT7KyzsaGyje21VnbUWdlRa2N8yVp+v+ZNHp18Kd/nTCQ9NoLMuAiSokzERBiJiTAS3X0bZtjTxrE7IdidHmwODzaHG6vDQ3uXm+q2Lqpbu6izOvD69lzbyRYTI9KiGZloIs9dS1RjEXXbNtBeXwdAbEoaWaPHkjVqLFmjxxEVnxDYkxQihBBYG+up3rGNym2bqdiyEWtjAwBRCYnkjB1P/oTJZI8dT3hk/7YneNvaKJx5AjHnnkPagw/263uHikGdyIXLRem8eXhtHeQvXIg+SnYX05rP7WbXGWfh9sGmP/6DqnYXVa1dVLbYae500m53Y3N6OJTLU1EgymQgJsJIekwEGXERZMSqt1lxZkakWUjcR88aIQRtdTWUbviRii2bqNq+GWenWtcfl5ZO+vBRpA8bQdrQESRkZgVsselg4nG7aSwvobZwF9U7t1GzYysdrS0AhEdGkTV6HNljjiJ77Hji0tL9Xh1Vc/c92L78kiHffjso/48HdSIHsK9fT/kllxJ36aWk3ndvQI8t/dzu1X8y/v4M0XPm7HMfr0/Q0V3SdvvUlW3Uy1W9Zs1hBizhBiLDDP2yxJ3P56WxvIzKrZuo3LqJmsKdOGxWAMIiIkgtGEZK/hCScvNJzsknLj19QCV3t8tJS1UljeWl1BUXUldcSGN5KT6vWi0WlZBI5ojRZIwYTcaIUSRmZge8N1jXxo2UXXQxqQ8+QNzFFwf02MFg0CdygLqH/kTr22+T+87bRIwfH/DjSyqf00nx3FMxJCWR++6CoG1U3F1iry3cSU3hTmoLd9BcWY7XoyY2Q5iJxOwcEjKyiM/IIj49k/iMTGKSU9H3ofeNvzntdtrqamipraa1pprmynIaK8tpq61BdC8FFxYRQUr+UFILhpI6ZBipBUOJTtS+s4AQgtLz5oEQ5H38UdBeO/4iEzng7eik5Iwz0Fss5H34wRGtPCL13e5Gq+xXXyFy6lStwzksXo+bluoqGspKaCgrobG8lNaaqp7qBgBFpyMqPoGYpBSiE5OITk7BkpCIOSaOyJhYzDGxmGNj+32RXyEEzs5OumztdNmsdLa2YmtuxNrchK25CVtzI+31ddjb2/a8SFGITU4lMTuXxOxckrJzSMzOJS41PWjHXrQueJe6Bx8kd8E7g65AJhN5N9tXX1F1400k3XYriTfeqEkMg5m3o5PiOXMwjRhOzquvah1Ov3Ha7bTWVNFSU0VrbTXWxgbaGxuwNjbQ0dLcU9LtzRBmIiwiApM5EpPZTFiEGUNYGHqDEZ3BgN5gQKc3IIQP4euekc/nw+fx4HY6cDkcuJ0O3A4HTnsnjg4bPq/358cxhmFJTMSSkEh0UgpxaRnEpaUTl5pOTGpaUK8avy/ejk6KZs7EMns26Y89qnU4AbW/RB683//8xHLiiVhOnUvTc//EcspcTPl5Woc0qLT8+994W1tJ/vWvtQ6lX5nMZrUKYsiwnz3n9XjobGvF3t6Gvb2NzvZW7G1tODo7cNo7cdntOLvsajLu7MDn8eDt3nweN+h06HQ6dW5snQ6dTo8xPIKw8HCi4uIxmsIJi4ggIjoGc3QMEZZo9X5MLJaERCIsA6tPvD4qkuizzqT9w49Iuedu9LGxWoekuUFXIgfwNDZSfPoZhA8bRvbrrwXtV8iBxtPaSvHsOUROm0bm35/ROhwphDl27qT07HNIvusuEq6+SutwAmZ/JfJBmcEMSUmk3PU77GvX0vbee1qHM2g0P/8Cvq4ukm6/TetQpBAXPnw4ERMm0LZgAVoURoPNoEzkADHnnUfktKk0/OVx3NXVWocz4Llra2l9+21izj4bU0GB1uFIA0DcxRfhKi/HvmqV1qFobtAmckVRSP3jQwig9v4H5Ke6nzU+/QwIQdLNv9I6FGmAsMydiz4mhtZ3Bu/cK7sN2kQOEJaZQfJvf0PnihW0f/SR1uEMWF1bttL+ySfEXXE5xowMrcORBgidyUTMvHnYli3DXd+gdTia6lMiVxTlAkVRtiqK4lMU5WcV8KEg7uKLMR99NPWPPoa7vl7rcAYcIQQNjz2GPi6OxBtu0DocaYCJu/gi8Hhoe/99rUPRVF9L5FuA84Bv+yEWTSg6HWl/egjhdlMnq1j6nW3pUuxr1pB06y3oLRatw5EGmLDsbCJnzKDt3XcRbrfW4WimT4lcCLFdCLGzv4LRSlhODsm/vp2Ob77B+umnWoczYAiXi4bHnyBsSAGxF1ygdTjSABV3yXw8jY3Yli3XOhTNDOo68t7iLruMiAkTqHv4z4O+vq2/tLz9Nu6KClLuuqtPK/9I0oFEHX88xvR0Wt95R+tQNHPQRK4oylJFUbbsYzusxRUVRblOUZS1iqKsbWxsPPKI/UTR60n788MIp5PaP9wnq1j6yNPaStNz/yRy+nSiZszQOhxpAFP0emLnX4x91SqcRUVah6OJgyZyIcRsIcSYfWz/OZwDCSFeEEJMFkJMTkpKOvKI/ciUl0fyb39L57ff0fbuu1qHE9Ka/vEcvo4Okn93p9ahSINA7Lx5KEbjoO2KKKtWfiLukvlETptG/WN/wVVWpnU4IclZUkLrggXEXngB4cN+PveIJPU3Q3w8llPn0v7JJ3g7OrUOJ+D62v3wXEVRqoCpwGeKoizun7C0o+h0pD3yZxSjkZq770F0zz0tHRohBHX/90d0EREk3XKL1uFIg0j8JZfg6+zEunDwdVjoa6+Vj4UQmUIIkxAiRQhxSn8FpiVjSgqp999P14YNNL/0stbhhBTrwoXYV60i+Td3YEiQ615KgRN+1FGYRo2k9a23B10bl6xa2Y/o00/DcupcGp99Fse2bVqHExK87e3UP/oY4ePGEXvhhVqHIw0yiqIQN38+zsJCun78UetwAkom8v1QFIW0Bx7AEBdH9e9+h6+rS+uQgl7DU0/hbW0l7cEH5NTAkiZizjgDXXQ0LW+9pXUoASX/2w5AHxtL2qOP4Coqpv7Pj2gdTlDr2rSJtgXvEnfZpYSPGqV1ONIgpYuIIHbePGxfLsFdV6d1OAEjE/lBRB13HAnXXkvb++9jXbRI63CCkvB4qH3wQQxJSSTdeqvW4UiDXNyll4DPR+uCwdMVUSbyQ5B06y1EjB9P7R/ux1VRoXU4Qaf17XdwbttOyu/vQR8VpXU40iAXlplJ1Ekn0fbue/icTq3DCQiZyA+BYjSS8dcnQK+n+o7fIFwurUMKGu7aWhqffprI6dOxnDIgOi1JA0D85ZfhbW3F+tng+BYtE/khMmZkkPanh3Bs2ULDX5/UOpygIHw+an7/e4QQpD5w/4Ba4FcKbeZjjsE0dAgtb74xKLoiykR+GKJPPpm4Sy6h5bXXsH31ldbhaK717Xew/7CSlLvuIiwrS+twJKmHoijEXXY5zm3bB0VXRJnID1PyXb/DNHIkNXfdjau8XOtwNOMsLaXhiSeIPH4GsRfKKWql4BNz1pnoYmJoeeNNrUPxO5nID5POZCLzmadRFIWqm2/G1zn45nUQHg81d9+NYjKR9tCfZJWKFJR0ERHEXXA+tqVLcdfUaB2OX8lEfgTCsrLI+NuTOItL1PlYBkEdXG/NL72MY+Mm0h64H2NKstbhSNJ+xc2fD0IM+FkRZSI/QpHTppF8553Yliyh+fnntQ4nYBzbttH47LNEn3Yq0aedpnU4knRAxowMLLNm0fbee/gcDq3D8RuZyPsg/qoriT7zTBqffgbb119rHY7f+ex2au66C0NcHCl/+IPW4UjSIYm7/DK87e20D+BlHGUi7wNFUUh76I+EjxxJzW/vxFlSqnVIfiOEoPb+B3AWFZP2yCMY4uK0DkmSDon56KMxjRxJy79fQ/h8WofjFzKR95EuPJzMZ/+OEhZG5fXX4wnCZez6Q+vrr2NduJCk224javpxWocjSYdMURQSrrkaV3Exnd99p3U4fiETeT8wpqeT9a9/4mlqouK66/F2dGgdUr/qXL2a+r88jmXObBKuv07rcCTpsEXPnYshNZXmV17VOhS/kIm8n0SMG0fmM0/jLCyk6uZb8A2QYfzuujqqf30HYTk5pD3yiOxqKIUkxWgk/oorsK9aRdeWrVqH0+9kIu9HUTNmkP7wn7CvXEnNXXeFfH2cz+Wi6tbbEA4Hmc/+XU6IJYW02AsvQBcVRcurA69ULhN5P4s5+2yS7/wtts+/oP7Pj4RsH3MhBHUP/h+OTZtIe/QRTPn5WockSX2ij4oi9oILsH7xxYAbICQTuR/EX3MN8VdeSeubb9L4t6dCLpkLIWj4y+O0f/QRiTfdRPScOVqHJEn9Iv6Ky0FRaHn9Da1D6VcykfuBoigk3/U7Yi+4gOYXXlBL5iFUzdL0j+doefVV4i67jMRbbtY6HEnqN8a0NKJPPZW2997Da7VqHU6/kYncTxSdjtQ//p9aMn/jDWr/8AeE16t1WAfV/Oq/aXr2WWLOPZeU398jGzelASfh6qvw2e20vf++1qH0G5nI/UhRFJLvvovEm26i/cOPqLnzToTbrXVY+9X67ns0PPYYlrlzSfvTQ3IBZWlACh81CvPUY2l5/Y0Bs0iM/E/1M0VRSLr1FpLvvBPros+puuVWfF1dWof1M20ffkTdgw8SNXMmGX95DEWv1zokSfKbhGuuwVNfT/vCz7QOpV/IRB4gCb+4htQHH6Djm28ou/AinKXBMZxfeL00/PWv1N57L5FTjyXj6adQwsK0DkuS/Cpy+nRMI0fS/MILIVHleTAykQdQ3MUXk/XCC3gaGyk7/wKsX36paTzejk6qbr6F5hdfIvbii8h6/nl04eGaxiRJgaAoCok33ICrrAzb4sVah9NnfUrkiqI8rijKDkVRNimK8rGiKLH9FNeAFTVjOnkffUhYQQHVt95G/WN/0aTe3FVVRfn8+XR8+y0p9/+BtAcfRDEaAx6HJGnFMmc2YQUFNP3zXyHVq2xf+loiXwKMEUKMA3YB9/Q9pIHPmJ5OzptvqOt/vvoq5VdciWPXroAcWwiBbdkyys6/AHd9PdkvvkD8JZcE5NiSFEwUnY7EG67HWVhIR4ivwdunRC6E+FII4en+cSWQ2feQBgddWBip9/+B9Mcfx1lSQuk551L7f/+Hp7XVb8d0lpZSed31VP3qZgxJSeS99y6R06b57XiSFOyiTz0VY3a2WioPsYF7vfVnHfk1wOf7e1JRlOsURVmrKMraxgE61euRiDnzDAq++Jy4+fNpe+99ik+ZS8vrr/drdYuvs5OGv/6VkrPOpmv9elLuuVut3snN7bdjSFIoUgwGEq79JY4tW+j83wqtwzliysE+hRRFWQqk7uOpe4UQ/+ne515gMnCeOISPtcmTJ4u1a9ceQbgDm7OwkPpHH6NzxQqM6elEn30WMWeeecTznDhLSrAu/Iy2Dz7A09BAzLnnkvybOzAkJvZz5JIUuoTLRdHJp6hVnm+9GdSD4BRF+VEIMflnj/f164SiKFcB1wOzhBD2Q3mNTOT7J4Sg45tvaH39DTpXrgSfj/AxY4g560wijzsOY2YmOpNp36/1enHX1mJb/CXtny3EuW07KArmY48h6dZbMU+YEODfRpJCQ8ubb1H/pz+R/fprRE6ZonU4++WXRK4oylzgSWCmEOKQ60tkIj807oYGrIsWYf3vpzi2bet53JCcjDEzE2NmBsLpwlNfj7u+Xl2dyKM2WYSPG0fM6adhmXuqXOlekg7C53BQNHsOpqFDyAniaW79lciLABPQ3P3QSiHEDQd7nUzkh89ZUoJjyxZcVVW4q6pxV1birq5GMZkwpKZgTE7BkJqKMTWFyGnTCMvJ0TpkSQopzS+/QsPjj5Pz5huYJ/8sVwYFv1WtHAmZyCVJCja+ri6KTz4FY1ZW0NaV7y+Ry5GdkiRJgC4igsRf3UTXunV0fPON1uEcFpnIJUmSusXOm4cxO5vGJ/8WUqM9ZSKXJEnqphiNJN12K85du7B+FjozI8pELkmS1Ev0qadiGjmSxqefCZn5ymUilyRJ6kXR6Ui+49e4q6po/eADrcM5JDKRS5Ik/UTk9OmYjz6apuf+ic9+SOMcNSUTuSRJ0k8oikLSHb/G29REy+tvaB3OQclELkmStA/mCROIOukkml96CU9Tk9bhHJBM5JIkSfuRfOdv8TmdNDzxV61DOSCZyCVJkvbDlJdHwtVX0/7JJ9h//FHrcPZLJnJJkqQDSLzhegxpadT98SGEx3PwF2hAJnJJkqQD0JnNpNx9N86dO2l9Z4HW4eyTTOSSJEkHYTl5DpHHHUfj008HZcOnTOSSJEkHoSgKKffdG7QNnzKRS5IkHYJgbviUiVySJOkQ9TR8PvggPqdT63B6yEQuSZJ0iHRmM2l//D+chUU0Pvmk1uH0kIlckiTpMETNmEHcpZfS8trrdPxvhdbhADKRS5IkHbbkO39LWEEBtffcg6e1VetwZCKXJEk6XLrwcDKeeBxPWxt19z+AFmsf7xWPpkeXJEkKUeEjR5J8+23Yliyh/aOPNY1FJnJJkqQjFH/11ZiPOYb6hx/GVVGhWRwykUuSJB0hRacj/dFHwGCg6uZb8NpsmsQhE7kkSVIfGNPSyHzqbzhLSqi69VZN1vmUiVySJKmPIqdNI/3hP2H/YSU1994X8MZPQ19erCjKQ8DZgA9oAK4SQtT0R2CSJEmhJObss3HX1tH41FMY09JIvuPXATt2X0vkjwshxgkhxgMLgfv7HpIkSVJoSrj+OmIvuojmF16g9Z13AnbcPpXIhRDWXj9GAtp2ppQkSdKQoiik/uE+PPX11D30J3RmMzFnn+334/a5jlxRlIcVRakELkWWyCVJGuQUg4GMJ/+KedIkau66m4YnnkB4vX495kETuaIoSxVF2bKP7WwAIcS9Qogs4C3g5gO8z3WKoqxVFGVtY2Nj//0GkiRJQUZnNpP98kvEXnwRzS+9TOVNN/m1a6LSX62riqJkA4uEEGMOtu/kyZPF2rVr++W4kiRJwaz1nXeoe/jPhGVnk/XcPwjLzT3i91IU5UchxOSfPt6nqhVFUYb2+vFsYEdf3k+SJGmgiZs/n+xXXsbb0kLphRdh90Mhtq915I92V7NsAk4GbuuHmCRJkgaUyClTyP3gAyLGjMGYkdHv79/XXivz+isQSZKkgSwsM4PsV172y3vLkZ2SJEkhTiZySZKkECcTuSRJUoiTiVySJCnEyUQuSZIU4mQilyRJCnEykUuSJIU4mcglSZJCXL/NtXJYB1WURqA84Ac+NIlAk9ZBHESwxyjj65tgjw+CP8aBGl+OECLppw9qksiDmaIoa/c1KU0wCfYYZXx9E+zxQfDHONjik1UrkiRJIU4mckmSpBAnE/nPvaB1AIcg2GOU8fVNsMcHwR/joIpP1pFLkiSFOFkilyRJCnEykUuSJIW4QZXIFUWZqyjKTkVRihRFuXsfz9+hKMo2RVE2KYqyTFGUnF7PeRVF2dC9/Vej+K5SFKWxVxy/7PXclYqiFHZvV2oU3996xbZLUZS2Xs8F4vy9oihKg6IoW/bzvKIoyjPd8W9SFGVir+cCcf4OFt+l3XFtVhTle0VRjur1XFn34xsURfHbgreHEOMJiqK09/pb3t/ruQNeHwGK785esW3pvu7iu5/z+zlUFCVLUZSvuvPIVkVRfrZqml+uQyHEoNgAPVAM5ANhwEZg1E/2OREwd9+/EXi313MdQRDfVcCz+3htPFDSfRvXfT8u0PH9ZP9bgFcCdf66j3E8MBHYsp/nTwM+BxTgWGBVoM7fIcY3bfdxgVN3x9f9cxmQGATn8ARgYV+vD3/F95N9zwSWB/IcAmnAxO77FmDXPv6P+/06HEwl8ilAkRCiRAjhAhagLhjdQwjxlRDC3v3jSiAzmOI7gFOAJUKIFiFEK7AEmKtxfPOBd/o5hgMSQnwLtBxgl7OB14VqJRCrKEoagTl/B41PCPF99/Eh8Nff7hgOdg73py/X7yE7zPi0uAZrhRDruu/bgO3ATxfp7PfrcDAl8gygstfPVfz8BPf2C9RPzd3CFUVZqyjKSkVRztEwvnndX8c+UBQl6zBfG4j46K6SygOW93rY3+fvUOzvdwjE+TtcP73+BPCloig/KopynUYx7TZVUZSNiqJ8rijK6O7HguocKopiRk2CH/Z6OKDnUFGUXGACsOonT/X7ddinxZcHKkVRLgMmAzN7PZwjhKhWFCUfWK4oymYhRHGAQ/sUeEcI4VQU5XrgNeCkAMdwKC4GPhBCeHs9FgznLyQoinIiaiKf3uvh6d3nLxlYoijKju7SaaCtQ/1bdiiKchrwCTBUgzgO5kxghRCid+k9YOdQUZQo1A+R24UQVn8co7fBVCKvBrJ6/ZzZ/dheFEWZDdwLnCWEcO5+XAhR3X1bAnyN+kkb0PiEEM29YnoJmHSorw1EfL1czE++0gbg/B2K/f0OgTh/h0RRlHGof9uzhRDNux/vdf4agI9RqzICTghhFUJ0dN9fBBgVRUkkiM5htwNdg349h4qiGFGT+FtCiI/2sUv/X4f+rPgPpg3120cJ6lf+3Y0xo3+yzwTUBpuhP3k8DjB1308ECunnhpxDjC+t1/1zgZViTyNJaXeccd334wMdX/d+I1AblZRAnr9ex8pl/w11p7N3I9PqQJ2/Q4wvGygCpv3k8UjA0uv+98Bcf8R3CDGm7v7boibCiu7zeUjXh7/j634+BrUePTLQ57D7XLwOPHWAffr9OvTLhRCsG2pr8S7UZH1v92N/RC19AywF6oEN3dt/ux+fBmzuvjg3A7/QKL5HgK3dcXwFjOj12mu6k0ARcLUW8XX//CDw6E9eF6jz9w5QC7hR6xd/AdwA3ND9vAL8ozv+zcDkAJ+/g8X3EtDa6/pb2/14fve529j997/Xj/8jB4vx5l7X4Ep6fejs6/oIdHzd+1wFLPjJ6wJyDlGrwwSwqdff8TR/X4dyiL4kSVKIG0x15JIkSQOSTOSSJEkhTiZySZKkECcTuSRJUoiTiVySJCnEyUQuSZIU4mQilyRJCnH/D6GY5UdruoU9AAAAAElFTkSuQmCC",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(ts, xs)"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [],
"source": [
"from scipy.integrate import solve_ivp"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [],
"source": [
"result = solve_ivp(eval_rhs, (t0, tf), x0, args=(p_vals,), t_eval=ts)"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x1e336ec4910>,\n",
" <matplotlib.lines.Line2D at 0x1e336ec5f30>,\n",
" <matplotlib.lines.Line2D at 0x1e336ec7c10>,\n",
" <matplotlib.lines.Line2D at 0x1e336ec50c0>,\n",
" <matplotlib.lines.Line2D at 0x1e336ec5870>,\n",
" <matplotlib.lines.Line2D at 0x1e336ec4eb0>]"
]
},
"execution_count": 49,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAAsTAAALEwEAmpwYAABfb0lEQVR4nO2dZ3gbVdqG71GzLRdZ7r2n9wIJSSCVVCD03mtgacsusJRvYQssS4el9957T0iHhJDeE/febbnIktXn+zGOk0BCmq2R7XPnmmuU0Wjm0Xj06Og973mPJMsyAoFAIOi5aNQWIBAIBIJjQxi5QCAQ9HCEkQsEAkEPRxi5QCAQ9HCEkQsEAkEPR6fGSWNiYuSMjAw1Ti0QCAQ9lg0bNjTIshz72+2qGHlGRgbr169X49QCgUDQY5EkqfRA20VoRSAQCHo4wsgFAoGghyOMXCAQCHo4wsgFAoGghyOMXCAQCHo4wsgFAoGghyOMXCAQCHo4wsgFAoGgu/F5ofQX+PHvYK3p8sOrMiBIIBAIej2OFihcCrk/QP4iaLeARgdpE2DA7C49lTBygUAg6ApkGWq3Q8FiyF8M5WvA54GQKOg3UzHv7GkQbOryUwsjFwgEgqPF1gBFy6FwmWLgbR1hk4RhMOEm6DcLUo8HjbZbZfQoI9+1uoqawhayRseRMsCMVidC/AKBwI+47EpLe49512xVtgdHQtYU6Hcy5MyA8AS/yupRRm5rdpK/oY6dq6oJMurIGBZD1qhY0oZEodN37zeeQCDog3hcULkBilcqS8Va8LpAo4fUcTDtXiVckjiy21vdf4SkxuTLY8eOlY+2+qHH7aViVxOFm+oo3tKA0+7BEKwla1Qs/Y9PIHmAGY1G6mLFAoGgT+BxKsZdsgpKfoLyteBpByRIHAGZJ0HmZEgbD0FhfpcnSdIGWZbH/nZ7j2qRA+j0WjKGx5AxPAav10dlbhP56+so2ljH7l9qMEYY6Dc2ngHjE4hNC1dbrkAgCGRcNqhYB6WrlaViHXgcgATxQ2HM5ZAxETImQYhZbbUHpce1yA+Gx+WldHsjeWtrKdnegM8jE5sWzuCJifQ7Lp4go75LzycQCHogdguUrYGy1Uped/VmJbNE0igdlOkdpp12Ahij1Fb7Ow7WIu81Rr4vDpub/HW17Pi5isaKNrR6DTmj4xg6OZmErK5P/REIBAFKc1mHcf+iGHf9LmW71gBJoyF9gmLeqcdDcIS6Wg+DPmXke5BlmfoyKztXVZO/tgaXw0tcRgQjpqWQPTpOZL0IBL0JnxfqdimmXbZGWVorlOcM4YpZp5+gDMhJHgP6YHX1HgXdZuSSJKUCbwHxgAy8JMvyU3/0Gn8Z+b64HB5y19SwdVkFzbV2Qk0Ghk5OYejkZIJDRdhFIOhxuB1Kx+Qe4y5fC84W5bnwRCU8knaC0jEZP0TVrJKuojuNPBFIlGV5oyRJ4cAG4HRZlnce7DVqGPkeZJ9M2U4LW5aWU77TgiFYy7ApKYyYkUpImEEVTQKB4DBob1LMunS1Yt5Vm5RUQIDYgfsbd2QaSL0ve63bslZkWa4GqjseWyVJ2gUkAwc1cjWRNBLpQ6NJHxpNQ0UbG74vYcPCUrYsq2DYScmMPDkNY4QwdIFAdVqr93ZKlv0CtTsAWcnhThoF4xYoMe7UcQHZMelPujRGLklSBrASGCrLcutvnrsWuBYgLS1tTGnpASeDVgVLlY3135dQsL4WrV7DyJPTGHVyGobgHpedKRD0XJpKO9IAf1bWliJluz60I749QWlxJ48Bg1FdrSrR7Z2dkiSFASuAB2RZ/uyP9lUztPJHNNfa+fXrIgrW12GMMHD8qZkMmpCIRis6RQXdi7e1FceOHaDRoouJRhcdjcZkQuqF4QFAKTDVVAIlP0PpKmXdUq48FxypZJKkn6CYd8II0IpGFXSzkUuSpAe+ARbKsvz4ofYPVCPfQ01RC6s+KaCmqIWopFAmnpVD2pBotWUJehHuykraVq2ifcsW2jdvwVVY+Pud9Hp0sTGET59B5DlnE9y/v/+FdiXNZVD8kzLUveTnvRklxhjFsDNOVNZxg0EjGk8Hojs7OyXgTcAiy/Kth/OaQDdyUFIXizbVs/rzQlrr28kZG8ekc/oRagpSW5qgB+OuraPhuedo/vRT8HjQRkYSMmIEIaNGEjxsGJJWi6ehEU9DPd7GRlwlpbQtX47sdhMyahSR55xDxJzZaEJC1H4rh8Za21GjZIUy3L2pRNlujFEG3WRMUsw7dkCv7JjsDrrTyCcBPwHbAF/H5rtlWf7uYK/pCUa+B6/bx8ZFpaz/vgS9QcsJZ2QzeGISkqjnIjgCvM3NNLz8Mk3vvIvs82E+91zMF12EITPjkOETT1MTLV98SfPHH+MqKkJjMhF/198wzZ8fWKEXR6sSJilaoVQH3DP4JjhSMe3MkxTjjhskjPso6ZMDgrqSphoby9/NpSq/mcQcE1MvHog5IVRtWYIAR5Zlmj/4gLrHn8DX1obptFOJuekmDCkpR3Ws9g0bqHviSdo3bCB81iwS7r8PnVmlGiBet5LHXbgMipZBxXqQvaALUVIAs6ZA1mRIGN4rcrgDAWHkXYAsy+xaXc3qTwvwuH1MPCuHoZOTA6tVJAgYZJeLmn8/QPNHHxE6cSJxd97RJXFu2eul8bXXqH/6f+giI0l88AHCTjyxCxQfBpYiKFiimHfxSnBZAUlJB8yeqph36jjQiRBkdyCMvAuxtThZ+tYuynZYSB8WzbRLBoncc8F+eJqaqLz5Fuzr1hF93XXE3nIzUhd34Dl27aLqjjtw5hdgvvhi4v92J5Kui7M7nFbFsAsWK/NP7olzR6Ypdbizpynhkj6ex+0vhJF3MbIss3VZBb98VoghRMv0ywaTPlRktgjAmZ9P+fU34KmrI/GBBzCdekq3ncvndFL32GM0vfU2EXPnkvTwf4/NzPedd7JgiTIQx+cBQ5gS495j3lFZIs6tAr2mHnmgIEkSI6alkjLAzKJXd/DNM1sYMSOVE87IRivyzvsstl9+oeJPNyKFGkl/521Chg/v1vNpgoJIuPtu9HFx1D36GLLHQ/JjjyLpj6B+kKNF6ZzMX6SYt7Va2R4/DE64UZm6LHUc6MSvzkBFtMi7AI/by+pPCti2opKkfpHMvHqISFPsg7Rv30HZpZeiT04m9eWX0Cf4d97GxjfeoO6h/xI2fTrJTzyOxnAQ45VlqN8NeQsh/8e9s70HmZQ4d7+TIXs6RCT6Vb/g0IjQih/I/bWG5e/sxmDUMfuaoSTmRKotSeAnXKWllFxwIZrgYNLffx99fJwqOizvvEvtv/9N2OTJJD/9FJqgjgaFu10ZjJO/EPIWQUuZsj1+KPSbqZh3ynGgFZVAAxlh5H6ioaKN71/cRlujgwln5zB8aorIaunleBoaKLngQnxWK+nvvUdQVqaqepo++JCa++8nbOJ4UhaciFTwoxI6cdtBb+yY7X2mspiSVdUqODJEjNxPxKSEce5dY1n8xi5+/iifhjIrUy4aiFYv4ua9EW+bjfJrr8NTX0/6G6+ra+KyDNVbMCcUI08Lo3bpGuotPxJ3khlGXQz9Z0H6pB45oYLgjxFG3g0EGfXMXTCMdd8Ws+7bEloa2pmzYJiod97LkF0uKm++GUduLinPPkPIyJH+F+F2KMPfc7+D3B/AWgVIRI05HqcnhcaVuwm+5r9EzJ3rf20CvyGMvJuQNBLHn5qFOSGUJW/t4pOH1jPvhhFEJYnRoL2Fuscew7Z6NYkPPED4lCn+O7HdonRU5n6nZJm4bUqp15xp0P9eJWQSFkvCJS6cl11O1T33YsjKInjgQP9pFPgVESP3AzXFLXz3/Da8Li8zrxlKuqik2OOxLl9OxYLrMV90EQn/d2/3n7CpBHZ/p5h36WplKHx4IgyYAwPmKoNyDhAycdfVUXL2OUgGAxkff6TecH5BlyA6O1XGanHw7XNbsVS2cdIFAxh6kuhk6qm4a+soPv10dHFxZHz04d7MkK5ElqFmK+z+Vllqtyvb4wYrxj1wLiSOOqxyr+1btlB68SUYjzuO1Jde7PrRnwK/IYw8AHA5PPz46g5KtjUyelY64+dniSqKPQzZ66Xsqqtp37KFzE8+Jig7u+sO7vMqIyl3fwu7v1Hqd0saSB0PA+cp5h2VdVSHbv70U6rvuZfoa68l7rY/d51mgV8RWSsBgCFYx5wFw1j5QR4bF5bS1uRg2qWD0OpERktPofHV17CvWUPiv//VNSbucSqpgbu+VsIm9kbQBikDc066QwmdhMYc82kizzqL9s2baXzlFcKmTsE4atQxH1MQOAgj9zMarYbJFw4gPDqYNV8UYWtxMue6YQQZxUCMQKd982bqn3qK8DmzMZ111tEfyGlVRlTu+lpZu6wQFKF0Ug46RRkSHxTedcI7iLvzTtpWraL6rrvJ/OJzNMEiDbG3IEIrKpL7aw1L39pFZLyRU28aQZhZfLACFW9bG8XzTwcg84vP0YYfodHaLZD7vWLehUvB64TQWCXePeg0yDzRL6VfbatXU3blVURdcQXxd97R7ecTdC0itBKADBiXQKjJwHcvbOPThzdw6k0jRXpigFL/1NO4q6pIf/fdwzdxa41i3Lu+VuaolL1gSoXjroJBpyqFqPw84ULohAlEXnA+ljfeIPzkGRhHj/br+QXdg2iRBwD15Va++d8WvB4f8/40gsRsk9qSBPvQvm07Jeedh/n880n4+//98c5NJXvNu3wtIEN0jmLcg05TJmBQuWSDz2aj6LT5oNOS9cUXPWP+TwEgslYCntaGdr56ejNtTU5mXT2EzBGxaksSALLHQ8m55+Gpryfru28P3Bqvz4VdX8HOr5SUQYCEYYpxDzotICcXtq35lbLLLyfqskuJv+suteUIDhMRWglwImJCOOv2MXzzzBa+f2EbUy4ayOBJSWrL6vM0vfsujp07SX7yib0mvifHe9fXink35CrbU46Hk/+ldFgeZZqgvwgdPw7zRRdheettwk8+GePY33mDoAchWuQBhsvhYeHLOyjb0cjxp2Yydu6hZ1kXdA/u6mqK5p1CyNgxpD7/PFLluo6wyVd7c7zTJ8Lg+Uqed0TP+uL12e0UnXIqmrAwMj/7VAwU6gGIFnkPwRCsY+4Nw1j+9m7Wfl2MrcXFSef3RyMGDvmdmn//G9nrJmGihPTEIGirBa0BsvbkeM+F0J5bbkFjNBL3tzupvPkWmj78kKiLLlJbkuAoEUYegGi1GqZdNgijKYiNC0uxtziZedUQdAb/Zjj0SVx2KFyC9dM3aFuyk9gRrRgqvoZ+M2DgqdB/JgT3ns7o8JNPxjh+PPVP/4+IuXNFLZYeigitBDhblpbz88f5JGabmHv9cIJDxcChLsdugbwflKHxBUvwORwUfh+PNiyczGfvQxowA/S9N7PDkZdH8RlnYj7vXBL+/ne15Qj+ABFa6aGMmJaKMcLA4jd28tkjGzjlphFERPdeU/EbTaXKkPjd3+6tJhiRDKMvwbJZg8f2FUnPP4c09Hi1lXY7wf37Y77gApree4/Ic88V5W57IKJF3kOozG3i+xe3odVpmPen4cSlR6gtqWfh80H1JmXyhdzv9lYTjB20tyBV0mg8FguFM2dhHD+e1GefUVezH/G2tFA4azZB/fqR9tabooM9QBEt8h5O8gAzZ/5VSU/8/LGNzLp6KBnDj72YUq/G3Q7FK/fOntNWs7ea4MwHDlhNsOHZZ/E5HMT95S8qiVYHrclE7K23UnP//Vh/+IGIOXPUliQ4AkSLvIdha3Hy7bNbaSi3ctL5/Rk6OUVtSYFFS2XHTPELoWgFeNrBEAY505Usk34zwRh1wJc6i4opOvVUIs89h8T77vOzcPWRvV6Kzz4Hb3Mz2d99K0Z8BiCiRd5LCDUFcfpto/jx1R2seD+P5rp2JpyZjUbbR0vhej1QuR7yFylLzTZle2QajL5EmXA44/AKUtU99hia4GBib7yxm0UHJpJWS8I9d1N68SVY3nyLmAXXqS1JcJgII++B7Klr/vMnBWxZUo6l2sbMq4b0nYyWtnooXKIYd8EScDSDpFWKUM34B/SffcTD4u3r19O2ZAmxt96CLrrn5oYfK8axYwmbOpXG117DfMH5aE29J9WyNyNCKz2cnT9XseL9XMKjg5l7/XCiEnth9USvGyrWQcFiZaneomwPjYN+JytL1lQIiTyqw8uyTMl55+OprSX7h+/7fEjBkZtL8fzTxWxCAYgIrfRSBk9Kwpxg5PsXt/HJf9cz86ohZAzrBZ2gliKlbnfhMqXD0tm6t9U97f+UmHfCiMOas/JQWL//HsfWrSQ++GCfN3GA4AEDiJg3D8vbbxN1ycXoYkUBt0BHtMh7CVaLg+9f2EZ9uZXj5mYwdm5Gz4qbtzcphl24DIqWKeVgAUxpyrRnOTMga3KXj6qU3W4K585DYzQq9Ua0YvQsgKu0lMK585TSvf93r9pyBB2IFnkvJzwqmDP+OpqV7+Wy7tsSKvOaOfnKwYE765DbARVrlfkqC5dB9WaQfUqGScaJMP5PSqs7KqtbS8C2fPUV7vJyUp57Tpj4PhjS04k86yyaPvqIqCuuwJCSrLYkwR8gWuS9kNxfa1jxXi4ancT0SwcFRm1zn1eJbRevUMy7bA14HEq4JOU4yJqitLyTx4DWP522sttN4Zy5aE0mMj75WAyC+Q3umhoKZ84iYt48kv7zoNpyBIgWeZ9iwLgE4jMiWPTqDr57fhvDpqZwwhnZ6P1ZdEuWoSFPyeUuXgElP4GjRXkubjCMuUIx7/QJEKzOKNWWL7/EXVFB/D13H5aJyz4Zt8uL2+HF5fDgcnhxO73IPhlZlpF9SscpgFavQafToNVr0Oo06IO0BBl1GIJ1SD2kkqU+IQHzhRdieestoq+6kqCcHLUlCQ6CaJH3YrxuH798XsiWpeVExIYw5cIBpA468GCYLqG5TIlzF61Q1m01ynZTGmSdBJlTIPMkCI/vPg2Hiex2Uzh7DtqoKDI++hBJknC7vDTX2mmutdNSZ8fa5MTe4sLW7MTe4sRudSP7ju3zIklgMOoIMuoJCdMTagrCaDIQajJgNAURZg4iPCqYsKhg/37xHgRPUxOF02cQOmkSKU8/pbacPo9okfdBtHoNk87tR8aIGJa/u5uvntrMwPEJTDy7H8FhXRC+sDVCyUolVFK0ApqKle2hsUqcO2syZE6GqMxjP1cXIssyFe9/RaUzFunEa9n+zBYsVTbampz77RcSrscYEURopIGYlDCMEQaCjHr0wVoMIVoMQTr0QVokrYQkSUgSSmtbBq/Hh9ftw9Oxdjk8OO0enHZ3x9pDu9VFc52dyvwmnDbP73SGhOsJjwomIiaEiNgQTLEhmGJCMMWFEGoK8kvLXmc2E3XFFTQ8+yyOnTsJHjy428/Zm/G4vN1SjrpLWuSSJL0GnALUybI89FD7ixa5//G4vKz/roRNi8owGHVMOjuH/scnHJkZuNuVSoFFy5WlZhsggyEcMiZ1GPdJSugkgOLNbqeX2pJWagqbqS5oobakFaddMU6NRsKcFEp0cijmeCOmOCPmBCOmWCP6IP+1iL1uH7ZWJ20WJ1aLY+/S6KC1vh1rowPfPr8GdHoNprgQTHFGIjvW5ngj5oTQrvmS3leb1UrBtOmETphAylNPdumx+woel5etyyvYuLCU024eedRF77p18mVJkk4C2oC3hJEHNo2VbSx7Zze1xa3EpIYx7rQs0odGHzhG7PMpVQILlyopgaW/gNcJGr2Sz501RTHvpNGgDZwfd852D9X5zVTmNVGV30x9eZsSEpEgKjEUM43oln1Ozs0Xk3LKJHR69UMYh8Ln9WG1OGmtb6eloZ3mOjsttXaa69ppbWjH5937OQ4O1WNOMBKZYCQqMVR5z4mhhJmDjrpDt+6pp2h84UWyvv5KxMqPAJ/Xx+5falj7TTG2ZidpQ6KZdE4O5oSjG7jXrUbecYIM4Bth5IGPzyeTv66WtV8X0drgIDHbxPjTs0jqZ1bCJUXLOkZRLgFbnfKiuMHK6MnsqUoHpSFwRpC6nV6qCpqp3N1EZV4T9WVWZBm0Og3xmREkZptIzIkkISsCg06mYPZsdLGxZHzwQa/IVPF5fbQ2Ojrj+021dppr7DTV2Gi3ujv30wdpiUoKJToplKjkMKKTQolODiMk3HDIc3iamiiYPoPwGdNJfvjh7nw7vQLZJ1O4qZ5fvyqiudZOfGYEJ5yRTXL/Y5uBSXUjlyTpWuBagLS0tDGlpaVdcl7B0eP1+Nj1cyXrvinA3iaTHF7MMO1HZAatRWOMVPK4s6cpBh6RqLbcTrxuH7UlLVTsbqIit4na4lZ8XhmNViI+M4Lk/maSB5hJyIz4XTyy6cOPqLnvPlJffomwE09U6R34j/Y2F03VNizVdixVNizVbTRW2nC07TV4o8lATEo4MalhxKQoS2Sc8Xdht9qHH8Hyxhtk//A9hrQ0f7+VHoHPJ1O4sY7135VgqbJhTjAy/vRsMkfEdEmjQXUj3xfRIlcZl00ZhJO/EPIW4bFa2Gafy1bn6bS5TISGSww6MY0hJyYHxIAir8dHbUkrlblNVOY1U1PUgtftQ5IgNi2clIGKcSfmRP5hpofs8ezNVPmwd7TGjwZZlrG3urBU2miobKOxso2G8jaaqm2dcXh9kJaY1DBi08KVJTWccJ2dopNPxjT/NBL/9S+V30Vg4fP6yF9fx4bvS2iqsWNOMDJmTgb9xsZ16QhrkbXS17E1Qt73ytRmhUuVwThBEZA9DV3/2YzqdzIjgqMo3d7I9pVVrP++lA3fl5Iy0Ez60BjSh0YTGW/0i1R7q4uaopbOpa7UitftAyA6JYyhJyaT1D+SpH6RR1TxsfWHhUre+F1/67MmDiBJEqGmIEJNQaQO3puO6nX7sNTYaCi3Ul/eRn2plZ0/V+FxKddeZ9BgmvIPjFs20m/hLhJHpmOKC+nT19Jpd7NrdTXbllfQ2uAgOjmUmVcPIXt0HBo/jhcQLfLejLUGdn4FO7+EstXKEPiIlL1Tm6VNAN2B46OtDe3s+LmK4s31NNXYATDFhpA2NJrELBNRSaFEJhjRHkNrw+eTaa1vV1qEFW2da2ujAwCNViI2LZyETJNi3DmRR52RIcsyxaefgezxkPX1V0hdUGyrL+DzyTTX2qkvs1JX2kptXiP1pS34tMp9E2TUEZcRQXxGROfaGHHomHtPp6nGxtZlFexeU4PH6SUxx8TIGWlkDo/p1rTQ7s5aeR+YAsQAtcB9siy/erD9hZF3I53m/YWSKoiszEs56FTFwBNHHHFqYEt9O2U7Gind3khFblNn61ijlYiMNxKVFIox3KAMdAnREWTUoQ/S7c2ldnvxuHw42z20NSkpdW0WJ23Nzs4BNpIEkfFGolPCiE0NJzHbRGx6eJdllLStXEn5tdeR+J//EHnG6V1yzL5K5b1/p2rxWoLve4oGi0xtcSuWKlvn3zIsKoj4DBPxmRHEZ0YQmxYeEIObjhVXu4fCTfXk/lpDZW4TGp1E/7HxDJ+WSmxauF80dHuM/EgQRt7FOFph19ew9UNlROUe8x5yOgw+HeK6blZ0r9tHU60dS5XSaWapasNSY8fR5sbV/vtBLfsiaSTCIoMIjw4mLCqIcHMwEbEhxKSEEZUY2i0DJfZQesmluCoqyFn4A5Kh97cYuxNXWRmFs+cQdfnlxN9xO6BkDtWXWaktaaWupJXa4lasFuWXlaSRiE4O7Wyxx2dEYE4w9ojqnF6vj/KdFvJ+raF4SwMet4+I2BAGnZDA4EnJfv/1IWLkvQ2vBwp+VMw793sl5m3OgJNuh6Fndal574tWr+nMbPgtPp+Mu2MEo9vpRdtRa0Rn0KDTa9HpNarUGbFv2oR93TolNi5M/JgxpKURcco8mj74gJhrr0EbGYk+SEtSP6XfYg+2Fqdi6h3mXrihjp0/VQHKgKbolP07U82JxoDI6W9rclC6vZGyHRbKd1twO7wEheoYOCFRqWOUGRFw/QKiRd7TqM+Dze/Alg+grRaM0TDkTBh+HqSMDagRlYFC+Z9upH39enKWLkETGjj57z0ZR24exfPnE3vrLcQsWHBYr5F9Mi317dSWtFJfZlWWcituhxdQbl1TnJHopFDMScpApsg4IxExwQQZu6ciptfto7GqTekDKLNSU9iCpcoGQJg5iLSh0WQMiyFtcBRanfq/IESLvCfjssOOz2DjW1D+q1L6tf9sGHWxMs2Zn8q+9kScBQW0LVlCzJ/+JEy8Cwke0J/QE0/E8s67RF1xBZqgQ09uLWmUPpXIeCMDxiUAHebe0E59mbUjz91GY5WNos317NvGDDLqMMWGEB4VTEiEAWOEgZBwZR0cqkOr0yq//vTKr0BZljv6Z5R+GrfLS3uri7YmZ8eilECwVNk6R8UGGXXEpoUzYHwC6UOjiUoMDbiW98EQRh7I1OfC+tdg8/vgbIGY/nDyv5TWdwBUEOwJNL7yKlJICOaLL1JbSq8j+qorKbv8Clq+/BLzuece1TEkjURknJHIuP1TWz0ur1J+oL6dlnqlDEFLQzuWahv23KbOWjlHQ5BRR2ikUmUybXB0Z3gnIia4xxj3bxFGHmh4PbD7a1j7CpT+rNQ1GTwfjrsK0k4QoZMjwF1VRcs332C+8AJ05mMbGi34PcZx4wgePBjL628QefbZXZrSqTNoD9oXA8ogsXari3arG6fd3dny9rh9eFxeJI2ETq/0zWgNSkvdGGEgNDIIQ3Dvs73e9456Ku1NsOFNWPsytFZAZBrMuB9GXgxhATDDTw/E8uZbAERffrm6QnopkiQRddWVVP3lr7QtW0b49Ol+O7dWpyHMHBwQI48DAWHkatNQAGuegy3vg9uu1PGe+7ASA9eo34PfU/FarTR/8gkRc+agT0pSW06vJWLWLOoff4LGV171q5EL9kcYuVpUbIBVT8Cub5TOymHnwvgFkDBMbWW9guZPPsVnsxF12WVqS+nVSDodUZddRu2DD2LfuAnj6FFqS+qTqJ9P05eQZcj/Ed44BV6ZpgzeOfE2+PMOOP1ZYeJdhOzxYHn7LYzHHUfI0CFqy+n1RJ51JhqTicbXDjqYW9DNiBa5P/D5IPdbWPmIMpN8RDLMehBGXwpB/hna25ew/vgjnqpqEu69V20pfQJNaCjmC86n8cWXcBYXE5QZWFP79QVEi7w78Xlh+6fwwkT48GJwWmH+s3DzZjjhT8LEuwFZlml8/Q306WmETZmitpw+Q9TFFyPp9Vhef0NtKX0SYeTdgc+nGPhz4+GTKxVDP/Nl+NM6ZRDPQSoOCo6d9k2bcWzdStRll4kKh35EFxODaf58Wr78Ek9Tk9py+hziTu9KZFnpvHxhkmLgkhbOeQNuWAPDzw2oeS17K5Y33kBjMhF5+ulqS+lzRF12KbLTSfMHH6gtpc8hjLwrkGXIXwwvTYEPL1ImKD7rVbh+FQw5A0TL0C+4ysuxLl6M+bzz0Bj9MwmGYC9BOTmETpqE5b338LlcasvpUwiHOVbK1ylZKO+eBe0WmP8c3PArDDtb5IH7Gcvbb4NWi/kiMRxfLaIuuwxvfQOt332ntpQ+hfitf7TU58GSf8DubyA0FuY+CqMvE/FvlfC2ttLyyaeY5s5BHx+ntpw+S+ikiRiys7G8+Ram+fN7bO2SnoZokR8p1hr46mZ4bhwUrYCp9yhZKMdfI0xcRZo/+RSf3S4GAKmMJElEXXopzl27sK9bp7acPoMw8sPFZYPl/4WnR8Pm9+D4a+GWzTD5Dgg6cGEfgX+QvV6a3n0X49ixBA8erLacPo9p/mloIyM7a90Iuh8RWjkUPq9i3Ev/DW01SiXC6fdBdLbaygQdtC1bhruykrg77lBbigDQBAcTef55NL74Eq7SUgzp6WpL6vUII/8jin+CH+6C2m2Qchyc+xakjVNbleA3WN5+B11SIuHTp6ktRdCB+cILaXz1NSxvv0PCvfcAylSALe1uGm0unB4vXp+Mxyfj7ViMBi2hQTrCg3SEBukwGrQixn6Y9Cgjf/zHPL7YVIkkgYQSj5MAg05DpFGPKURPZIgBk1FPQkQwaVFG0qKNpJqNhBzJxL6WYvjx/5QJjU2pcPZrynRq4qYKOBy5udh//ZW4v/4FSdejbudehyzL1Lc5KayzUdTQTvTwCbg+/JgF+jFUenRYbC68vsOfWlKnkUiMDCY5MoQUs5HkyBDSo40MToogOzYMfQ+YvNlf9Kg7P9Ucwui0SABklPRtGXC4vbTY3ZQ02Glub6bZ7sbp8e332tjwIAYmhDM02cTQJBNDkyNIizLu/43vtMJPj8Evz4JGB1PvhQk3gj7Eb+9RcGQ0vfMOUnAwkWefrbaUPofd5WFLeQsby5rYVNbM5vImGtr25o8PMo/hcfcKZpf+Ssn0M4gOMxAdGkRUqIFgvRa9VkKrkdBpNGgksLu82FwerA4PNqeHJrubquZ2Kpvb+Tm/gVqro3P6N4NOw4D4cIYkRTA8JZITsqPJiDb22RZ8r5x8WZZlmuxuyix2yix2yi12ihts7KpuJa/Wirtjjr7wYB3HZ0RxQlYUs+WVJK97CKmtBkZcoMTBIxK7TaPg2PE0NVEwZSqm+fNJ/Oc/1JbT6/H5ZLZVtrAir57luXVsqWjpbGFnxYQyKs3M0GSltZwdF0ZiRDDll16Ku6qK7EULj/kXk8vjo7TRxs7qVnZWtbKjqpUdVS002d0AJEQEc0J2NCdkRXNS/1gSTL1v0ok+NfmyJElEhRqICjUwMjVyv+dcHh95tVZ2VLWwubyFhvx1jCh6nhRNHtvJ5tuUZ8lJm8J0XTSRBzy6IFBo/vgTZKdTzMfZjTjcXlbk1fPD9hpW5tXTaHMhSTA82cT1k7MZk25mZGok5tADp96aL72EyptvwbpsGREnn3xMWgw6Df3iw+kXH878kcmA0mgrarDxS2EjvxQ1sjKvns83VQIwIsXEzCEJzBqSQE5c784s65Ut8sPCboGl/4L1r+MNiWJz/1v40H0SKwoaqW11otVIjMuMYubgeGYOSSApUoRXAgnZ46FgxskYMjJIf+N1teX0KtxeH6sKGvh6SzWLdtRgdXowG/VM7h/LlAFxnNgvhuiwoMM6luz1UnjyTPTJyaS/3f3piLIsk1trZcmuOhbtqGFLRQsAWbGhnDYiibNGp5Aa1XPLNxysRd73jNzng83vwI/3gaNFGcgz5S4Iiex4Wvn5uGhnDYt21JJf14YkwfjMaM4YncycoQmEB+vV0S7opPWHhVTeeispzz1L+DSRrdIVFNa38cHaMj7dWInF5iI8SMesoQmcNiKJCdnR6I6yc7Hx1deoe+QRMj//jOBBg7pY9R9T1dzOjztr+X57Nb8WW5BlOD4zirNHpzBnWM/7LPcKIy/ZspG6kiJACZ/QkbUiaTRotFo0Wi2SRllrtVq0ej0anR6dTqesW0vQ/vI0uobtaJNHoZv5d3RJw9AZDGj1+gN2lBTVt/H1lmo+31RBSaOdYL2GkwcncM6YFCblxKDR9M3OFbUpuehiPLW1ZC/8AUkratocLQ63l4U7anjv1zJ+Lbag00icPDieM0YlM3lALEG6Y7+23pYW8qdMJWLOHJIefKALVB8dlc3tfL6xgk83VlLcYCNYr2H+iGQuOSGdockm1XQdCb3CyBe/8hxbfuymYjyShE5vQGcwoAsKQm8I+t3a6tVQbvVQ1OSmzashNCyU0TkJHN8/kciIMIKMYQSFhhIcqqyDQoyiJnY34Ni5k+IzzyLuzjuJvuJyteX0SOqsDt75pZR3fi3DYnORFmXkguPTOHtMCrHhhxc2ORKq77+fls8+J2f5MnRRUV1+/CNBlmU2lTfz8fpyvthURbvby5h0M5eekM6coYkYdIH7me0VRu71uPF5vSCDjAyyrKQgyj58Xi+yT1n7vB68Hi9etwvf7oV4Vz+Lt70Vb7+5eAafjUfS43W78biceFwda7cLj8vVsc2F2+nE43Iqa6cTt9OhLA4HLoeyVpIfD44kaQgOD8cYYcJoiiQkwkSoKZKwqGjCo2M6lljCoqLRihzow6bqnnto/e57+q1YjjYiQm05PYrcGiuv/FTEl5urcPt8zBgUz2UnZDAhO7pbf106CwspmncKsbfeQsyCBd12niOlxe7m4w3lvLOmlJJGOzFhQVw+IZ1LxmdgMgZe2KVXGPkR0VAA394GxSsgaTSc8gQkjeyyw8uyjMfpZEdpHZ/+Wsjy7eX4nA5GxBuYlhlORrgGp62Ndmsr9pYW7K3N2FtasDU34Wq373csSdIQEReHOSGJyIQkzIlJRCWlEJueSWikucs078En+3B5Xbh9bjw+T+fa4/Pgk32di1f27qNRYs8/jUaDXtKj1WjRSlp0Gh0GrYFgbTA6ja5bc3k7Uw5PP53Ef9zfbefpbawvsfC/pQWsyKsnWK/hnDGpXDExg6xY/2VzlF11Nc78fHKWLEbSB5ZJ+nwyPxU08OrPxazMq8do0HLB8WlcNSkzoBId+o6Rux3w8+Pw8xOgC4Hp/wdjr+z22uCtDjfv/1rGa6uKqW11MiA+nGtOyuK0EUm/+6nmardjbWzE2liPtbGB1vpamqqraKqpormmCld7e+e+RlMk0alphCcnYkyJJzg1Dk+4njZPG22uNtrcbdjddtrcbdjcNuxuO3aPnXZP+36L0+vE6XHi9Dpx+bqv6L+ERJA2iCBdECG6EIw6o7LolXW4IZwwQxjhhnAiDBGEG8KJDIrsXExBJiKDItEe5O/V+Mor1D36GJlffUlw//7d9j56A7Iss6bIwtNL8vmlqJGoUANXTszgonHpB00X7E7aVqyg/LoFJD36KKZT5vn9/IfLzqpWXlpZyNdbq5GA+SOTuXFaDpkxoWpL6yNGXrQCvvkzWAph2Dkw8wEIj+/68/wBLo+Pr7ZU8fLKInJrrSRHhnD9lGzOGZuCjJtmZzMtzhaanc2djzsXVwstjmZsLc14G1rR1rcT3OQlokWDuc2A1qe0dNsNXurMTuojndREO2iMcBFsCCFUH7qfaYboQjqXIF2QYrAdi0FrQK/Ro9fo0Wl0nWuNpNlvkZCQkZFlmT3/fD6lte7xefDISkve5XXh8rqULwqvC4fXQbunvfOLZc/a6rJidVlpc7fhk30HvIYSEuZgM1HBUUQHRxMVEkVsSCxxQTGMvvFlpOQEIl58kvjQeIK0XR/P7enIssyqgkaeWpLHupImYsODuO6kLC4cl4bRoF4IT/b5KJozF21kJBkfBv50cBVNdl75qZgP1pXh8vg4vcPQ/fkr5rf0biO3NcCie2HL+2DOhFMeh+zuTUnz+rw0O5tpcjTR5GxS1h2Pm53NWNotFFnqKG6qw+GzotHaQXPwlnCILoQIQwSmIBOmIBMRhoi9S1AEYdpQDE1uqGrBUVZHW0kltvoGAAzGUNKGDCN92CjSho3EnJgU8EOVfbIPu9tOi0v5UmtxKOs919LisNDY3kijo5HG9kYa2hsYtsvO7Z/6ePRMDWsHKL9yooKjSAhNIMGYQGJYIkmhSSSHJZMUlkRSWBKmoJ6RjdBVbCi18MjCXNYUWUg0BXP9lGzOHZtKsD4wMnssb79D7QMPkPHRh4QMH662nMOi3urkpZWFvL2mFJfHx/yRydykkqH3TiOXZaXE7KJ7wdkKE2+Fk/56VLVRPD6PYsAOS6eR7Fn2mLTFYek0mhZni9LhegDC9eFEBiuhAnOQGZcrhNwqLzVNWsL0EcwZlM3codnEGqM6wwlH07K0tzRTvnMbpVs3UbptM631dQCYE5PIHjuenONOILFffzS9YMo5WZYpvvwynCUlNL3zAHXORmpttVTbqqmx13Q+trlt+70uXB9OSngKKeEppIankhqeSlp4GukR6cQZ4wL+C+9w2VHVwmOL8li6u46YMAM3Ts3hgnFpXZI+2JV422wUTJlC2NSpJD/ysNpyjoh6q5OXfyrirV9KcHtlzhmTws3T+/k1ht77jLyxUAmjFK+A1HFw6lMQt3ewgSzLWN1WLO2W/Uy50dH4u20Wh+WgxqyRNJ2GbA42d/7kNwebMQftfRwZFElUsGLMeu3vO3L2xCufWJzH2mILyZEh3DK9H2eOTj7qgRa/PX5LbQ0lWzZSuOFXyrZvxef1YDRFkj3meAZMOInUIcN6rKk7CwooOuVUYm+7jZhrrzngPrIs0+pqpbKtkqq2KirbKim3llPRVkGFtYLKtko8Pk/n/iG6EFLDU0mPSCcjIoNMUyYZERlkmDIIN4T7660dE+UWO48uyuXLzVVEBOtYMCWbyydkqBpCORQ1Dz5I0/sfkLNkMfq4njctX73VybPLCnjv1zKQ4NLx6dwwNYeow+h3cLuc6PSGo25A9Aojd7mdNBRuo7nycyyb3qJRp8cyYBaWmGwszqZOk250NGJxWPb70O6LKcikGHCQmeiQaKKCo/aac7BZict2/N9kMB204+1okGWZn/IbeHRRLlsrWsiMCeXPJ/fnlGGJXZr+5bTbKN60noL1v1K8aR2u9nZCzVEMnHAigyZNJS4zu0e1Rmv++U+aP/mUnBXL0ZmPLpPH6/NSY6+hrLWMstYySlpLKLOWUdpaSoW1Yr8snZiQGDJNmWSZspQlUlnHhsQGxHVrsrl4ZlkBb/9SikYDV07M5LrJ2ZhCAisb5EC4SkspnD2HmBtuIPamG9WWc9RUNNl5anE+n26swGjQce1JWVx9YubvvkR9Xi9l27ewe9UK8teu5ux7/k1ivwFHdc5eYeTfXj4N045qbrxei1e798Nk0BiIDonu7Bjr7CQLjvrd/yODI9Fr1L/ZZVlm0c5aHl+UR26tlaHJEfxt9iAm9Yvp8nO5XU6KNqxj96rlFG1cj8/rwZyUwrCpJzP4pGndkuLYlXitVvInTyFi1iyS/vNgt5zD7XVT3lZOSUsJxS3FytJaTHFzMVa3tXO/cEM42aZssiOzyYnMITsym37mfkQHR/vF4B1uL2+sLuHZZQXYnB7OGZPKn0/u3+Mq/ZVft4D2HTvIWboEjaFnz3VbUGfl0YV5/LCjhtjwIG47uT9nj06mobiAXT8vJ/eXn7C3NBNkDKXfuAmMPeVMolNSj+pcvcLINzx0OsY3cmm6ahLGy/6kGHRINEZdz61D7PXJfLWlkkcX5lHZ3M6J/WK4c/bAbhsy7GhrI+/XVexYsYSq3J1otFqyx4xj2LSZpI8YFZChF8tbb1P74INkfPIJIUOH+PXcsizT0N5AYUshRc1FFDYXUtBcQGFLIS3Ols79zEFmcsw55EQqS39zf3IicwgzdE2HmCzLfLO1moe+301lczvTBsZx5+yBDEjoGSGg39L28yrKr76apIf/i+m009SW0yVsKG3isc/W4Mpdx7D2fMIcTWj1erJHH8/ASZPJHDkW3TF+aXWrkUuSNBt4CtACr8iy/NAf7X+0Ri47Wik++0KQJDK/+qrHmveBcLi9vLOmlGeWFdBsdzN/ZBK3zxpAirn7KrU1VpSzbdkidq5YQru1lfCYWEacPJdh02ZijAiMbI/OlDWzmYwP3ldbTieyLNPoaKSguYCCpgIKmgvIb86noKkAu2fvgK+k0CT6mfvR39y/c0mLSEOnOfwY9qayJv71zU42ljUzKDGCe+cNYmJO1/9y8yeyLFM07xQ0oaFkfvyR2nKOCbfTQf6vq9mxYjFlO7aBLNMYnsKmoBzMg8dw1/zRDEvpms9Ttxm5JElaIA84GagA1gEXyLK882CvOZbOzuYvvqD6b3eR+uILhE2efFTHCGRa2t28sKKQ134uRgaumJjBDVNyujX26fW4KVj3K1t+/I7yHVvR6Q0MmHgSo2afSnymupNMt/30M+XXXEPSI49gOvUUVbUcDrIsU2WrIr8pf+/SnE9JSwkeWemzMWgMnSGZfQ0+OiR6v2NVNbfz3x928+XmKmLDg7h95gDOGpOCtpcUarO89x61//wXGR9+QMiIEWrLOSJkWaa2qIDtyxax6+cVuNrtmOITGHziNCVcGRPHB2vLeGJxPhabizNHJXP77AEkmo4tw6U7jfwE4H5Zlmd1/P8uAFmW/3Ow1xyLkcsuFwUzZ2FISyP9rTeP6hg9garmdh5dlMvnmyqJDNFzy/R+XDQ+vdvnKWwoK2HTwm/Y+dMyPE4nyQOHMPbUM8kefZwqBcDKr1tA+84d9FuyBKkHx1JdXhfFLcXkNeXttzS0N3TuExMSQ39zfzIjciirMbF0qw7ZGcs1J/ZnwZRswoICNxPlaPDZbORPnkLYlCkkP/qI2nIOC6fdxs6VS9m2dBH1pcXo9Ab6j5/I0GkzSRk09HdRglaHm+eWFfLaqmI0ElxzYhbXTT76v2V3GvnZwGxZlq/u+P8lwDhZlm/8zX7XAtcCpKWljSktLT3qcza+/gZ1//1vjxpUcLRsr2zhP9/vYlVBI5kxodw5eyCzhsR3e1jJ0dbG9uU/svH7r7A21GNOSmHsvNMZdNJU9Ab/jKZ0lZVROGt2j89u+CMa2xvJb84nz5JHriWXdVU7qG4vAUnJoNFJOrIiszpb7QPMA+gf1d9vnavdTe1//oPl3ffIWbIEfXxgpiLKskxtYT5bFn/P7tUr8TidxGVmM2zaLAZOPIng0EP3g5Rb7DyyMJevtlTx0iVjmDkk4ai0qG7k+3KseeTeNhsFU6cSOmECKU89edTH6SnIsszy3Hoe/G4X+XVtHJ8Zxb3zBjE8JbLbz+3zeslb8zPrvv6MuuJCQiJMjJ59KiNnnUJwWPeObKt96L9Y3nknoD/kXcWG0ib++c1OtpQ3MzwljKumhqMPqSG3KZf8pnxym3Kps9d17h8VHEU/cz/6RXaEZ6L6k23KJljXs7JXOlMRr7+e2JtvUlvOfrgdDnatWs6WRd9TV1KIPiiYgZMmM2LGHOKzco7qmLtrWhkQHx54eeT+Dq3soe6xx2l89VWyv/8OQ3r6MR2rp+Dx+vhwfTlP/JhHQ5uL00cmcfvsgST7YWSZLMtU7NzGuq8+pXjzBgwhIQyfMYcx804nzNz19aV9djv5U6YSNmkiyY8/3uXHDxQqmuz894dcvt5SRXxEEHfMGsgZo5IPOKag2dFMfnM+uZbcztBMQXMBTq8TUAavpYWnKQZv7kf/yP70M/cjOSy5S8dCdDXlC66nfds2cpYtDYhUREtVJVt+/I4dyxfjtNuISctgxMlzGTRpCkFGdaeJ604j16F0dk4HKlE6Oy+UZXnHwV7TFUburqujcPoMTGefReJ99x3TsXoaVofSIfrKT/7rEN2XupIi1n7xMXlrVqHRaRk6ZQbHnXYWprij+7l4IJo+/Iia++4j/b13MY4e3WXHDRTanB6eX17AKz8VI0lw7UnZLJicdcQjMr0+L+XWcvKa8jpDNPnN+VRYKzpHKgdrgzvz3vuZ+3Xmv8cbuz9Edzi0rVpF+VVXk/jQf4g8/XRVNPh8Xoo3rWfTD99QunUTGq2WfuMmMnLWPJIHDA6I6wTdn344F3gSJf3wNVmW/3A+p66qtVL9f/9Hy1dfk7NsqeqzjqjBbztEb57ej4vGpftthpOmmirWffUpO5YvQZZ9DD5xGuPOOAdzYvIxHVeWZYrnnw5aLZmffRowH6KuwOuT+Xh9OY/9mEe91cnpI5O4Y/bALq/XYXfbKWop6sya2ZNB0+ho7NwnXB9OVmQWOZE5ZJmyyI7MVsXgZVmm6JRT0QQHk/HJx349t8PWxo7li9m08BtaamsIi4pm+IzZDJ8+OyAHyvWKAUG/xVlURNHcecTccD2xN9/cBcp6Jvt2iGZEG7l91kDmDkvw2wfC2tjAuq8/ZdvihXg9HgZMOJFxZ5xLTOrRhbxsa9dSdullJD7wbyLPOquL1arHirx6Hvx2F7m1Vsakm7ln3iBGp/nXLJocTcqApo6BTXseNzubO/cJ1YeSGZFJVmQWmabMziU1PLXbRkU3ffABNff/w2+/wBorytn0w9fsXLkUt9NB8sDBjJp9GjnHjQ/o2bp6pZEDlP/pRtrXr1fiayrHr9RElmWW59Xz0He7ya21MiLFxJ1zBjIh238DR2zNTaz/5nO2LPoOt9NBv3ETGH/m+cRlZB3RcSpuuRX7mjXkrFiOJrhndd4diF3VrTz43S5+ym8gLcrI3+YMZM5Q/33RHg4Wh4XC5o7Rqy2FFLcUU9RStF8Hq1bSkhKeQkZEhlJozJRBeng6aRFpxBnj0EhH/0twT59I6IQJpDz5RFe8pd8h+3wUb9nAxu++onTrJrR6PQMnTg6I8RKHS681cvvGTZReeCHxd99N1KWXdMkxezJen8znmyp5fFEuVS0OJveP5c7ZAxmc5L+5LdutrWz87ks2fv81rnY72WPHMf7M80nI7nfI17qrqymYcTJRl19G/O23+0Ft91FusfPEj3l8vrmSiGA9N03L4ZIT0gOutOwfYXPb9taeaSmmpLVEKTbWWtbZyQpKHD4lPIW08LT9SganhqeSGJZ4WC352ocfwfLmm+Qs/hF9YmKXvQdXu53ty5eweeHXNFVXEWaOYsTMeQyfMTtgRjAfLr3WyAFKLroYd3UVOQsXBtxcgGrhcHt5+xdlyH+rw80pw5O4dUY/sv1YDN9ha2PT91+z8bsvcdjayBg5hvFnnEfywMEHfU3dk0/S+OJLZP/4I4aUY4u1q0Vjm5NnlhXw7poyJAkun5jB9ZOziTSqn5HRVfhkHzW2GkpbSym3llPaWkpZa1ln2eB9TV4jaYgzxv1u0o8EYwIJYcqkIEa9EXdlJQUnzyT6qquI+8ttx6yxuaaaTQu/YfuyH3G120nMGcDouafRb9zEgA6f/BG92sitS5dRccMNJD3yMKZTT+2y4/YGWtrdvLSykNdXleBwezlzdAq3TO9HapT/wlBOu53NC79hw7df0G5tJXXwMMadeR5pQ0fsF17wOZ0UTJ1GyMiRpD73rN/0dRUt7W5e+7mYV38uxu7ycO7YVG6Z0e+Yh2X3NHyyj4b2Bsqt5ZRbyztrw++pE19rr/3dNH+mIBNxxjiueLee1IJWVj97OTHmZGJCYjqX6JDoQ07AIssyZdu2sPGHryjauA6NRkP/8ZMYPee0oy4dG0j0aiOXfT6KTjsNSasj84vPAyr2GCg0tDl5YXkhb60pxeeTOfe4VG6Ykt2tRbl+i9vhYOuSH1j39WfYmiwk5gxg3JnnkjX6eCRJovnzL6i+6y7SXn+N0BNO8JuuY8XqcPP6qhJe+amIVoeH2UMS+Ous/uTE9czKhN2N2+em3l5Pta1ameHJVkONrYZaey0h24q44vkiXpyjZcnI33+O98y+ZQ42ExXUMalLcCThGNHurMe+Ng9nXROGsFCyp5zE0OkziYtLJUgb1Ct8oVcbOdBpAqkvv0TYiSd26bF7EzUtDv63NJ+P1pcjy3DGqGRumOrfGcI9Lhc7Vixm7Zef0FpfR0xqOsfNPxvDE08judxkff11j/jQWR1u3vqllJd/KqLZ7mbGoHhundGv20oQ9wVkWab4jDORvR6M77+ExWnpnLO1ob2BxvbG/ebJddY3k5DvJas8BINHQ0OEk10ZVkoSbXj36YrQaXSdE5HvOzG5QWvonIzcoDVg0BjQaXRoJS06jW6/CckllHtSkiQkJHyyT5mMXPYhy3LnRORen/f3k5L7nLi9blxeF7cfdzvDY4+utEivN/K+Ukyrq6hqbuellUW8v7YMt9fHqSOSuGFKjl/rW3s9HnJXr2Ttl5/QWFFGiNPNiOMncvxf7kAfFLjZKg1tTl5fVcxbv5RidXiYNjCOW2f080vJhL5A86efUX3PPaS98Tqh48f/7nmf10vhhl/ZvOg7yrZtRqPVkT1uPFlTT0KfEk2rq5U2VxttbmWxuW1YXVbaPe2/W/YYrcvrwul14vK68MgevLJXMWSf8lhGRpbl/dZaSYskSWjQIEnSfuavk5T1vl8Qex7fPOpmhsUOO6pr0+uNHPpWMa2uos7q4NWfinl7TSl2l5fJ/WO5+sRMJuXE+K1VLPt8rP3TdWyvKqM5WE9IeAQjZ81j5Mx5GE2RftFwOJRb7Ly0soiP1pfj8vqYPSSB66dkCwPvYnxOJwVTphIyatR+fSVtlka2LV3E1iU/0GZpJDw6luEzZjNs2syAHLzTHfQJI/e22SiYNo3Q8eNJefqpLj9+b6bJ5uLdX0t585dS6q1OBsSHc9WJmZw2Iolgffemy7lraiiYPgPzpZfiOXUu6775jKINa9HpDQw+aRpjTjmdqKSUbtVwMPZMmv3m6hIW7axBq5E4c1QK107O8msGUF9jT/ZS5nffUmmpZ9vShRRtWIcs+8gYMZoRM+eRNWosGm3PSeXsCvqEkQPUPfEkjS+9RNZ33xKUmdkt5+jNOD1evt5SzSs/FbG7xorZqOfsMSlccHwaWd1kXHVPPEnjyy+TvWhRZ8phY2U5G7/9kh0rl+B1u8kcNZZRs08lY/gov9RFt7s8fLGpijdXl5BbayXSqOe841K5fEJGn8tCUYOGnTtYdeMCKpNiaXe7MJoiGTplBsOmzSIyoetyzHsafcbIPQ0NFEybjmn+fBL/9c9uOUdfQJZlVhU08u6vpfy4sxaPT2Z8VhQXjktn5uD4Lmul+xwO5Wf02DGkPvPM7563tzSzedF3bF38PbbmJsyJyYycNY8hk2d0eSU6WZbZWNbMJxsq+GZLFVanh8GJEVw+IYPTRnb/L5O+jsvRTt6aVexcsYTyndsAiLU5GPfXu8iZNLnH5n53JX3GyAGq77+flk8/I3vx4l5fx9of1FkdfLy+gg/WlVFuaSc8SMesoQmcOiKJidnR6I5h1qK9HVtvEDp+3EH383rc5K1ZxaYfvqY6Pxd9cAiDJk1m+PTZR10beg9Vze18vqmSTzdUUNRgI0SvZc7QBC4Yl8bYdHOPyKDpqfi8Xsp2bGX3z8vJW7MKt9NBZEIiQ06aTlZCCo1XXUPcHXcQfeUVaksNCPqUkbvKyymcNZuoyy4j/s47uu08fQ2fT2Z1YSNfbK5k4fYarE4P0aEG5g5L5OTB8RyfGXVErVZZlik+8yzweMj86svDNsyawnw2L/yG3F9+xuNyEp+Vw/Dpsxk48SQMIYfXSi+sb2PhjhoWbq9hS0ULAMdnRHH2mBTmDk/sddOqBRKyLFNTkMeuVcvJXf0T9pZmDCFGBpwwiSGTZ5A0YFDnvVB68SW4qirJWbQISbTI+5aRA1T+9Xbali4lZ+kStJGR3XquvojD7WV5bj1fb61iya5aHG4fIXotE3OimTIgjikDYg852Mi+fj2lF19Cwj//gfncc49cg62NXT8vZ+viH2goK0EfFEy/409g8EnTSR06DM0+kynYnB7WlVj4pbCRpbvryK9rA2BEiolZQxOYNyyR9Gj/5dL3NfaYd96vq8hfu5qW2hq0ej1Zo45j4KTJZI06Dt0BJpWwLl5MxY03kfzkE0TMnq2C8sCizxm5IzeP4vnzibn5JmJvuKFbz9XXaXd5WVPUyLLcOpburqOiqR2AJFMwo9PNjOlYBiVG7Dd5dMXNt2D79Vf6LV+GJuToOxD3mMS2pQvJW7MKp91GSKSZ8MHjqIsfwmqLgS0VLXh8MnqtxNj0KGYNiWfmkIQurwMu2IvP66Uqdxf5a1eTv/YXrI31aLRa0oaOYMAJJ5Jz/AmHnO9S9nopnD0HXXQ0GR+87yflgUufM3LomIF9yxZyli7p0yVu/YksyxTW2/gpv54NpU1sLG2iqsUBgEGrITMmlJy4MIZqbEz+5wJ8515E1J//TExY0BFPiGF3eahsaqeiuZ3KpnbKLXZ2VVhozdtCSsMO0trL0eLDHhJFcP/RDJ90EieOG45RhE26DUdbG8VbNlC0YS0lmzfgsLWh1evJGDGafsdPIHvMuCOe69Xy1tvUPvggGR9+QMiIEd2kvGfQJ43cvnEjpRdeRPzddxF16aXdfj7BgalqbmdjWRPbKlooqGujoL6N2cvf55SiVVwx824aQiIBMIXoiQ0PIixIh14roddq0Gk16DUSLq+PNqcHm9ODzenF6nDT6vDsdx6DVkP/hDCGJJoYnBRBvwjQlm2jbMMvlO/Yjiz7MCcmkXP8BLJGjSWp/6A+l4fc1fi8XmoK8yjdupnSbZuoytuN7PMREh5B1ujjyBp9HBkjRh9238WB8LbZKJgyhbCTTuzV87ceDn3SyAFKLr4Yd2UVOQt/QAqAiV0F4G1rI3/yFOTxkyi/4W80tLloaHNSb1UWm8uDxyvj9vpw+2Q8Xh8GnYawIB2hBh3GIC3hQTriTcEkR4aQYg4hKTKEuPBgtAeYtBiUNMaCdWvIXfMzFTu34fN6CQ4NI33EaLJHH0f6iNE9rja1Gvh8XhrKSqnYtYPyHVsp37EVp90GkkR8ZjYZI0aTNfo4EnL679dHcazU/vdhLG+9Rc6Pi9AnJXXZcXsafdbI21aupPza60h88EEizzzDL+cU/DGNb7xB3UP/JePjjwkZNtTv53fabZRs2UTxpnUUb96AvaUZgJi0DNKGDCd16AhSBg05ZPy2L+B2OKgtKqAqfzeVu3dQuXunYtxARGwc6cNGkj58FKlDhnfrF+GeWuVRV1ze4yccORb6rJF3VlNzucj65mu/jAoUHBzZ46Fw5ix0SYlkvPOO2nKQfT5qivIp27aFsh1bqdq9E4/bhSRpiElLJ7HfAJL6DyKx3wDMicm9Oqfc43LRWFFGXUkR1QW51BTk0VBeiuxTaoebk1JIGTSElEFDSRk4hIhY/47RqLj1z9hWrSJn2TK0YX0zw+hgRt7re30kSSL6mqup+stfsf64mIhZM9WW1KexLl6Cu6qKuLv+prYUACSNhsScAUpt9DPOxeN2U52/m7LtW6nO383uVSvZuvgHAIJDw4jNyCIuI5PY9Cxi0zOJSk5F18NmpfJ6PLTU1WKpqsBSWU59aTH1pcVYqio6TTsoNJSE7P6MG3M8CTkDSMzpr3oBs+grLsf6ww80f/Ix0ZdfrqqWQKPXt8hBSWEqmjsPKdRI5qef9upWVaBTcsGFeBoayP7he6Qe0NEo+3xYqiqoyt9NdX4u9SVFNJSX4XEpU5lJGg2muHjMCUmYE5OJTEzCHJ9IeEwsYVExXV5G4LA0yzLt1lbaLI201NfSWldHS30NrfV1NFVX0VxTjc+7t6M4LDqGuPS9X05xGZlExicG5K/X0osvwVVZSc6ivjmtY59tkQNIWi3R115L9T330LZ8OeFTp6otqU/SvmUL7Zs2EX/3XT3CxEEx6uiUNKJT0hg2Vfk15/N5aaquor60mIayUpqqK2mqqaJi1w7cTsd+rzeEhBAWFUNopJmQsHCCw8MJCY8gOCycIGMoOoMBXVAQekMQOoNhnywaCUmSkGUZn8eDx+3C43bjdbtwO5247DYcNhtOuw2nzUZ7Wyu2piZsLU3Ym5v3M2oAfXAIprh4opKSyTluPFFJKUQlpWBOSu5RfQFRV11JxfU30Pr995hOO01tOQFDn2iRA8huN4Wz56CNiSbjgw9Eq1wFKm+7jbaVP5GzfHmvjHHKsoytyUJzbTVtlkaslkasjfW0NTZia26ivc2Ko2PZE8I4VgwhIQQZwwgOCyPUHEWoyUyo2UxopJmwqGhMsfFExMUTHBrWK+75vj6tY59ukQNIej3R11xDzf33Y1u9mrCJE9WW1KdwV1XRunARUZde2itNHJT+mLCoaMKiov9wP9nnw9lux9XejsflxONSWtlupwN8PmSAfRpYWr0erU6PzmBAq9ej0xsICg0lyBja5/LgJY2G6CuupPqee7CtWk3YJPE5hj5k5ACmM8+g4YUXaHj+eWHkfsbytpKhEnXxRSorUR9JoyE4NKxHhTQCiYhTT6H+ySdpfPUVYeQdBF5vRjeiMRiIvuoq2tdvwLZ2rdpy+gze1laaP/yQiNmz0Scnqy1H0MPRGAyYL70E+y9raN+xQ205AUGfMnKAyHPORhsTQ8Pzz6stpc/Q9MGH+Ox2oq+6Um0pgl6C+bzz0ISGYnntdbWlBAR9zsg1wcFEX3kl9l/WYN+0SW05vR6f04nl7bcInTiR4MGD1ZYj6CVoIyKIPPdcWn/4AVdFpdpyVKfPGTmA+fzz0EZGila5H2j56iu89Q1EX32V2lIEvYyoSy8BScLyxhtqS1GdPmnkGqORqMsvx7byJ9q3bVNbTq9F9vmwvPoawYMHYxw/Xm05gl6GPjER07x5NH/yCR6LRW05qtInjRzAfPFFaE0m6p/+n9pSei3WJUtwlZQQfc3VfS7fV+Afoq+7FtnpxPLGm2pLUZU+a+TasDCirr4K208/Yd8oYuVdjSzLNL7yCvrUVMJPPlltOYJeSlBWFuGzZtH07rt4W1rUlqMafdbIAaIuughtdDT1Tz+ttpReR/uGDTi2bCXqisvFpLmCbiVmwXX4bDYs776rthTV6NNGrjEaib7mauxr1mD7VeSVdyWNr7yKNiqKyDPPVFuKoJcTPHAgYVOn0vTmW3jbbGrLUYU+beQA5vPPRxcXR/3TT6NG3ZneiDM/n7blyzFffBGa4GC15Qj6ADELrsPb0kLzhx+oLUUV+ryRa4KDib7uWto3bMC2arXacnoFDS+8iMZoxHzBBWpLEfQRQkaMIHTCCTS+/gY+h+PQL+hl9HkjB4g85xx0iYmiVd4FOIuKaP3uO8wXXYTObFZbjqAPEb1gAd6GBpo/+VRtKX7nmIxckqRzJEnaIUmST5Kk35VW7CloDAZirl+AY+tW2pYvV1tOj6bxxReRgoOJuuJytaUI+hjG444jZPRoGl95BdnlUluOXznWFvl24ExgZRdoUZXIM85An5pK/f/+12W1ovsarpISWr7+BvMFF6CLilJbjqCPIUkSMdcvwFNTQ/MXX6gtx68ck5HLsrxLluXcrhKjJpJeT8yfbsC5cxfWH35QW06PpOHFl5S671deobYUQR8ldNIkgocOpfGFF/H1oVa532LkkiRdK0nSekmS1tfX1/vrtEeE6dRTCRowgLrHn+hTN0FX4Covp+WrrzCffx66mBi15Qj6KJIkEXvLzbirqmj+5BO15fiNQxq5JEmLJUnafoBl/pGcSJbll2RZHivL8tjY2NijV9yNSFotcX/9C+6KCpo/6JtpTEdL40svIWm1RF0pimMJ1CV00iRCxoyh8fkX8LW3qy3HLxzSyGVZniHL8tADLF/6Q6C/CZ00CeMJ42l47nm8VqvacnoE7spKmj//gshzz0UfH6e2HEEfR5Ik4m69BU99PU3vva+2HL8g0g9/gyRJxP31r3ibm2l8+RW15fQIGl5+GUmSRKlaQcBgPO44QidOpPHll/G2taktp9s51vTDMyRJqgBOAL6VJGlh18hSl5AhQ4g49VQsb76Ju6ZGbTkBjbuqipZPP8N09lnoExLUliMQdBJ76y14m5uxvNn7KyMea9bK57Isp8iyHCTLcrwsy7O6SpjaxN5yC/h8osztIah/+n8gScRce63aUgSC/QgZNoywGdOxvP4G3uZmteV0KyK0chAMKcmYL7qIli++wJGbp7acgMSRm0fLl19ivuRi9ImJassRCH5H7M0347PZaHz1NbWldCvCyP+AmAXXoQkLo+6xR9WWEpDUP/44mvBwYq65Rm0pAsEBCe7fn4h587C88w6eAE177gqEkf8B2shIYhYswLbyJ6zLlqktJ6CwrV1L24oVxFx7DdrISLXlCAQHJfbGPyG7XNQ/95zaUroNYeSHIOriizBkZVH7n4fwOZ1qywkIZFmm7rHH0MXHY774YrXlCAR/iCEjA/P559P84Uc48npnmFQY+SGQDAbi77kbd1kZltdfV1tOQGD98UccW7YSe/NNot64oEcQc+OflDDpQ//tlRVOhZEfBmETJxI+cyYNL7yIu6pKbTmqIns81D/xJIacbEzzj2hwr0CgGjqzmdgb/4Rt9WraVqxQW06XI4z8MIm/8w4Aav/7sMpK1KX5s89wFRcTd9ttYi5OQY/CfMEFGDIylFa52622nC5FGPlhok9OJua6a7EuXIhtdd+cScjbZqPhf88QMno0YVOnqi1HIDgiJL2euDvvwFVSQtP7vWvovjDyIyDqyivRp6ZS8+8H+lzheoCGZ57B09BA/J13IEmS2nIEgiMmbMoUQidMoP7Z5/A0Naktp8sQRn4EaIKCiL/7LlxFRVjefkdtOX7FkZuH5e23iTznHEJGjFBbjkBwVEiSRNzf7sRntdLwbO9JRxRGfoSET51K2NSp1D/zDK6yMrXl+AVZlqn55z/RhocT++db1ZYjEBwTwf37E3neuTS9/z7O/Hy15XQJwsiPgoT7/o6k1VJ9z719Ylq4li++pH3DBuL++hcxobKgVxB7001ow8Ko/r+/94rPsDDyo0CfkED8XX/Dvm4dTb18AgpvSwt1jzxCyMiRmM48U205AkGXoIuKIu6uv9G+eXOvqFkujPwoMZ15JqGTJlH36GO4KirUltNt1D35JN7mZuVXiEbcLoLeg2n+fEInTqT+8cd7/PgQ8ck8SiRJIvGf/0CSJKrv/b9eOVqsfdt2mj/4EPOFFxI8aJDacgSCLkWSJBL+8Q9kWaa6Y91TEUZ+DOiTkoi74w7sa9bQ/OFHasvpUmSXi5r77kMbHU3sLTerLUcg6BYMKcnE3XoLthUraf32O7XlHDXCyI+RyHPPwXjCeOoefhh3ZaXacrqM+meexbFzJwl//z+04eFqyxEIug3zxRcTPHw4tQ880GNzy4WRHyOSJJH4r38DUHnnncgej8qKjh37unU0vvwyprPOJGLmTLXlCATdiqTVkvivf+G1Wql76CG15RwVwsi7AENKMgn/+Aft6zdQ/+STass5JrxWK5V33ok+NZWEu+9WW45A4BeCB/Qn+pqrafnyK1p/6HlTDwsj7yJMp55C5Pnn0fjKq1iXLFFbzlFT889/4amtI/mRh9GEhqotRyDwG7HXX0/w8OFU33tvjxvsJ4y8C4m/6y6Chwyh6m934SovV1vOEdPyzbe0fv01MTdcL4bhC/ocksFA8uOPg0ZD5a1/xteD6ikJI+9CNEFBJD/1JEgSlbfc2qNmFHJXVlLzj38QMmoUMdddp7YcgUAVDCnJJP3nQRw7d1LXg0pWCyPvYgwpKSQ99BCOnTupffA/ass5LHzt7VTc+mfwekl6+L+izrigTxM+fTpRl11G07vv9ph4uTDybiB82lSir7ma5g8/pCnA88tlr5eqO+7AsX07SY88jCE1VW1JAoHqxP3lth4VLxdG3k3E3nILoSedSM0//hHQ3+p1jzyK9cfFxP/tTsKnT1dbjkAQEOwbL6+45Va8bTa1Jf0hwsi7CUmnI+WppwgZOZLK22+nbdUqtSX9Dst772F54w3MF12E+dJL1ZYjEAQUhpRkkh95GGdeHhU33RjQnZ/CyLsRTUgIqS88T1B2NhU33oR90ya1JXXStmIFtf9+gLApU4i/+y4x449AcADCJk8m8YF/Y/9lDVV/vR3Z61Vb0gERRt7NaCMiSHv5JXSxsZQvuB5HXp7akmjfto3KP99G0MABJD/2KJJWq7YkgSBgiTz9dOLv+hvWRYuouT8wi2sJI/cDuthY0l57FU1QEOVXXY0jVz0zb/vpJ0ovuxxtZCSpz78gBv0IBIdB1GWXEX3ddTR//DH1Tz6ltpzfIYzcTxhSUkh79RUASi+4AOvSpX7X0PzZ55QvuB5DejrpH7yPPj7O7xoEgp5K7K23EHnuuTS++CINL70cUC1zYeR+JKhfPzI++RhDVhYVf7qRhhdf8svNIMsy9c89R/XddxM6bhzpb7+FPk6YuEBwJEiSRMJ9fydi7hzqH3+cmr//HTlAOkCFkfsZfXw86e+8TcScOdQ/8QRVt9+Bz+HotvMpdcXvp+Hp/2GafxqpLzyPNiys284nEPRmJK2WpEcfJXrBdTR//AmlV1yJp7FRbVnCyNVAExxM0mOPEnvrrbR+8w2lF15E+9atXX4e25pfKTrjTJo/+ojoa68l8aGHkAyGLj+PQNCXkDQa4m69laTHHsWxYwfFZ5+DY9cuVTUJI1cJSZKIWXAdKc8+g7u2lpJzz6Pytr/gqjj2ySncdXVU/uWvlF1+ObLTScrzzxF3259FiqFA0IWY5s0j/Z13QJYpueBCmj78SLX5CCQ1AvZjx46V169f7/fzBireNhuNr76C5fU3wOvFfMklxFx3LVqT6YiO47FYaPniSxqefRbZ5SL6mmuIvvYaNMHB3SNcIBDgqa+n8ra/YF+3DkN2NnF/vpWw6dO7peEkSdIGWZbH/m67MPLAwV1TQ/1TT9PyxRdIOh0hY8YQOmECoRMnEDxo0AFnsXdVVNK2ZDHWHxdj37gRfD5CTzyRhHvvwZCersK7EAj6HrIs07ZkCXWPP4GrqIiQkSOJ++tfMI79neceE8LIexCO3btp+fIrbKtW4ewYQKQ1m9GnpYLbg+x2I3s8+JwOPFXVAAT170/4jOmEz5hB0KBBIowiEKiA7PHQ/PnnNPzvGTx1dejT0zAedxyhxx+P8bjj0CcmHtPxhZH3UDz19dh++QXbqtV4GhqUErN6HZJOj6TXEzxwIOEzpovWt0AQQPja22n+9DNsq1Zh37ABX2srAPqUFBL//S9Cx48/quN2i5FLkvQIcCrgAgqBK2RZbj7U64SRCwSCvoLs9eLMy8O+bh32deuI/fNtBGVlHtWxusvIZwJLZVn2SJL0XwBZlu881OuEkQsEAsGRczAjP6b0Q1mWF8myvCffZg2QcizHEwgEAsGR05V55FcC3x/sSUmSrpUkab0kSevr6+u78LQCgUDQtznk5IySJC0GEg7w1D2yLH/Zsc89gAd492DHkWX5JeAlUEIrR6VWIBAIBL/jkEYuy/KMP3pekqTLgVOA6XIglQMTCASCPsIxTZcuSdJs4A5gsizL9q6RJBAIBIIj4Vhj5M8A4cCPkiRtliTphS7QJBAIBIIj4Jha5LIs53SVEIFAIBAcHaL6oUAgEPRwVBmiL0lSPVDq9xMfHjFAg9oiDkGgaxT6jo1A1weBr7G36kuXZTn2txtVMfJARpKk9QcaORVIBLpGoe/YCHR9EPga+5o+EVoRCASCHo4wcoFAIOjhCCP/PS+pLeAwCHSNQt+xEej6IPA19il9IkYuEAgEPRzRIhcIBIIejjBygUAg6OH0KSOXJGm2JEm5kiQVSJL0twM8f5skSTslSdoqSdISSZLS93nO21GGYLMkSV+ppO9ySZLq99Fx9T7PXSZJUn7HcplK+p7YR1ueJEnN+zznj+v3miRJdZIkbT/I85IkSU936N8qSdLofZ7zx/U7lL6LOnRtkyRptSRJI/Z5rqRj+2ZJkrptVpbD0DhFkqSWff6Wf9/nuT+8P/yk7/Z9tG3vuO+iOp7r9msoSVKqJEnLOnxkhyRJtxxgn66/D2VZ7hMLoEWZji4LMABbgMG/2WcqYOx4fD3w4T7PtQWAvsuBZw7w2iigqGNt7nhs9re+3+x/E/Cav65fxzlOAkYD2w/y/FyUmvkSMB741V/X7zD1TdhzXmDOHn0d/y8BYgLgGk4BvjnW+6O79P1m31NRZjDz2zUEEoHRHY/DgbwDfI67/D7sSy3y44ECWZaLZFl2AR8A8/fdQZblZfLeKo7+nvHokPr+gFnAj7IsW2RZbgJ+BGarrO8C4P0u1vCHyLK8ErD8wS7zgbdkhTVApCRJifjn+h1SnyzLqzvODyrNuHUY1/BgHMv9e9gcoT417sFqWZY3djy2AruA5N/s1uX3YV8y8mSgfJ//V/D7C7wvV7H/jEfBHTMcrZEk6XQV9Z3V8XPsE0mSUo/wtf7QR0dIKhNYus/m7r5+h8PB3oM/rt+R8tv7TwYWSZK0QZKka1XStIcTJEnaIknS95IkDenYFlDXUJIkI4oJfrrPZr9eQ0mSMoBRwK+/earL78Njqn7YW5Ek6WJgLDB5n83psixXSpKUBSyVJGmbLMuFfpb2NfC+LMtOSZKuA94EpvlZw+FwPvCJLMvefbYFwvXrEUiSNBXFyCfts3lSx/WLQykbvbujdepvNqL8LdskSZoLfAH0U0HHoTgVWCXL8r6td79dQ0mSwlC+RG6VZbm1O86xL32pRV4JpO7z/5SObfshSdIM4B7gNFmWnXu2y7Jc2bEuApajfNP6VZ8sy437aHoFGHO4r/WHvn04n9/8pPXD9TscDvYe/HH9DgtJkoaj/G3ny7LcuGf7PtevDvgcJZThd2RZbpVlua3j8XeAXpKkGALoGnbwR/dgt15DSZL0KCb+rizLnx1gl66/D7sz8B9IC8qvjyKUn/x7OmOG/GafUSgdNv1+s90MBHU8jgHy6eKOnMPUl7jP4zOANfLeTpLiDp3mjsdR/tbXsd9AlE4lyZ/Xb59zZXDwjrp57N/JtNZf1+8w9aUBBcCE32wPBcL3ebwamN0d+g5DY8Kevy2KEZZ1XM/Duj+6W1/H8yaUOHqov69hx7V4C3jyD/bp8vuwW26EQF1QeovzUMz6no5t/0RpfQMsBmqBzR3LVx3bJwDbOm7ObcBVKun7D7CjQ8cyYOA+r72ywwQKgCvU0Nfx//uBh37zOn9dv/eBasCNEl+8ClgALOh4XgKe7dC/DRjr5+t3KH2vAE373H/rO7ZndVy7LR1//3u68TNyKI037nMPrmGfL50D3R/+1texz+XAB795nV+uIUo4TAa27vN3nNvd96EYoi8QCAQ9nL4UIxcIBIJeiTBygUAg6OEIIxcIBIIejjBygUAg6OEIIxcIBIIejjBygUAg6OEIIxcIBIIezv8DhtPCOU5T/c0AAAAASUVORK5CYII=",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(result.t, np.transpose(result.y))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "MinerU",
"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.10.16"
}
},
"nbformat": 4,
"nbformat_minor": 2
}