-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathsalpeter.py
More file actions
313 lines (313 loc) · 62.9 KB
/
salpeter.py
File metadata and controls
313 lines (313 loc) · 62.9 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
{
"cells": [
{
"cell_type": "markdown",
"id": "5f6bfe51",
"metadata": {},
"source": [
"## Random Numbers"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "7fdedd38",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from scipy.integrate import quad\n",
"import matplotlib.pyplot as plt\n",
"import time "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1ef7e344",
"metadata": {},
"outputs": [],
"source": [
"def randomgen(x,n,a= float(125648),c =float(158756),m =(245018664)):\n",
" output= []\n",
" for i in range(n):\n",
" x = (a*x + c)%m\n",
" output.append(x)\n",
" return output\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f468ed0f",
"metadata": {},
"outputs": [],
"source": [
"print(randomgen(2,10))"
]
},
{
"cell_type": "markdown",
"id": "f87a2272",
"metadata": {},
"source": [
"### Tried Implementing in a different way"
]
},
{
"cell_type": "code",
"execution_count": 212,
"id": "f857792f",
"metadata": {},
"outputs": [],
"source": [
"def pstar(x,a=2.3):\n",
" f = x**(-a)\n",
"# print('function',f)\n",
" return f\n",
"\n",
"def p(x):\n",
" z = quad(pstar,0.1,150)\n",
"# print('what is z',z)\n",
" return pstar(x)/z[0] "
]
},
{
"cell_type": "code",
"execution_count": 213,
"id": "e37c4dbe",
"metadata": {},
"outputs": [],
"source": [
"# x=np.linspace(0.1,150,10)\n",
"# z = quad(pstar,0.1,150)\n",
"# print(z[0])\n",
"# plt.loglog(x,p(x))\n",
"# print(pstar(x)/z[0])"
]
},
{
"cell_type": "code",
"execution_count": 214,
"id": "73f40466",
"metadata": {},
"outputs": [],
"source": [
"def cdf(Niteration = 10):\n",
" x = np.linspace(0.1,10000,Niteration)\n",
" f =np.cumsum(p(x))\n",
" f/= np.max(f)\n",
" return x,f\n",
"\n",
"#plt.plot(cdf(100)[0],cdf(100)[1])\n",
"\n",
"\n",
"def inverse_sampling(Niteration = 10):\n",
" sample =[]\n",
" a,b = cdf(Niteration)\n",
" print(\"length of a is:\",len(a))\n",
" for i in range(len(a)):\n",
" u = np.random.uniform(0,1)\n",
" z = a[np.argmin(np.abs(cdf(Niteration)[1]-u))]\n",
" sample.append(z)\n",
" return sample\n",
"\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 221,
"id": "cfae7053",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"length of a is: 100000\n"
]
},
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7f32a6a22ee0>"
]
},
"execution_count": 221,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGhCAYAAACphlRxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABHcUlEQVR4nO3dd3gU5d7/8femN0IISOhIEemgoSOGGgxFQFFUelEwICJWxAYHRQ+C4EMRRASsURSOIC0ovRMJghQBQRCCkRICIZCQzO+POeZ3IiXZsMnsbj6v68p1PTs7O/PdOCd8npn7vr82wzAMRERERFyEh9UFiIiIiNhD4UVERERcisKLiIiIuBSFFxEREXEpCi8iIiLiUhReRERExKUovIiIiIhL8bK6AEfLzMzk5MmTFClSBJvNZnU5IiIikguGYXDhwgXKlCmDh8fN7624XXg5efIk5cuXt7oMERERyYPjx49Trly5m+7jduGlSJEigPnlg4ODLa5GREREciM5OZny5ctn/Tt+M24XXv5+VBQcHKzwIiIi4mJyM+RDA3ZFRETEpSi8iIiIiEtReBERERGX4nZjXkREBDIyMkhPT7e6DJEs3t7eeHp6OuRYCi8iIm7EMAxOnTpFUlKS1aWIXCMkJIRSpUrd8jpsCi8iIm7k7+BSsmRJAgICtFinOAXDMLh06RKJiYkAlC5d+paOp/AiIuImMjIysoJL8eLFrS5HJBt/f38AEhMTKVmy5C09QtKAXRERN/H3GJeAgACLKxG5vr+vzVsdj6XwIiLiZvSoSJyVo65NhRcRERFxKRrzIiJSCJxISuVcSlqBna9YoA9lQ/xzvX/Lli2pX78+kydPzr+iJJs33niDRYsWER8fD0C/fv1ISkpi0aJFltaVGwovIiJu7kRSKm0nriU1PaPAzunv7cmqZyNyHWC+/fZbvL2987kquZkpU6ZgGIbVZeSKwktuZWbC3kVQswt4OGaRHRGRgnAuJY3U9Awm96hP1ZJB+X6+Q4kXGRETz7mUtFyHl9DQ0HyuKmcZGRnYbDY8PArniIqiRYtaXUKuFc7/Qnmxcz4s6A8fRcKpPVZXIyJit6olg6hdtmi+/+QlILVs2ZIRI0Zkvb799tt56623GDBgAEWKFKFChQrMmjUr6/2mTZvy0ksvZTvGX3/9hbe3N6tXrwYgLS2NF154gbJlyxIYGEjjxo1Zs2ZN1v5z584lJCSEJUuWULNmTXx9ffn9999Zs2YNjRo1IjAwkJCQEJo3b87vv/+e9bnFixcTHh6On58flStXZsyYMVy9evWG3+1mxzt8+DBdunQhLCyMoKAgGjZsyKpVq7J9/vbbb2fcuHH06dOHoKAgKlasyH/+8x/++usvunTpQlBQEHXq1GHHjh3XfLdFixZRrVo1/Pz8aNeuHcePH79hnf369aNr167Z/psMHz6cF154gdDQUEqVKsUbb7yR7TP79+/nnnvuwc/Pj5o1a7Jq1SpsNlu+P3pSeMktDy/wKQIndsCsCFj1BqSnWl2ViIjbmjhxIg0aNGDnzp1ER0fz5JNPsn//fgB69uzJF198ke0xR0xMDGFhYURERADQv39/Nm7cyJdffsnPP//MQw89xH333cfBgwezPnPp0iXGjx/P7Nmz+eWXXwgNDaVr165ERETw888/s3nzZp544omsWTIrVqygV69eDB8+nL179zJz5kzmzp3Lm2++ed3vcPXq1Zse7+LFi3To0IFVq1axc+dO2rdvT+fOnTl27Fi247z33ns0b96cnTt30rFjR3r37k2fPn3o1asXP/30E1WrVqVPnz7Zfh+XLl3izTffZN68eWzcuJHk5GQeeeQRu/4bzJs3j8DAQLZu3cq///1vxo4dS2xsLACZmZl07dqVgIAAtm7dyqxZsxg9erRdx88zw82cP3/eAIzz58/nw8FPGMaXPQ3j9WDzZ3Jdwzj0g+PPIyKSB6mpqcbevXuN1NTUbNt3/5FkVHxxibH7j6QCqSMv54uIiDCefvrprNcVK1Y0evXqlfU6MzPTKFmypDFjxgzDMAwjMTHR8PLyMtatW5e1T9OmTY3nn3/eMAzDOHTokGGz2YwTJ05kO0+bNm2MUaNGGYZhGB9//LEBGPHx8VnvnzlzxgCMNWvWXLfOFi1aGG+99Va2bZ988olRunTp6+6f0/Gup2bNmsb//d//Zb3+5+8iISHBAIxXX301a9vmzZsNwEhISMj23bZs2ZK1z759+wzA2Lp1q2EYhvH6668b9erVy3q/b9++RpcuXbJeR0REGPfcc0+22ho2bGi8+OKLhmEYxrJlywwvL6+scxqGYcTGxhqAsXDhwut+txtdo4Zh37/fuvNij+Ay0ONTeORzKFIGzh2FT7rBt09AymmrqxMRcSt169bN+r9tNhulSpXKWl7+tttuo127dnz22WcAHDlyhM2bN9OzZ08AfvrpJwzDoFq1agQFBWX9rF27lsOHD2cd18fHJ9t5QkND6devX9YdkClTppCQkJD1flxcHGPHjs12zMcff5yEhAQuXbp0zXfI6XgpKSm88MIL1KxZk5CQEIKCgti/f/81d17+t8awsDAA6tSpc822v38/AF5eXjRo0CDrdfXq1QkJCWHfvn03/qX/w/+eF8xl/f8+x4EDByhfvjylSpXKer9Ro0a5PvatUHjJi+odYdg2aDwEsMHPMTC1Aez8DFxkpLaIiLP75+wjm81GZmZm1uuePXuyYMEC0tPT+fzzz6lVqxb16tUDzEcanp6exMXFER8fn/Wzb98+pkyZknUMf3//axZO+/jjj9m8eTPNmjUjJiaGatWqsWXLlqzjjhkzJtsxd+/ezcGDB/Hz87vu97jZ8Z5//nm++eYb3nzzTdavX098fDx16tQhLS37tPb//V38Xe/1tv3v7+d/t+e07UZu9t/AMAzLFkRUeMkr3yIQ9Q4M+gHC6kDqOfhPNMzrDKcPWV2diIjb69q1K5cvX2b58uV8/vnn9OrVK+u9u+66i4yMDBITE6latWq2n/+9U3Ajd911F6NGjWLTpk3Url2bzz//HIC7776bAwcOXHPMqlWr3nSW0o2Ot379evr160e3bt2oU6cOpUqV4ujRo7f2i/mvq1evZhvEe+DAAZKSkqhevbpDjl+9enWOHTvGn3/+mbVt+/btDjl2ThReblW5cHhiNbQdA17+cHQ9zGgGayfA1YJbEEpEpLAJDAykS5cuvPrqq+zbt4/HHnss671q1arRs2dP+vTpw7fffsuRI0fYvn0777zzDkuXLr3hMY8cOcKoUaPYvHkzv//+OytXruTXX3+lRo0aALz22mvMnz+fN954g19++YV9+/YRExPDK6+8kqfjVa1alW+//Zb4+Hh27drFY489ds3dk7zy9vbmqaeeYuvWrfz000/079+fJk2aOOzRTrt27ahSpQp9+/bl559/ZuPGjVkDdvP7jozWeXEET2+4Z4S5Bsz3I+Hwj7B6HOxZAJ2nQIUmVlcoIsKhxItudR4wHx117NiRe++9lwoVKmR77+OPP2bcuHE8++yznDhxguLFi9O0aVM6dOhww+MFBASwf/9+5s2bx5kzZyhdujTDhg1j8ODBALRv354lS5YwduxY/v3vf+Pt7U316tUZNGhQno733nvvMWDAAJo1a0aJEiV48cUXSU5OdsjvJiAggBdffJHHHnuMP/74g3vuuYc5c+Y45NgAnp6eLFq0iEGDBtGwYUMqV67MhAkT6Ny58w0foTmKzTDca5BGcnIyRYsW5fz58wQHBxd8AYYBuxfA8pfg0n8H8TYYAG1eB/+Qgq9HRAqNy5cvc+TIESpVqpTtHw9XWGFXHGvu3LmMGDGCpKSkAj3vxo0bueeeezh06BBVqlS55v0bXaNg37/fuvPiaDYb1H0IqraB2Fdh56ewYw7s/94cI1Ozq7mPiEgBKRviz6pnI5y6t5G4poULFxIUFMQdd9zBoUOHePrpp2nevPl1g4sjKbzkl4BQ6DIN6j4CS0bAmUPwdT+odh90eBdCyltdoYgUImVD/BUmxOEuXLjACy+8wPHjxylRogRt27Zl4sSJ+X5ePTYqCOmXYcMkWD8JMtPBOxBavwKNB6tPkog4zM1uyYs4A0c9NtJso4Lg7QetXoYhG6B8E0hPgRWj4MPWkLDL6upERERcisJLQSpZHfovg06TwbcoJMTDrFawYjSkpVhdnYiIiEtQeCloHh7QoL+5Qm+tbmBkwOapMK0JHIy1ujoRERGnp/BilSKl4KG58NhXULQ8nD8Gn3WHBQPgYmKOHxcRESmsFF6sVq09RG+BpsPA5gF7vjH7JMXNAwetsigiIuJOFF6cgW8QtH8THv8RSteDy+dh8XCY2xH+OmB1dSIiIk5F4cWZlLkLBv0I7d8C7wA4tglmNIfV4+HqFaurExGx3MGDBwkLCyMgIICNGzdaXY5YROHF2Xh6QdOhMHQr3NHeXBdm7dtmiDm6werqREQsc/LkSSIjI7nnnnsYOHAgnTp1Yvfu3Tl+7uzZszz11FPceeedBAQEUKFCBYYPH8758+dv+rkZM2ZQt25dgoODCQ4OpmnTpixbtsxRX0dugVbYdVYhFeCxGPhlISx7Ec4cNB8j3dUb2o01V/AVESkkzp07lxVc5s6di6enJ0WKFKF9+/Zs2LCBypUr3/CzJ0+e5OTJk7z77rvUrFmT33//nSFDhnDy5EkWLFhww8+VK1eOt99+m6pVqwIwb948unTpws6dO6lVq5bDv6PYwXAz58+fNwDj/PnzVpfiOJfOGcZ3TxvG68Hmz7+rGMbPXxtGZqbVlYmIE0lNTTX27t1rpKamWl2KXRITE42wsDDjzTffzNq2ZcsWw9vb21ixYoWRkpJiNGnSxBg8eLCRkZGR7bPjxo0zqlSpYiQkJNh1zq+++srw8fEx0tPT7fpcsWLFjNmzZ9v1Gfn/bnaN2vPvt+68uAL/EOg8Ger2gMVPw+kD8M1AiP8cOk2CYrdbXKCIOC3DgPRL1pzbOyBXjWhvu+025syZQ9euXYmMjKR69er06tWL6OhoIiMjAdi8efN1Pzt69GhGjx5td2l/L0Hv5ZW7fwYzMjL4+uuvSUlJoWnTpnafTxxL4cWVVGwKQ9bDximwbgIc/sFc3K7VKGgy1BwvIyLyv9IvwVtlrDn3yyfBJzBXu3bo0IHHH3+cnj170rBhQ/z8/Hj77bfzpawzZ87wr3/9i8GDB+e47+7du2natCmXL18mKCiIhQsXUrNmzXypS3JPA3ZdjZcvRLwAT26G21vA1VSIfQ0+bAkn4qyuTkQkz959912uXr3KV199xWeffZan5pJvvfUWQUFBWT/Hjh3L9n5ycjIdO3akZs2avP766zke78477yQ+Pp4tW7bw5JNP0rdvX/bu3Wt3XeJY6irtygwD4j+Dla9A6jlzkbtGT5gdq32LWF2diBSw63bsdYHHRn/75ZdfaNCgAenp6SxcuJDOnTvbfcqzZ89y9uzZrNe333571qOhCxcu0L59ewICAliyZEmewlHbtm2pUqUKM2fOtPuz4riu0nrOYKcTSamcS0nLcb9igT6UDfHP32JsNrirlzmlesXLsPsr2PoB7FsMHd6F6h3y9/wi4vxstlw/urFSWloaPXv2pEePHlSvXp2BAweye/duwsLC7DpOaGgooaHXzsZMTk6mffv2+Pr68t133+UpuAAYhsGVK1p3y2oKL3Y4kZRK24lrSU3PyHFff29PVj0bkf8BBiDoNnjwQ6j3CCx5BpJ+hy8fhRr3Q9S/Ibh0/tcgInILRo8ezfnz53n//fcJCgpi2bJlDBw4kCVLltzysS9cuEBkZCSXLl3i008/JTk5meTkZMAcLOzp6QlAmzZt6NatG8OGDQPg5ZdfJioqivLly3PhwgW+/PJL1qxZw/Lly2+5Jrk1Ci92OJeSRmp6BpN71KdqyaAb7nco8SIjYuI5l5JWMOHlb1XbmH2S1r4Dm/4P9n0Hv62Btq9D+ACzo7WIiJNZs2YNkydPZvXq1VmPCz755BPq1q3LjBkzePLJJ2/p+HFxcWzduhUga82Wvx05coTbb78dgMOHD3P69Oms9/7880969+5NQkICRYsWpW7duixfvpx27drdUj1y6xRe8qBqySBqly1qdRnX5xMA7cZAne7mtOoTcfD9s7ArBjpPgTCNkhcR59KyZUvS09OzbatQoQJJSUkOO35uhncePXo02+uPPvrIIecXx9P/K+6uStWBgbEQNQF8guCPbTCzBfwwFtJTra5OREQkzxRe3JmHJzR+AoZugzs7QuZVWD8RZjSD39ZaXZ2IiEieKLwUBkXLwqOfQ49PoUhpOPsbzL8fFj4JKWesrk5ERMQuCi+FSY3OZrfqho8DNtj1OUxtALu+NNeCEBERcQEKL4WNX1Ho+K45HqZkTUg9CwsHw/wucOaw1dWJiAO42dqj4kYcdW0qvBRW5RvC4HXQ5nXw8oMja82xMOsnQkZ6zp8XEafj7e0NwKVLFq2oK5KDv6/Nv6/VvNJU6cLM0xtajIRaXc3F7X5bY85G2r3AnFZdvpHVFYqIHTw9PQkJCSExMRGAgIAAbHYszy+SXwzD4NKlSyQmJhISEpK1MGBeKbwIhFaG3ovg569gxShI3AsfRULDgdDmNfNRk4i4hFKlSgFkBRgRZxISEpJ1jd4KhRcx2WxQrwdUbQuxr5oNH7fPhv3fmy0GanS2q8GaiFjDZrNRunRpSpYsec3CbyJW8vb2vuU7Ln9TeJHsAotD1+lQtwcsGWFOq/6qN9zZATpMgKLlrK5QRHLB09PTYf9QiDgbDdiV66scAU9ughbPgYcXHFgK0xrDlg8gM+fGlCIiIvlF4UVuzNsf2rwKQzZA+caQdhGWvwiz28Kp3VZXJyIihZTCi+SsZA3ovxw6TgLfYDj5E8yMgJWvQpqmZIqISMFSeJHc8fAwZx8N3QY1u4CRAZveh+lN4NAqq6sTEZFCROFF7BNcGh6eD49+CcHlIOl3+PRB+GYQXPzL6upERKQQUHiRvLkzyuyT1CQabB6w+2uzT9JP89UnSURE8pVThpdu3bpRrFgxunfvbnUpcjO+QXDfeBj0A5SqA5eT4LunYG4nOH3Q6upERMRNOWV4GT58OPPnz7e6DMmtsnfD42sgchx4B8DvG8w+SWvegatXrK5ORETcjFOGl1atWlGkSBGryxB7eHpBs6cgegtUbQcZabDmLfigBfy+2erqRETEjTg8vKxbt47OnTtTpkwZbDYbixYtumaf6dOnU6lSJfz8/AgPD2f9+vWOLkOsUqwi9Pwaus+BwNvg9AH4+D5Y/DSknrO6OhERcQMODy8pKSnUq1ePqVOnXvf9mJgYRowYwejRo9m5cyctWrQgKiqKY8eO5el8V65cITk5OduPWMxmg9oPwrDtcHdfc1vcXJjaCPZ8owG9IiJySxweXqKiohg3bhwPPPDAdd+fNGkSAwcOZNCgQdSoUYPJkydTvnx5ZsyYkafzjR8/nqJFi2b9lC9f/lbKF0fyLwb3vw/9lkLxOyAlERYMgM8fhqS8hVUREZECHfOSlpZGXFwckZGR2bZHRkayadOmPB1z1KhRnD9/Puvn+PHjjihVHOn25vDkRmg5Cjx94OBKs0/SpqmQcdXq6kRExMUUaHg5ffo0GRkZhIWFZdseFhbGqVOnsl63b9+ehx56iKVLl1KuXDm2b99+w2P6+voSHByc7UeckJcvtHwJhmyEis0h/RKsHA2zW8PJnVZXJyIiLsTLipPabLZsrw3DyLZtxYoVBV2SFJTbqkHfJRD/Kax8BRJ2wYetofGT0Oplc+0YERGRmyjQOy8lSpTA09Mz210WgMTExGvuxogb8/CAu/vAsB1QuzsYmbBlmtkn6VcFVxERubkCDS8+Pj6Eh4cTGxubbXtsbCzNmjUryFLEGQSVhO4fQc9vIKQCnD9uDub9qi9cOJXz50VEpFByeHi5ePEi8fHxxMfHA3DkyBHi4+OzpkKPHDmS2bNnM2fOHPbt28czzzzDsWPHGDJkiKNLEVdxR1tzcbtmT4HNE/YuMqdV75gDmZlWVyciIk7G4WNeduzYQatWrbJejxw5EoC+ffsyd+5cevTowZkzZxg7diwJCQnUrl2bpUuXUrFiRUeXIq7EJ9BsL1DnYVg83BzEu+QZ2BUDnSdDyRpWVygiIk7C4eGlZcuWGDksQhYdHU10dLSjTy3uoHRds9Hjtlnww7/g+BazxcA9I6DFc+DtZ3WFIiJiMafsbSSFnIcnNHkShm6FalGQmQ7rJpjNHo+ss7o6ERGxmMKLOK+Q8vDoF/DwfAgqBWcPw7zOsGgoXDprdXUiImIRhRdxbjYb1OwCw7ZBg4GAzVwjZmoDczyM+iSJiBQ6Ci/iGvyKQqdJMGAF3FYDLp2BhU/AJ93g7G9WVyciIgVI4UVcS4XGMHgdtH4VPH3ht9UwvSlseA8y0q2uTkRECoDCi7geLx+49zmI3gyV7oWrl2HVGzCrJfwRZ3V1IiKSzxRexHUVrwJ9voOuH4B/KPy5B2a3gaUvwOVkq6sTEZF8ovAirs1mg/qPmn2S6j0KGLBtJkxrDPuWWF2diIjkA4UXcQ+BxaHbB9B7ERSrBBdOQkxP+LInJJ+0ujoREXEghRdxL1VamWNhWjwLHl6wf4nZJ2nbh5CZYXV1IiLiAAov4n68/aHNa+aspHINIe0CLH0OPoqEU3usrk5ERG6Rwou4r7Ba5rowHd4FnyJwYgfMijBnJqWnWl2diIjkkcMbM4q1TiSlci4lLcf9igX6UDbEvwAqspiHJzR6HKp3hGUvwL7F5powvyyETu9BldZWVygiInZSeHEjJ5JSaTtxLanpOY/t8Pf2ZNWzEYUjwAAEl4Een8L+7+H75+DcUXN13ro9oP1bEFjC6gpFRCSXFF7cyLmUNFLTM5jcoz5VSwbdcL9DiRcZERPPuZS0whNe/la9o7mw3Y/jYOtM+DkGDq6EyDeh/mPm1GsREXFqCi9uqGrJIGqXLWp1Gc7LtwhEvQN1HobFT8Ofu+E/0bDrC+g0GUpUtbpCERG5CQ3YlcKrXDg8sRrajgEvfzi6HmY0g7UT4GrO44ZERMQaCi9SuHl6wz0jzLVhqrSGjCuwehzMbAHHtlhdnYiIXIfCiwhAaCXo9S08MBsCSsBf+2FOe1jyDKQmWV2diIj8D4UXkb/ZbFD3IRi2He7qZW7bMQemNTKnVhuGtfWJiAig8CJyrYBQ6DIN+i6B4lXh4p/wdT/44hFIOm51dSIihZ7Ci8iNVGoBQzZCxIvg4Q2/Lje7VW+erj5JIiIWUngRuRlvP2j1MgzZABWaQnoKrBgFs9tAwi6rqxMRKZQUXkRyo2R16LcUOk8B36JwcifMagUrRkNaitXViYgUKgovIrnl4QHh/cwBvbUeACMDNk+FaU3gYKzV1YmIFBoKLyL2KhIGD30Mj30NRSvA+WPwWXdYMAAuJlpdnYiI21N4EcmrapEwdAs0HQY2D9jzDUxtAHHzIDPT6upERNyWwovIrfAJhPZvwuOroXR9uHweFg+HuR3hrwNWVyci4pYUXkQcoUx9GPQDtH8LvAPh2CaY0RxWj4erV6yuTkTErSi8iDiKpxc0HWo+SrqjPWSmw9q3zRBzdIPV1YmIuA2FFxFHC6kAj8XAQ3MhsCScOWg+RvrPMLh01urqRERcnsKLSH6w2aBWN3NadXh/c9vOT8w+SbsXqE+SiMgt8LK6AHF+J5JSOZeSluN+xQJ9KBviXwAVuRD/EOg8Ger2gMVPw+kD8M1AiP8cOk2CYrdbXKCIiOtReJGbOpGUStuJa0lNz7mXj7+3J6uejVCAuZ6KTWHIetg4BdZNgMM/mIvbtRoFTYaa42VERCRX9BdTbupcShqp6RlM7lGfqiWDbrjfocSLjIiJ51xKmsLLjXj5QsQL5uq8S0bA0fUQ+xrs/tpsO1A23OoKRURcgsKL5ErVkkHULlvU6jLcQ4mq0HcxxH8GK1+BU7thdlto9AS0fgV8i1hdoYiIU9OAXREr2GxwVy8Yuh3qPAxGJmz9AKY1hv1Lra5ORMSpKbyIWCnoNnjwQ+j1LYRUhOQT8OWjENMbkhOsrk5ExCkpvIg4g6ptIHoLNB8BNk/Y9505rXr7bPVJEhH5B4UXEWfhEwDtxsDgtebg3SvJ8P2z8PF98Odeq6sTEXEaCi8izqZUHRgYC1ETwCcIjm+FmS3gh7GQnmp1dSIillN4EXFGHp7Q+AkYug2qd4LMq7B+IsxoBr+ttbo6ERFLKbyIOLOiZeGRz6DHp1CkNJz9DebfDwufhJQzVlcnImIJhRcRV1Cjs3kXptETgA12fQ5TG8CuL9UnSUQKHYUXEVfhFwwdJpjjYUrWgtSzsHAwzO8CZw5bXZ2ISIFReBFxNeUbmjOS2rwOXn5wZK05Fmb9RMhIt7o6EZF8p/Ai4oo8vaHFSIjeDJVbwdXL5mykmffC8W1WVycikq/U20gK3ImkVM6lpOW4X7FAHzV5zEloZei9EH7+ClaMgsS98FEkNBwIbV4DP/WjEhH3o/AiBepEUiptJ64lNT0jx339vT1Z9WyEAkxObDao1wPuaGc2eoz/zFyZd//3EPVvc7CvzWZ1lSIiDqPwIgXqXEoaqekZTO5Rn6olg26436HEi4yIiedcSprCS24FhELX6VC3BywZYU6r/qo33NnBHOhbtJzVFYqIOITCi1iiaskgapfVI418UTkCntwE696FjZPhwFI4sg5avwqNHjcXwBMRcWEasCvijrz9oc2rMGQDlG8MaRdh+Yswuy2c2m11dSIit0ThRcSdlawB/ZdDx0ngGwwnf4KZEbDyVUi7ZHV1IiJ5ovAi4u48PMzZR0O3Qc0uYGTApvdhehM4tMrq6kRE7KbwIlJYBJeGh+fDo19CcDlI+h0+fRC+GQQX/7K6OhGRXFN4ESls7oyCoVuhSTTYPGD312afpJ8+UZ8kEXEJCi8ihZFvENw3Hgb9AKXqwOUk+G4YzOsMpw9aXZ2IyE0pvIgUZmXvhsfXQOQ48A6Ao+vNPklr3oGrV6yuTkTkuhReRAo7Ty9o9hREb4Gq7SAjDda8BR+0gN83W12diMg1tEidOLVDiRdz3Ec9kBykWEXo+TX88i0sexFOH4CP74PwftD2DfAvZnWFIiKAwos4qWKBPvh7ezIiJj7HfdUDyYFsNqj9IFRpDbGvw0/zIG4u7F8KUW9DrQfUJ0lELKfwIk6pbIg/q56NyLH7tHog5RP/YnD/+1DvEVj8NJz+FRYMgF1fQseJEFLB6gpFpBBTeBGnVTbEX4HEahWbmS0GNrwH6yfCwZUwrTG0Gg2Nh5jjZURECpgG7IrIzXn5QsuXYMhGqNgc0i/BytEwuzWc3Gl1dSJSCDlleOnWrRvFihWje/fuVpciIn+7rRr0XQL3/x/4hUDCLviwNSx/Ga7kPLBaRMRRnDK8DB8+nPnz51tdhoj8k4cH3N0Hhm2H2t3ByIQt08w+Sb+usLo6ESkknDK8tGrViiJFilhdhojcSFBJ6P4R9PzGHLx7/jh8/jB81RcunLK6OhFxc3aHl3Xr1tG5c2fKlCmDzWZj0aJF1+wzffp0KlWqhJ+fH+Hh4axfv94RtYrc0KHEi+w5cf6mPyeSUq0u0/3c0dZc3K7ZcLB5wt5FMLUR7JgDmZlWVycibsruqQIpKSnUq1eP/v378+CDD17zfkxMDCNGjGD69Ok0b96cmTNnEhUVxd69e6lQwZxeGR4ezpUr1y49vnLlSsqUKWNXPVeuXMl2rOTkZDu/kbgyrQfjBHwCIfJfUOchWDzcHMS75BnYFQOdJ0PJGlZXKCJuxu7wEhUVRVRU1A3fnzRpEgMHDmTQoEEATJ48mRUrVjBjxgzGjx8PQFxcXB7Lvdb48eMZM2aMw44nrkXrwTiR0nXNRo/bPoQfxsLxLWaLgXtGQIvnwNvP6gpFxE04dJGGtLQ04uLieOmll7Jtj4yMZNOmTY48VZZRo0YxcuTIrNfJycmUL18+X84lzknrwTgRD09oMgSqd4Slz8Ovy2DdBNjzrXkXptK9VlcoIm7AoQN2T58+TUZGBmFhYdm2h4WFcepU7gfxtW/fnoceeoilS5dSrlw5tm/ffsN9fX19CQ4OzvYjIhYLKQ+PfgEPz4egUnD2MMzrDIuGwqWzVlcnIi4uX5bHtP2j94lhGNdsu5kVKzTlUsTl2WxQswtUbgmrxpiDeOM/Ne/G3Pe2OUZGfZJEJA8ceuelRIkSeHp6XnOXJTEx8Zq7MSJSSPgVhU6TYMAKuK0GXDoD3z4Onz4AZ49YXZ2IuCCHhhcfHx/Cw8OJjY3Ntj02NpZmzZo58lQi4moqNIbB66D1q+DpC4d/NBe32/AeZKRbXZ2IuBC7w8vFixeJj48nPj4egCNHjhAfH8+xY8cAGDlyJLNnz2bOnDns27ePZ555hmPHjjFkyBCHFi4iLsjLB+59DqI3m4N3r16GVW/ArJbwh+NmIYqIe7N7zMuOHTto1apV1uu/Z/r07duXuXPn0qNHD86cOcPYsWNJSEigdu3aLF26lIoVKzquapE8OpSYux48xQJ9NIMpPxWvAn2+g11fwoqX4c89MLsNNHoCWr8Cfhp4LyI3Znd4admyJYZh3HSf6OhooqOj81yUiKPZs5gdaEG7AmGzQf1H4Y5Is0v1ri9g20zYtxg6TIAanayuUEScVL7MNhJxNrldzA60oF2BCywO3T6Auj3MlXnPHYGYnlC9kxligu1bdVtE3J/CixQaWszOyVVpZY6FWTcBNk6B/Uvgt7XQ9nVoMMBcAE9EBCftKi0ihZS3P7R5zZyVVK4hpF2Apc/BR5Fwao/V1YmIk1B4ERHnE1YLBqyEDu+CbzCc2AGzIsyZSenqDi5S2Cm8iIhz8vCARo/D0K1QozNkXjXXhJnexFwjRkQKLY15EbkFJ5JSczUIWFOvb0FwGejxKexfaj5COncUPulmDvBt/xYElrC6QhEpYAovIjeQ05owZ1LSGPJJHKnpGTkeS1OvHaB6B6jUAn4cB1tnws8xcHAlRL4J9R9TnySRQkThReQf7FkTxt/bk3kDGlE80OeG+2jqtQP5FoGod6DOw7D4afhzN/wn2lwjptNkKFHV6gpFpAAovIj8gz1rwuhxkEXKhcMTq2HzNFjzNhxdDzOawb3PQ/OnzTYEIuK2FF5ErkNrwrgAT2+4ZwTU7ALfjzQH8a4eB3sWQOf3zUaQIuKWNNtIRFxbaCXo9S08MBsCSsBf+2FOpLlab2qS1dWJSD5QeBER12ezQd2HYNh2uKuXuW3HHJjWGH5ZBDn0YxMR16LwIiLuIyAUukyDvkugeFW4eAq+7gtfPAJJx62uTkQcROFFRNxPpRYwZCNEvAge3vDrcvMuzObpkJnz1HYRcW4KLyLinrz9oNXLMGQDVGgK6SmwYhTMbgMJu6yuTkRugcKLiLi3ktWh31LoPAV8i8LJnTCrFawYDWkpVlcnInmg8CIi7s/DA8L7mQN6az0ARgZsngrTmsDBWKurExE7KbyISOFRJAwe+hge+xqKVoDzx+Cz7rBgAFxMtLo6EcklhRcRKXyqRcLQLdB0GNg8YM83MLUBxM2DzEyrqxORHCi8iEjh5BMI7d+Ex1dD6fpw+TwsHg5zO8JfB6yuTkRuQuFFRAq3MvVh0A/Q/i3wDoRjm2BGc1g9Hq5esbo6EbkO9TYScSInklLVENIKnl7QdCjU6AzfPwcHV8Dat83HSZ0nw+33WF2hiPwPhReRAnIo8eJN3z+TksaQT+JITc95ETV/b09WPRuhAONoIRXgsRjYuwiWvgBnDpqPke7qDe3Gmiv4iojlFF5E8lmxQB/8vT0ZEROf477+3p7MG9CI4oE+N9znUOJFRsTEs/3IWc6VDMrx3Ao4drLZoFY3qNwKVr0BcR/Dzk/MVXrvextqP2juIyKWUXgRyWdlQ/xZ9WyEwx4H2RuGdIcmj/xDzEdGdXvA4qfh9AH4ZiDEfw6dJkGx2y0uUKTwUngRKQBlQ/wdFiByG4b+vkNzLiVN4eVWVGwKQ9bDximwbgIc/sFc3K7VKGgy1BwvIyIFSv+rE3FBjgxDkgtevhDxgrk675IRcHQ9xL4Gu7822w6UDbe6QpFCRVOlRURyq0RV6LsYukwD/2JwajfMbgvLXoQrF6yuTqTQUHgREbGHzQZ39YKh26HOw2BkwtYPYFpj2L/U6upECgWFFxGRvAi6DR78EHp9CyEVIfkEfPkoxPSG5ASrqxNxawovIiK3omobiN4CzUeAzRP2fQfTGsH22eqTJJJPFF5ERG6VTwC0GwOD15qDd68kw/fPwsf3wZ97ra5OxO0ovIiIOEqpOjAwFqImgE8QHN8KM1vAD2MhPdXq6kTchqZKi7i5nNoS/E2r8TqIhyc0fgKqd4RlL8D+JbB+IvyyEDpNhsoRVlco4vIUXkTclD0r8YJW43W4omXhkc9g32JY+jyc/Q3m3w/1HoPIcRBY3OoKRVyWwouIm7KnLYFW481HNTpDpQj48V+w7UPY9fl/+ySNN1sPqE+SiN0UXkTcmFbidRJ+wdBhgrkuzOKnIfEXWDj4v32S3oPiVayuUMSlaMCuiEhBKd/QnJHU5nXw8oMja2FGM3NMTEa61dWJuAyFFxGRguTpDS1GQvRmqNwKrl42ZyPNvBeOb7O6OhGXoPAiImKF0MrQeyF0mwUBxSFxL3wUaa4Pc/m81dWJODWFFxERq9hsUK8HDNsB9XsChrky77TGsPc7MAyrKxRxSgovIiJWCwiFrtOhz3cQWgUuJMBXveHLx+D8H1ZXJ+J0NNtIRLLkZkE7LWaXjypHwJObYP27sGEyHFgKR9ZB61eh0ePmAngiovAiIvYtaOfv7ckHvcMpHuiT4zEVcvLA2w9avwK1HzSnVR/fCstfhJ9j4P73zRYEIoWcwouI5HpBuzMpaQz5JI6+c3KeFaMVe29RyRrQfzn8NBdi34CTP8HMCGg6FFqOMptBihRSCi8iAuR+QbvchByt2OsgHh7QYADc2QGWvQh7F8Gm92Hvf6DTJKja1uoKRSyh8CIidtGqvRYoUgoengcHlsH3z0HS7/Dpg1DnIWg/HoJus7pCkQKl2UYiIq7izigYuhWaRIPNA3Z/DVMbwE+faFq1FCoKLyIirsQ3yGzqOOgHc/Du5ST4bhjM6wynD1pdnUiBUHgREXFFZe+Gx9dA5DjwDoCj680+SWvegatXrK5OJF8pvIiIuCpPL2j2FERvgartICMN1rwFH7SA3zdbXZ1IvlF4ERFxdcUqQs+vofscCLwNTh+Aj+8z14lJPWd1dSIOp/AiIuIObDZzYbth2+Huvua2uLkwtRHs+UYDesWtKLyIiLgT/2LmSrz9l0GJapCSCAsGwOcPQ9Ixq6sTcQit8yIiljqRlJrjonegdgN2q9gMhmyADe/B+olwcKXZrbrVaGg8xBwvI+KidPWKiGVOJKXSduJaUtMzctxX7QbywMsXWr4EtR6AJSPg942wcjTs/go6T4Eyd1ldoUieKLyISL7JqUv1ocSLpKZnMLlHfaqWDLrpfmo3cAtuqwZ9l0D8p7DyFUjYBR+2hsZPQquXzbVjRFyIwouIOJy9XaobVgpVKMlvHh5wdx+odh8sHwV7FsCWabDvO+jwLtx5n9UViuSawouIOFxuu1SDxrIUuKCS0P0jqPcofP+MOYj3ix5QsytEvWP2URJxcgovIpIv1MDRyd3R1lzcbs142Dzd7Fh9eDW0ewPu7mfeqRFxUro6RUQKK59As73AE2vMwbtXzsOSZ+DjKEjcZ3V1Ijek8CIiUtiVrms2erzvbfAOhONbzBYDP46D9MtWVydyDYUXEREBD09o8iQM3QrVoiAzHdZNMJs9HllndXUi2Si8iIjI/xdSHh79Ah6eD0Gl4OxhmNcZFg2FS2etrk4EcMLwcuHCBRo2bEj9+vWpU6cOH374odUliYgULjYb1OwCw7ZBg4GAzVwjZmoD2BWjPkliOaebbRQQEMDatWsJCAjg0qVL1K5dmwceeIDixYtbXZqIWCynRe9AU68dyq8odJoEdXuYHar/2gcLn4BdX5jbQytbXaEUUk4XXjw9PQkICADg8uXLZGRkYCjlixRq9i56pzYCDlahMQxeB5veh7X/ht9Ww/SmZuuBpsPA09vqCqWQsTu8rFu3jgkTJhAXF0dCQgILFy6ka9eu2faZPn06EyZMICEhgVq1ajF58mRatGiR63MkJSURERHBwYMHmTBhAiVKlLC3TBFxI7ld9E5tBPKRlw/c+xzU6mb2STqyDla9AbsXQOf3oVy41RVKIWJ3eElJSaFevXr079+fBx988Jr3Y2JiGDFiBNOnT6d58+bMnDmTqKgo9u7dS4UKFQAIDw/nypUr13x25cqVlClThpCQEHbt2sWff/7JAw88QPfu3QkLC7tuPVeuXMl2rOTkZHu/koi4AC165ySKV4E+38GuL2HFy/DnHpjdBho9Aa1fAb9gqyuUQsDu8BIVFUVUVNQN3580aRIDBw5k0KBBAEyePJkVK1YwY8YMxo8fD0BcXFyuzhUWFkbdunVZt24dDz300HX3GT9+PGPGjLHzW4iISJ7ZbFD/Ubgj0uxSvesL2DYT9i2GDhOgRierKxQ359DZRmlpacTFxREZGZlte2RkJJs2bcrVMf7888+suyfJycmsW7eOO++884b7jxo1ivPnz2f9HD9+PO9fQEREci+wOHT7AHovgmKV4MJJiOkJX/aE5JNWVyduzKEDdk+fPk1GRsY1j3jCwsI4depUro7xxx9/MHDgQAzDwDAMhg0bRt26dW+4v6+vL76+vrdUt4iI3IIqrSB6szmYd9P7sH8J/LYW2r4ODQaYC+CJOFC+zDay2WzZXhuGcc22GwkPDyc+Pj4fqhKRwiI3U6pB06odytvfDCt1upvTqv/YDkufM8fGdJ4CpWpbXaG4EYeGlxIlSuDp6XnNXZbExMQbDrgVEXEUe6ZUg6ZV54uwWjBgBeyYA6vGwIkdMCsCmj0FES+aIUfkFjk0vPj4+BAeHk5sbCzdunXL2h4bG0uXLl0ceSoRkWvkdko1aFp1vvLwhEaPQ/WOsOwFcyDvhvfgl4XQ6T2o0trqCsXF2R1eLl68yKFDh7JeHzlyhPj4eEJDQ6lQoQIjR46kd+/eNGjQgKZNmzJr1iyOHTvGkCFDHFq4iMj1aEq1EwkuAz0+hf3fw9Ln4dxR+KSbuWJv+7cgUGt4Sd7YHV527NhBq1atsl6PHDkSgL59+zJ37lx69OjBmTNnGDt2LAkJCdSuXZulS5dSsWJFx1UtIuIgajlQAKp3hEr3wo/jYOtM+DkGDq6EyDeh/mPm1GsRO9gdXlq2bJnjcv3R0dFER0fnuSgRkfymlgMFzLcIRL0DdR+G756GP3fDf6L/2ydpMpSoanWF4kKcrreRiEhBUMsBi5QNhydWw5bpsHo8HF0PM5rBvc9D86fNNgQiOVB4EZFCy57xMXq85ECe3mZQqdkFloyEwz/A6nGwZ4E5rbpCE6srFCen8CIichN6vJSPit0Ovb4xmzsufwn+2g9z2kN4f2j7BviHWFygOCuFFxGRm9DjpXxms0Hdh6BqG4h9FXZ+CnEfw4Gl5hiZml01oFeuofAiIpIDTb8uAAGh0GUa1H0EloyAM4fg635Q7T7o8C6ElLe6QnEiDm3MKCIicksqtYAhG83VeD284dflMK0xbJ4OmRlWVydOQuFFRESci7cftHoZhmyACk0hPQVWjIIPW0PCLqurEyeg8CIiIs6pZHXot9ScgeRbFBLiYVYrWDEa0lKsrk4spPAiIiLOy8MDwvvBsO1QqxsYGbB5KkxrAgdjra5OLKLwIiIizq9IGDw0Fx77GopWgPPH4LPusGAAXEy0ujopYJptJCLiQFrMLp9Vi4SKm2HNeHOV3j3fwKFV0O5fcFdv806NuD2FFxERB9BidgXINwjavwl1usPip81BvIuHw64vofNkuO1OqyuUfKbwIiLiAFrMzgJl7oJBP8K2mWbH6mObYEZzaPEstBgJXr5WVyj5ROFFRMRBtJidBTy9oOlQqNEZvn8ODq6AtW+bj5M6T4bb77G6QskHejgoIiKuL6QCPBYD3T+GwJJw5iDM7Qj/GQaXzlpdnTiYwouIiLgHmw1qP2BOqw7vb27b+QlMa2Q2fzQMa+sTh9FjIxERC2hWUj7yDzEfGdXtYfZJ+ms/fDMQ4j+HTpPMbtbi0hReREQKkGYlFaCKTWHwetg4BdZNgMM/mIvbtRoFTYaa42XEJem/nIhIAdKspALm5QMRz5ur8y4ZAUfXQ+xrsPtrs+1A2XCrK5Q8UHgRESlgmpVkgRJVoe9i89HRytFwajfMbguNnoDWr4BvEasrFDsovIiIODGNjXEgmw3u6gnV2sOKl+HnGNj6AexbDB3eheodrK5QcknhRUTECWlsTD4KLAEPzDIH9H4/Es4dhS8fhRr3Q9S/Ibi01RVKDhReRESckMbGFICqbeDJzbD2Hdj0f7DvO/htDbR5DRoMVJ8kJ6bwIiLipDQ2pgD4BEC7MVDnIbNP0okdsPQ5+Pkrc0BvWE2rK5TrUKwUEREpVRsGroSoCeBTBP7YBjNbwKoxkJ5qdXXyDwovIiIiAB6e0PgJGLoVqneCzKuwYRJMb2o+ThKnofAiIiLyv4qWhUc+gx6fQZHScO4IzO8C3w6GlNNWVycovIiIiFxfjU4wdBs0fBywwc9fwtSGEP+F+iRZTOFFRETkRvyCoeO7MDAWStaE1LOwaIh5J+bMYaurK7QUXkRERHJSviEMXgdtXgcvPziyFmY0g3XvwtWbT2cXx1N4ERERyQ1Pb2gxEqI3Q+WWcPUy/PgvmBUBx7dZXV2hovAiIiJij9DK0HsRdJsFAcUhcS98FAnfPwuXz1tdXaGg8CIiImIvmw3q9YBhO6B+T8CA7bNhWmPY+x8N6M1nWmFXRMQN5KaBI6iJo8MFhELX6WafpCXPwNnD8FUfuLMDdJgARctZXaFbUngREXFh9jRwBLOJ4we9wyke6JPjcRVy7FA5Ap7cBOvfhQ3vwYGlcGQdtH4VGj1uLoAnDqPwIiLiwnLbwBHgTEoaQz6Jo++cnAeXqlN1Hnj7QetXoNYDsGQEHN8Ky1+En2PMPkml61pdodtQeBERcXH2NHBUp+oCEFYT+i+Hn+ZC7Btw8ieY1RKaDoWWo8xmkHJLFF5ERAoRdaouIB4e0GCAOfZl2YuwdxFset8czNtpElRta3WFLk2zjURERPJLkVLw8Dx4NAaCy0HS7/Dpg7BgIFxMtLo6l6XwIiIikt/uvM/sVt1kKNg8YM8Cs0/ST/M1rToPFF5EREQKgm8Q3PcWDPoBStWBy0nw3VMwtxOcPmh1dS5F4UVERKQglb0bHl8DkePAOwB+32D2SVrzDly9YnV1LkEDdkVE5Lpys/Cd1oPJI08vaPYU1LjfbCtwKBbWvGU+Tuo8BSo2s7pCp6bwIiIi2diz8J3Wg7lFxSpCz6/hl29h2Utw+lf4OAru7gvtxoB/MasrdEoKLyIikk1uF77TejAOYrNB7QehSmuIfR1+mmf+HFgGUW+bi97ZbFZX6VQUXkRE5BpaD8YC/sXg/vf/2ydphHkXZsEA2PUldJwIIRWsrtBpaMCuiIiIM7m9OQzZYK7G6+kDB1ea3ao3TYWMq1ZX5xQUXkRERJyNly+0fAmGbISKzSH9EqwcDbNbw8mdVldnOYUXERERZ3VbNei7BO7/P/ArCgm74MPWsPxluJLzbDB3pTEvIiJyS3IzpTq3NPX6Ojw84O4+UO0+WD7KnE69ZRrs+w46vGuu3lvIKLyIiEie2DOlOrc09fomgkpC94+g3qPw/TOQdAy+6AE1u0LUO2YfpUJC4UVERPIkt1Oqc0tTr3PpjrYQvQXWvA2bp5kdqw+vhnZvwN39zDs1bk7hRURE8kxTqi3iEwiR/4I6D8Hi4eYg3iXPwK4Y6DwZStawusJ85f7xTERExF2Vrms2erzvbfAOhONb4IMW8OM4SL9sdXX5RuFFRETElXl4QpMnYehWqBYFmemwbgLMaApH1lldXb5QeBEREXEHIeXh0S/g4fkQVArO/gbzOsOioXDprNXVOZTCi4iIiLuw2aBmFxi2DRoOAmwQ/ylMbWiOhzEMqyt0CIUXERERd+NX1OyHNGAF3FYDLp2GhU/Apw+Yd2RcnMKLiIiIu6rQGAavg9avgqcvHP4RpjeFDe9BRrrV1eWZpkqLiIjLOZGUmqv1ZbRiL+DlA/c+B7W6md2qj6yDVW/Az19D5ylQvqHVFdpN4UVERFzKiaRU2k5cS2p6Ro77asXe/1G8CvT5DnZ9AStGQ+Iv8FE7aPS4eWfGL9jqCnNN4UVERFzKuZQ0UtMzmNyjPlVLBt1wP63Yex02G9R/DO5ob3ap3vUFbJsF+5ZAhwlQo5PVFeaKwouIiLikqiWDqF22qNVluKbA4tDtA6jbw1yZ99wRiOkJ1TtB1L+haFmrK7wpDdgVEREprKq0gujN0OJZ8PCC/UtgWmPYOgsyc34sZxWnDC9eXl7Ur1+f+vXrM2jQIKvLERERcV/e/tDmNXNWUrmGkHYBlj0PH0XCqT1WV3ddTvnYKCQkhPj4eKvLEBERCxxKvHhL70sehdWCASthx0ewagyc2AGzIqDZUxDxohlynIRThhcRESl8igX64O/tyYiY+Bz39ff2pFigT/4XVdh4eJizj6p3hGUvwL7F5powvyyETu9BldZWVwjkIbysW7eOCRMmEBcXR0JCAgsXLqRr167Z9pk+fToTJkwgISGBWrVqMXnyZFq0aJHrcyQnJxMeHo6/vz9vvvkmERER9pYpIiIupmyIP6uejdD6Lc4guAz0+BT2L4Wlz8G5o/BJN6jzMLR/C4Jus7Q8u8NLSkoK9erVo3///jz44IPXvB8TE8OIESOYPn06zZs3Z+bMmURFRbF3714qVKgAQHh4OFeuXLnmsytXrqRMmTIcPXqUMmXKsGfPHjp27Mju3bsJDr7+/PMrV65kO1ZycrK9X0lERJxE2RB/hRJnUr0DVGoBP46DrTNh91dwKNYMMPUfs6wsu8NLVFQUUVFRN3x/0qRJDBw4MGug7eTJk1mxYgUzZsxg/PjxAMTFxd30HGXKlAGgdu3a1KxZk19//ZUGDRpcd9/x48czZswYe7+GiIiI5IZvEYh6B+o+DN89DX/uhpM7LQ0vDp1tlJaWRlxcHJGRkdm2R0ZGsmnTplwd49y5c1l3Uv744w/27t1L5cqVb7j/qFGjOH/+fNbP8ePH8/4FRERE5PrKhsMTq+G+t80VeS3k0AG7p0+fJiMjg7CwsGzbw8LCOHXqVK6OsW/fPgYPHoyHhwc2m40pU6YQGhp6w/19fX3x9fW9pbpFREQkFzy9ocmTVleRP7ONbDZbtteGYVyz7UaaNWvG7t2786MsERERcQMOfWxUokQJPD09r7nLkpiYeM3dGBEREZG8cGh48fHxITw8nNjY2GzbY2NjadasmSNPJSIiIoWU3Y+NLl68yKFDh7JeHzlyhPj4eEJDQ6lQoQIjR46kd+/eNGjQgKZNmzJr1iyOHTvGkCFDHFq4iIiIFE52h5cdO3bQqlWrrNcjR44EoG/fvsydO5cePXpw5swZxo4dS0JCArVr12bp0qVUrFjRcVWLiIg42ImkVC2Q5yLsDi8tW7bEMIyb7hMdHU10dHSeixIRESlIJ5JSaTtxLanpOXdS9vf2ZNWzEQowFlJvIxERKfTOpaSRmp7B5B71qVoy6Ib7HUq8yIiYeM6lpCm8WEjhRURE5L+qlgyidtmiVpchOXDobCMRERGR/KbwIiIiIi5F4UVERERcisKLiIiIuBQN2BUREbGQ1pexn8KLiIiIRbS+TN4ovIiIiFhE68vkjcKLiIiIxbS+jH00YFdERERcisKLiIiIuBQ9NhIREckHuZlFdCjxYgFV414UXkRERBzM3llExQJ9CqAq96HwIiIi4mC5nUUEWr8lLxReRERE8olmEeUPDdgVERERl6LwIiIiIi5Fj41ERERcRG5mJxWGMTQKLyIiIk6uWKAP/t6ejIiJz3HfwtADSeFFRETEyZUN8WfVsxG5WjemMPRAUngRERFxAWVD/N06kNhDA3ZFRETEpSi8iIiIiEtReBERERGXovAiIiIiLkXhRURERFyKwouIiIi4FE2VFhERsVNOK93mZiVcyTuFFxERkVyyd6XbYoE++V9UIaTwIiIikku5XekWCkePIasovIiIiNhBK91aTwN2RURExKUovIiIiIhLUXgRERERl6LwIiIiIi5F4UVERERcisKLiIiIuBSFFxEREXEpWudFRETEzeS2PYGrLqSn8CIiIuIm7GlfAGYLg1XPRrhcgFF4ERERcRP2tC84lHiRETHxnEtJU3gRERER69jbviA3j5ic7fGSwouIiEghZG+HbGd6vKTwIiIiUgjl9hGTMz5eUngREREppFy1Q7bWeRERERGXovAiIiIiLkXhRURERFyKwouIiIi4FIUXERERcSkKLyIiIuJSFF5ERETEpSi8iIiIiEtReBERERGXovAiIiIiLkXhRURERFyKwouIiIi4FIUXERERcSlu11XaMAwAkpOTHX7sixeSybxyiYsXkklOtt3yfq5Qn6OPadXvRkQKn/z4W1cYFdTv5u9/t//+d/xmbEZu9nIhf/zxB+XLl7e6DBEREcmD48ePU65cuZvu43bhJTMzk5MnT9K6dWt27Nhxw/0aNmzI9u3b7XovOTmZ8uXLc/z4cYKDgx1Wc3642fdztnPk9Tj2fC63++a0n64b5zmHrhvnoesmb/vqusn+HQzD4MKFC5QpUwYPj5uPanG7x0YeHh6UK1cOLy+vm/6H8/T0vOH7N3sPIDg42Okvipy+gzOdI6/Hsedzud03p/103TjPOXTdOA9dN3nbV9fNtd+haNGiufqs2w7YHTp0aJ7fz+mzrqAgvoOjzpHX49jzudzuq+tG101e9tV1o+smL/vqusn7d3C7x0b5KTk5maJFi3L+/HmnT7TiPHTdSF7oupG8KCzXjdveeckPvr6+vP766/j6+lpdirgQXTeSF7puJC8Ky3WjOy8iIiLiUnTnRURERFyKwouIiIi4FIUXERERcSkKLyIiIuJSFF5ERETEpSi85JNu3bpRrFgxunfvbnUp4sSWLFnCnXfeyR133MHs2bOtLkdchP6+iL2OHz9Oy5YtqVmzJnXr1uXrr7+2uqRboqnS+WT16tVcvHiRefPmsWDBAqvLESd09epVatasyerVqwkODubuu+9m69athIaGWl2aODn9fRF7JSQk8Oeff1K/fn0SExO5++67OXDgAIGBgVaXlie685JPWrVqRZEiRawuQ5zYtm3bqFWrFmXLlqVIkSJ06NCBFStWWF2WuAD9fRF7lS5dmvr16wNQsmRJQkNDOXv2rLVF3YJCGV7WrVtH586dKVOmDDabjUWLFl2zz/Tp06lUqRJ+fn6Eh4ezfv36gi9UnNqtXkcnT56kbNmyWa/LlSvHiRMnCqJ0sZD+/kheOPK62bFjB5mZmZQvXz6fq84/hTK8pKSkUK9ePaZOnXrd92NiYhgxYgSjR49m586dtGjRgqioKI4dO5a1T3h4OLVr177m5+TJkwX1NcRit3odXe+Jrc1my9eaxXqO+PsjhY+jrpszZ87Qp08fZs2aVRBl5x+jkAOMhQsXZtvWqFEjY8iQIdm2Va9e3XjppZfsOvbq1auNBx988FZLFBeQl+to48aNRteuXbPeGz58uPHZZ5/le63iPG7l74/+vhReeb1uLl++bLRo0cKYP39+QZSZrwrlnZebSUtLIy4ujsjIyGzbIyMj2bRpk0VViavJzXXUqFEj9uzZw4kTJ7hw4QJLly6lffv2VpQrTkJ/fyQvcnPdGIZBv379aN26Nb1797aiTIfysroAZ3P69GkyMjIICwvLtj0sLIxTp07l+jjt27fnp59+IiUlhXLlyrFw4UIaNmzo6HLFSeXmOvLy8mLixIm0atWKzMxMXnjhBYoXL25FueIkcvv3R39f5H/l5rrZuHEjMTEx1K1bN2u8zCeffEKdOnUKulyHUHi5gX+OPTAMw67xCJo1IpDzdXT//fdz//33F3RZ4uRyum7090Wu52bXzT333ENmZqYVZeULPTb6hxIlSuDp6XnNXZbExMRrUq3Ijeg6krzQdSN5URivG4WXf/Dx8SE8PJzY2Nhs22NjY2nWrJlFVYmr0XUkeaHrRvKiMF43hfKx0cWLFzl06FDW6yNHjhAfH09oaCgVKlRg5MiR9O7dmwYNGtC0aVNmzZrFsWPHGDJkiIVVi7PRdSR5oetG8kLXzT9YOtfJIqtXrzaAa3769u2btc+0adOMihUrGj4+Psbdd99trF271rqCxSnpOpK80HUjeaHrJjv1NhIRERGXojEvIiIi4lIUXkRERMSlKLyIiIiIS1F4EREREZei8CIiIiIuReFFREREXIrCi4iIiLgUhRcRERFxKQovIiIi4lIUXkRERMSlKLyIiIiIS1F4EREREZfy/wBoy1j0D03LoAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"### Plot Histogram ####\n",
"\n",
"mybins= np.logspace(-1,2.15,45)\n",
"plt.hist(inverse_sampling(100000),bins = mybins,density=True,histtype ='step',log = True,label = 'inverse sampling')\n",
"plt.loglog(x,pstar(x),label='x^-2.3')\n",
"plt.legend()\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "65d84894",
"metadata": {},
"source": [
"## Inverse sampling (Salpeter Mass Function)"
]
},
{
"cell_type": "code",
"execution_count": 227,
"id": "8fd17bf7",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"15.347031041410887\n"
]
},
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7f32ab07bc10>"
]
},
"execution_count": 227,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAG9CAYAAAAcFdw9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABe4ElEQVR4nO3de1yP9/8/8Me784EiEVIyOZQUKlbEGstyzGyzmeQ4kUNyytjBZksSOVSY445iDjOHWXydypxSG4uRlRwiod46qNT1+8PP+7P3Kt7Vu6734XG/3brddr2u631dz+xST6/X8/V6SQRBEEBERESkgXTEDoCIiIiorjDRISIiIo3FRIeIiIg0FhMdIiIi0lhMdIiIiEhjMdEhIiIijcVEh4iIiDQWEx0iIiLSWEx0iIiISGMx0SEiIiKNxUSHiIiINJbaJzp///03unTpIvsyNjbGnj17xA6LiIiIVIBEkzb1zM/Ph52dHW7cuAFTU1OxwyEiIiKRqX2Pzr/t3bsXffv2ZZJDREREAFQg0Tlx4gQGDx6Mli1bQiKRVDrsFBMTgzZt2sDIyAiurq44efJkpffavn07RowYUccRExERkboQPdEpKCiAi4sL1qxZU+n5uLg4BAcHY8GCBUhOToaXlxd8fX2RmZkpd51UKkViYiIGDBhQH2ETERGRGlCpGh2JRILdu3fDz89P1tajRw9069YNsbGxsjYHBwf4+fkhLCxM1vbtt9/i0KFD+O677174jOLiYhQXF8uOy8vL8fDhQzRp0gQSiUR53wwRERHVGUEQ8PjxY7Rs2RI6OlX32+jVY0zVVlJSgqSkJISGhsq1+/j44NSpU3Jt27dvx4cffvjSe4aFhWHRokVKjZOIiIjEcfPmTbRq1arK8yqd6OTk5KCsrAxWVlZy7VZWVrh7967sOC8vD2fPnsXOnTtfes/58+cjJCRE7rO2tra4efMmzMzMlBc8ERER1RmpVAobGxs0bNjwhdepdKLz3H+HlARBkGszNzfHvXv3FLqXoaEhDA0NER0djejoaJSVlQEAzMzMmOgQERGpmZeVnYhejPwilpaW0NXVleu9AYDs7OwKvTzVFRQUhNTUVJw7d65W9yEiIiLVpdKJjoGBAVxdXREfHy/XHh8fD09PT5GiIiIiInUh+tBVfn4+0tLSZMfp6elISUmBhYUFbG1tERISAn9/f7i5ucHDwwPr169HZmYmAgMDa/Xc/w5dERGJoby8HCUlJWKHQaRy9PX1oaurW+v7iD69/NixY/D29q7QHhAQgC1btgB4tmDg0qVLkZWVBScnJ6xYsQK9e/dWyvOlUinMzc2Rl5fHGh0iqlclJSVIT09HeXm52KEQqaRGjRqhefPmldbhKPr7W/RER2xMdIhIDIIgIDMzE6WlpS9dB4RI2wiCgMLCQmRnZ6NRo0Zo0aJFhWsU/f0t+tCVWDh0RURievr0KQoLC9GyZUuYmJiIHQ6RyjE2NgbwbAJSs2bNajyMpbX/hOCsKyIS0/N/ZBkYGIgcCZHqev6PgNLS0hrfQ2sTHSIiVcCtZ4iqpoy/H1qb6ERHR8PR0RHu7u5ih0JERER1RGtrdIKCghAUFCQrZiIiUgW3c4vwqKD+pps3NjWAdSPjenuepqtsc2pNkJGRgTZt2iA5ORldunSRzZh+9OgRGjVqJHZ4L6S1iQ4Rkaq5nVuEfpHHUVRaf5MkjPV1cXhWH4WTnTFjxmDr1q2YNGkS1q5dK3duypQpiI2NlVseRBk+++wz7NmzBykpKUq539GjRxEREYEzZ86gqKgIdnZ28PX1RUhICKytrZXyDE3n6emJrKwstegoYKJDRKQiHhWUoKi0DFEjusC+WYM6f15adj6C41LwqKCkWr06NjY22LZtG1asWCGbGfPkyRP8+OOPsLW1ratwlWLdunWYMmUKAgICsHPnTtjZ2SEzMxPffPMNIiMjsXz58hrdt6SkRKsKyw0MDNC8eXOxw1CMoKXWrFkjODg4CO3btxcACHl5eWKHRERapKioSEhNTRWKiopkbRdv5Qqt5+0TLt7KrZcYavK8gIAAYejQoULnzp2F7777Ttb+/fffC507dxaGDh0qBAQEyNqfPHkiTJs2TWjatKlgaGgo9OzZUzh79qzs/NGjRwUAwuHDhwVXV1fB2NhY8PDwEK5cuSIIgiBs3rxZACD3tXnzZkEQBCE3N1eYOHGi0LRpU6Fhw4aCt7e3kJKSUmXsN2/eFAwMDITg4OBKzz969EgQBEHIyckR3nvvPcHa2lowNjYWnJychB9++EHu2j59+ghBQUHCzJkzhSZNmgi9e/cWBEEQAAi7d++WXffnn38K3t7egpGRkWBhYSFMnDhRePz4cZUxPnz4UBg5cqRgaWkpGBkZCfb29sKmTZtk5+fOnSu0a9dOMDY2Ftq0aSMsXLhQKCkpkZ3/9NNPBRcXF2Hjxo2CjY2NYGpqKgQGBgpPnz4VwsPDBSsrK6Fp06bC4sWL5Z4LQIiJiRHefPNNwcjISLCzsxO2b98uO5+eni4AEJKTkwVB+N//t+d/Zps3bxbMzc2FX3/9VejYsaNgamoq9O/fX7hz547sHqWlpcK0adMEc3NzwcLCQpg7d64wevRoYejQoVX+eVT29+S5vLw8hX5/a20xMqeXExHV3NixY7F582bZ8aZNmzBu3LgK182dOxc7d+7E1q1bceHCBdjb26N///54+PCh3HULFixAZGQkzp8/Dz09Pdm9RowYgVmzZqFTp07IyspCVlYWRowYAUEQMHDgQNy9excHDhxAUlISunXrhr59+1a493M7duxASUkJ5s6dW+n557UmT548gaurK/bt24dLly7hww8/hL+/P86cOSN3/datW6Gnp4fExESsW7euwv0KCwvx5ptvonHjxjh37hx27NiBw4cPY+rUqVX+uX788cdITU3FwYMHcfnyZcTGxsLS0lJ2vmHDhtiyZQtSU1OxcuVKfP3111ixYoXcPa5fv46DBw/i119/xY8//ohNmzZh4MCBuHXrFo4fP47w8HAsXLgQp0+frvDs4cOH448//sCoUaPw/vvv4/Lly1XGWtn3u2zZMnz77bc4ceIEMjMzMXv2bNn58PBwfP/999i8eTMSExMhlUqxZ88ehe9fYy9Mg7SAohkhEZEyqXuPzv379wVDQ0MhPT1dyMjIEIyMjIT79+/L9ejk5+cL+vr6wvfffy/7fElJidCyZUth6dKlgiDI9+g8t3//fgGA7M/meS/Fvx05ckQwMzMTnjx5Itfetm1bYd26dZXGPnnyZMHMzEzh7/XfBgwYIMyaNUt23KdPH6FLly4VrsO/enTWr18vNG7cWMjPz5ed379/v6CjoyPcvXu30ucMHjxYGDt2rMJxLV26VHB1dZUdf/rpp4KJiYkglUplbf379xfs7OyEsrIyWVuHDh2EsLAwubgDAwPl7t2jRw9h8uTJgiAo1qMDQEhLS5N9Pjo6WrCyspIdW1lZCREREbLjp0+fCra2tnXeo8MaHSIiqjZLS0sMHDgQW7dulfWu/LvnAXjWs1BaWoqePXvK2vT19dG9e/cKPQXOzs6y/36+3H92dnaVNT9JSUnIz89HkyZN5NqLiopw/fr1Sj8jCIJC67KUlZVhyZIliIuLw+3bt1FcXIzi4mKYmprKXefm5vbC+1y+fBkuLi5yn+vZsyfKy8vx999/w8rKqsJnJk+ejOHDh+PChQvw8fGBn58fPD09Zed/+uknREVFIS0tDfn5+Xj69GmF7Q/s7OzQsGFD2bGVlRV0dXXlthmxsrJCdna23Oc8PDwqHFenANzExARt27aVHbdo0UL2jLy8PNy7dw/du3eXndfV1YWrq2ud7/XGRIeIiGpk3LhxsmGY6OjoCueF/7+V4n+Ti8oSDn19fdl/Pz/3ol+A5eXlaNGiBY4dO1bhXFXTndu3b4+8vDxkZWVVunfSc5GRkVixYgWioqLQuXNnmJqaIjg4uMIu8/9NfP7rRYlVVe2+vr64ceMG9u/fj8OHD6Nv374ICgrCsmXLcPr0abz33ntYtGgR+vfvD3Nzc2zbtg2RkZFy9/j3n+XzZ1XWpkiCUZ0F+yp7hvCf7TQrexfqmtbW6BARUe28+eabKCkpQUlJCfr371/hvL29PQwMDJCQkCBrKy0txfnz5+Hg4KDwcwwMDCrsS9itWzfcvXsXenp6sLe3l/v6b8/Sc2+//TYMDAywdOnSSs/n5uYCAE6ePImhQ4di1KhRcHFxwSuvvIJr164pHO9zjo6OSElJQUFBgawtMTEROjo6aN++fZWfa9q0KcaMGYPvvvsOUVFRWL9+veyzrVu3xoIFC+Dm5oZ27drhxo0b1Y6rKv+t2Tl9+jQ6duyolHubm5vDysoKZ8+elbWVlZUhOTlZKfd/Ea3t0eGmnkREtaOrqysbgqpsw0VTU1NMnjwZc+bMgYWFBWxtbbF06VIUFhZi/PjxCj/Hzs4O6enpSElJQatWrdCwYUP069cPHh4e8PPzQ3h4ODp06IA7d+7gwIED8PPzq3RYycbGBitWrMDUqVMhlUoxevRo2NnZ4datW/jmm2/QoEEDREZGwt7eHjt37sSpU6fQuHFjLF++HHfv3q1WcgYAH3zwAT799FMEBATgs88+w/379zFt2jT4+/tXOmwFAJ988glcXV3RqVMnFBcXY9++fbLn2tvbIzMzE9u2bYO7uzv279+P3bt3VyumF9mxYwfc3NzQq1cvfP/99zh79iw2btyotPtPmzYNYWFhsLe3R8eOHbF69Wo8evSozrdB0dpEhysjE5GqSsvOV5vn/Lc+5L+WLFmC8vJy+Pv74/Hjx3Bzc8OhQ4fQuHFjhZ8xfPhw7Nq1C97e3sjNzcXmzZsxZswYHDhwAAsWLMC4ceNw//59NG/eHL17964yiQCeLWrYvn17LFu2DMOGDZMtGDho0CCEhIQAeDb7KD09Hf3794eJiQk+/PBD+Pn5IS8vT+GYgWc1K4cOHcKMGTPg7u4OExMTDB8+/IVr9RgYGGD+/PnIyMiAsbExvLy8sG3bNgDA0KFDMXPmTEydOhXFxcUYOHAgPv74Y3z22WfViqsqixYtwrZt2zBlyhQ0b94c33//PRwdHZVybwCYN28e7t69i9GjR0NXVxcffvgh+vfvX+NdyRUlEepjgEyFPU908vLyXvoXlohIWZ48eYL09HS0adMGRkZGANRjZWTSTGJsXVFeXg4HBwe8++67+OKLLyq9prK/J88p+vtba3t0iIhUjXUjYxye1Yd7XZFGunHjBn777Tf06dMHxcXFWLNmDdLT0zFy5Mg6fS4THSIiFWLdyJiJB2kkHR0dbNmyBbNnz4YgCHBycsLhw4erXftUXUx0iIiItFx9VLHY2NggMTGxzp/zX1o7vTw6OhqOjo5wd3cXOxQiIiKqI1qb6HCvKyIiIs2ntYkOERERaT4mOkRERKSxmOgQERGRxmKiQ0RERBqLiQ4REdWJa9euwcrKCiYmJqJMKyYCmOgQEVEduHPnDnx8fNCrVy+MHz8egwYNwsWLF1/6uYcPH2LatGno0KEDTExMYGtri+nTp790n6nY2Fg4OzvDzMwMZmZm8PDwwMGDB5X17ZAa44KBtXD4n8PwtvOGrk7dbkhGRKROHj16JEtytmzZAl1dXTRs2BD9+/dHQkICXnnllSo/e+fOHdy5cwfLli2Do6Mjbty4gcDAQNy5cwc//fRTlZ9r1aoVlixZAnt7ewDA1q1bMXToUCQnJ6NTp05K/x5JjQhaas2aNYKDg4PQvn17AYCQl5dXrc+fvnlawGcQ2q5sK6w9t1YoKi2qo0iJSBMVFRUJqampQlGRev3syM7OFqysrIQvv/xS1nb69GlBX19fOHTokFBQUCC8+uqrwqRJk4SysjK5zy5evFho27atkJWVVa1nbt++XTAwMBBKS0ur9bnGjRsLGzZsqNZnSLW86O9JXl6eQr+/uXt5DXcv3/HXDkzePxkPih4AAJo3aI6Zr85EoFsgzAy5CzoRvVhVuzIXlBRU+RldHV0Y6Sl2rY5EB8b6xi+91tTAtDphAwAOHDgAPz8/nDp1Ch07dkTXrl0xcOBAREVFVfteitiwYQPmz5+P+/fvK3R9WVkZduzYgYCAACQnJ8PR0bFO4qK6p4zdy5no1DDRAZ794NhwYQOW/b4Mt6S3AADmhuYIcg/CfK/5aGDQoC5CJiINUNUPcMkiSZWfGdBuAPaP3C87Nv3KFIWlhZVe26d1Hxwbc0x23DSiKXIKcypcJ3xas18BQUFBOHz4MNzd3fHHH3/g3LlzFX4RKcODBw/QrVs3+Pv7Y/HixS+89uLFi/Dw8MCTJ0/QoEED/PDDDxgwYIDSY6L6o4xEh8XItWBqYIoZr87A9enXsWXoFjhYOiCvOA/fXfwOhrqGYodHRFRnli1bhqdPn2L79u34/vvva5TkfPXVV2jQoIHsKzMzU+68VCrFwIED4ejoiE8//fSl9+vQoQNSUlJw+vRpTJ48GQEBAUhNTa12XKRZ2KNTix6d/yoXyrH3770oKy/DcMfhAIDSslLM/m02JrpOhFMzJ2WETEQaQJ2HrgDgr7/+gpubG0pLS7F7924MHjy42vd4+PAhHj58KDu2s7ODnt6zOTKPHz9G//79YWJign379tUokerXrx/atm2LdevWVfuzpBqU0aPDWVdKpCPRgV9HP7m2uL/isOrsKqw6uwqD2g/C/F7z4WnjKU6ARKTyqpN41NW1L1NSUoIPPvgAI0aMQMeOHTF+/HhcvHgRVlZW1bqPhYUFLCwsKrRLpVL0798fhoaG2Lt3b42HxARBQHFxcY0+S5qDQ1d1zNnKGW87vg0JJNh3dR96buqJ3pt74+C1g9DyzjQiUlMLFixAXl4eVq1ahblz58LBwQHjx49Xyr0fP34MHx8fFBQUYOPGjZBKpbh79y7u3r2LsrIy2XV9+/bFmjVrZMcfffQRTp48iYyMDFy8eBELFizAsWPH8MEHHyglLlJf7NGpY85Wztjxzg5cfXAVEYkR2PrHVpzMPImTP5yEs5UzTow5AXMjc7HDJCJSyLFjxxAVFYWjR4/Khgu+/fZbODs7IzY2FpMnT67V/ZOSknDmzBkAkK2J81x6ejrs7OwAANevX0dOzv+Kq+/duwd/f39kZWXB3Nwczs7O+PXXX/HGG2/UKh5Sf6zRUWKNjiJuS29jxekVWJe0Dq4tXOVmRZQL5dCRsJONSBu8qPaAiJ7hrCs1ZG1mjWU+y3Aj+Aa+Hvy1rP1+wX3Yr7JHeEI48p68eKlzIiIiUoxGJDrp6enw9vaGo6MjOnfujIKCqmciqAoLYwu0a9JOdrzhwgak56Yj9EgobKNs8dGRj3Av/56IERIREak/jUh0xowZg88//xypqak4fvw4DA3Vbw2b2Z6zsdVvKxybOkJaLEVYQhjsVtohaH8Q0h+lix0eERGRWlL7ROevv/6Cvr4+vLy8ADybrvh8HQZ1oq+rj9Euo3Fx8kX8/N7PeLXVq3jy9AlizsfAKdYJ0mKp2CESERGpHdETnRMnTmDw4MFo2bIlJBIJ9uzZU+GamJgYWSGSq6srTp48KTt37do1NGjQAEOGDEG3bt3w1Vdf1WP0yqcj0cGQDkNwatwpHAs4hv5t+2NU51Fy+2ddyr4kYoRERETqQ/REp6CgAC4uLnLrIfxbXFwcgoODsWDBAiQnJ8PLywu+vr6ypcJLS0tx8uRJREdH4/fff0d8fDzi4+Pr81uoExKJBH3s+uDXUb8iemC0rP1C1gV0ju0Mr81e2H91P9fiIVJz/DtMVDVl/P0QPdHx9fXF4sWL8dZbb1V6fvny5Rg/fjwmTJgABwcHREVFwcbGBrGxsQCAVq1awd3dHTY2NjA0NMSAAQOQkpJS5fOKi4shlUrlvlSdns7/huIuZF2Aga4BEjITMOjHQXBZ64IfLv6Ap+VPRYyQiKpLV1cXwLNVhomocoWFzzat1dfXr/E9VLqYpaSkBElJSQgNDZVr9/HxwalTpwAA7u7uuHfvHh49egRzc3OcOHECkyZNqvKeYWFhWLRoUZ3GXZcmdJuAAe0GIOp0FGLPx+Ji9kV8sOsDLPi/BZjjOQfju46HoZ76FWMTaRs9PT2YmJjg/v370NfXh46O6P/uJFIZgiCgsLAQ2dnZaNSokewfBjWh0olOTk4OysrKKuyfYmVlhbt37wJ49sPiq6++Qu/evSEIAnx8fDBo0KAq7zl//nyEhITIjqVSKWxsbOrmG6gjLRu2xNI3lmJ+r/mIOReDlWdWIiM3A1+e/BLjuypnGXYiqlsSiQQtWrRAeno6bty4IXY4RCqpUaNGaN68ea3uodKJznMSiUTuWBAEuTZfX1/4+voqdC9DQ0O1nH5emcbGjbGg9wLM9JiJzcmbYWpgKuvNeVr+FEsTl2Jc13Fo3qB2LwkR1Q0DAwO0a9eOw1dEldDX169VT85zKp3oWFpaQldXV9Z781x2dna1d8n9r+joaERHR8ttEqeuTPRNENQ9SK5tZ+pOLPi/Bfj8+OcY02UM5njOQVuLtiJFSERV0dHR4RYQRHVIpQeFDQwM4OrqWmEWVXx8PDw9PWt176CgIKSmpuLcuXO1uo+qat6gOTxaeaC4rBjrktah/Zr2eH/n+0i5myJ2aERERPVG9EQnPz8fKSkpsplS6enpSElJkU0fDwkJwYYNG7Bp0yZcvnwZM2fORGZmJgIDA2v13OjoaDg6OsLd3b2234JK6mPXB4njEnF8zHH42vuiXCjHtkvb0HVdVwz4fgDyS/LFDpGIiKjOib57+bFjx+Dt7V2hPSAgAFu2bAHwbMHApUuXIisrC05OTlixYgV69+6tlOfX9+7lYkm5m4LwxHBs/2s7Xm31KhLGJlSofSIiIlIXiv7+Fj3REZu2JDrPXX94HY9LHqNL8y4AgIdFDzHwh4EIcg/CiE4joK9b87UKiIiI6ouiv79FH7qi+tXWoq0syQGA2HOxOH3rNPx3+6Pd6nZYc3YNCksLxQuQiIhIibQ20dH0Gh1FBXUPwpevf4mmJk1xI+8Gph2cBrsoO3x54kvkPskVOzwiIqJa4dCVlg1dVaWotAibUzYj4lQEMnIzAAAWxhbIDM6EqYGpuMERERH9B4euqFqM9Y0xxX0Krk27hu+GfQenZk4Y0mGIXJKT9ThLxAiJiIiqT2sTHQ5dVU5PRw8fOH+APwP/xGrf1bL2i/cuwmaFDd776T0kZyWLGCEREZHiOHTFoSuFLDu1DHPi58iO+7ftj/m95qN3696cpk5ERPWOQ1ekVLM9ZyNlUgred3ofOhIdHLp+CK9tfQ2emzyx9++9KBfKxQ6RiIioAvbosEen2v559A8iEiOwOWUzisuK0cy0GTJmZMBY31js0IiISEuwR+clWKNTc680fgWxg2KREZyB0J6hWOC1QJbklAvl2JS8iWvxEBGRSmCPDnt0lGr35d14a/tbsDSxxIweMxDkHoTGxo3FDouIiDQMe3RIFBKJBG0atUFOYQ4+PvoxbKNsMee3Objz+I7YoRERkRZijw57dJTuaflT7PhrB5YkLsGf9/4EABjoGmC082is8l3FWh4iIqo19uiQaPR09PB+5/eRMikF+0fuh5etF0rKSpByLwVGekZih0dERFpET+wAxBIdHY3o6GiUlZWJHYrGkkgkGNBuAAa0G4DEzERZGwDkPcnDmJ/HYHr36XjN7jWuxUNERHWCQ1ccuhLFkoQlmH9kPgCgh3UPhPYKxZAOQ6AjYScjERG9HIeuSKWN6DQCQe5BMNIzwpnbZzAsbhicYpywNWUrSstKxQ6PiIg0BHt02KMjqnv597DqzCpEn4tGXnEegGfr9Pw15S/W8xARUZXYo0NqwaqBFb7s+yUyZ2YivF84mjdoDo9WHnJJDhcfJCKimmKPDnt0VMqTp0/wuPgxmpo2BQBcybmC7l93x8RuExHiEQJrM2uRIyQiIlXAHh1SS0Z6RrIkBwC+//N7PC55jOWnl6PNyjaYsHcC/s75W8QIiYhInWhtosO9rtTD596f4+AHB9G7dW+UlpdiY/JGOEQ74O3tb+P8nfNih0dERCqOQ1cculIbp26ewpKEJfjl6i8AgCbGTXAr5BaLlomItJCiv7+1dsFAUj+eNp7Y+/5eXMq+hKWJS9HRsqMsyREEAYeuH4JPWx+uxUNERDLs0WGPjkbYd3UfBv84GB0tO2Jez3kY2XkkDHQNxA6LiIjqCIuRSatkF2SjkVEjXMm5grE/j0XbVW0RdToK+SX5YodGREQiYo8Oe3Q0hrRYivVJ67H89+XIys8CAFgYW2Ba92n4yOsj9vAQEWkQ9uiQ1jEzNMNsz9n4Z8Y/WD9oPewt7PGw6CH2XNkDfR19scMjIiIRMNEhjWOkZ4SJrhNxJegK4t6Ow9I3lsp2R39c/BhT9k/BlZwrIkdJRET1gYkOaSxdHV282+ld+LT1kbV9feFrxJ6PhWO0I4ZvH45zt8+JGCEREdU1rU10uGCgdurdujeGdhgKAQJ2Xd6F7hu6o983/XD4n8PQ8nI1IiKNxGJkFiNrpdT7qQhPDMcPF3/A0/KnAIAe1j1wYuwJFi0TEakBFiMTvYBjU0ds9duKtGlpmNZ9Goz1jNG6UWu5JKesvEzECImISBnYo8MeHQJwv+A+ip4WwdbcFgCQ9jANfb/pixk9ZuBD1w/RwKCByBESEdG/sUeHqBqamjaVJTkAsPb8WmTmZWLWb7Ngu8IWnx79FDmFOSJGSERENcFEh6gSX77+Jb4e/DXaWbTDoyeP8PmJz9E6qjWCfw3GzbybYodHREQK4tAVh67oBcrKy7Dr8i6EJYQh+W4yAKCZaTPcmnkL+rpchJCISCwcuiJSAl0dXbzT6R0kfZiEQ6MOwdvOGx92+1CW5AiCgJS7KeIGSUREVdITOwAidSCRSODT1gc+bX3kZmPF/xOP/t/1h7edN+b3mo9+r/STrcJMRETiY48OUTXp6ujK/vvivYvQ09HD0Yyj8PnOB+5fu+On1J84NZ2ISEVoRKKjp6eHLl26oEuXLpgwYYLY4ZAWmeU5C9enX8eMHjNgom+CpKwkvLPjHThEO2DDhQ2yxQiJiEgcGlGMbGlpiZycmk39ZTEyKUtOYQ5Wn1mN1WdX49GTR+jUtBP+nPwndCQa8e8JIiKVwmJkonpmaWKJRd6LkDkzE8t9lmPx64tlSU5haSG+OP4F7hfcFzlKIiLtInqic+LECQwePBgtW7aERCLBnj17KlwTExODNm3awMjICK6urjh58qTcealUCldXV/Tq1QvHjx+vp8iJKtfAoAFmesyEX0c/WdvGCxvxybFP0DqqNWYcnIHMvEzxAiQi0iKiJzoFBQVwcXHBmjVrKj0fFxeH4OBgLFiwAMnJyfDy8oKvry8yM//3iyIjIwNJSUlYu3YtRo8eDalUWl/hEymkXZN2cG3hiqKnRVh1dhXarmqLMXvGIPV+qtihERFpNJWq0ZFIJNi9ezf8/PxkbT169EC3bt0QGxsra3NwcICfnx/CwsIq3MPX1xdffPEF3NzcKn1GcXExiouLZcdSqRQ2Njas0aE6JwgCjqQfwZKEJTiSfkTWPqzjMGx/Zzv0dLjaAxGRojSiRqekpARJSUnw8fGRa/fx8cGpU6cAAI8ePZIlLrdu3UJqaipeeeWVKu8ZFhYGc3Nz2ZeNjU3dfQNE/yKRSNDvlX44PPowzkw4g7cc3oIEz9bcYZJDRFQ3VPqna05ODsrKymBlZSXXbmVlhbt37wIALl++jEmTJkFHRwcSiQQrV66EhYVFlfecP38+QkJCZMfPe3SI6lN36+7Y+e5OXMm5Ikt2ACAjNwMjfhqBWR6zMNxhuNyaPUREVH0qneg899+VZgVBkLV5enri4sWLCt/L0NAQhoaGSo2PqKY6WnaUO446HYWzt89ixE8jYG9hj7meczHaZTQM9fjOEhHVhEoPXVlaWkJXV1fWe/NcdnZ2hV6e6oqOjoajoyPc3d1rdR8iZfq498f4rM9nsDC2QNrDNHy470O0WdkGEYkRkBazyJ6IqLpUOtExMDCAq6sr4uPj5drj4+Ph6elZq3sHBQUhNTUV586dq9V9iJSpiUkTfPrap8gMzsSK/ivQyqwVsvKzMPfwXDjHOnOlZSKiahJ96Co/Px9paWmy4/T0dKSkpMDCwgK2trYICQmBv78/3Nzc4OHhgfXr1yMzMxOBgYG1em50dDSio6NRVlbzPYlu5xbhUUHJS69rbGoA60bGNX4OaR9TA1MEvxqMKe5T8P2f3yM8MRx+Hf1kRcuCIODO4zuwNrMWOVIiItUm+vTyY8eOwdvbu0J7QEAAtmzZAuDZgoFLly5FVlYWnJycsGLFCvTu3Vspz6/pFhC3c4vQL/I4ikpfnigZ6+vi8Kw+THaoxsqFcjx5+gQm+iYAgKPpR/HGt29gZOeRmNdzHjo16yRyhERE9UvR39+iJzpiq2mic+l2HgatTkDUiC6wb9agyuvSsvMRHJeCfdN6wcnaXBkhE2HBkQX4KuEr2fGQDkMQ2jMUHjYeIkZFRFR/FP39LfrQlViUMXQFAPbNGjCBoXr3Zd8vMcxhGJYkLMGuy7uw9++92Pv3XvRu3RuhPUPxpv2bFWYrEhFpI5UuRq5LLEYmdefW0g0/vfsTLgddxviu46Gvo48TN05g5qGZEKDVHbVERDJam+gQaYoOlh2wYcgG/DPjH4S8GoKPe38s2zX9ydMn2HBhA548fSJylERE4tDaoav6lpadr9B1nKFFNdXKrBUi+0fKtW1N2YrA/YH4+OjHmPnqTAS6BcLMkHu6EZH20NpER1k1Oi/T2NQAxvq6CI5LUeh6ztAiZWpg0AA2Zja4Kb2JeYfn4auTXyHIPQgzXp2BZqbNxA6PiKjOcdZVLWddKTKbStH1djhDi+pCSVkJfrz4I8ITw3E55zIAwEjPCOO7jsfKN1dyPy0iUkucdaVCrBsZs4eGRGOga4CALgHwd/HH3r/3IiwhDGdvn0VGbgaTHCLSeEx0iLSEjkQHfh39MLTDUBzLOIbGxo1l525Jb2HqgamY4zkHPW17ihglEZFyae2sK27qSdpKIpHAu403ujTvImtb8fsK/Pz3z+i1uRe8NnvhwLUD0PJRbSLSEFrboxMUFISgoCDZGJ8qUWSGFmdnkTIFugVCWizF1j+2IiEzAQN/GAhnK2eE9gzFO53eke2xRUSkbliMXA/FyIri/lkkttvS21hxegXWJa1DfsmzhNvZyhnJk5Jla/MQEakCFiOrIetGxjg8q89LZ2g9n531qKCEiQ4plbWZNZb5LMNHXh8h5lwMVp5Zib5t+solOfkl+WhgUPX+bkREqoSJjorhDC1SBRbGFljYeyFmvjoTpeWlsvbnw1pT3KZgxqsz0LxBcxGjJCJ6Oa3ti2YxMtHLmRqYopFRI9nxDxd/gLRYiiWJS2AXZYcp+6fgn0f/iBcgEdFLaG2iw009iapvzYA12DNiD3pY90BxWTFiz8ei3ep2+GDXB/jz3p9ih0dEVAGHrtQY98+i+qYj0cHQjkMxpMMQnLhxAksSl+DXtF/xw8UfcObWGVyddpVFy0SkUpjoqCHun0Vik0gk6GPXB33s+iA5KxnhieFyRcslZSU48s8RvGn/JiQSicjREpE2Y6KjhhSdnQVwhhbVva4tumLb29vk2r778zuM3zsenZt1xrye8zDCaQTX4iEiUfAnj5ri7CxSZY+LH6OBQQNczL6IUbtHYeHRhZjjOQdju4yFsT7fWyKqPxxMJyKlm/HqDGQGZ2Kx92I0NWmKjNwMBB0Igt1KO4SdDEO5UC52iESkJbQ20eH0cqK61di4MRb0XoCM4Ays9l2N1uatkV2QjaMZR1mwTET1RmuHrlR5r6u6wP2zSCwm+iaY2n0qJrlOQtxfcbC3sJedy3qchS9OfIFZHrPQ1qKtiFESkabS2kRHW1RnhhZnZ1Fd0tfVxyjnUXJtK8+sROz5WKxLWod3O72L0J6hcGnuIlKERKSJmOhoOO6fRapscPvBuJh9EQeuHcC2S9uw7dI2+Nr7IrRXKLxsvTg1nYhqjYmOFuAMLVJVPW17Yv/I/fjj7h8ITwxH3F9xOJh2EAfTDqJvm76I949nskNEtcKKQCISnUtzF/ww/AdcnXoVga6BMNQ1hIOlg1ySU1ZeJmKERKSumOgQkcpoa9EWsYNikT4jHQt7L5S1n7l1Bvar7RF9NhpFpUUiRkhE6oZDVySH+2eRKmjRsIXcccz5GGTkZmDqwalYdHwRgl8NxhT3KXI7qxMRVYaJDgHg/lmk2tYOXIse1j0QcSoCGbkZWPB/C7AkYQkC3QIx89WZFRIjIqLntDbRiY6ORnR0NMrKOO4PcP8sUm3G+saY4j4FH7p+iLhLcViSuASXsi8h4lQE9lzZg7+n/s2iZSKqlNYmOtq2YKAiODuLVJ2ejh4+cP4AIzuPxP5r+xGWEIb3Or0nS3JKy0rx1/2/0KV5F3EDJSKVwWJkIlI7EokEg9oPQuK4RAR1D5K1b7u0DV3XdcWb372JYxnHIAiCiFESkSpgokNEau3f+2ZdfXAVOhIdHLp+CN5bveG5yRM/X/mZm4gSaTGtHbqi2uP+WaRqvnj9C4ztOhbLTi3DpuRNOH3rNPzi/ODY1BHzes6Dv7M/a3mItAwTHao27p9FquyVxq8gZmAMPunzCVaeXomY8zFIvZ+KzSmbMdpltNjhEVE9Y6JD1cb9s0gdNG/QHGH9whDaKxSx52PRw7qH7Nz9gvtYl7QOQe5BaGzcWMQoiaiuMdGhGuEMLVIX5kbmCO0VKte26swqLD65GOGJ4ZjkOgkhHiFo2bClSBESUV1iMTIRaZ1uLbrB2coZ+SX5iPw9Em1WtsHEvRNx7cE1sUMjIiVjokNEWmeYwzCkTErB/pH70cu2F0rKSrAheQM6rOkA/93+nJZOpEE0ZuiqsLAQDg4OeOedd7Bs2TKxw6F/4ewsUkUSiQQD2g3AgHYDkJCZgCUJS7D/2n6Y6JlwZhaRBtGYROfLL79Ejx49Xn4h1RvOziJ10cu2F/aN3Ic/7/0pt1HohawLmLJ/CkJ7hWJIhyFya/YQkXrQiETn2rVruHLlCgYPHoxLly6JHQ79f5ydRerG2cpZ7njZqWU4c/sMhsUNg4OlA+b1nIeRnUdCX1dfpAiJqLqqleiMGzdOoes2bdqk8D1PnDiBiIgIJCUlISsrC7t374afn5/cNTExMYiIiEBWVhY6deqEqKgoeHl5yc7Pnj0bEREROHXqlMLPpfrB2Vmkzlb0X4E2jdog+lw0Ludcxpifx+Djox9jtudsjO86HqYGpmKHSEQvUa1+2C1btuDo0aPIzc3Fo0ePqvyqjoKCAri4uGDNmjWVno+Li0NwcDAWLFiA5ORkeHl5wdfXF5mZmQCAn3/+Ge3bt0f79u2r9VwiopexamCFL/t+iRvBNxDeLxzNGzTHTelNzPh1Bnpu6smiZSI1UK0encDAQGzbtg3//PMPxo0bh1GjRsHCwqJWAfj6+sLX17fK88uXL8f48eMxYcIEAEBUVBQOHTqE2NhYhIWF4fTp09i2bRt27NiB/Px8lJaWwszMDJ988kml9ysuLkZxcbHsWCqV1ip+Uh5FipYBFi5T/TM3MsfcnnMxvcd0fPPHN1iauBSjXUbLipaflj/Fvfx7sDazFjlSIvoviVDNf5IUFxdj165d2LRpE06dOoWBAwdi/Pjx8PHxqfVMBYlEIjd0VVJSAhMTE+zYsQPDhg2TXTdjxgykpKTg+PHjcp/fsmULLl269MJZV5999hkWLVpUoT0vLw9mZmYKx3rpdh4GrU7Avmm94GRtrvDnqKLbuUXoF3kcRaVlCl3PwmUS29PypygrL4OhniGAZ7umj949GqNdRmOO5xx0sOwgcoREmk8qlcLc3Pylv7+rXYxsaGiI999/H++//z5u3LiBLVu2YMqUKSgtLUVqaioaNGhQq8D/LScnB2VlZbCyspJrt7Kywt27d2t0z/nz5yMkJER2LJVKYWNjU6s4qXYULVoGWLhMqkFPRw96Ov/78Xks4xhKy0uxMXkjNiVvwlsOb2F+r/lwbekqYpREBNRy1pVEIoFEIoEgCCgvL1dWTJU+598EQai092jMmDEvvZehoSEMDQ0RHR2N6OholJUp1otAdYtFy6TO1g5ai9Euo7EkYQl+ufoLdl7eiZ2Xd6LfK/0wv9d8eNt5c20eIpFUe1GI4uJi/Pjjj3jjjTfQoUMHXLx4EWvWrEFmZqZSe3MAwNLSErq6uhV6b7Kzsyv08lRXUFAQUlNTce7cuVrdh4gIADxtPLH3/b24OPki/J39oSvRxeF/DuPjox+LHRqRVqtWojNlyhS0aNEC4eHhGDRoEG7duoUdO3ZgwIAB0NFR/kJaBgYGcHV1RXx8vFx7fHw8PD09lf48IqLacmrmhG+GfYO06WmY1n0aPu79saw352HRQ2xJ2YKSspcP0xKRclRr6Grt2rWwtbVFmzZtcPz48QrFwM/t2rVL4Xvm5+cjLS1Ndpyeno6UlBRYWFjA1tYWISEh8Pf3h5ubGzw8PLB+/XpkZmYiMDCwOqFXwKEr9cZtJUjV2TWywyrfVXJta86uwafHPsXHRz/GLI9ZmNBtAhoYKLcnnIjkVSvRGT16tNLHmc+fPw9vb2/Z8fNC4YCAAGzZsgUjRozAgwcP8PnnnyMrKwtOTk44cOAAWrduXavnBgUFISgoSFa1TeqB20qQOrMytUKLBi1wS3oLMw/NxBcnvsD07tMxtftUNDFpInZ4RBqp2tPLNY2i09P+i9PLxXM7t0jhbSX4/4dUTfHT4mdr8ZxairSHz3qzTfRNEOQehPB+4SxaJlJQnU0v1xQculJfnKFF6sxQzxATXSdiXNdx2Hl5J5YkLEHy3WTkFOYwySGqAzVKdIYNG1bpX0iJRAIjIyPY29tj5MiR6NBBdRfN4tAVEYlJV0cX73Z6F+84voP4f+LxSuNXZOcu3ruIz45/htCeoXC3dhcxSiL1V6NEx9zcHHv27EGjRo3g6uoKQRCQnJyM3Nxc+Pj4IC4uDuHh4Thy5Ah69uyp7JiJFMZtJUjVSSQS+LT1kWsLTwzHrsu7sOvyLvRt0xehvULRt01f9vgQ1UCNEp3mzZtj5MiRWLNmjWxaeXl5OWbMmIGGDRti27ZtCAwMxLx585CQkKDUgIkUUZ2iZYCFy6RaPvL6CLo6uvjh4g84kn4ER9KPwK2lG0J7hsKvox90dXTFDpFIbdSoGLlp06ZITEyssGP41atX4enpiZycHFy8eBFeXl7Izc1VVqxK9e8anatXr7IYWQMpUrQMsHCZVNeN3BtY/vtyfH3haxQ9LQIA9G7dG8fHVL60B5E2qdNi5KdPn+LKlSsVEp0rV67IinuNjIxUupuVNTqaj0XLpO5aN2qNlb4rsbD3Qqw+uxqrz66Gr72v7Hy5UI7C0kKuxUP0AjVKdPz9/TF+/Hh89NFHcHd3h0QiwdmzZ/HVV19h9OjRAIDjx4+jU6dOSg2WiEgbNTVtis+9P8cczznQkfxvFfo9V/Zgwt4JmNZ9Gqb1mAZLE0sRoyRSTTVKdFasWAErKyssXboU9+7dA/BsR/GZM2di3rx5AAAfHx+8+eabyouUiEjLNTRsKHcc91ccHj15hM9PfI5lvy/DxG4TEeIRAltzW5EiJFI9tV4wUCqVAkC16ltUAWt06Lnn/y+jRnSBfbMXDwFwdhapkrLyMuy+shthCWG4kHUBAKCno4cPOn+AeT3nwaGpg8gREtUdRWt0apXo3L9/H3///TckEgk6dOgAS0v16zblysh0O7cI/SKPo6j05YtHcnYWqSJBEHD4n8MISwjD0YyjAIDu1t1xZsIZkSMjqjt1WoxcUFCAadOm4ZtvvkF5eTkAQFdXF6NHj8bq1athYmJSs6iJRGDdyBiHZ/VReFuJRwUlTHRIpUgkErzR9g280fYNnLl1BuGJ4fB39pedz3uShzO3z+CNV95Q6UkiRHWhRolOSEgIjh8/jl9++UW2IGBCQgKmT5+OWbNmITY2VqlBEtU1ztAiTdGjVQ/sGrFLri32fCzmH5mPbi26IbRnKN5yeItr8ZDW0Hn5JRXt3LkTGzduhK+vL8zMzGBmZoYBAwbg66+/xk8//aTsGImIqBZKy0phom+CC1kX8O5P78Ih2gEbLmxA8dNisUMjqnM1SnQKCwthZWVVob1Zs2YoLCysdVD1ITo6Go6OjnB35z4yVD1p2fm4dDvvpV+3c4vEDpUIAPBxn49xI/gGPun9CRobNca1h9cw8ZeJaLOyDaJOR4kdHlGdqlExct++fdGkSRN88803MDIyAgAUFRUhICAADx8+xOHDh5UeaF1hMTIpqjpFywALl0k15Zfk4+ukrxH5eyRuP76NDzp/gO/e+k7ssIiqrU6LkVeuXIk333wTrVq1gouLCyQSCVJSUmBkZIRDhw7VOGgiVaZo0TLAwmVSXQ0MGmCmx0xMcZ+C7y9+j1dbvSo7dyXnCmLOxWC252yuxUMao0aJjpOTE65du4bvvvsOV65cgSAIeO+99/DBBx/A2Jg/1ElzsWiZNIWhniHGdR0n1xaeGI4tKVsQez4WIzuPxFzPuejUjCvck3qrUaIDAMbGxpg4caIyYyEiIhH5O/vjZt5NHEk/gm/++Abf/PENhnYYitBeoXI9P0TqROFEZ+/evQrfdMiQITUKhkjTpGXnv/QarrZMquL1Nq/j9Tav49ztcwhPDMeuy7vw898/4+e/f8bbjm9jxzs7xA6RqNoUTnT8/PwUuk4ikch2MCfSVo1NDWCsr4vguJSXXsuiZVI17tbu+Ondn3Al5woiEiPw7Z/fwrmZs+y8IAgoF8q5Fg+pBYUTnecrIGuKf+91RaRsXG2ZNEFHy47YOHQjFnkvQgOD/+0Dt+/qPoT8FoK5nnMx2mU0DPUMRYyS6MVqXKOj7oKCghAUFCSbnkakbCxcJk3RyqyV3PHapLVIe5iGD/d9iE+PfYqZr87EJLdJMDNUr82dSTsonOisWrVK4ZtOnz69RsEQEZHqi3s7DhsubEDk75G4Jb2FuYfn4quErxDkHoTpPaajmWkzsUMkklE40VmxYoXc8f3791FYWIhGjRoBAHJzc2FiYoJmzZox0SGqJkWKlgEWLpNqaGDQAMGvBmOK+xT8cPEHhCeG40rOFXx58kv8X/r/4dT4U2KHSCSjcKKTnp4u++8ffvgBMTEx2LhxIzp06AAA+PvvvzFx4kRMmjRJ+VESaajqFC0DLFwm1WKga4AxXcZgtMto/HzlZ4QlhGFGjxmy8/kl+cjIzYBTMycRoyRtV6MtINq2bYuffvoJXbt2lWtPSkrC22+/LZcUqTpuAUFiu51bVK3VlvnOkap6/utEIpEAACJPRWJ2/GwMbj8Y83vNh4eNh5jhkYap0y0gsrKyUFpaWqG9rKwM9+7dq8ktibQWi5ZJUzxPcJ7LyM2ABBL8cvUX/HL1F/Ru3Rvze81H/7b9K1xLVFdqtHt53759MXHiRJw/f16WwZ8/fx6TJk1Cv379lBogERGpp9UDVuNy0GWM7zoe+jr6OHHjBHy/90XXdV2x/a/tYodHWqJGic6mTZtgbW2N7t27w8jICIaGhujRowdatGiBDRs2KDtGIiJSUx0sO2DDkA1In5GOkFdDYKpvij/u/YGdl3eKHRppiRoNXTVt2hQHDhzA1atXZZt6Ojg4oH379sqOr85wwUBSV9xWgtSRtZk1IvtHYkHvBYg+G43BHQbLzqU9TMPO1J0IdAuEuRHrz0i5qlWM7OHhAT8/PwwZMgQODg51GVe9YTEyqYvbuUXoF3kcRaUvT845O4vUyaRfJmH9hfUwMzRDkHsQZvSYAasGVmKHRSquToqRAwMDsXfvXixevBgtWrTA0KFDMWTIEPTq1YuFZUR1jNtKkKbqY9cHJzNP4nLOZYQlhGHF6RUY12UcZnvORpvGbcQOj9RctWp0AgICsHPnTuTk5CAqKgpSqRQjRoxAs2bNMGbMGOzevRuFhYV1FSuR1rNuZAwna/MXftk3a/DyGxGpkJGdR+LSlEvYM2IPelj3wJOnTxBzPgbtVrfD1ANTxQ6P1FyNipENDQ0xYMAArFu3Dnfu3MG+fftgbW2NTz75BJaWlhg0aBASExOVHSsREWkoHYkOhnYcit/H/46jAUfh09YHZUIZGhs1ll1Tg2XfiJSzqWePHj3Qo0cPfPnll7h+/Tr27t2LrKwsZdyaiGqIRcukjiQSCV6zew2v2b2GC1kXYGtuKzv32/XfsPjkYoT2DMWAdgNYMkEKUfru5W3btsXMmTOVfVsiUlB1tpVg0TKpsm4tuskdR/4eiYTMBAzKHITOzTojtFco3u30LvR0lP6rjDRItd6ONm3avDSDlkgkuH79eq2CIqKaY9Eyaaotfluw4vcVWJu0FhezL+KDXR9g4f8txBzPORjTZQyM9fkeU0XVSnSCg4OrPJeRkYF169ahuLi4tjERUS1xWwnSRC0btkSETwQ+8voI0eeisfLMSqTnpmPKgSnYkboD/xfwf2KHSCqoWonOjBkzKrQ9fPgQX3zxBWJjY9GjRw+Eh4crLTgiIqL/amzcGAt7L0SIRwg2XtiIZb8vw9guY2XnC0sLIS2WonmD5iJGSaqixgObRUVFWL58OSIiImBnZ4ddu3ZhwIAByoxNIY8fP8brr7+O0tJSlJWVYfr06Zg4cWK9x0GkrhQpWgZYuEyqx0TfBNN6TEOgW6BcWcXGCxsxJ34OxnV9thbPK41fETFKElu1E52ysjJ8/fXXWLRoEYyMjLB69WqMGjVKtOp3ExMTHD9+HCYmJigsLISTkxPeeustNGnSRJR4iNRFdYqWARYuk+rS19WXOz6ZeRLFZcWIPR+LdUnr8J7Te5jXcx6crZxFipDEVK1EZ/v27Vi4cCHy8vLw0UcfYfLkyTAwMKir2BSiq6sLExMTAMCTJ09QVlbGtRaIFKBo0TLAwmVSL3FvxyHoRhCWJC7Br2m/4oeLP+CHiz9gQLsBCO0ZCq/WXmKHSPWoWonOe++9B2NjY7z//vu4ceMGQkNDK71u+fLlCt/zxIkTiIiIQFJSErKysrB79274+fnJXRMTE4OIiAhkZWWhU6dOiIqKgpfX/17U3Nxc9OnTB9euXUNERAQsLS2r820RaS0WLZMmkkgk6GPXB33s+iA5KxnhieHYkboDB64dgIGuARMdLVOtRKd3795Knz5eUFAAFxcXjB07FsOHD69wPi4uDsHBwYiJiUHPnj2xbt06+Pr6IjU1Fba2zxaSatSoEf744w/cu3cPb731Ft5++21YWXFDOCJl4yKEpG66tuiKbW9vw+KHixGRGIHx3cbLzt3IvYGEzASMcBrBtXg0WLV2L69rEomkQo9Ojx490K1bN8TGxsraHBwc4Ofnh7CwsAr3mDx5Ml5//XW88847lT6juLhYbgq8VCqFjY0Ndy8negHunE6aaOqBqYg+Fw27RnaY4zkHY7uM5Vo8akTpu5eHhIQodJ1EIkFkZKSit32hkpISJCUlVRgi8/HxwalTpwAA9+7dg7GxMczMzCCVSnHixAlMnjy5ynuGhYVh0aJFSomPSFtwEULSRG0atUFTk6bIyM1A0IEgLDq+CME9gjHZfTIaGTUSOzxSEoUTneTkZLnjpKQklJWVoUOHDgCAq1evQldXF66urkoLLicnB2VlZRWGoaysrHD37l0AwK1btzB+/HgIggBBEDB16lQ4O1ddWT9//ny5pO15jw4RvRjreUjTzPKchcnuk7E5eTMiTkXgRt4NfPR/HyEsIQxze87Fwt4LxQ6RlEDhROfo0aOy/16+fDkaNmyIrVu3onHjZzvLPnr0CGPHjpUrElaW/05dFwRB1ubq6oqUlBSF72VoaAhDQ0NER0cjOjoaZWUv74onIiLNZKJvgqDuQfjQ9UPE/RWHJQlL8Nf9v5Bfotj6UqT6alR9FRkZid9++02W5ABA48aNsXjxYvj4+GDWrFlKCc7S0hK6urqy3pvnsrOza11sHBQUhKCgINkYHxEpDxchJHWjr6uPUc6jMLLzSBy4dgBuLd1k545lHMO6pHWY13MeujTvIl6QVCM1SnSkUinu3buHTp06ybVnZ2fj8ePHSgkMAAwMDODq6or4+HgMGzZM1h4fH4+hQ4cq7TlEpBxchJDUnY5EB4PaD5JrC0sIw2/Xf8O2S9vga++L0F6h8LL1Em2hXKqeGiU6w4YNw9ixYxEZGYlXX30VAHD69GnMmTMHb731VrXulZ+fj7S0NNlxeno6UlJSYGFhAVtbW4SEhMDf3x9ubm7w8PDA+vXrkZmZicDAwJqELsOhKyLl4yKEpInC+4XDwtgC2//ajoNpB3Ew7SA8WnkgtFcoBrUfBB2Jjtgh0osINVBQUCBMnjxZMDQ0FHR0dAQdHR3BwMBAmDx5spCfn1+tex09elQAUOErICBAdk10dLTQunVrwcDAQOjWrZtw/PjxmoRdqby8PAGAkJeXV63PXbyVK7Set0+4eCtXabEQaRP+HSJ1k/YgTQj8JVAw/MJQwGcQ8BmEYduGiR2W1lL093et1tEpKCjA9evXIQgC7O3tYWpqqpTkqz4pOg//v7iODlHt8O8Qqau7+XcRdToKMediEPVmFMZ1HQcAKH5ajDKhDCb6JiJHqB0U/f1dq/42U1NTODs7w8XFRe2SnOjoaDg6OsLd3V3sUIiISI00b9AcS/otQebMTIxyHiVr35yyGXZRdvjyxJd4VPRIxAjp37R2zWvOuiJSDdxWgtTVfxcV/Cn1J9wvvI+FRxciPDEck1wnIcQjBC0athAnQAKgxYkOEYmrOjO0ODuL1MGvo37F9r+2Y0nCElzMvohlvy/DqrOrMMZlDOb0nAN7C3uxQ9RKWpvocNYVkbi4rQRpGj0dPYzsPBLvO72P/df2Y0nCEiTeTMT6C+tx+/Ft7Bu5T+wQtZLWJjocuiISH7eVIE0kkUgwqP0gDGo/CAmZCQhLCENor//t2Xjn8R1cfXAVfVr34Vo89UBrEx0iUi9cbZnUUS/bXtg/cr9c27JTy7Di9Aq82upVhPYMxeAOg7kWTx1iokNEKo2rLZOm0ZXowlDXEKdvnYZfnB8cmzpiXs95eN/pfejr6osdnsap1To6moDr6BCpvtu5RdVabZl/L0nV3cu/92wtnvMxkBZLAQC25rZY6LUQE10nihydelD097fW9uiwGJlIfbCWhzSNVQMrhPV7Vruz9vxarDi9Apl5mbicc1ns0DSO1g4KBgUFITU1FefOnRM7FCIi0lLmRuaY12se0mekI3ZgLEI8QmTnTt08hdm/zcZt6W0RI1R/WpvoEBERqQpjfWMEugWilVkrWduXJ79E5O+ReGXVK5i4dyKuPbgmYoTqS2uHrohIc3G1ZdIEU92n4nHxY5zMPIkNyRuwMXkj3nZ8G/N6zoNrS1exw1MbTHSISGNwtWXSJL7tfOHbzheJmYlYkrgE+67uw47UHdiRugMTuk7A10O+FjtEtaC1iQ6LkYk0D1dbJk3U07YnfrH9BRfvXcTSU0vx48Uf5Xp0npY/hY5Eh2vxVEFrEx2ujEykmThDizRVZ6vO+HbYt/j8tc/lNgr95o9vEHEqAvN6zsPIziNhoGsgYpSqh+kfERGRGmnTuA2M9Ixkx+uT1uNKzhWM/Xks7FfZY+XplSgoKRAxQtWitT06RETcVoI0wW/+v2Hd+XVYfno5bkpvIvhQML448QWm95iOqd2nwsLYQuwQRcVEh4i0DreVIE1iZmiGOT3nYFqPadiashURpyJw/dF1fHrsU5y/cx57398rdoiiYqJDRFpH0aJl4H+Fy+fSH+JRswYvvJY9PyQmIz0jTHKbhPHdxmNn6k6EJYRh5qszZefvF9zHw6KH6GDZQcQo65/WJjqcdUWk3RQtWuaUdVI3ejp6GOE0Au92eleufdmpZYg4FYG3HN5CaK9QuLV0EynC+qW1iQ5nXRGRIjhlndSVRCKRO76TfwcCBOy8vBM7L+9Ev1f6IbRnKF5v83qFazWJ1iY6RESK4pR10gTfDvsW83rOQ3hiOH68+CMO/3MYh/85DPeW7ljYeyGGdBgidoh1gtPLiYiItIRTMyd8O+xbpE1Pw7Tu02CsZ4xzd87ht+u/iR1anWGiQ0REpGXsGtlhle8qZARnYIHXAsz2nC07l3QnCSt+X4H8EsWWX1B1THSIiIi0VDPTZlj8+mLYNbKTtS0+uRghv4WgdVRrLDq2CA8KH4gXoBIw0SEiIiKZge0Gwt7CHg+LHuKz45/BNsoWM3+diVvSW2KHViNMdIiIiEhmQrcJuBJ0BXFvx6Fr864oLC1E1JkovLLyFSw4skDs8KqNs66IiJSI20qQJtDV0cW7nd7FO47vIP6feCxJWIKjGUfRvEFz2TWCIKjFtHQmOkRESlCTbSXW+ruiiemLd5pmQkRikkgk8GnrA5+2Pjhz6ww6W3WWnfvuz++w9Y+tCO0Vir5t+qps0qO1iQ5XRiYiZarOthIPCkoQ+G0SAjadfem1XG2ZVEWPVj3kjpefXo6Uuyk4kn4Eri1cEdorFMM6DoOujq5IEVZOIgiCIHYQYnq+MnJeXh7MzMwU/tyl23kYtDoB+6b1gpM1V1Ymouq5nVuk8GrL/DlDqigzLxORpyLx9YWvUfS0CADQvkl7zPWci1HOo2CoZ1inz1f097fW9ugQEYmJqy2TurM1t8VK35VY2HshVp9djTVn1+Dqg6uY8MsEHLp+CNvf2S52iAA464qIiIhqoalpU3zu/TluBN9ApE8kWjZsifFdx8vO5z7JRU5hjmjxMdEhIiKiWmto2BAhHiH4Z/o/8GnrI2uPPBWJ1lGtMePgDGTmZdZ7XBy6IiJScYpMWefsLFIV/67NEQQBp2+fRmFpIVadXYWY8zGY13MeFr++uN7iYaJDRKSiqjNlnbOzSBVJJBL8Nuo3HEk/grCEMPxf+v/BuqF1vcbARIeISEUpOmX9+eysRwUlTHRI5UgkEvR7pR/6vdIPZ2+fRedmnV/+ISViokNEpMI4O4s0SXfr7vX+TBYjExERkcZS+0Tn5s2beO211+Do6AhnZ2fs2LFD7JCIiIhIRaj90JWenh6ioqLQpUsXZGdno1u3bhgwYABMTU3FDo2IiIhEpvaJTosWLdCiRQsAQLNmzWBhYYGHDx8y0SEircOd04kqEj3ROXHiBCIiIpCUlISsrCzs3r0bfn5+ctfExMQgIiICWVlZ6NSpE6KiouDl5VXhXufPn0d5eTlsbGzqKXoiIvHVZOd0TkUnbSF6olNQUAAXFxeMHTsWw4cPr3A+Li4OwcHBiImJQc+ePbFu3Tr4+voiNTUVtra2susePHiA0aNHY8OGDfUZPhGR6KqzczqnopO2ET3R8fX1ha+vb5Xnly9fjvHjx2PChAkAgKioKBw6dAixsbEICwsDABQXF2PYsGGYP38+PD09X/i84uJiFBcXy46lUqkSvgsiInFxGjpR5URPdF6kpKQESUlJCA0NlWv38fHBqVOnADxbXnrMmDF4/fXX4e/v/9J7hoWFYdGiRXUSLxGRuuC2EqQtVDrRycnJQVlZGaysrOTarayscPfuXQBAYmIi4uLi4OzsjD179gAAvv32W3TuXPnKi/Pnz0dISIjsWCqVsqaHiLQGt5UgbaPSic5zEolE7lgQBFlbr169UF5ervC9DA0NYWhoiOjoaERHR6OsrEypsRIRqTJuK0HaRqUTHUtLS+jq6sp6b57Lzs6u0MtTXUFBQQgKCoJUKoW5uXmt7kVEpE5Yz0PaRKVXRjYwMICrqyvi4+Pl2uPj419adExEREQkeo9Ofn4+0tLSZMfp6elISUmBhYUFbG1tERISAn9/f7i5ucHDwwPr169HZmYmAgMDa/VcDl0RERFpPtETnfPnz8Pb21t2/LxQOCAgAFu2bMGIESPw4MEDfP7558jKyoKTkxMOHDiA1q1b1+q5HLoiIno5rrZM6k70ROe1116DIAgvvGbKlCmYMmVKPUVERERcbZk0heiJjlg4dEVEVLWarLZ8Lv0hHjVr8MJr2fND9U1rEx0OXRERvZiis7O4Ng+pMq1NdIiISDm4Ng+pMq1NdDh0RUSkPFybh1SVSq+jU5eCgoKQmpqKc+fOiR0KERER1RGtTXSIiIhI8zHRISIiIo2ltYlOdHQ0HB0d4e7uLnYoREREVEe0NtFhjQ4REZHm09pZV0REJA5uK0H1iYkOERHVC24rQWJgokNERPWiJttKcHFBqi0mOkREVG+4sCDVN61NdLgyMhGR6lOknoe1PPQiWpvocFNPIiLVxY1CSVm0NtEhIiLVxY1CSVmY6BARkUqqTj0Pp6xTVZjoEBGR2uKUdXoZJjpERKS2OGWdXkZrEx3OuiIi0gycsk4vorWJDmddERFRVW7nFinUS8SaH9WntYkOERFpp5cVLj8oKEHgt0koKn15jz9rflQfEx0iItIK1V2bZ+u47mhialDlNaz5UQ9MdIiISCtUp3CZQ1Kag4kOERFpDRYuax8mOkRERLXAxQpVGxMdIiKiGuBiheqBiQ4REVENcLFC9aC1iQ4XDCQiotpizY/q0xE7ALEEBQUhNTUV586dEzsUIiIiqiNam+gQERGR5mOiQ0RERBqLiQ4RERFpLCY6REREpLGY6BAREZHGYqJDREREGouJDhEREWksJjpERESksTQi0Rk2bBgaN26Mt99+W+xQiIiISIVoxBYQ06dPx7hx47B161axQyEiIqqSIjudc5dz5dKIRMfb2xvHjh0TOwwiIqJKVWenc+5yrlyiJzonTpxAREQEkpKSkJWVhd27d8PPz0/umpiYGERERCArKwudOnVCVFQUvLy8xAmYiIiomhTd6byudjm/nVuk0C7rgOb1KIme6BQUFMDFxQVjx47F8OHDK5yPi4tDcHAwYmJi0LNnT6xbtw6+vr5ITU2Fra2tCBETERFVn1g7nd/OLUK/yOMoKi1T6HpN61ESPdHx9fWFr69vleeXL1+O8ePHY8KECQCAqKgoHDp0CLGxsQgLC6v284qLi1FcXCw7lkql1Q+aiIhITTwqKEFRaRmiRnSBfbMGL7y2rnqUxCR6ovMiJSUlSEpKQmhoqFy7j48PTp06VaN7hoWFYdGiRcoIj4iISG3YN2sAJ2tzscOodyo9vTwnJwdlZWWwsrKSa7eyssLdu3dlx/3798c777yDAwcOoFWrVjh37lyV95w/fz7y8vJkXzdv3qyz+ImIiEhcKt2j85xEIpE7FgRBru3QoUMK38vQ0BCGhoZKi42IiIhUl0r36FhaWkJXV1eu9wYAsrOzK/TyVFd0dDQcHR3h7u5eq/sQERGR6lLpRMfAwACurq6Ij4+Xa4+Pj4enp2et7h0UFITU1NQXDnMRERGRehN96Co/Px9paWmy4/T0dKSkpMDCwgK2trYICQmBv78/3Nzc4OHhgfXr1yMzMxOBgYG1em50dDSio6NRVqbYdDsiIiJSP6InOufPn4e3t7fsOCQkBAAQEBCALVu2YMSIEXjw4AE+//xzZGVlwcnJCQcOHEDr1q1r9dygoCAEBQVBKpXC3Fz7qtCJiIi0geiJzmuvvQZBEF54zZQpUzBlypR6ioiIiIg0heiJDhEREclTZPNPoO62a9CkzUe1NtFhjQ4REama6mz+CSh/uwZN3HxUaxMd1ugQEZGqUXTzT6ButmsQe/PRuqC1iQ4REZEqqu7mny8bZlJ0GKymz1d1WpvocOiKiIjUWXWHmRqbGtR9UCpIaxMdDl0REZE6q84wl7oUDtcFrU10iIiI1J2mDTPVBZXeAoKIiIioNrQ20eGmnkRERJpPaxMdbupJRESk+bQ20SEiIiLNx0SHiIiINBYTHSIiItJYWpvosBiZiIhI82ltosNiZCIiIs2ntYkOERERaT4mOkRERKSxmOgQERGRxmKiQ0RERBqLiQ4RERFpLK1NdDi9nIiISPNpbaLD6eVERESaT2sTHSIiItJ8THSIiIhIYzHRISIiIo3FRIeIiIg0FhMdIiIi0lhMdIiIiEhj6YkdABEREamntOx8ha5rbGoA60bGdRxN5bQ20YmOjkZ0dDTKysrEDoWIiEitNDY1gLG+LoLjUhS63lhfF4dn9REl2dHaRCcoKAhBQUGQSqUwNzcXOxwiIiK1Yd3IGIdn9cGjgpKXXpuWnY/guBQ8KihhokNERETqwbqRsWjDUdXBYmQiIiLSWEx0iIiISGMx0SEiIiKNxUSHiIiINBYTHSIiItJYTHSIiIhIYzHRISIiIo2lEYnOvn370KFDB7Rr1w4bNmwQOxwiIiJSEWq/YODTp08REhKCo0ePwszMDN26dcNbb70FCwsLsUMjIiIikal9j87Zs2fRqVMnWFtbo2HDhhgwYAAOHTokdlhERESkAkRPdE6cOIHBgwejZcuWkEgk2LNnT4VrYmJi0KZNGxgZGcHV1RUnT56Unbtz5w6sra1lx61atcLt27frI3QiIiJScaInOgUFBXBxccGaNWsqPR8XF4fg4GAsWLAAycnJ8PLygq+vLzIzMwEAgiBU+IxEIqnTmImIiEg9iF6j4+vrC19f3yrPL1++HOPHj8eECRMAAFFRUTh06BBiY2MRFhYGa2truR6cW7duoUePHlXer7i4GMXFxbLjvLw8AIBUKq1W3PmPpSgvLkT+YymkUiZWRERElamr35fPf29X1uEhR1AhAITdu3fLjouLiwVdXV1h165dctdNnz5d6N27tyAIglBaWirY29sLt27dEqRSqWBvby/k5ORU+YxPP/1UAMAvfvGLX/ziF7804OvmzZsvzC1E79F5kZycHJSVlcHKykqu3crKCnfv3gUA6OnpITIyEt7e3igvL8fcuXPRpEmTKu85f/58hISEyI7Ly8vh6uqKCxcuVDrk5e7ujnPnzlVol0qlsLGxwc2bN2FmZlbTb7FeVfW9qOozanqv6nxO0WsVue5F1/A9Eu8ZfI9UE9+jml/L9+gZQRDg6uqKli1bvvA6lU50nvtvAiIIglzbkCFDMGTIEIXuZWhoCENDwwpt5ubmlV6vq6v7wv/hZmZmavFCAC//XlTtGTW9V3U+p+i1ilz3omv4Hon3DL5HqonvUc2v5Xv0PwYGBtDReXG5sejFyC9iaWkJXV1dWe/Nc9nZ2RV6eWojKCioRufUTX18L8p8Rk3vVZ3PKXqtItfxPVLNZ/A9Uk18j2p+Ld+j/1Hke5H8/9oYlSCRSLB79274+fnJ2nr06AFXV1fExMTI2hwdHTF06FCEhYWJEOUzUqkU5ubmyMvLU5vMl1QP3yNSBr5HpAya+h6JPnSVn5+PtLQ02XF6ejpSUlJgYWEBW1tbhISEwN/fH25ubvDw8MD69euRmZmJwMBAEaN+Ntz16aefVhgGI6oOvkekDHyPSBk09T0SvUfn2LFj8Pb2rtAeEBCALVu2AHi2YODSpUuRlZUFJycnrFixAr17967nSImIiEjdiJ7oEBEREdUVlS5GJiIiIqoNJjpERESksZjoEBERkcZiokNEREQai4lOPRg2bBgaN26Mt99+W+xQSI3s27cPHTp0QLt27bBhwwaxwyE1xZ8/VFs3b97Ea6+9BkdHRzg7O2PHjh1ih1QtnHVVD44ePYr8/Hxs3boVP/30k9jhkBp4+vQpHB0dcfToUZiZmaFbt244c+YMLCwsxA6N1Ax//lBtZWVl4d69e+jSpQuys7PRrVs3/P333zA1NRU7NIWwR6ceeHt7o2HDhmKHQWrk7Nmz6NSpE6ytrdGwYUMMGDAAhw4dEjssUkP8+UO11aJFC3Tp0gUA0KxZM1hYWODhw4fiBlUNWp/onDhxAoMHD0bLli0hkUiwZ8+eCtfExMSgTZs2MDIygqurK06ePFn/gZJaqe17defOHVhbW8uOW7Vqhdu3b9dH6KRC+POJlEGZ79H58+dRXl4OGxubOo5aebQ+0SkoKICLiwvWrFlT6fm4uDgEBwdjwYIFSE5OhpeXF3x9fZGZmSm7xtXVFU5OThW+7ty5U1/fBqmY2r5XlY0oSySSOo2ZVI8yfj4RKes9evDgAUaPHo3169fXR9jKI5AMAGH37t1ybd27dxcCAwPl2jp27CiEhoZW695Hjx4Vhg8fXtsQSQ3V5L1KTEwU/Pz8ZOemT58ufP/993UeK6mu2vx84s8feq6m79GTJ08ELy8v4ZtvvqmPMJVK63t0XqSkpARJSUnw8fGRa/fx8cGpU6dEiorUnSLvVffu3XHp0iXcvn0bjx8/xoEDB9C/f38xwiUVxZ9PpAyKvEeCIGDMmDF4/fXX4e/vL0aYtSL67uWqLCcnB2VlZbCyspJrt7Kywt27dxW+T//+/XHhwgUUFBSgVatW2L17N9zd3ZUdLqkJRd4rPT09REZGwtvbG+Xl5Zg7dy6aNGkiRrikohT9+cSfP/QiirxHiYmJiIuLg7Ozs6y+59tvv0Xnzp3rO9waYaKjgP/WRgiCUK16Cc6Wocq87L0aMmQIhgwZUt9hkZp52XvEnz+kiBe9R7169UJ5ebkYYSkFh65ewNLSErq6uhV6b7Kzsytkv0SK4ntFysD3iJRBG94jJjovYGBgAFdXV8THx8u1x8fHw9PTU6SoSN3xvSJl4HtEyqAN75HWD13l5+cjLS1Ndpyeno6UlBRYWFjA1tYWISEh8Pf3h5ubGzw8PLB+/XpkZmYiMDBQxKhJ1fG9ImXge0TKoPXvkahzvlTA0aNHBQAVvgICAmTXREdHC61btxYMDAyEbt26CcePHxcvYFILfK9IGfgekTJo+3vEva6IiIhIY7FGh4iIiDQWEx0iIiLSWEx0iIiISGMx0SEiIiKNxUSHiIiINBYTHSIiItJYTHSIiIhIYzHRISIiIo3FRIeIiIg0FhMdIlJZDx48QLNmzZCRkSF2KFWaPXs2pk+fLnYYRFQFJjpEVG/GjBkDiURS6WaBU6ZMgUQiwZgxY2RtYWFhGDx4MOzs7AAAGRkZkEgk0NPTw+3bt+U+n5WVBT09PUgkknpNjObOnYvNmzcjPT293p5JRIpjokNE9crGxgbbtm1DUVGRrO3Jkyf48ccfYWtrK2srKirCxo0bMWHChAr3aNmyJb755hu5tq1bt8La2rruAq9Cs2bN4OPjg7Vr19b7s4no5ZjoEFG96tatG2xtbbFr1y5Z265du2BjY4OuXbvK2g4ePAg9PT14eHhUuEdAQAA2b94s17ZlyxYEBATItT169AgffPABmjZtCmNjY7Rr1072uWPHjkEikSA3N1d2fUpKilyP0JYtW9CoUSMcOnQIDg4OaNCgAd58801kZWXJPWfIkCH48ccfa/TnQUR1i4kOEdW7sWPHyiUqmzZtwrhx4+SuOXHiBNzc3Cr9/JAhQ/Do0SMkJCQAABISEvDw4UMMHjxY7rqPP/4YqampOHjwIC5fvozY2FhYWlpWK9bCwkIsW7YM3377LU6cOIHMzEzMnj1b7pru3bvj5s2buHHjRrXuTUR1j4kOEdU7f39/JCQkICMjAzdu3EBiYiJGjRold01GRgZatmxZ6ef19fUxatQobNq0CcCzRGnUqFHQ19eXuy4zMxNdu3aFm5sb7Ozs0K9fvwrJ0MuUlpZi7dq1cHNzQ7du3TB16lQcOXJE7prnQ2aqXDRNpK30xA6AiLSPpaUlBg4ciK1bt0IQBAwcOLBCT0tRURGMjIyqvMf48ePh4eGBr776Cjt27MDvv/+Op0+fyl0zefJkDB8+HBcuXICPjw/8/Pzg6elZrVhNTEzQtm1b2XGLFi2QnZ0td42xsTGAZ70/RKRa2KNDRKIYN24ctmzZgq1bt1YYtgKeJUOPHj2q8vNOTk7o2LEj3n//fTg4OMDJyanCNb6+vrhx4waCg4Nx584d9O3bVzbspKPz7MefIAiy60tLSyvc47+9RBKJRO4zAPDw4UMAQNOmTauMl4jEwUSHiETx5ptvoqSkBCUlJejfv3+F8127dkVqauoL7zFu3DgcO3as0kTpuaZNm2LMmDH47rvvEBUVhfXr18vaAcgVFqekpNTgOwEuXboEfX19dOrUqUafJ6K6w0SHiEShq6uLy5cv4/Lly9DV1a1wvn///vjrr79e2KszceJE3L9/v9Ip6ADwySef4Oeff0ZaWhr++usv7Nu3Dw4ODgAAe3t72NjY4LPPPsPVq1exf/9+REZG1uh7OXnyJLy8vGRDWESkOpjoEJFozMzMYGZmVum5zp07w83NDdu3b6/y83p6erC0tISeXuXlhgYGBpg/fz6cnZ3Ru3dv6OrqYtu2bQCeDUn9+OOPuHLlClxcXBAeHo7FixfX6Pv48ccfMXHixBp9lojqlkT472AzEZGKOHDgAGbPno1Lly7JampUzf79+zFnzhz8+eefVSZcRCQe/q0kIpU1YMAAXLt2Dbdv34aNjY3Y4VSqoKAAmzdvZpJDpKLYo0NEREQaSzX7gomIiIiUgIkOERERaSwmOkRERKSxmOgQERGRxmKiQ0RERBqLiQ4RERFpLCY6REREpLGY6BAREZHGYqJDREREGuv/AerpyE02qyQVAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"integrand = lambda x: x**(-2.3)\n",
"norm = quad(integrand,0.1,150)\n",
"print (norm[0])\n",
"\n",
"def inverse(y,a = 2.3):\n",
" m = ((a-1)*norm[0]*y + (0.1)**(1-a))**(-1/1.3) \n",
" return m\n",
"\n",
"def salpeter(x,a=2.3):\n",
" f = 10**5*x**(-a)\n",
"# print('function',f)\n",
" return f\n",
"\n",
"\n",
"N = 10**6\n",
"mass = []\n",
"for i in range(N):\n",
" y=np.random.normal(0,1)\n",
" m=inverse(y)\n",
" mass.append(m)\n",
" \n",
"\n",
"x=np.linspace(0.1,150)\n",
"plt.hist(mass,bins=np.logspace(-1,2.15,num=45),histtype='step',label='Monte Carlo sampling')\n",
"plt.xlim(8*10**-2,2*10**2)\n",
"plt.ylim(10**0,10**7)\n",
"plt.xscale('log')\n",
"plt.yscale('log')\n",
"plt.xlabel('M(Msun)')\n",
"plt.plot(x,salpeter(x,a = 1.3),'g--',label='x^-2.3')\n",
"plt.ylabel('dN/dlogM')\n",
"plt.legend()"
]
},
{
"cell_type": "markdown",
"id": "150abceb",
"metadata": {},
"source": [
"## Box_Muller Method\n",
"If we have two random numbers X and Y obtained from a Guassian Distribution with thr same Std deviation $\\sigma$. the probability that thr point with position vector (x,y) falls in some small element $dxdy$ is\n",
"\n",
"\\begin{equation}\n",
" p(x)dxp(y)dy = \\frac{1}{2\\pi\\sigma^{2}}dxexp\\left(\\frac{-(x^2+y^2)}{2\\sigma^2}\\right)dxdy\n",
"\\end{equation}\n",
"Changing the coordinate to $(r,\\theta)$ we cann integrate the fuctiions seperately and invert it,so the probabiity we get will be a normal distribution between (0,1).\n",
"### The algorithm.\n",
"So, we need to chose two uniform random number z1 and z2.<br> \n",
"The sampling for $r$ and $\\theta$ will respectivly be:<br>\n",
"\\begin{equation}\n",
" r = \\sqrt{-2\\sigma ln(1-z_1)}\\\\\n",
" \\theta = 2\\pi z_2\n",
"\\end{equation}\n",
"and then $x = r cos(\\theta) $ and $y = r sin(\\theta)$"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "940344ce",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.5360599421957188\n"
]
}
],
"source": [
"for i in range(100):\n",
"z1 = np.random.normal(loc=0.0,scale =2.0)\n",
"print (z1)\n",
"z2 = np.random.normal(0,1)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0016f191",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}