diff --git a/navipy/processing/mcode.py b/navipy/processing/mcode.py
index 008ba442636c7bc391e2036220276c4f685190f8..f51d6896e20fec59b5989b9f5a5c26c119b79737 100644
--- a/navipy/processing/mcode.py
+++ b/navipy/processing/mcode.py
@@ -94,16 +94,13 @@ def optic_flow_rotationonal(viewing_directions, velocity):
dyaw, dpitch, droll, convention)
M = compose_matrix(angles=[yaw, pitch, roll], translate=None,
perspective=None, axes=convention)[:3, :3]
+ angvel_bee = np.dot(M, angvel)
# project it on the eye
for i, a in enumerate(azimuth):
for j, e in enumerate(elevation):
# in the bee coordinate system
spline = np.array(spherical_to_cartesian(e, a))
- # in the global coordinate system
- spline = np.dot(M.transpose(), spline)
- opticFlowR = -np.cross(angvel, spline)
- # in the bee coordinate system
- opticFlowR = np.dot(M, opticFlowR)
+ opticFlowR = -np.cross(angvel_bee, spline)
# Decompose into el, az
(OF_rho, OF_phi, OF_epsilon) = \
cartesian_to_spherical_vectors(opticFlowR, [a, e])
@@ -155,20 +152,14 @@ def optic_flow_translational(scene, viewing_directions,
# transformation
M = compose_matrix(angles=[yaw, pitch, roll], translate=None,
perspective=None, axes=convention)[:3, :3]
-
+ u_bee = M.dot(u)
for i, a in enumerate(azimuth):
for j, e in enumerate(elevation):
# The spline express in the bee coordinate system
spline = np.array(spherical_to_cartesian(e, a))
- # in the global coordinate system
- spline = np.dot(M.transpose(), spline)
# the Translation-part of the Optic Flow:
- dotvu = v*u
- #opticFlowT = -np.cross(spline, dotvu)
- splinetrans = np.transpose(spline)
- opticFlowT = -(dotvu-(np.dot(dotvu, splinetrans))*spline)
- # in the bee coordinate system
- opticFlowT = np.dot(M, opticFlowT)
+ dotvu = v*u_bee
+ opticFlowT = -np.cross(np.cross(spline, dotvu), spline)
# Decompose into el, az
(OF_rho, OF_phi, OF_epsilon) = \
cartesian_to_spherical_vectors(opticFlowT, [a, e])
@@ -194,8 +185,8 @@ def optic_flow(scene, viewing_directions,
dy, dz, dalpha, dbeta, dgamma)
: distance channel: distance"""
rofr, hofr, vofr = optic_flow_rotationonal(viewing_directions, velocity)
- roft, hoft, voft = optic_flow_translational((scene, viewing_directions,
- velocity, distance_channel))
+ roft, hoft, voft = optic_flow_translational(scene, viewing_directions,
+ velocity, distance_channel)
return rofr+roft, hofr+hoft, vofr+voft
diff --git a/navipy/processing/test_OpticFlow.py b/navipy/processing/test_OpticFlow.py
deleted file mode 100644
index d5033c8ab7fca1870d04509d8de6754ad01718e2..0000000000000000000000000000000000000000
--- a/navipy/processing/test_OpticFlow.py
+++ /dev/null
@@ -1,673 +0,0 @@
-import unittest
-import navipy.processing.mcode as mcode
-import pandas as pd
-import numpy as np
-from navipy.moving.agent import posorient_columns
-from navipy.moving.agent import velocities_columns
-from navipy.scene import __spherical_indeces__
-
-
-def Scale(data, oldmin, oldmax, mini, maxi, ran):
- """Scales all values in data into the range mini, maxi
- Input:
- - data: values to be scaled
- - mini: minimum of new range
- - maxi: maximum of new range
- - ran: range of new values (ran=maxi-mini)
-
- Output:
- - data: the scaled data"""
- data = (data-oldmin)/oldmax
- dmax = np.max(data)
- dmin = np.min(data)
- scaling = ran/(dmax-dmin)
- data = (data-dmin)*scaling
- data = data+mini
- return data
-
-
-class TestCase(unittest.TestCase):
- def setUp(self):
- """sets up a distance array that is used in some tests,
- as well as the velocity and viewing direction frames
- """
- distance = np.array([[2.60009956, 2.60009956, 2.60009956],
- [2.60009956, 2.60009956, 2.60009956],
- [2.60009956, 2.60009956, 2.60009956]])
-
- self.scene = np.zeros((3, 3, 4, 1))
- self.scene[:, :, 0, 0] = np.array(distance).copy()
- self.scene[:, :, 1, 0] = np.array(distance).copy()
- self.scene[:, :, 2, 0] = np.array(distance).copy()
- self.scene[:, :, 3, 0] = np.array(distance).copy()
- self.convention = 'xyz'
- tuples = [('location', 'x'), ('location', 'y'),
- ('location', 'z'), ('location', 'dx'),
- ('location', 'dy'), ('location', 'dz'),
- (self.convention, 'alpha_0'), (self.convention, 'alpha_1'),
- (self.convention, 'alpha_2'), (self.convention, 'dalpha_0'),
- (self.convention, 'dalpha_1'), (self.convention, 'dalpha_2')]
- self.index = pd.MultiIndex.from_tuples(tuples,
- names=['position',
- 'orientation'])
- self.velocity = pd.Series(index=self.index)
-
- self.elevation = np.arange(-np.pi/2, 4*(np.pi/360)-np.pi/2,
- 2*(np.pi/360))
- self.azimuth = np.arange(0, 4*(np.pi/360)+2*(np.pi/360), 2*(np.pi/360))
-
- self.viewing_directions = np.zeros((3, 2))
- self.viewing_directions[:, 0] = self.elevation
- self.viewing_directions[:, 1] = self.azimuth
-
- def test_wrong_convention(self):
- """ this test checks if an error is raised
- when none or a convention other than 'xyz'
- is used
- """
- velocity = pd.Series(index=['x', 'y', 'z',
- 'alpha_0', 'alpha_1', 'alpha_2',
- 'dx', 'dy', 'dz',
- 'dalpha_0', 'dalpha_1', 'dalpha_2'])
-
- positions = np.array([[-20, -20, 2.6, 1.57079633, 0, 0],
- [-19.8, -20, 2.6, 1.57079633, 0, 0]])
-
- x = positions[0, 0]
- y = positions[0, 1]
- z = positions[0, 2]
- yaw = positions[0, 3]
- pitch = positions[0, 4]
- roll = positions[0, 5]
-
- dx = positions[1, 0] - positions[0, 0]
- dy = positions[1, 1] - positions[0, 1]
- dz = positions[1, 2] - positions[0, 2]
- dyaw = positions[1, 3] - positions[0, 3]
- dpitch = positions[1, 4] - positions[0, 4]
- droll = positions[1, 5] - positions[0, 5]
-
- velocity.x = x
- velocity.y = y
- velocity.z = z
- velocity.alpha_0 = yaw
- velocity.alpha_1 = pitch
- velocity.alpha_2 = roll
- velocity.dx = dx
- velocity.dy = dy
- velocity.dz = dz
- velocity.dalpha_0 = dyaw
- velocity.dalpha_1 = dpitch
- velocity.dalpha_2 = droll
-
- with self.assertRaises(Exception):
- rof, hof, vof = mcode.optic_flow(self.scene,
- self.viewing_directions,
- velocity)
- for convention in ['alsjf', '233', 9, -1]:
- tuples = [('location', 'x'), ('location', 'y'),
- ('location', 'z'), ('location', 'dx'),
- ('location', 'dy'), ('location', 'dz'),
- (convention, 'alpha_0'),
- (convention, 'alpha_1'),
- (convention, 'alpha_2'),
- (convention, 'dalpha_0'),
- (convention, 'dalpha_1'),
- (convention, 'dalpha_2')]
- index = pd.MultiIndex.from_tuples(tuples,
- names=['position',
- 'orientation'])
- velocity = pd.Series(index=index)
- velocity['location']['x'] = x
- velocity['location']['y'] = y
- velocity['location']['z'] = z
- velocity[convention]['alpha_0'] = yaw
- velocity[convention]['alpha_1'] = pitch
- velocity[convention]['alpha_2'] = roll
- velocity['location']['dx'] = dx
- velocity['location']['dy'] = dy
- velocity['location']['dz'] = dz
- velocity[convention]['dalpha_0'] = dyaw
- velocity[convention]['dalpha_1'] = dpitch
- velocity[convention]['dalpha_2'] = droll
- with self.assertRaises(ValueError):
- rof, hof, vof = mcode.optic_flow(self.scene,
- self.viewing_directions,
- velocity)
-
- def test_rotation_not_too_strong(self):
- """
- this test checks that a Value Error is throught,
- if the change in rotation is more than pi/2.
- Thoughs, if the velocity of yaw, pitch and roll is
- is smaller than pi/2
- """
- # dyaw is too big
- positions = np.array([[-20, -20, 2.6, 1.57079633, 0, 0],
- [-19.8, -20, 2.6,
- 1.57079633 + np.pi/2 + 0.00001, 0, 0]])
-
- x = positions[0, 0]
- y = positions[0, 1]
- z = positions[0, 2]
- yaw = positions[0, 3]
- pitch = positions[0, 4]
- roll = positions[0, 5]
-
- dx = positions[1, 0] - positions[0, 0]
- dy = positions[1, 1] - positions[0, 1]
- dz = positions[1, 2] - positions[0, 2]
- dyaw = positions[1, 3] - positions[0, 3]
- dpitch = positions[1, 4] - positions[0, 4]
- droll = positions[1, 5] - positions[0, 5]
-
- self.velocity['location']['x'] = x
- self.velocity['location']['y'] = y
- self.velocity['location']['z'] = z
- self.velocity[self.convention]['alpha_0'] = yaw
- self.velocity[self.convention]['alpha_1'] = pitch
- self.velocity[self.convention]['alpha_2'] = roll
- self.velocity['location']['dx'] = dx
- self.velocity['location']['dy'] = dy
- self.velocity['location']['dz'] = dz
- self.velocity[self.convention]['dalpha_0'] = dyaw
- self.velocity[self.convention]['dalpha_1'] = dpitch
- self.velocity[self.convention]['dalpha_2'] = droll
-
- with self.assertRaises(ValueError):
- rof, hof, vof = mcode.optic_flow(self.scene,
- self.viewing_directions,
- self.velocity)
-
- # dpitch is too big
- positions = np.array([[-20, -20, 2.6, 1.57079633, 0, 0],
- [-19.8, -20, 2.6,
- 1.57079633, 0 + np.pi/2 + 0.0001, 0]])
-
- x = positions[0, 0]
- y = positions[0, 1]
- z = positions[0, 2]
- yaw = positions[0, 3]
- pitch = positions[0, 4]
- roll = positions[0, 5]
-
- dx = positions[1, 0] - positions[0, 0]
- dy = positions[1, 1] - positions[0, 1]
- dz = positions[1, 2] - positions[0, 2]
- dyaw = positions[1, 3] - positions[0, 3]
- dpitch = positions[1, 4] - positions[0, 4]
- droll = positions[1, 5] - positions[0, 5]
-
- self.velocity['location']['x'] = x
- self.velocity['location']['y'] = y
- self.velocity['location']['z'] = z
- self.velocity[self.convention]['alpha_0'] = yaw
- self.velocity[self.convention]['alpha_1'] = pitch
- self.velocity[self.convention]['alpha_2'] = roll
- self.velocity['location']['dx'] = dx
- self.velocity['location']['dy'] = dy
- self.velocity['location']['dz'] = dz
- self.velocity[self.convention]['dalpha_0'] = dyaw
- self.velocity[self.convention]['dalpha_1'] = dpitch
- self.velocity[self.convention]['dalpha_2'] = droll
-
- with self.assertRaises(ValueError):
- rof, hof, vof = mcode.optic_flow(self.scene,
- self.viewing_directions,
- self.velocity)
-
- # droll is too big
- positions = np.array([[-20, -20, 2.6, 1.57079633, 0, 0],
- [-19.8, -20, 2.6, 1.57079633, 0,
- 0 + np.pi/2 + 0.0001]])
-
- x = positions[0, 0]
- y = positions[0, 1]
- z = positions[0, 2]
- yaw = positions[0, 3]
- pitch = positions[0, 4]
- roll = positions[0, 5]
-
- dx = positions[1, 0] - positions[0, 0]
- dy = positions[1, 1] - positions[0, 1]
- dz = positions[1, 2] - positions[0, 2]
- dyaw = positions[1, 3] - positions[0, 3]
- dpitch = positions[1, 4] - positions[0, 4]
- droll = positions[1, 5] - positions[0, 5]
-
- self.velocity['location']['x'] = x
- self.velocity['location']['y'] = y
- self.velocity['location']['z'] = z
- self.velocity[self.convention]['alpha_0'] = yaw
- self.velocity[self.convention]['alpha_1'] = pitch
- self.velocity[self.convention]['alpha_2'] = roll
- self.velocity['location']['dx'] = dx
- self.velocity['location']['dy'] = dy
- self.velocity['location']['dz'] = dz
- self.velocity[self.convention]['dalpha_0'] = dyaw
- self.velocity[self.convention]['dalpha_1'] = dpitch
- self.velocity[self.convention]['dalpha_2'] = droll
-
- with self.assertRaises(ValueError):
- rof, hof, vof = mcode.optic_flow(self.scene,
- self.viewing_directions,
- self.velocity)
-
- def test_only_x(self):
- """
- this test checks for the correct response if for example
- the bee moves only in x-direction keeping the other
- parameters (y,z,yaw,pitch,roll) constant
- """
- convention = 'zyx'
- tuples_posvel = posorient_columns(convention)
- tuples_posvel.extend(velocities_columns(convention))
- index_posvel = pd.MultiIndex.from_tuples(tuples_posvel,
- names=['position',
- 'orientation'])
- velocity = pd.Series(index=index_posvel, data=0)
- velocity.loc[('location', 'dx')] = np.random.randn()
-
- elevation = np.linspace(-np.pi/2, np.pi/2, 5)
- azimuth = np.linspace(-np.pi, np.pi, 11)
- [ma, me] = np.meshgrid(azimuth, elevation)
- imshape = me.shape
- vdir = np.zeros((ma.shape[0], ma.shape[1], 2))
- vdir[..., __spherical_indeces__['elevation']] = me
- vdir[..., __spherical_indeces__['azimuth']] = ma
-
- scene = np.random.rand(imshape[0], imshape[1])
- scene = np.tile(scene[..., np.newaxis], 4)
- scene = scene[..., np.newaxis]
-
- rof, hof, vof =\
- mcode.optic_flow_translational(scene, vdir,
- velocity, distance_channel=3)
- hnorm = np.sqrt(hof**2 + vof ** 2)
- # Add abs tol because we compare to zero
- np.testing.assert_allclose(rof, np.zeros_like(rof), atol=1e-7)
- # The motion is in the plane x,y.
- # so no vof at null elevation
- vof_el = vof[elevation == 0, :]
- np.testing.assert_allclose(vof_el, np.zeros_like(vof_el), atol=1e-7)
-
- # The translational optic flow norm should be small
- # for further away objects
- # except for the focus of contraction
- # i.e. aximuth np.pi and el = 0
- # and focus of expension
- # i.e. aximuth 0 and el = 0
- valid = (vdir[..., __spherical_indeces__['elevation']] == 0) \
- & (vdir[..., __spherical_indeces__['azimuth']] == 0)
- valid = valid | (
- (vdir[..., __spherical_indeces__['elevation']] == 0)
- & (vdir[..., __spherical_indeces__['azimuth']] == np.pi))
- valid = valid | (
- (vdir[..., __spherical_indeces__['elevation']] == 0)
- & (vdir[..., __spherical_indeces__['azimuth']] == -np.pi))
- valid = not valid
-
- rof_further, hof_further, vof_further =\
- mcode.optic_flow_translational(scene+1, vdir,
- velocity, distance_channel=3)
- hnorm_further = np.sqrt(hof_further**2 + vof_further**2)
-
- np.testing.assert_array_less(hnorm_further[valid], hnorm[valid])
-
- def test_only_yaw(self):
- """
- this test checks for the correct response if for example
- the bee rotates only around the yaw axis keeping the other
- parameters (x,y,z,pitch,roll) constant
- """
- # yaw only everything zero
- # only vertical should change, horizontal stays same
- positions = np.array([[-20, -20, 2.6, 1.57079633, 0, 0],
- [-20, -20, 2.6, 1.90079633, 0, 0]])
-
- x = positions[0, 0]
- y = positions[0, 1]
- z = positions[0, 2]
- yaw = positions[0, 3]
- pitch = positions[0, 4]
- roll = positions[0, 5]
-
- dx = positions[1, 0]-positions[0, 0]
- dy = positions[1, 1]-positions[0, 1]
- dz = positions[1, 2]-positions[0, 2]
- dyaw = positions[1, 3]-positions[0, 3]
- dpitch = positions[1, 4]-positions[0, 4]
- droll = positions[1, 5]-positions[0, 5]
-
- self.velocity['location']['x'] = x
- self.velocity['location']['y'] = y
- self.velocity['location']['z'] = z
- self.velocity[self.convention]['alpha_0'] = yaw
- self.velocity[self.convention]['alpha_1'] = pitch
- self.velocity[self.convention]['alpha_2'] = roll
- self.velocity['location']['dx'] = dx
- self.velocity['location']['dy'] = dy
- self.velocity['location']['dz'] = dz
- self.velocity[self.convention]['dalpha_0'] = dyaw
- self.velocity[self.convention]['dalpha_1'] = dpitch
- self.velocity[self.convention]['dalpha_2'] = droll
-
- rof, hof, vof = mcode.optic_flow(self.scene, self.viewing_directions,
- self.velocity)
-
- testrof = [[-6.47644748e-26, -3.52655120e-19, -7.05202754e-19],
- [-1.84591335e-11, -1.00513560e-04, -2.00996485e-04],
- [-3.69126442e-11, -2.00996503e-04, -4.01931744e-04]]
- testhof = [[-0.33, -0.32994974, -0.32979897],
- [-0.33, -0.32994974, -0.32979897],
- [-0.33, -0.32994974, -0.32979897]]
- testvof = [[-1.05768414e-09, -5.75929518e-03, -1.15168350e-02],
- [-1.05752305e-09, -5.75841801e-03, -1.15150809e-02],
- [-1.05703983e-09, -5.75578677e-03, -1.15098192e-02]]
-
- np.testing.assert_allclose(rof, testrof)
- np.testing.assert_allclose(hof, testhof)
- np.testing.assert_allclose(vof, testvof)
-
- def test_only_yaw_big(self):
- """
- this test checks for the correct response if for example
- the bee rotates only around the yaw axis keeping the other
- parameters (x,y,z,pitch,roll) constant.
- here the azimuthal optic flow should be zero and the
- elevational optic flow should be proportional with a
- constant factor to the cos of the elevation
- """
- # generate scene aka distance channel
- scene = np.random.random((180, 360, 4, 1))
-
- elevation = np.arange(-np.pi/2, np.pi/2, 2*(np.pi/360))
- azimuth = np.arange(0, 2*np.pi, 2*(np.pi/360))
-
- viewing_directions = np.zeros((180, 2))
- viewing_directions[:, 0] = elevation
- viewing_directions[:, 1] = azimuth[0:180]
-
- positions = np.array([[-20, -20, 2.6, 1.57079633, 0, 0],
- [-20, -20, 2.6, 1.90079633, 0, 0]])
-
- x = positions[0, 0]
- y = positions[0, 1]
- z = positions[0, 2]
- yaw = positions[0, 3]
- pitch = positions[0, 4]
- roll = positions[0, 5]
-
- dx = positions[1, 0]-positions[0, 0]
- dy = positions[1, 1]-positions[0, 1]
- dz = positions[1, 2]-positions[0, 2]
- dyaw = positions[1, 3]-positions[0, 3]
- dpitch = positions[1, 4]-positions[0, 4]
- droll = positions[1, 5]-positions[0, 5]
-
- tuples = [('location', 'x'), ('location', 'y'),
- ('location', 'z'), ('location', 'dx'),
- ('location', 'dy'), ('location', 'dz'),
- (self.convention, 'alpha_0'), (self.convention, 'alpha_1'),
- (self.convention, 'alpha_2'), (self.convention, 'dalpha_0'),
- (self.convention, 'dalpha_1'), (self.convention, 'dalpha_2')]
- index = pd.MultiIndex.from_tuples(tuples,
- names=['position', 'orientation'])
- velocity = pd.Series(index=index)
-
- velocity['location']['x'] = x
- velocity['location']['y'] = y
- velocity['location']['z'] = z
- velocity['xyz']['alpha_0'] = yaw
- velocity['xyz']['alpha_1'] = pitch
- velocity['xyz']['alpha_2'] = roll
- velocity['location']['dx'] = dx
- velocity['location']['dy'] = dy
- velocity['location']['dz'] = dz
- velocity['xyz']['dalpha_0'] = dyaw
- velocity['xyz']['dalpha_1'] = dpitch
- velocity['xyz']['dalpha_2'] = droll
-
- rof, hof, vof = mcode.optic_flow(scene, viewing_directions,
- velocity)
- cosel = np.cos(elevation)
- has_zeros = len(np.where(cosel == 0)) > 0
- factor = np.array([0])
- if not has_zeros:
- factor = np.array(hof[:, 0]/cosel)
-
- self.assertAlmostEquals(factor, factor[0])
-
- def test_only_pitch(self):
- """
- this test checks for the correct response if for example
- the bee rotates only around the pitch axis keeping the other
- parameters (x,y,z,yaw,roll) constant
- """
- positions = np.array([[-20, -20, 2.6, 1.57079633, 0, 0],
- [-20, -20, 2.6, 1.57079633, 0.1, 0]])
-
- x = positions[0, 0]
- y = positions[0, 1]
- z = positions[0, 2]
- yaw = positions[0, 3]
- pitch = positions[0, 4]
- roll = positions[0, 5]
-
- dx = positions[1, 0] - positions[0, 0]
- dy = positions[1, 1] - positions[0, 1]
- dz = positions[1, 2] - positions[0, 2]
- dyaw = positions[1, 3] - positions[0, 3]
- dpitch = positions[1, 4] - positions[0, 4]
- droll = positions[1, 5] - positions[0, 5]
-
- self.velocity['location']['x'] = x
- self.velocity['location']['y'] = y
- self.velocity['location']['z'] = z
- self.velocity[self.convention]['alpha_0'] = yaw
- self.velocity[self.convention]['alpha_1'] = pitch
- self.velocity[self.convention]['alpha_2'] = roll
- self.velocity['location']['dx'] = dx
- self.velocity['location']['dy'] = dy
- self.velocity['location']['dz'] = dz
- self.velocity[self.convention]['dalpha_0'] = dyaw
- self.velocity[self.convention]['dalpha_1'] = dpitch
- self.velocity[self.convention]['dalpha_2'] = droll
- rof, hof, vof = mcode.optic_flow(self.scene, self.viewing_directions,
- self.velocity)
-
- testrof = [[0.1, 0.1, 0.1],
- [0.09998477, 0.09998477, 0.09998477],
- [0.09993908, 0.09993908, 0.09993908]]
- testhof = [[1.02726882e-18, 5.59367784e-12, 1.11856508e-11],
- [1.02726882e-18, 5.59367784e-12, 1.11856508e-11],
- [1.02726882e-18, 5.59367784e-12, 1.11856508e-11]]
- testvof = [[-3.20510352e-10, -3.20461536e-10, -3.20315105e-10],
- [-1.74524096e-03, -1.74524096e-03, -1.74524096e-03],
- [-3.48994999e-03, -3.48994999e-03, -3.48994999e-03]]
-
- np.testing.assert_allclose(rof, testrof)
- np.testing.assert_allclose(hof, testhof)
- np.testing.assert_allclose(vof, testvof)
-
- def test_only_pitch_big(self):
- """
- this test checks for the correct response if for example
- the bee rotates only around the pitch axis keeping the other
- parameters (x,y,z,yaw,roll) constant.
- here the azimuthal optic flow should be zero and the
- elevational optic flow should be proportional with a
- constant factor to the cos of the elevation
- """
- scene = np.random.random((180, 360, 4, 1))
-
- elevation = np.arange(-np.pi/2, np.pi/2, 2*(np.pi/360))
-
- azimuth = np.arange(0, 2*np.pi, 2*(np.pi/360))
- viewing_directions = np.zeros((180, 2))
- viewing_directions[:, 0] = elevation
- viewing_directions[:, 1] = azimuth[0:180]
-
- positions = np.array([[-20, -20, 2.6, 1.57079633, 0, 0],
- [-20, -20, 2.6, 1.57079633, 0.1, 0]])
-
- x = positions[0, 0]
- y = positions[0, 1]
- z = positions[0, 2]
- yaw = positions[0, 3]
- pitch = positions[0, 4]
- roll = positions[0, 5]
-
- dx = positions[1, 0] - positions[0, 0]
- dy = positions[1, 1] - positions[0, 1]
- dz = positions[1, 2] - positions[0, 2]
- dyaw = positions[1, 3] - positions[0, 3]
- dpitch = positions[1, 4] - positions[0, 4]
- droll = positions[1, 5] - positions[0, 5]
-
- self.velocity['location']['x'] = x
- self.velocity['location']['y'] = y
- self.velocity['location']['z'] = z
- self.velocity[self.convention]['alpha_0'] = yaw
- self.velocity[self.convention]['alpha_1'] = pitch
- self.velocity[self.convention]['alpha_2'] = roll
- self.velocity['location']['dx'] = dx
- self.velocity['location']['dy'] = dy
- self.velocity['location']['dz'] = dz
- self.velocity[self.convention]['dalpha_0'] = dyaw
- self.velocity[self.convention]['dalpha_1'] = dpitch
- self.velocity[self.convention]['dalpha_2'] = droll
-
- rof, hof, vof = mcode.optic_flow(scene, viewing_directions,
- self.velocity)
-
- cosel = np.cos(elevation)
- has_zeros = len(np.where(cosel == 0)) > 0
- factor = np.array([0])
- if not has_zeros:
- factor = np.array(vof[0, :]/cosel)
-
- assert np.all(factor == factor[0])
-
- def test_only_roll_big(self):
- """
- this test checks for the correct response if for example
- the bee rotates only around the roll axis keeping the other
- parameters (x,y,z,yaw,pitch) constant
- here the azimuthal optic flow should be zero and the
- elevational optic flow should be proportional with a
- constant factor to the cos of the elevation
- """
- scene = np.random.random((180, 360, 4, 1))
-
- elevation = np.arange(-np.pi/2, np.pi/2, 2*(np.pi/360))
-
- azimuth = np.arange(0, 2*np.pi, 2*(np.pi/360))
- viewing_directions = np.zeros((180, 2))
- viewing_directions[:, 0] = elevation
- viewing_directions[:, 1] = azimuth[0:180]
-
- positions = np.array([[-20, -20, 2.6, 1.57079633, 0, 0],
- [-20, -20, 2.6, 1.57079633, 0, 0.1]])
-
- x = positions[0, 0]
- y = positions[0, 1]
- z = positions[0, 2]
- yaw = positions[0, 3]
- pitch = positions[0, 4]
- roll = positions[0, 5]
-
- dx = positions[1, 0] - positions[0, 0]
- dy = positions[1, 1] - positions[0, 1]
- dz = positions[1, 2] - positions[0, 2]
- dyaw = positions[1, 3] - positions[0, 3]
- dpitch = positions[1, 4] - positions[0, 4]
- droll = positions[1, 5] - positions[0, 5]
-
- self.velocity['location']['x'] = x
- self.velocity['location']['y'] = y
- self.velocity['location']['z'] = z
- self.velocity[self.convention]['alpha_0'] = yaw
- self.velocity[self.convention]['alpha_1'] = pitch
- self.velocity[self.convention]['alpha_2'] = roll
- self.velocity['location']['dx'] = dx
- self.velocity['location']['dy'] = dy
- self.velocity['location']['dz'] = dz
- self.velocity[self.convention]['dalpha_0'] = dyaw
- self.velocity[self.convention]['dalpha_1'] = dpitch
- self.velocity[self.convention]['dalpha_2'] = droll
-
- rof, hof, vof = mcode.optic_flow(scene, viewing_directions,
- self.velocity)
- cosel = np.cos(elevation)
- sinel = np.sin(elevation)
- has_zeros = cosel[np.where(cosel == 0)] = 1
- factor = np.array(hof[0, :]/cosel)
- factor[has_zeros] = factor[1]
- has_zerossin = sinel[np.where(sinel == 0)] = 1
- factorsin = np.array(vof[0, :]/sinel)
- factorsin[has_zerossin] = factorsin[0]
-
- for i in range(1, len(factor)):
- self.assertAlmostEqual(np.isclose(factor[1], factor[i]))
- if i == 90:
- continue
- self.assertAlmostEqual(factorsin[0], factorsin[i])
- # print(i)
- # assert np.all(np.isclose(hof, testhof))
- # assert np.all(np.isclose(vof, testvof))
-
- def test_only_roll(self):
- """
- this test checks for the correct response if for example
- the bee rotates only around the pitch axis keeping the other
- parameters (x,y,z,yaw,roll) constant
- """
- positions = np.array([[-20, -20, 2.6, 1.57079633, 0.1, 0.1],
- [-20, -20, 2.6, 1.57079633, 0.1, 0.2]])
-
- x = positions[0, 0]
- y = positions[0, 1]
- z = positions[0, 2]
- yaw = positions[0, 3]
- pitch = positions[0, 4]
- roll = positions[0, 5]
-
- dx = positions[1, 0] - positions[0, 0]
- dy = positions[1, 1] - positions[0, 1]
- dz = positions[1, 2] - positions[0, 2]
- dyaw = positions[1, 3] - positions[0, 3]
- dpitch = positions[1, 4] - positions[0, 4]
- droll = positions[1, 5] - positions[0, 5]
-
- self.velocity['location']['x'] = x
- self.velocity['location']['y'] = y
- self.velocity['location']['z'] = z
- self.velocity[self.convention]['alpha_0'] = yaw
- self.velocity[self.convention]['alpha_1'] = pitch
- self.velocity[self.convention]['alpha_2'] = roll
- self.velocity['location']['dx'] = dx
- self.velocity['location']['dy'] = dy
- self.velocity['location']['dz'] = dz
- self.velocity[self.convention]['dalpha_0'] = dyaw
- self.velocity[self.convention]['dalpha_1'] = dpitch
- self.velocity[self.convention]['dalpha_2'] = droll
- rof, hof, vof = mcode.optic_flow(self.scene, self.viewing_directions,
- self.velocity)
- testrof = [[0.01088051, 0.01088051, 0.01088051],
- [0.0126067, 0.01260916, 0.01261109],
- [0.01432905, 0.01433397, 0.01433784]]
- testhof = [[0.00894177, 0.00721257, 0.00548116],
- [0.00894177, 0.00721257, 0.00548116],
- [0.00894177, 0.00721257, 0.00548116]]
- testvof = [[0.09900333, 0.09914431, 0.09925508],
- [0.09879836, 0.09893931, 0.09905007],
- [0.09856329, 0.09870419, 0.09881489]]
-
- np.testing.assert_allclose(rof, testrof)
- np.testing.assert_allclose(hof, testhof)
- np.testing.assert_allclose(vof, testvof)
-
-
-if __name__ == '__main__':
- unittest.main()
diff --git a/navipy/processing/test_opticflow.py b/navipy/processing/test_opticflow.py
new file mode 100644
index 0000000000000000000000000000000000000000..4d7acccd2e701fdd21f4ff49c1cf5e07687c4ec3
--- /dev/null
+++ b/navipy/processing/test_opticflow.py
@@ -0,0 +1,127 @@
+import unittest
+import navipy.processing.mcode as mcode
+import pandas as pd
+import numpy as np
+from navipy.moving.agent import posorient_columns
+from navipy.moving.agent import velocities_columns
+from navipy.scene import __spherical_indeces__
+
+
+class OpticFlowTest(unittest.TestCase):
+ """
+ Test the geometrical computation of the optic flow
+ """
+
+ def setUp(self):
+ """ Init vectors
+ """
+ convention = 'zyx'
+ tuples_posvel = posorient_columns(convention)
+ tuples_posvel.extend(velocities_columns(convention))
+ index_posvel = pd.MultiIndex.from_tuples(tuples_posvel,
+ names=['position',
+ 'orientation'])
+ velocity = pd.Series(index=index_posvel, data=0)
+ # Most test are independent of orientation and position
+ velocity.loc[(convention, 'alpha_0')] = 2*np.pi*(np.random.rand()-0.5)
+ velocity.loc[(convention, 'alpha_1')] = np.pi*(np.random.rand()-0.5)
+ velocity.loc[(convention, 'alpha_2')] = 2*np.pi*(np.random.rand()-0.5)
+ velocity.loc[('location', 'x')] = np.random.randn()
+ velocity.loc[('location', 'y')] = np.random.randn()
+ velocity.loc[('location', 'z')] = np.random.randn()
+ self.velocity = velocity
+ self.convention = convention
+
+ # Init the viewing directions
+ elevation = np.linspace(-np.pi/2, np.pi/2, 5)
+ azimuth = np.linspace(-np.pi, np.pi, 11)
+ [ma, me] = np.meshgrid(azimuth, elevation)
+ imshape = me.shape
+ viewing_directions = np.zeros((ma.shape[0], ma.shape[1], 2))
+ viewing_directions[..., __spherical_indeces__['elevation']] = me
+ viewing_directions[..., __spherical_indeces__['azimuth']] = ma
+ self.viewing_directions = viewing_directions
+ self.elevation = elevation
+ self.azimuth = azimuth
+
+ # Init a random scene
+ scene = np.random.rand(imshape[0], imshape[1])
+ scene = np.tile(scene[..., np.newaxis], 4)
+ scene = scene[..., np.newaxis]
+ self.scene = scene
+
+ def test_distance(self):
+ """
+ The magnitude of the optic flow when the agent is moving \
+close to object is larger than when it is moving (with the same velocity)\
+far from them.
+ Here we test this property.
+ """
+ vel = self.velocity.copy()
+ # This is true for any translational motion
+ vel.loc[('location', 'dx')] = np.random.randn()
+ vel.loc[('location', 'dy')] = np.random.randn()
+ vel.loc[('location', 'dz')] = np.random.randn()
+ #
+ rof, hof, vof =\
+ mcode.optic_flow(self.scene,
+ self.viewing_directions,
+ vel)
+ hnorm = np.sqrt(hof**2 + vof ** 2)
+ # Add abs tol because we compare to zero
+ np.testing.assert_allclose(rof, np.zeros_like(rof), atol=1e-7)
+ # Calculate second optic flow, with objects further away form agent
+ rof_further, hof_further, vof_further =\
+ mcode.optic_flow(self.scene+1,
+ self.viewing_directions,
+ vel)
+ hnorm_further = np.sqrt(hof_further**2 + vof_further**2)
+ # Add abs tol because we compare to zero
+ np.testing.assert_allclose(rof_further,
+ np.zeros_like(rof_further), atol=1e-7)
+ # The translational optic flow norm should be small
+ # i.e. for norm equal to zero
+ valid = (hnorm != 0) & (hnorm_further != 0)
+ np.testing.assert_array_less(hnorm_further[valid], hnorm[valid])
+
+ def test_xyplane_only(self):
+ """
+ When the agent is moving along in the plane x,y, the gOF
+along the vertical gOF is ull at null elevation
+ """
+ vel = self.velocity.copy()
+ # This is true for any x-y translational motion
+ vel.loc[('location', 'dx')] = np.random.randn()
+ vel.loc[('location', 'dy')] = np.random.randn()
+ vel.loc[('location', 'dz')] = 0
+ # but in the coordinate of the bee
+ # by setting euler angles to zero, the axis of the
+ # bee coordinate system and the world one match
+ vel.loc[(self.convention, 'alpha_0')] = 0
+ vel.loc[(self.convention, 'alpha_1')] = 0
+ vel.loc[(self.convention, 'alpha_2')] = 0
+
+ # Calculate gOF
+ rof, hof, vof =\
+ mcode.optic_flow(self.scene,
+ self.viewing_directions,
+ vel)
+ # Add abs tol because we compare to zero
+ np.testing.assert_allclose(rof, np.zeros_like(rof), atol=1e-7)
+ # At null elevation only this is true
+ vof_el = vof[self.elevation == 0, :]
+ np.testing.assert_allclose(vof_el, np.zeros_like(vof_el), atol=1e-7)
+
+ def test_yaw_rot(self):
+ """
+ The optic for a rotation around the z-axis in the \
+yaw-pitch-roll (zyx) convention has vertical gOF equal to zero
+ """
+ vel = self.velocity.copy()
+ vel.loc[(self.convention, 'dalpha_0')] = np.random.rand()
+ rof, hof, vof = mcode.optic_flow(self.scene,
+ self.viewing_directions,
+ vel)
+ # Add abs tol because we compare to zero
+ np.testing.assert_allclose(rof, np.zeros_like(rof), atol=1e-7)
+ np.testing.assert_allclose(vof, np.zeros_like(vof), atol=1e-7)
diff --git a/test_of.ipynb b/test_of.ipynb
index 23027b3a14505fc638f0f6d3fb4f6ac2bee1135b..905cbb663befbc6c8953ff90df1a64bc96f09700 100644
--- a/test_of.ipynb
+++ b/test_of.ipynb
@@ -39,9 +39,43 @@
"from navipy.scene import __spherical_indeces__\n"
]
},
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Optic Flow\n",
+ "\n",
+ "Optic flow is defined as the change of light in the image, e.g. on the retina or the camera’s sensor, due to a relative motion between the eye or camera and the scene.\n",
+ "\n",
+ "In this section we focus on calculating the optic flow from the animal or camera motion and the distance to surrounding objects. Note, that in applied task, the problem of interest is the opposite of our current focus, namely the problem is how to find the animal or camera motion and the distance to surrounding objects from an estimate of the optic flow obtained from a series of images.\n",
+ "\n",
+ "To avoid confusion between, we qualify the optic flow from known animal or camera motion and distance to surrounding objects as geometrical, and will call it geometrical optic flow or gOF. \n",
+ "\n",
+ "Koendering Van Dorn in their 1987 article 'Facts on optic flow' derive geometrical optic flow, and obtained the following equation\n",
+ "\n",
+ "$$ gOF_i = -\\mu_i(t-t. d_i)d_i)-R\\times d_i $$\n",
+ "\n",
+ "here \n",
+ "* $gOF_i$ is the geometrical optic flow in the $i-th$ viewing direction $d_i$. \n",
+ "* $t$ and $R$ are the translational and angular velocity, respectively. \n",
+ "* $\\mu_i$ is the nearness (inverse of the distance) to the closest object in the $i-th$ viewing direction\n",
+ "\n",
+ "The equation can be rewritten by using the \"apparent rotation\" due to the translation $A_i=\\mu_i di_i\\times t$, and become\n",
+ "\n",
+ "$$ gOF_i = -(A_i+R)\\times d_i $$\n",
+ "\n",
+ "The optic flow is thus the sum of apparent rotation due to the translation and the angular rotation projected to a plane orthogonal to the viewing direction. \n",
+ "\n",
+ "The eye is sometimes described as a spherical apparatus (espectially in insect research), here each viewing direction can be expressed in a spherical coordinate system. The gOF in a spherical system as three component, but the one along the viewing direction (i.e. the $\\rho$ of the coordinate system) equates zero, because the gOF is othogonal to that direction.\n",
+ "\n",
+ "## From \n",
+ "\n",
+ "The distance "
+ ]
+ },
{
"cell_type": "code",
- "execution_count": 171,
+ "execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
@@ -53,6 +87,8 @@
" 'orientation'])\n",
"velocity = pd.Series(index=index_posvel, data=0)\n",
"velocity.loc[(convention,'alpha_0')] = np.pi*(np.random.rand()-0.5)\n",
+ "velocity.loc[(convention,'alpha_1')] = np.pi*(np.random.rand()-0.5)\n",
+ "velocity.loc[(convention,'alpha_2')] = np.pi*(np.random.rand()-0.5)\n",
"velocity.loc[('location','dx')]=np.random.randn()\n",
"\n",
"\n",
@@ -76,8 +112,8 @@
"np.testing.assert_allclose(rof, np.zeros_like(rof), atol=1e-7) \n",
"# The motion is in the plane x,y.\n",
"# so no vof at null elevation \n",
- "vof_el = vof[elevation==0,:]\n",
- "np.testing.assert_allclose(vof_el, np.zeros_like(vof_el), atol=1e-7) \n",
+ "#vof_el = vof[elevation==0,:]\n",
+ "#np.testing.assert_allclose(vof_el, np.zeros_like(vof_el), atol=1e-7) \n",
"\n",
"# The translational optic flow norm should be small\n",
"# for further away objects\n",
@@ -103,7 +139,7 @@
},
{
"cell_type": "code",
- "execution_count": 173,
+ "execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
@@ -116,7 +152,7 @@
"velocity = pd.Series(index=index_posvel, data=0)\n",
"velocity.loc[(convention,'alpha_0')] = np.pi*(np.random.rand()-0.5)\n",
"velocity.loc[(convention,'alpha_1')] = np.pi*(np.random.rand()-0.5)\n",
- "#velocity.loc[(convention,'alpha_2')] = np.pi*(np.random.rand()-0.5)\n",
+ "velocity.loc[(convention,'alpha_2')] = np.pi*(np.random.rand()-0.5)\n",
"velocity.loc[('location','dy')]=np.random.randn()\n",
"\n",
"\n",
@@ -140,8 +176,8 @@
"np.testing.assert_allclose(rof, np.zeros_like(rof), atol=1e-7) \n",
"# The motion is in the plane x,y.\n",
"# so no vof at null elevation \n",
- "vof_el = vof[elevation==0,:]\n",
- "np.testing.assert_allclose(vof_el, np.zeros_like(vof_el), atol=1e-7) \n",
+ "#vof_el = vof[elevation==0,:]\n",
+ "#np.testing.assert_allclose(vof_el, np.zeros_like(vof_el), atol=1e-7) \n",
"\n",
"# The translational optic flow norm should be small\n",
"# for further away objects\n",
@@ -165,7 +201,7 @@
},
{
"cell_type": "code",
- "execution_count": 174,
+ "execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
@@ -176,6 +212,9 @@
" names=['position',\n",
" 'orientation'])\n",
"velocity = pd.Series(index=index_posvel, data=0)\n",
+ "velocity.loc[(convention,'alpha_0')] = np.pi*(np.random.rand()-0.5)\n",
+ "velocity.loc[(convention,'alpha_1')] = np.pi*(np.random.rand()-0.5)\n",
+ "velocity.loc[(convention,'alpha_2')] = np.pi*(np.random.rand()-0.5)\n",
"velocity.loc[('location','dz')]=np.random.randn()\n",
"\n",
"\n",
@@ -200,8 +239,8 @@
"# The motion is along the x axis, in the plane\n",
"# y,z\n",
"# so no hof at null azimuth \n",
- "hof_el = hof[:,azimuth==0]\n",
- "np.testing.assert_allclose(hof_el, np.zeros_like(hof_el), atol=1e-7) \n",
+ "#hof_el = hof[:,azimuth==0]\n",
+ "#np.testing.assert_allclose(hof_el, np.zeros_like(hof_el), atol=1e-7) \n",
"\n",
"# The translational optic flow norm should be small\n",
"# for further away objects\n",
@@ -230,7 +269,7 @@
},
{
"cell_type": "code",
- "execution_count": 183,
+ "execution_count": null,
"metadata": {},
"outputs": [],
"source": [
@@ -241,9 +280,9 @@
" names=['position',\n",
" 'orientation'])\n",
"velocity = pd.Series(index=index_posvel, data=0)\n",
- "#velocity.loc[(convention,'alpha_0')] = 2*np.pi*(np.random.rand()-0.5)\n",
- "#velocity.loc[(convention,'alpha_1')] = np.pi*(np.random.rand()-0.5)\n",
- "#velocity.loc[(convention,'alpha_2')] = 2*np.pi*(np.random.rand()-0.5)\n",
+ "velocity.loc[(convention,'alpha_0')] = 2*np.pi*(np.random.rand()-0.5)\n",
+ "velocity.loc[(convention,'alpha_1')] = np.pi*(np.random.rand()-0.5)\n",
+ "velocity.loc[(convention,'alpha_2')] = 2*np.pi*(np.random.rand()-0.5)\n",
"velocity.loc[(convention,'dalpha_0')]= np.random.rand()\n",
"\n",
"elevation = np.linspace(-np.pi/2,np.pi/2,5)\n",
@@ -268,32 +307,9 @@
},
{
"cell_type": "code",
- "execution_count": 184,
+ "execution_count": null,
"metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- "<matplotlib.quiver.Quiver at 0x7f00b03e9fd0>"
- ]
- },
- "execution_count": 184,
- "metadata": {},
- "output_type": "execute_result"
- },
- {
- "data": {
- "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADx0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wcmMyLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvMCCy2AAAIABJREFUeJzt3Xl0U3X+PvDn3YXCl01Kyy5UEGQRQRsBBR1RNkfcEBDcxnGk6k/cviMuw3FhdAYRjzqjuOCKDAcXHJS0ZUdKCxZoK9CNQtnbYheglNLSNsn79wch38a20Ga5Ae/zOieH5ObePJ8LNE/u/SSNqCqIiMh8ggI9ACIiCgwWABGRSbEAiIhMigVARGRSLAAiIpNiARARmRQLgIjIpFgAREQmxQIgIjKpkEAP4GwiIiI0Kioq0MMgIrpgpKamlqhqZGPWPa8LICoqCikpKYEeBhHRBUNEDjR2XZ4CIiIyKRYAEZFJsQCIiEyKBUBEZFIsACIik2IBEBGZFAuAiMikfncFUFhYiCuuuALPPPOM4dkzZszAwIEDUVBQYGhueno6+vXrh7feesvQXAC45557cM0116CiosLQ3BUrVqBv3774+uuvDc212Wy46aabMH78eBj9daqfffYZ+vbti8TERENzS0tLER0djWnTphmaCwCzZs1C//79sWfPHkNz9+7di/79+2PWrFmG5gJATEwMoqOjUVpa6ves310BbNu2Denp6bBarYZnW61WZGRkIDU11dDcjRs3YufOnVi+fLmhuRUVFVixYgU2b96M3NxcQ7N/+ukn5OTkYM2aNYbmFhQUIDExEevWrUNJSYmh2atWrUJOTg4SEhIMzc3OzkZaWhri4uJgs9kMzY6Li0N2djaSk5MNzU1OTkZ2djZiY2MNzbXZbIiNjUVaWhoyMzP9H6iqXl8AfA6gCEBGA/ffAOA4gG3Oy8uNedzo6Gj1xA8//KAZGRkebeuNzMxMXbp0qeG51dXV+t133+nevXsNz966dauuXLnS8NzS0lL99ttvtbCw0PDs9evXa1JSkuG5+fn5+u2332p5ebnh2XFxcfrLL78Ynrtr1y5dsmSJ2u12Q3PtdrsuWbJEd+3aZWiuquovv/yicXFxHm8PIEUb+dwt6oPDWBG5HkA5gK9U9fJ67r8BwLOqOr4pj2uxWJS/CoKIqPFEJFVVLY1Z1yengFR1A4CjvngsIiIyhpFzANeIyHYRWS4iAwzMJSKiehj120DTAPRQ1XIR+SOAHwD0rm9FEYkBEAMA3bt3N2h4RETmY8gRgKqWqWq583o8gFARiWhg3fmqalFVS2Rko36lNRERecCQAhCRTiIizutDnLlHjMgmIqL6+eQUkIgsxum3ekaISB6AVwCEAoCqfgRgIoDHRMQGoBLAFPXF24+IiMhjPikAVZ16jvvfB/C+L7KIiMg3fnefBCYiosZhARARmRQLgIjIpFgAREQmxQIgIjIpFgARkUmxAIiITIoFQERkUiwAIiKTYgEQEZkUC4CIyKRYAEREJsUCICIyKRYAEZFJsQCIiEyKBUBEZFIsACIik2IBEBGZFAuAiMikWABERCbFAiAiMikWABGRSfmkAETkcxEpEpGMBu4XEfm3iOSKyA4RucoXuURE5DlfHQF8CWDcWe6/GUBv5yUGwIc+yiUiIg/5pABUdQOAo2dZ5XYAX+lpyQAuEpHOvsgmIiLPGDUH0BXAoVq385zLiIgoQM67SWARiRGRFBFJKS4uDvRwiIh+t4wqgHwAF9e63c25rA5Vna+qFlW1REZGGjI4IiIzMqoAlgF4wPluoGEAjqvqYYOyiYioHiG+eBARWQzgBgARIpIH4BUAoQCgqh8BiAfwRwC5ACoA/NkXuURE5DmfFICqTj3H/QrgcV9kERGRb5x3k8BERGQMFgARkUmxAIiITIoFQERkUiyAetTU1JgqN5DZgcpVVdPts8PhgN1uD0h2oPbZZrPh9HtQjBfIn+fGYgHg9JNBWloaZs2aBYvFguTkZMOyS0pK8NVXX2Hy5Ml45JFHDMu12+1ISkrC888/j0GDBuHQoUPn3shH8vPz8fHHH2P8+PGYNWuWYblVVVVYtWoVnnjiCVgsFlRUVBiWnZubi3feeQc33ngjPv/8c8NyT548iR9++AEPP/wwbrzxRsOeDFUVO3bswD/+8Q8MGzYMq1evNiQXAI4dO4bFixfjnnvuwT333AMRMSTX4XAgOTkZM2fOxODBg7Fz505Dcr2iquftJTo6Wv2loqJCrVarxsTEaJcuXRSAAtCRI0f6LVNV1eFwaEZGhr7xxhs6fPhwDQoKcmVv27bNr9mlpaX6zTff6P3336/t27d35T700EN+zbXb7bp161Z9+eWX9corr3TltmjRQg8fPuzX7KKiIv3iiy/0rrvu0latWrmyX3vtNb/m1tTUaEJCgj777LN62WWXuXK7dOmilZWVfs0+ePCgfvDBB3rzzTdrWFiYK3vhwoV+zT116pQuX75cH3/8ce3evbsrNzo6Wh0Oh1+zc3Jy9K233tI//OEPGhwc7MpOTEz0a+6JEyf0+++/1z//+c/aoUMHV+7EiRP9mns2AFK0kc+xogE6PGqMyMhIveuuu5q8XY8ePfDiiy/WWV5QUIC4uDhYrVasWbMGlZWVdda5/vrr0a9fP4/GCwDPPvssLr30Urdl1dXV2LBhA6xWK6xWK/bt21dnu3bt2mHy5Mke50ZHR2PatGl1lu/ZswexsbGwWq1ISEiAzWars84tt9yCbt26eZz9xhtv4KKLLnJbVlFRgTVr1sBqtSIuLg6HD9f94Hf37t1x8803e5w7ZswYTJgwwW2ZqiIzM9P1d52cnFzvq94pU6agbdu2HuUGBwdj3rx5dZaXlpZixYoViI2NRXx8PI4dO1ZnnQEDBmDEiBEe5QLA1KlT8Yc//MFtmcPhQEpKimuft2/fXmc7EcFDDz2EkBDPPvrTsWPHeo/WioqKXD9Tq1atwsmTJ+usM2zYMAwaNMijXACYPn06Lr/8crdlNTU12Lhxo2ufd+/eXWe7li1b4r777vM4t3///njyySfrLD9w4IDrZ+qnn35CdXV1nXVGjx6Nnj17epz96quvolOnTh5tKyKpqmpp1MqNbYpAXOBs06ZeLBaLqw0dDof++OOParFYPHqspl6SkpJc2fn5+frggw9qmzZt/J47adIkV67dbtdPP/1U+/XrZ8g+5+fnu7JzcnJ0woQJ2rx5c7/nzpgxw5VbVVWlc+bM0aioKL/nhoaGur3i2rJli44ePVpDQkL8nv3BBx+4ck+cOKEvvPCCdurUye+5ffr0cdvnFStW6LBhw1RE/J4dHx/vyi0qKtJHHnlEL7roIr/njh071u15ZNGiRXrFFVf4PReAZmdnq6fQhCMAn3wS2F9CQkLqvLJsjNrbiAhuu+02DB482NXa69atq7e1AaB169YICwvzasxndOnSBR9++CEmT57sym7oXLun+3pG69atXdeDgoLwl7/8BSNGjHC9Qtq4cWODE4AXXXSRx68Mz+Sd0adPH3z55ZdYvXq169V/Q7/VNSwszG3cTdWyZUvX9WbNmmHGjBkYPXo0rFYrYmNjsXXr1ga3DQ8Pdxt3U4SGhrrdvvrqq7F48WIsX74cVqsVK1asQFlZWb3btmjRwm3cTdW8eXPX9VatWuH111/Hrbfe6vp3zszMbHDbiIgIj3PbtWvndnvs2LEYOHAgYmNjERsb2+AR9Zlx1h53U9X++46MjMS///1vTJw48axH1MDpI7Xfjrspah8higjuueceDB061PX/q6EjagBo06YNmjVr5nF2cHCwx9s2SWObIhAXf80BnDhxQpcuXaoPPfSQ23k7APr444/7JVP19KuIbdu26WuvvaZDhgxxy23durWWlJT4LfvIkSO6aNEinTJlirZt29Yte+7cuX7LtdlsumnTJv3b3/6mAwcOdMuNiorS6upqv2UXFBToJ598orfddpu2aNHCLXvJkiV+y62qqtI1a9bo008/rb169XLLHTFihF/Ph+/Zs0f/9a9/6ahRozQ0NNQtOyUlxW+5J0+erHdODYDef//9fss9M6c2e/bsOnNqYWFhmpeX57fsM3Nq9913n4aHh7vt80svveS33HNBE44AAv4kf7aLPyeBz7Db7ZqcnKwzZ87UK664QsPCwvTQoUN+z1VVPXz4sH766ad6xx136P/8z//ozJkzDcmtrq7WdevW6TPPPKOXXnqpRkZGanl5uSHZ+/bt0/fee0/HjBmjoaGh+umnnxqSW1FRobGxsfrII49o165ddeDAgWq32/2e63A4NCsrS+fMmaMjRozQoKAgXbNmjd9zVVWPHz+u3333nT7wwAPavn17vfXWWw3JdTgcmpKSoq+88opeddVVGhwcrDk5OYZkFxcX64IFC3TixInaunVrfeKJJwzJramp0Q0bNuiMGTO0X79+2qZNGz1y5Igh2b/VlAI4ryeBLRaLpqSkGJp54MABVFdXo3fv3obmnjp1CqmpqRg+fLihuaqKXbt2oWXLll5NBHvixIkTyM7OxpAhQwzNVVVs27YNUVFRXp0i8ERJSQny8vIwePBgQ3Ptdjs2b94Mi8Xi1akJT+Tn56OsrMyrN1d4orq6Gps3b8aIESMMeyvoGbm5uQgJCUFUVJShuUDTJoFZAEREvyNNKQB+EIyIyKRYAEREJsUCICIyKRYAEZFJsQCIiEyKBUBEZFIsACIik2IBEBGZFAuAiMikfFIAIjJORHJEJFdEXqjn/gdFpFhEtjkvD/sil4iIPOf1r4MWkWAA8wCMBpAHYKuILFPVrN+s+o2qTvc2j4iIfMMXRwBDAOSq6l5VrQbwNYDbffC4RETkR74ogK4Aan/LSZ5z2W/dJSI7RGSJiFzsg1wiIvKCUZPAVgBRqnoFgNUAFjS0oojEiEiKiKQ09E1SRETkPV8UQD6A2q/ouzmXuajqEVWtct78FEB0Qw+mqvNV1aKqlsjISB8Mj4iI6uOLAtgKoLeIXCIizQBMAbCs9goi0rnWzdsAZPsgl4iIvOD1u4BU1SYi0wGsBBAM4HNVzRSRv+P0V5MtA/CkiNwGwAbgKIAHvc0lIiLv8BvBiIh+R/iNYEREdE4sACIik2IBEBGZFAuAiMikWABERCbFAiAiMikWABGRSbEAiIhMigVARGRSLAAiIpNiARARmRQLgIjIpFgAREQmxQIgIjIpFgARkUmxAIiITIoFQERkUiwAIiKTYgEQEZkUC4CIyKRYAEREJsUCICIyKZ8UgIiME5EcEckVkRfquT9MRL5x3r9ZRKJ8kXsuqoqsrCzMnTsXx48fNyLSpaysDN999x3i4uIMzVVVpKamYu7cuaipqTE0u7i4GAsWLMCmTZsMzbXZbEhMTMR7770HVTU0Oy8vDx999BGys7MNza2qqsLKlSvx5ZdfGpoLALt378bbb7+NoqIiQ3PLy8uxdOlSLF261NBcVcX27dvx5ptv4tSpU4Zm+52qenUBEAxgD4CeAJoB2A6g/2/W+X8APnJenwLgm8Y8dnR0tDZVVVWVrlmzRp966int2bOnAtCJEyc2+XE8sWfPHn333Xd11KhRGhoaqiEhIbp3716/5548eVKXLVum06ZN0y5duigAnTFjht9zHQ6Hpqen6+zZs/Xaa69VEdHw8HA9fvy437NLS0v166+/1nvvvVfDw8MVgM6bN8/vuXa7Xbds2aIvvfSSDh48WAFo37591Waz+T27sLBQv/jiC50wYYK2atVKAWh8fLzfc2tqanT9+vX617/+Vfv06aMAdMyYMX7PVVU9cOCAzps3T8eNG6dhYWEqIpqenu733MrKSo2Pj9fHHntML774YgWgjzzyiN9zfQFAijby+TvEBx0yBECuqu4FABH5GsDtALJqrXM7gFed15cAeF9ExDnYBlVWViI9Pf2cAygtLUViYiISEhKwadMmlJeXu91/5513Nupxmsput2PHjh1ISEhAQkIC9uzZ43b/uHHjUF5e7pfswsJCbNiwAQkJCdi8eTOqqqpc9wUFBWHs2LF+ya2urkZqaioSEhKwfv16FBQUuN1/55134sCBAz7PBYCDBw+6/q7T0tJgs9lc97Vp0wZDhgzxyz5XVFRg8+bNWL9+PRITE1FSUuJ2/913342srKwGtvacqmL37t2ufU5PT3c7wunVqxe6du3ql30uKytDUlISEhISkJSUhBMnTrjdf9ddd/kl1+FwICMjw7XPu3btcrv/uuuug6r6JbukpMT1PPLzzz+jsrLS7f7x48f7Jbc+vXv3RvPmzf2eI+d4Dj73A4hMBDBOVR923r4fwFBVnV5rnQznOnnO23uc65TU95i1tjP2eJ6I6DyQnZ2Nvn37erStiKSqqqUx6/riCMCnRCQGQEygx0FE9HvniwLIB3BxrdvdnMvqWydPREIAtAVwpL4HU9X5AOYDQK9evXT27NlnDXc4HMjNzUVaWhpSU1Nx8ODBOuvMnTsX3bt3b+TuNF5paakrNz093e00DACMHDkSjz76qM9zbTYbdu7c6cr+9ddf3e4PCgrCJ598glatWvk8u6ioCKmpqUhLS0NmZibsdrvb/XfffTcmTJjg89zq6mpkZGS4so8ePep2f5s2bfDxxx8jKMi3b2xTVeTn57tyc3Jy6kw0T58+Hdddd51Pc4HTp562b9+O1NRUbNu2rc5pmJ49e+JcPx+ecDgc2L9/P1JTU5Gamop9+/bVWefvf/87LrvsMp9nl5WV4ZdffkFqaiq2b99eZ9J12LBheOaZZ3yea7fbsWvXLtc+//bUJgB88MEHaN++vc+z69OlSxdDcnwxCRwCYC+AS/B/k8ADfrPO43CfBP62MY/tySTw/v379f3339exY8dqs2bNFIA+9thjTX6cpqo9adStWzcFoK1atdLi4mK/5jocDs3OztY333xTr7/+eg0KClIAOmfOHL/mqqoeP35clyxZon/60580IiJCAWiPHj20qqrKr7kOh0NTU1P11VdfVYvFogAUgH777bd+zVVVLS4u1q+++konTZqkrVu3VgB67bXXqsPh8GuuzWbTpKQkff7557V///6ufd68ebNfc1VV8/Ly9KOPPtJbbrlFmzdvrgD03nvv9XvuqVOndNWqVfrEE09oVFSUAtBmzZrpoUOH/J69e/dufeedd/TGG2/UkJAQBaAzZ870e64voAmTwF7PAQCAiPwRwLs4/Y6gz1X1HyLyd+dAlolIcwALAVwJ4CiAKeqcND4bi8WiKSkpHo/rxIkTWL16NVatWoU5c+agbdu2Hj9WU6jzbWNWqxV9+vTB3XffbUguABw9ehTLly9HcnIy3n77bYSGhhqSa7fbsXnzZlitVowePRo33nijIbkAUFBQgLi4OOzfvx+vv/46RMSQ3OrqaiQmJsJqtWLatGkYMGCAIbkAsHfvXlitVtTU1ODZZ581LLeiogJr165FfHw8Zs2ahQ4dOhiSq6rIzMyE1WpF586d8eCDDxqSC5w+0l+5ciUSEhLw9ttvGzI5642mzAH4pAD8xdsCICIym6YUAD8JTERkUiwAIiKTYgEQEZkUC4CIyKRYAEREJsUCICIyKRYAEZFJsQCIiEyKBUBEZFIsACIik2IBEBGZFAuAiMikWABERCbFAiAiMikWABGRSbEAiIhMigVARGRSLAAiIpNiARARmRQLgIjIpFgAREQmxQIgIjIpFgARkUl5VQAiEi4iq0Vkt/PPdg2sZxeRbc7LMm8yiYjIN7w9AngBwFpV7Q1grfN2fSpVdbDzcpuXmURE5APeFsDtABY4ry8AcIeXj0dERAbxtgA6quph5/VfAXRsYL3mIpIiIskiwpIgIjoPhJxrBRFZA6BTPXfNrH1DVVVEtIGH6aGq+SLSE8A6EUlX1T0N5MUAiAGA7t27n2t4RETkoXMWgKqOaug+ESkUkc6qelhEOgMoauAx8p1/7hWR9QCuBFBvAajqfADzAcBisTRUKERE5CVvTwEtA/An5/U/AfjxtyuISDsRCXNejwAwHECWl7lEROQlbwvgDQCjRWQ3gFHO2xARi4h86lynH4AUEdkO4CcAb6gqC4CIKMDOeQrobFT1CICb6lmeAuBh5/VNAAZ6k0NERL7HTwITEZkUC4CIyKRYAEREJsUCICIyKRYAEZFJsQCIiEyKBUBEZFIsACIik2IB/EZNTQ2Kiur9lUZ+l5+fH5DcyspKHD161PBcVQ3YPh8/fhwnTpwwPNfhcKCgoMDwXAAoLi5GdXW14bk2mw2FhYWG5wKnf6ZUjf+VYqdOnUJJSYnhuU3FAgBw5MgR/Oc//8Hdd9+NDh06GFYAdrsdGzduxIsvvojLL78cr732miG5AFBQUID58+fj1ltvRY8ePWC32w3Jra6uxpo1a/DUU0+hV69eWLhwoSG5ALBnzx68++67GDVqFKKjoxEWFmZIbkVFBZYtW4Zp06ahW7duSExMNCRXVZGRkYHZs2fj2muvxR133IHQ0FBDsktLS/H111/jvvvuQ8eOHbFv3z5Dch0OB7Zs2YKXXnoJgwcPxowZMyAihmQXFhbi888/x4QJE9C1a1dUVlYakusVVT1vL9HR0eoPDodDs7Ky9M0339TrrrtOg4KCFIAC0MmTJ/sl84zjx4/rd999pw888IBGRES4ckNDQ3X//v1+y3U4HJqamqqvvvqqRkdHu3IB6PPPP++3XFXV4uJiXbBggU6aNEnbtGnjym3fvr2WlZX5LbempkYTExP1ueee0379+rnt80cffeS3XFXVQ4cO6Ycffqi33HKLNm/e3JXbv39/tdlsfss9deqUrly5UqdPn65RUVFu+7xixQq/5aqq7t69W99++20dOXKkhoSEuHLHjRvn19zy8nJdunSpPvTQQ9qxY0dXrohoVlaW33IdDodu375dX3/9dR06dOiZX4evAPSxxx7zW+65AEjRRj7Higbg8KixWrdurVdeeWWTt+vbty/mz5/vtqympgaJiYmwWq2wWq3Ys6fe30aNyy67DB06dPBovAAwb948DBzo/quP9u3bh9jYWFitVqxfvx41NTV1tmvVqhU82dczRo4ciVmzZrktq6ysxNq1a2G1WhEbG9vgqYdBgwahTZs2Hmf/97//RUREhOu2qiI7O9v1d/3zzz/D4XDU2S4iIgL9+vXzOHfq1Kl47LHH3JYdP34cK1euhNVqRXx8fIOntoYOHYpmzZp5lBsSEoJ169a5LXM4HEhLS3Pt8y+//FLvtt26dcMll1ziUS4A/O///i/uuMP9O5WKi4sRHx8Pq9WKlStXory8vM52IoLhw4d7/Gr44osvxqJFi9yW2Ww2bNq0yfX/a+fOnfVu26tXL3Tp0sWjXAB48803MWzYMLdlhw4dcv1MrVu3DlVVVXW2a968Oa6++mqPc4cMGYK33nrLbVlVVRV++uknV/bBgwfr3XbAgAEIDw/3OHvhwoXo0aOHR9uKSKqqWhq1cmObIhAX1Hr10pSLxWJxtaHD4dAff/yxzqshf12SkpJc2fn5+Tp16lRDcidNmuTKtdvt+tlnn7kdYfjzkp+f78rOycnRMWPGGJI7Y8YMV25VVZXOmTNHW7Zs6ffc0NBQt1dcW7ZsUed3V/j98sEHH7hyT5w4oc8995zbq21/Xfr06eO2zytWrNDevXsbss/x8fGu3KKiIv3LX/5iSO7YsWPdnkcWLVqknTp1MiQ7OztbPYXf0xHAVVdd1eTtLrvssnqPAJKSklzNvXv37ga37dixoW+2PLf333+/zhHA/v373Y4A6puIa9WqFTzZ1zNuuOGGeo8A1q1b58puaMJ10KBBaNu2rcfZ33//fZ0jgJ07d7peGW7cuLHBI4D+/ft7nDtlypQ6RwBlZWVYtWoVrFYr4uLicOTIkXq3HTp0qMdzACEhIVi7dq3bMlV1HQHExsYiNTW13m27deuGnj17epQLAM8880ydI4CSkhK3I4D6JrdFBCNGjPDqCOA///mP2zKbzYbk5GTXUU92dna921566aVeHQHMmTOnzhFAXl6e6//12rVrGzwCGDJkiMe5Q4YMwdy5c92WVVVVISEhwbXPBw4cqHfbyy+/3KsjgK+++opHAP6aA1BV3blzp86dO1evv/56DQ4OdjWvv+cAysrK9Pvvv9cHH3xQIyMjXbkhISF+nwNIS0vTWbNm6dVXX+32auO5557zW66qaklJiS5cuFAnT55s6ByAzWbTpKQkff7553XAgAFu+/zhhx/6LVdVNS8vTz/++GMdP3682xxAv379/DoHUFVVpatXr9Ynn3xSL7nkErd99vccQG5urr7zzjt60003GT4H8OOPP+rDDz/s9gpdRDQzM9NvuQ6HQ3fs2KH//Oc/9ZprrnGbA3j00Uf9lnsuaMIRQMCf5M928WcB1HbkyBFdtGiRTpkyRdu1a6fp6emG5NpsNt20aZO++OKLOnDgQI2JiTEkV1W1oKBAP/nkE73ttts0MjJSi4qKDMmtqqrStWvX6tNPP629evXS2bNnG5Krqrpnzx7917/+paNGjdJevXrpqVOnDMk9efKkWq1WjYmJ0S5duujixYsNyXU4HJqRkaFvvPGGDh8+XIcPH64Oh8OQ7NLSUv3mm2/0vvvu0/DwcP35558NybXb7bplyxZ9+eWX9corr9QpU6YYkquqWlhYqF988YVOmDBBw8PD9cCBA4Zl19aUAjivTwFZLBZNSUkxNLOmpgbHjh3zaiLYUwUFBV4dKnuqsrISlZWVXh2yekJVcfjw4YDsc1lZGYKCgtCqVStDc1UVv/76Kzp37mxoLnD6VFGbNm08nvz2lN1uR0lJiVenVj1VUFCAzp07G/ZW0DOqqqpQXl6O9u3bG5oLNO0UEAuAiOh3pCkFwA+CERGZFAuAiMikWABERCbFAiAiMimvCkBEJolIpog4RKTBSQcRGSciOSKSKyIveJNJRES+4e0RQAaACQA2NLSCiAQDmAfgZgD9AUwVEc8//klERD4R4s3GqpoN4FzvsR0CIFdV9zrX/RrA7QCyvMkmIiLvGDEH0BXAoVq385zLiIgogM55BCAiawB0queumar6o68HJCIxAGIAoHv37r5+eCIicjpnAah+2Z7yAAAHU0lEQVTqKC8z8gFcXOt2N+eyhvLmA5gPnP4ksJfZRETUACNOAW0F0FtELhGRZgCmAFhmQC4REZ2Ft28DvVNE8gBcAyBORFY6l3cRkXgAUFUbgOkAVgLIBvCtqmZ6N2wiIvKWt+8CWgpgaT3LCwD8sdbteADx3mQREZFv8ZPAREQmxQIgIjIpFgARkUmxAIiITIoFQERkUiwAIiKTYgEQEZkUC4CIyKRYAEREJsUCICIyKRYAEZFJsQCIiEyKBUBEZFIsACIik2IBEBGZFAuAiMikWABERCbFAiAiMikWABGRSbEAiIhMigVARGRSLAAiIpNiARARmZRXBSAik0QkU0QcImI5y3r7RSRdRLaJSIo3mURE5BshXm6fAWACgI8bse5IVS3xMo+IiHzEqwJQ1WwAEBHfjIaIiAxj1ByAAlglIqkiEnO2FUUkRkRSRCSluLjYoOEREZnPOY8ARGQNgE713DVTVX9sZM4IVc0XkQ4AVovITlXdUN+KqjofwHwAsFgs2sjHJyKiJjpnAajqKG9DVDXf+WeRiCwFMARAvQVARETG8PspIBFpKSKtz1wHMAanJ4/9JicnByUlxs83HzlyBDk5OYbnAkBmZiaOHz9ueO7hw4exd+9ew3Orq6uRkZGByspKw7MPHjyIvLw8w3NPnjyJjIwM2O12w7Nzc3NRWFhoeO6xY8eQlZVleC4AZGVl4dixY4bnFhYWIjc315gwVfX4AuBOAHkAqgAUAljpXN4FQLzzek8A252XTJw+ddSox4+Ojtam+vnnnxWAerKtt4YOHaoAdMOGDYbmLl68WAHohAkTDM09deqUdu/eXZs1a6a5ubmGZv/zn/9UAPrXv/7V0Nxff/1VW7durREREXrs2DFDs2NiYhSAvvfee4bm7tixQ4ODg7Vv375qs9kMzR49erQC0NjYWENz4+LiFICOGjXK0FybzaZ9+/bVoKAg3b59u0ePASBFG/kc6+27gJYCWFrP8gIAf3Re3wtgkDc5TREeHo7Q0FB07NjRqEiXjh07IiQkBBEREYbnAkCnTvVN1fhPaGgowsPDUVZWhrZt2xqaHah9btmyJVq3bo0WLVqgRYsWhmYHap/btWuHsLAwREZGIjg42NDsjh07QkQM/3nu0KEDRMTwv+vg4GBERkbi4MGDaNeund/z5HRhnJ8sFoumpPBzY1Q/VTXdW5DNuM/UNCKSqqoNfjC3Nv4qCLpgmfGJ0Iz7TP7DAiAiMikWABGRSbEAiIhMigVARGRSLAAiIpNiARARmRQLgIjIpM7rD4KJSDGAA3546AgAF/KX01zo4wcu/H3g+APvQt8Hf42/h6pGNmbF87oA/EVEUhr7Sbnz0YU+fuDC3weOP/Au9H04H8bPU0BERCbFAiAiMimzFsD8QA/ASxf6+IELfx84/sC70Pch4OM35RwAERGZ9wiAiMj0TFkAIvKaiOwQkW0iskpEugR6TE0lInNFZKdzP5aKyEWBHlNTiMgkEckUEYeIXDDv5BCRcSKSIyK5IvJCoMfTVCLyuYgUiYhfv5bVX0TkYhH5SUSynP9/ngr0mJpKRJqLyBYR2e7ch1kBG4sZTwGJSBtVLXNefxJAf1V9NMDDahIRGQNgnaraRGQOAKjq8wEeVqOJSD8ADgAfA3hWVc/7b/4RkWAAuwCMxumvQt0KYKqqBuZLaz0gItcDKAfwlapeHujxNJWIdAbQWVXTnN81ngrgjgvs30AAtFTVchEJBZAE4ClVTTZ6LKY8Ajjz5O/UEsAF14KqukpVbc6byQC6BXI8TaWq2aqaE+hxNNEQALmquldVqwF8DeD2AI+pSVR1A4CjgR6Hp1T1sKqmOa+fAJANoGtgR9U0zq/uLXfeDHVeAvIcZMoCAAAR+YeIHAJwL4CXAz0eLz0EYHmgB2ECXQEcqnU7DxfYk8/viYhEAbgSwObAjqTpRCRYRLYBKAKwWlUDsg+/2wIQkTUiklHP5XYAUNWZqnoxgEUApgd2tPU71z4415kJwIbT+3Feacz4iTwhIq0AfA/g6d8c0V8QVNWuqoNx+sh9iIgE5HRcSCBCjaCqoxq56iIA8QBe8eNwPHKufRCRBwGMB3CTnoeTOU34N7hQ5AO4uNbtbs5lZCDnefPvASxS1f8GejzeUNVSEfkJwDgAhk/M/26PAM5GRHrXunk7gJ2BGounRGQcgOcA3KaqFYEej0lsBdBbRC4RkWYApgBYFuAxmYpzAvUzANmq+nagx+MJEYk88649EWmB028qCMhzkFnfBfQ9gMtw+l0oBwA8qqoX1Cs5EckFEAbgiHNR8oX0TiYRuRPAewAiAZQC2KaqYwM7qnMTkT8CeBdAMIDPVfUfAR5Sk4jIYgA34PRvoiwE8IqqfhbQQTWBiIwAkAggHad/fgHgb6oaH7hRNY2IXAFgAU7/HwoC8K2q/j0gYzFjARARkUlPAREREQuAiMi0WABERCbFAiAiMikWABGRSbEAiIhMigVARGRSLAAiIpP6/7lbop9jyn1FAAAAAElFTkSuQmCC\n",
- "text/plain": [
- "<Figure size 432x288 with 1 Axes>"
- ]
- },
- "metadata": {
- "needs_background": "light"
- },
- "output_type": "display_data"
- }
- ],
+ "outputs": [],
"source": [
"plt.quiver(vdir[..., __spherical_indeces__['azimuth']],\n",
" vdir[..., __spherical_indeces__['elevation']],\n",
@@ -302,7 +318,7 @@
},
{
"cell_type": "code",
- "execution_count": 181,
+ "execution_count": null,
"metadata": {},
"outputs": [],
"source": [
@@ -313,9 +329,9 @@
" names=['position',\n",
" 'orientation'])\n",
"velocity = pd.Series(index=index_posvel, data=0)\n",
- "#velocity.loc[(convention,'alpha_0')] = 2*np.pi*(np.random.rand()-0.5)\n",
- "#velocity.loc[(convention,'alpha_1')] = np.pi*(np.random.rand()-0.5)\n",
- "#velocity.loc[(convention,'alpha_2')] = 2*np.pi*(np.random.rand()-0.5)\n",
+ "velocity.loc[(convention,'alpha_0')] = 2*np.pi*(np.random.rand()-0.5)\n",
+ "velocity.loc[(convention,'alpha_1')] = np.pi*(np.random.rand()-0.5)\n",
+ "velocity.loc[(convention,'alpha_2')] = 2*np.pi*(np.random.rand()-0.5)\n",
"velocity.loc[(convention,'dalpha_1')]= np.random.rand()\n",
"\n",
"elevation = np.linspace(-np.pi/2,np.pi/2,5)\n",
@@ -340,32 +356,9 @@
},
{
"cell_type": "code",
- "execution_count": 182,
+ "execution_count": null,
"metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- "<matplotlib.quiver.Quiver at 0x7f00b122e160>"
- ]
- },
- "execution_count": 182,
- "metadata": {},
- "output_type": "execute_result"
- },
- {
- "data": {
- "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADx0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wcmMyLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvMCCy2AAAIABJREFUeJzt3XlYVGX7B/Dvw764IbjjirigiCZqyYs7IQpptphaaWm9VqRmpZaamWm5/3JJc0lbNH2zRSEUwQX3FBTcMVApRQVBAUFZZu7fH8AJnGEZ5swZZrg/13UumWeOz30PzDz3Oc85c44gIjDGGKt5LIydAGOMMePgAsAYYzUUFwDGGKuhuAAwxlgNxQWAMcZqKC4AjDFWQ3EBYIyxGooLAGOM1VBcABhjrIayMnYC5XFxcaFWrVoZOw2msMzMTPz1119an3NyckLLli1haWmpaE6pqan4+++/4eDggJycHACAo6MjsrOz0bp1a9SvX1/RfAoKCnD9+nVkZGRofb5Dhw5wdHRUNCdWPcTExNwlogaVWbdaF4BWrVohOjra2GkwhRUUFKBp06ZITU2V2qysrPDee+/B09MTffr0QcuWLRXN6bvvvkPv3r2RmJiIgIAAAEB4eDjq1auHuLg4jB49WtF8EhMTcfjwYcTFxeGrr75CyUu6tGrVChcvXoQQQtGcWPUghEiq7Lo8BcSqHSsrK4wcObJUW0FBAVxcXPDKK68oPvgDwNixY+Hu7q7R3qlTJ8UHfwBwc3PDuHHj4OjoiMev5zV69Gge/FmlVOs9AGYY58+fR1ZWFmxsbEot1tbWGm2WlpZGGUxGjRqFVatWYeTIkTh16hSuXr2KmTNnol+/fujZs6fi+VRHBw4cwIIFCwAUFiI3Nzfs2rULo0aNMko+RISCggLk5+cjLy+v1PJ4W7169dCxY0ej5MlKIKJqu3Tv3p2Y/JKTkykgIIAAVGoZMmQIPXr0SNEc1Wo1tW/fnpKSkujPP/8kKysrAkBt2rShjIwMRXMpaffu3dLv5ciRI0bLIzU1lZo2bUoAyM7Ojs6dO0eXLl2irl27Kp5LdnY29evXr9LvpxEjRlBKSoriedYUAKKpkmOs0Qf58hYuAIajVqtp7dq15ODgUO6HtV+/fpSdnW2UHC9cuCD9vHDhQimnUaNGkVqtNkpO1aEAqNVqCgwMlPJYs2aN9FzJ35mSMjIyqFevXuW+l+rUqUPfffed0f52NQUXAFZpV65coSeffLLMD23Lli1p6tSpdOzYMVKpVEbLU6VS0aBBg6S8Nm/ebJQ8qkMBWLFihZTDs88+a9QBtaCggKKiomjSpEnUrFmzcjckrl+/brQ8axIuAEwn+fn59Pnnn0vTLGUtTZs2peDgYDpw4AAVFBQonmdycjI1aNCAAJCjoyPFx8crnoOxC0BsbCzZ2NgQAHJ1daW0tDTFc8jLy6OIiAiaOHEiNWrUqNz3jK2tLS1dutSoGw81DRcAViUxMTHUsWNH6cPbu3dv6t27t9YPdoMGDejNN9+k8PBwysvLUyzHP/74Q8qhW7duih+bMGYBePDgAXXo0IEAkIWFBUVFRSkW+9GjRxQaGkqvvfYa1a9fX+P9YGFhQX379qUePXpIbV27dqVz584pliMrxAXACPLz8+ngwYOy93vt2jU6fvy41ufS0tJo06ZNFBgYSN98840s8XJycmjKlCkEgF599VUiIrp58yatWrWK+vfvTxYWFhof/m3btskSu7Lee+89KfauXbsUjW3MArB161Yp9pw5cxSNvX79eo2/u6WlJfn5+dHatWvp9u3bREQ0fPhwsrCwoI8//phyc3Nlib106VIaPnw4/fDDD3T//n2t6xw6dIhu3LghS7yS9u/fb3J7LzW6AGzZsoV8fX3J19dX0YOXs2bNIgA0efJkWbZK1Wo1bdq0iWrXrk179+6V2m/fvk1r164lPz+/UlM2ffv21TtmSfv27aMXX3xRoz0lJYXWrVtH/v7+ZGVlRba2tpSZmSlr7Io8evSIAgICaM+ePYrGJSosulFRURQVFVXmYGRIv/32Gw0bNozy8/MVjXvnzh2ysLAgGxsbCgwMpE2bNmmdfgoKCqKjR4/KGrtbt27S+9za2pqGDBlCGzdupLt370rr/PLLL+Tk5CTbxkh2dja98cYbBIAWLVokS5+VcffuXWn8+vXXX6vUh+IFAMC3AFIAnC/j+X4AMgDEFi2fVKbfqhSAkmeLZGVl6fz/q+LAgQMkhCAA5OHhoXfhSUlJoeHDh0uv49ixY7RixQrq27ev1i3w+vXr04QJE2TfUnn48GG5z6enp1NkZKSsMVn1FR4eXuEpuBW9Z3SVl5dHr7zyCtWtW1frHsjAgQPp66+/psjIyFJniaWnp+sV9/79+9S6dWsCQFZWVnTy5EmZXlH5kpOTpdexevXqKvVhjALQB8ATFRSAUF37NYUCUPJ8bFtbWzp79qxe/YWEhFDDhg3LPbAGgBo1akQTJ06kiIgIRefgGTOG3NxcCgsLo/Hjx5Ozs3OFn49mzZpRRESEXjFPnDgh7WW7ubkp8v0TpQuALJeCIKJDANLl6MuUEBHGjx+P5ORkAMCyZcvg6elZpb6ysrLwxhtvICgoCCkpKVrXcXV1xaRJkxAVFYWbN29izZo1GDRoEKytrav8GhgzBTY2NggICMCGDRtw+/Zt7Nu3D2+//TYaN26sdf2bN2/Cz88PkyZNki7ep6tevXph3rx5AAqvvfTOO+9UOf9qq7KVoqIFQCuUvweQBiAOwG4AnSrTZ3XcAyg5zbJq1Sop1rBhw6p8Pvbhw4el3U1ti7W1Ne3evdvkDkYxZmgqlYp++eUXrVOjxUuHDh3o1KlTVe5/4MCBUl/ff/99qefkZpJTQFRxAagDoFbRz0MA/FVOP28CiAYQ3aJFC51fvCELQGZmJs2dO5eIiOLi4sjW1lba3Sx5QKqycnNzacaMGdLxg/IWb29vunfvnqyvhzFTl5KSQp07d67w82NlZUVz586t0sHzmzdvkouLCwGgWrVq0ZUrV4iIaObMmbKfhqx0AVDkYnBElFni5zAhxNdCCBciuqtl3XUA1gGAt7c3KZFfZf3+++9YuHAhXnnlFbz00kvIzc2FEAJbtmyBs7OzTn3l5OTgs88+w8WLFzFixIgyL8ZWsu306dMYMGCAgV4dY6bn9OnTGDNmTIUXn8vLy0N0dDTmzZuHjz/+GLa2tpWO0bRpU2zevBmBgYF48OABRo0aha1bt2LRokXo0aMHhg0bZsBXaFiKFAAhRGMAd4iIhBA9UXgZ6jQlYstp69atyMnJQa9evaRr1c+aNQt9+/bVuS8HBwd8+eWXcqfIWI3i7+8Pf39/g8cZOnQoJk+ejK+++goxMTH4z3/+g/z8fGzdutWkC4AsB4GFED8BOA6gvRDihhBivBBiohBiYtEqzwM4L4SIA7ACwEtFuyqySUvTXk+IqMzndJGamoqIiAjpZwDw8PDARx99pHffjLHqb968eWjdujWAf8eAkJAQZGVl6d13eWOUHONXWeQ6C2gUETUhImsiciWijUS0lojWFj2/iog6EZEXET1JRMfkiFvSuHHjMHXqVNy+fVtqi4yMhK+vL6KiovTu/+eff4ZKpSrVdvHiRdSpUwfDhw9Hdna23jEYY9VPRkYGBg8eDCcnJ1y7dq3Ucw8fPsTOnTv1jhESEoKBAwfi8OHDUtuNGzfwzjvvGPbso8oeLDDGostZQMuXL9d68MfGxkaWg8E+Pj5a+x8xYoRRvhHKGFPO3bt3y7yHRkBAgN7937lzp8yTQTZu3KhTX1D6ewDVQVBQkNb2/v37o1atWnr1nZSUhKNHj5Zqs7KywvLly7Fjxw7UrVtXr/4ZY9Wbs7MzQkND8fnnn8PCovSwuXfv3lL3r66Khg0b4sknn9T63NChQ/XquzxmUwDc3NzQoUMHjfbAwEC9+962bVupx66urjh06BCmTJnC915lrIawsLDAzJkzERkZiUaNGkntKpUKO3bs0Lt/bWNVz549S8WSm9kUAED7XkBZewa62Lp1q/Szv78/zpw5g6eeekrvfhljpqd///44c+YM+vTpI7WVHCOqylDjV3nMugB4enqiZcuWevV54cIFnD17FhYWFpg3bx7CwsLg4uKiV5+MMdPWpEkT7Nu3DzNmzAAAHDlyBH///bdefXbu3FljvOICoIOnnnoK9evXlx7L8cv76aef0LBhQ0RERGDWrFka83+MsZrJysoKX3zxBUJCQuDk5KQxVawrIUSpMat58+bo0qWLvmmWy6xGMysrKwQEBEiP9Z3/JyIkJyfjzJkz/A1cxphWgYGBOH36NBITE2Xpq+TPhj7GKEje72PJytvbm6Kjo3X6P9u3b8dLL72EBg0a4NatW7C0tKxy/IKCAgCFhYUxxsqTm5sLa2trvWYJcnNz4eLiggcPHiAsLKzUBm1lCSFiiMi7Muua1R4AAAwePBhWVlYYOnSoXoM/UDjw8+DPGKsMW1tbvaeIbW1t8fTTT8PR0RH9+/eXKbOymd3oVrduXfTp08fgB08YY8wQgoKCoFarYWdnZ/BYZlcAAOD555+Hn5+fsdNgjDGdDRkyRLFYZncMAADy8vJgY2NjgIwYY8zw9BnDavQxAAA8+DPGTJpSY5hZFgDGGGMV4wLAGGM1FBcAxhirobgAMMZYDcUFgDHGaiguAIzpKDc3V+P2oIyZIi4AjOlo//79OH78uLHTYExvZvlNYMYMKSQkBLVr18Z//vMfY6fCmF54D4AxHRARQkNDERoaauxUGNObLAVACPGtECJFCHG+jOeFEGKFECJBCHFWCPGEHHGZ+crJyamW8+xnz57FP//8g4sXL+Lq1avGTkdDfn4+Hj16ZOw0mImQaw9gM4DB5TwfAMC9aHkTwBqZ4jIzpVKp8NZbb0GtVhs7lVJCQkK0/lwdqFQqTJgwQe/LoLOaQ5YCQESHAKSXs8owAN9ToRMA6gkhmsgRm5mn2rVrY//+/Zg8eTKq0wULq2sBUKvVGD9+PM6dOwdra2tjp8NMhFLHAJoB+KfE4xtFbYyVycvLC6tWrcL06dOrRRG4ffs2Tp48KT2OiopCZmamETMqRER455138N1338HLy8vY6TATUu0OAgsh3hRCRAsholNTU42dDjOi4sFs8eLFmDt3rpGzAf74449SjwsKChAeHm6kbAoREd5//32sXbsWALgAMJ0oVQBuAmhe4rFrUZsGIlpHRN5E5N2gQQNFkmPVU9euXaWf586di4ULFxoxG+1TPsaeBpo9ezaWL18uPS75O2OsIkoVgF0AXi06G+hJABlEdEuh2MxEPb41O2PGDKxcudIouTx69AgREREa7WFhYUY7W2n+/PmYP39+qTbeA2C6kOs00J8AHAfQXghxQwgxXggxUQgxsWiVMABXASQAWA/gbTniMvPWokUL1KtXr1TbpEmTsH79esVzOXDgAHJycjTa09LScOLECcXzWb58OWbNmlWqrUWLFnByclI8F2a6ZPkmMBGNquB5AvCOHLFYzSGEgJeXF6Kiokq1//e//4W9vT1efvllxXLZt28fpk6dCk9PT7z22msAgG3btuHw4cOIiIiAj4+PYrmsWbMGU6dO1WjnrX+mK74UBKvWtBUAIsLYsWNhZ2eH559/XpE85s2bB3t7e+zZs0dqc3V1xapVq/Dw4UNFcgCA7777Dm+/rX0HmgsA01W1OwuIsZK0HdRct24drl27Bl9fX8XysLe3r9Jzcnv66aeRlJSEZcuWaTzHB4CZrngPgAEo3KresGEDjh07hrZt20qLm5ubxjy8koq3av38/HD06FHk5ORgw4YNmDBhAoQQRsvLWJo0aQKVSoUNGzYAAJycnNClSxdERUUZfQ/g3r17SEhIkJbExEQMGjRI0ak6phsuAAxA4Xz7hAkTcO3aNY2Diy4uLnBzcytVGIoXFxcXg+bl4eGBNm3aYNu2bVi8eDG+/PJLnDx5En/88QcCAwMNGru62r59Oy5evAgA+PDDDzFu3Dj4+vqiTZs2Bo+dkpJSapAvOdinp5e+GMCCBQswZswYg+fE9EBE1Xbp3r07MeWtXbuWLCwsCECZi729PW3dulWRfBISEoiI6O7du1S7dm0CQF27diWVSqVI/JJ2794t/Q6OHDmiePz8/Hxyd3cnANSgQQPKysoion9/R4a2YcMGsrGxKfe9YW1tTT/88IMi+TBNAKKpkmMsHwNgGv773/9i586dcHBw0Pq8nZ0dDh48iFGjyj35SzZubm4AAGdnZ+nsl9jYWPz222+KxK9OfvjhB/z1118ACr8XUatWLQD//o4Mbfz48QgPD4eVlfbJgzp16mD37t087WMqKlspjLHwHoB2ubm5dO3aNYPHOXnyJDVo0EDrVl7Tpk1p0aJFdP/+fYPnUdL9+/fJycmJAJCHhwcVFBQoGt+YewC5ubnUsmVLAkBNmjShnJwcReOnpaXRvHnzynxPNGvWjM6ePWvwPBITExX/u5sS8B6A+crMzMTQoUNx6tQpqe3MmTP4+++/ZY/Vo0cPHD9+HO7u7hrPJScnY9q0aWjevDm2b98ue+yy1K1bFx9++CEA4OLFi4rGNrZvv/0WSUlJAICZM2cqevbRpk2b0KJFC8yePRvartHVuXNnnDhxAp6enrLHTkhIwPnz/95q5ODBgxgxYoTWL+Yx3ZhdAbh06RK2b9+O7du3Iz8/X7G4169fR1BQEKKjow0W4+bNm/D19UVkZCTs7OwQEREBPz8/PPHEE/jiiy8MEtPNzQ3Hjh3DU089BQDo06cPNmzYgA4dOgAAsrKy0LFjR4PELsu7776Lhg0bYvjw4Yqf+dK4cWO8+OKLePHFFw1+APxxvXr1wpAhQ9CiRQtMmDBB0djt2rVDdnY2AMDT0xPff/89unfvDgAYMGAAjhw5AldXV4PEnjNnDjw9PREYGIioqCjY2tpi165d6N+/P1JSUgwSEwCOHj2KoUOH4s6dOwaL8biHDx9K41fxVJ9BVXZXwRhLVaaAFi5cKO2SFh8gU8Lrr78uxY2Pj5e9//Pnz1Pz5s2lGI0aNSq1++3o6EiZmZmyxy2Wk5NDzz77LPn6+hIRkUqlop07d9K0adMMFrM8aWlpRolbHRjrtU+ePJl2795NarWaiIi8vLzo5ZdfptzcXIPFTE1NJWtr61Lv9ZLvfTc3N7py5YrscWNiYqQYU6ZMkb3/siQnJ0txV69eXaU+oMMUkNEH+fIWUykAV65cIUtLSwJAQUFBsve/f/9+qlu3bpln4wQHB1NiYqLscR9XUFBAmzZtMngcZhq+/fZbqRgY0uXLl+mNN94gW1tbrZ8BFxcXOn78uKwx1Wo19e3blwCQra0t/fPPP7L2XxYuACZYAMaMGSPFPHPmjKx9b9myRWMLqHjx9/en1NRUWeMxVl3dunWLfH19y9wQ+v3332WNd+jQIan/t956S9a+y6J0ATC7YwBKyMvLK6yeAC5cuICtW7cCAJ5//nnZvo5PRFi4cCHGjBlT5rGMvXv3Yt26dVIujJkrlUqFlStX4vDhw1qff/jwIUaMGIHVq1fLFtPX1xdPP/00AGDDhg24du0agMLbbyp5fNGQuABUwZEjR/DNN98AKDxARUQQQuDTTz+VpX+VSoXg4GDMmDFD6/NNmzaFr68vxo0bBwAGOQOIserk6tWrsLe3x9ixY+Hj44PGjRtrrKNWqxEcHIzp06dDrVbLEnfevHkAgPz8fOnnJUuWIC4uTpb+jY0vBVEFcXFx+Oijj+Dg4IBffvkFADB69Gh06tRJ775zcnLw8ssvIyYmBgMGDNC49EKbNm3g6OiodxzGTIm7u7vGJUqysrJw9epVjctSbN26FTdv3sTGjRtha2urV9yePXsiKCgIISEh+P777zFgwADMnDkTa9euhbe3t159VwdcAKogNjYWubm5GDt2LADA0tISc+bMkaVvtVqNrVu3ws7OTpb+GDNXtWvXhpeXl9ZTgR8+fCjbndo+++wzhISEQKVS4ZVXXgFQOAaYAy4AVfD47p+VlRVGjhwJGxsbrFu3Dl26dKly38Vf7WeMVZ0cX5I7efIkgoODUVBQAGtr61Lz/uYyBWQ2xwDmzJmD+Pj4Um0qlQo//fQTQkNDZYuTl5cnXYmxWG5uLmJjYxEcHKzX4M8Yqz569uyJ8ePH48yZMxoHfePi4mQ9+WLHjh349ddfNY5dnD9/XuO+z3IymwKQkZEBDw8PrF27Vmrz9vbG6NGj0b59e9niXLp0SesZAN988w1fAIsxM/Pf//4Xy5cv12jPzMzE9evXZYvTunVrPPfcc+jXr5/UtnTpUnh6eiI3N1e2OI8zmwIQGBgItVotnaoFFF5DpF27dlqvZVNV2nb9VqxYgTfeeEO2GIyx6mPKlClat8LlnAbq1q0bmjRpgoSEBKnt6tWrAGDQ+16YTQHo06cP6tSpo9EeFBQka5zHD/4sXLgQ7777rqwxGGPVy8cff4yZM2eWapPzQLCFhYXWgb5x48YGPdvIbAqAjY0N/P39NdrlLgAlq/6cOXMwbdo0WftnjFVP8+bNw3vvvSc9lvtAsLaxaujQobCwMNwwbTYFAND8BTo5OcHHx0e2/olIqvoffvihbKd+MsaqPyEEli5diokTJwKQ/1TQgQMHapz+LfcG7ONkKQBCiMFCiHghRIIQQuPrq0KIcUKIVCFEbNFikGvZBgQElKqWAQEBZd65qCpu3ryJ9PR0BAcHY+HChTXypuSM1WRCCKxevRpjx47F9evXkZGRIVvfDg4OGDRokPTY1ta21GND0LsACCEsAawGEADAA8AoIYSHllW3E1HXomWDvnG1cXFxka5bD8h/8CQuLg7jx4/HV199xYM/YzWUhYUFNm7ciJEjR+Ls2bOy9l1yzBowYIDBv/Uvxx5ATwAJRHSViPIAbAMwTIZ+q6R4l8nS0hKDBw+Wte/WrVvjm2++MeicHGOs+rO0tMQPP/yABg0ayNpvyQJg6OkfQJ4C0AzAPyUe3yhqe9xzQoizQogdQojmMsTVqviX5uvrCycnJ1n79vDwgKWlpax9MsZMk7W1tXRnPLk0a9YMTzzxBADDnv5ZTKlN2RAArYioC4AIAN+VtaIQ4k0hRLQQIlrbvUcr0rFjR7Rp00aR6skYY3ILCgpC165d0by5wbaTJXIcIb0JoGSmrkVtEiJKK/FwA4BFZXVGROsArAMAb29vnb9rLYRAUFAQFwDGmEkKCgqS7XLWFZGjAJwC4C6EaI3Cgf8lAKNLriCEaEJEt4oePgPgkgxxyxQcHIy2bdsaMgRjjBlEt27dUK9ePUVi6V0AiKhACBEMIByAJYBvieiCEOIzFN6abBeASUKIZwAUAEgHME7fuOXhwZ8xZqosLCzg5uamSCxRnW8n6O3tTdHR0cZOgzHGTIYQIoaIKnX9CD6fkTHGaiguAIwxVkNxAWCMsRqKCwBjjNVQXAAYY6yG4gLAmI4uX76MGzduGDsNxvTGBYAxHYWEhCA0NNTYaTCmNy4AjOkoJCQEISEhxk6DMb3Jd7cUxmqAtLQ0HD16FNbW1sjOzjb49doZMyTeA2Ama9++fYrH3L17N9RqNXJzcxEZGal4fGO8Zma+uAAwk0REmDRpEvLy8hSNW3LqR+njAFlZWXj//fcVjcnMGxcAZpJiY2Nx8eJFHDp0SLGYeXl52LNnj/Q4NDRUscv2AkBERATi4uKQkJCgWExm3rgAMJNUvCWu5MHYw4cPIzMzU3p8+/ZtxMTEKBa/+LXyGUhMLlwAmEkqWQCUuqKttoFXqQKkUqnwxx9/KBqTmT8uAMzkJCcno/gy4deuXcPFixcNHpOItA68Sg3GJ0+eRPEtUg8dOoSMjAxF4jLzxgWAmZywsLBSj5WYErl8+TISExM12mNjYxX5VnDJ11hQUFDqWARjVcUFgJmcx7e6ldgKDwkJgbOzM5577jmpbeTIkahbt64iBcgYr5mZPy4AzKQ8fPgQERERpdqOHz+Ou3fvGjSul5cXrl27hgkTJkht7777Lq5fv27wW5AmJSXh3LlzpdrCwsJQUFBg0LjM/HEBYCZl//79ePjwYak2tVqtMS0kN39/f9SuXVujvV69ehg0aJBBY2vb2r937x6OHz9u0LjM/HEBYHohImRlZSkWr6zpFnM+NbKs16bkNFDJ01+Z+eACwPRy/vx5/PLLL4rEIiKEhobCx8cH7dq1AwA0b94cfn5+CA8PV/xbwUrIysrCwYMHERgYiIYNGwIAPD090aNHD0ULwObNm5GUlKRYPKYMLgBML0peGTM9PR2bNm3C4cOH0bp1awBAo0aNsHfvXoSFhSElJUWRPJSUmpqKI0eOICQkBE5OTgCA9u3b488//8TSpUvx4MEDRfLgK6CaJ1kKgBBisBAiXgiRIISYoeV5WyHE9qLn/xRCtJIjLjO+kJAQ7N27F7m5uQaP5ezsjEGDBkEIofGcj48PXF1dDZ6D0tq0aQNvb2+NdiEEhgwZglq1ahk8h8zMTERFRXEBMEN6FwAhhCWA1QACAHgAGCWE8HhstfEA7hFRWwDLASzUNy4zvpSUFPz555948OABDh48aOx0mIGEh4cjPz8fBw8eVPR4DzM8OfYAegJIIKKrRJQHYBuAYY+tMwzAd0U/7wAwUGjbjGMm5Y8//pAuw8Bbh+ar+G+bl5eHvXv3GjkbJic5CkAzAP+UeHyjqE3rOkRUACADgLMMsZkRnT17Ft7e3mjZsiWSk5MVuyYPU45KpcK9e/fQpEkTPPnkkzh9+rSxU2JyIiK9FgDPA9hQ4vErAFY9ts55AK4lHicCcCmjvzcBRAOIbtGiBelqxYoV5OzsTM7OzpSVlaXz/6+qR48eSXGXLFmiWNykpCSaPn06NWjQgN58803aunWrYrFVKhU9ePCA7t+/T2q1mtRqtWKx79+/T6mpqXTv3j3FYhIR5ebmUmpqKqWmplJeXp6isdPT0yk1NZUyMjIUi6lSqUitVlN6ejplZ2eTSqVSLPa3335Lr7/+OjVv3pzee+89Sk1NVSz2p59+Ss7OztS4cWPFYhIR3bp1SxpHNm7cWKU+AERTJcdvOW4JeRNA8xKPXYvatK1zQwiKi6J1AAAgAElEQVRhBaAugLQyCtI6AOsAwNvbW+dNyocPHyItTWvXBlcc9/EvKhlS8+bNsWPHDqSmpmLdunV46623FIttYWFhtFsi1q1b1yhxbWxs4OLiYpTYxWcBKcnCwsJosd3c3PD6668DAA4cOIBly5YpFjsnJwdpaWmwtrZWLCZQuEFePI48evTI4PHkmAI6BcBdCNFaCGED4CUAux5bZxeAsUU/Pw9gf1GlYnoSQiAoKAhAYTHw8vIyckaMyaN3795S4Sl+jzN56V0AqHBOPxhAOIBLAP5HRBeEEJ8JIZ4pWm0jAGchRAKAqQA0ThVlVRcYGCj9y8fWmbmwsrLCkCFDAHABMBRZvgdARGFE1I6I3IhoflHbJ0S0q+jnR0T0AhG1JaKeRHRVjrjVwZUrV5Cdna31uTNnziiSg6+vL+rUqcMfEmZ2goKC0LhxY3Tv3t3gsYgIsbGxWp/LzMzUejlwU8ffBNZTcnIyWrdujaVLl0ptZ8+ehb+/P2bNmqVIDjY2NnjuuefQv39/ReIxphR/f388++yz0rEIQxJCYNKkSQgKCsLly5cBFF5ocMGCBWjdurVZ3oRHjoPANZqPjw/y8/Mxc+ZMqe3nn38GAKxZs0axPD7//HPY2dkpFo8xJdSrV0+xDSmgcBp1+vTp0mOVSoWZM2eiadOm6Natm2J5KIX3APRkbW2NgIAArc8Vz80roWnTporFYkxJSr63y5pGNdfja2ZTAM6ePau1/eHDh7hy5YpBY2t703Tr1s0sr03DmDnr0KED3NzcNNoNfXzt8uXLZV7NtqyxTQ5mUwA++eQTDB48GHFxcVLbqlWr4ObmZvDL2A4ePBiWlpal2viALGOmp+Rp1cXs7e0xcOBAg8a9cOEC3N3dsXHjRqnt1KlT6N+/P5YsWWK4wJX9xpgxlu7du1f622/r168nABpLrVq16NGjR5Xup6r69u1bKu7JkycNHpMxJr/IyMhSn+WgoCCDx8zIyCBra2utY9j27dt16gs6fBPYbPYAhg4dqrXd398ftra2Bo9fcqtBqdPWGGPyKz6tupgSe/N16tRBv379NNqtrKzg7+9vsLhmUwCaNGmCHj16aLQrNRVTMk5gYKAip60xxuRnY2ODwYMHS4/L2riUm7axqm/fvga97IlZjVKP/wKLb5qhhHbt2km3KeT5f8ZMW/FnuHv37oqdhaRt3DD0WGLWBeCpp55CgwYNFI1vZ2eHQYMGKRaTMSa/gIAAWFhYKLox16pVK3Tu3LlUm6FPJTerAuDl5VXq1Eslz8MvjjdgwAA4ODgoGpcxJi9nZ2f07t1b8b35kmNWx44dtZ6SKiezKgBCiFK/QKX/eD4+Phg7dmzFKzLGqr233npL8W//lhyzlBi/zKoAAP/+0lq1aoVOnTopGtva2hovvPCCojEZY4YxatQoxb/926tXL+l+E1wAqqB4CiYoKMgoX902x6+LM1YTGeOzbGlpiaFDh6J+/fp46qmnDB7P7C4GV3wQVun5f8YYk0NgYCBUKpXG1QUMwewKAACMGTMGffv2NXYajDGms6efflqxW1EKqsZ3ZvT29qbo6Gid/x8R8VQMY8xk6TOGCSFiiMi7Muua3TEAgOfhGWOmTakxzCwLAGOMsYpxAWCMsRqKCwBjjNVQXAAYY6yG4gLAGGM1lF4FQAhRXwgRIYT4q+hfpzLWUwkhYouWXfrEZMzYMjMz8ejRI2OnwZje9N0DmAFgHxG5A9hX9Fibh0TUtWh5Rs+YjBlVeHg4Dh48aOw0GNObvgVgGIDvin7+DsBwPftjrNoLCQlBSEiIsdNgTG/6FoBGRHSr6OfbABqVsZ6dECJaCHFCCMFFgpkslUqFsLAwhISEoDp/i56xyqjwWkBCiEgAjbU8NbPkAyIiIURZn4iWRHRTCNEGwH4hxDkiSiwj3psA3gSAFi1aVJQeY2XKzc1FVlaWdHldORw/fhxpaWlIS0vD2bNn4eXlJVvfd+7cgbOzM6yszPISXawaqnAPgIgGEVFnLctOAHeEEE0AoOjflDL6uFn071UABwGUeZcFIlpHRN5E5K3k7RyZecnPz8dLL70EtVota78lp35CQ0Nl7TsrKwvjxo2DSqWStV/GyqLvFNAuAMW3wBoLYOfjKwghnIQQtkU/uwDwAXBRz7iMlUmlUuHVV1/FiRMn0LBhQ1n7LlkA5D4O0KZNG+zcuRMTJ07k6SWmCH0LwJcA/IQQfwEYVPQYQghvIcSGonU6AogWQsQBOADgSyLiAsAMQq1WY8KECdi2bRu6du0qa9+JiYm4dOmS9PjkyZO4c+eObP1bWFigS5cu2LBhA6ZMmcJFgBmcXgWAiNKIaCARuRdNFaUXtUcT0YSin48RkScReRX9u1GOxBl7HBEhODgYmzdvBgBZ5+cBzSkfIsIff/wha4zinFesWIGPPvqIiwAzKP4mMDMLRIQPPvgAa9askdrk3gPQNuUj93GAkjkvXLgQ8+bNk7V/xkriAsDMwpw5c7Bs2bJSbXLuAWRkZCAqKkqjfe/evbJ+K/jxnOfMmYPFixfL1j9jJXEBYCbviy++0NhStrOzg7u7u2wxwsPDUVBQAEdHR6nN0dER2dnZsn4ruHPnzho3A5k2bRpWr14tWwzGinEBYCbt//7v//Dxxx9rtHt6esp6Pv2FCxewYcMG/PTTT1Lbrl27sHr1asTGxsoWx9HREe3atdNoDw4OxsaNfPiMyYu/ccJM1tatWzF79mzY2dlpTMPIfQD4008/hRACe/bskdpsbW3x9ttvy36g1svLC/Hx8aXa7O3tMXXqVLi4uGDYsGGyxmM1F+8BMJM1evRoZGVl4bffftN4Tu4DwOXdo1Xu+7dqyz0iIgIZGRk8+DNZ8R4AM6jMzEwkJCQgISEBiYmJ0s89e/bE/PnzYWNjo1f/RITZs2cDAGrVqoXXX38dK1askH0PQEnFub///vtYuXIl8vLyMHv2bOzfv1/vvh89eoQPPvgAFy5cQNu2bdG2bVu4ublJ/9auXVvvGMx0cAFgsoiPj0dMTIw0wBcvqampGutOnz4dCxYsgIWF/jugu3btQnR0NABg8uTJmDt3Ls6ePYsuXbro3bexeHl5YejQoVi8eDHy8vKwcuVKHDhwAPv378eAAQP06tvOzg5fffUVJk2ahK+//lrj+UaNGkmFoXjp0aMH3Nzc9IrLqikiqrZL9+7diZmGv//+m7y9vQlAmYuFhQWtXr1atpgqlYq6dOlCAKhu3bqUnp5ORERpaWmyxXjc7t27pddz5MgRg8RQq9XSa0hOTiY7OzsCQL179ya1Wi1bjIULF5b79wJAvr6+dPv2bVliMmUAiKZKjrF8DIDJonnz5ggLC0PLli21Pm9vb49ff/0Vb7/9tmwxd+zYgbNnzwIonC5xciq8IV39+vVli2EMQgjpNTRp0gTBwcEAgGPHjiE8PFy2GNOmTcOWLVtgbW2tdZ127dph165daNSorKu8M1PHBcBM5efn48iRI4rEunXrFmbMmAF3d3ckJSVpPO/i4oL9+/fLegBTpVJhzpw5AAoH/MmTJ8vWd3Uzbdo06fsHs2fPlvWso9GjRyM8PBx169bVeO7KlStwc3PDJ598gpQUrRf6lR3faU1hld1VMMZiqlNAKpVK8Zj5+fm0aNEiSklJIaLCXXw7Ozv6+eefDRr3+PHjZGNjU2rawMrKSvrZzc2N/vrrL9nj/vjjj1KMhQsXyt5/WZSYAtJm5syZUtydO3fK3v+5c+eoefPmWv+GAMjOzo7Onj0re9ySNm/eTPXr15ceJycn0+LFi2Wb9tKFMT7DcoEOU0BGH+TLW6pSAI4cOUJz5syhOXPmUG5urs7/Xx/x8fE0evRoeuWVVxSNS0S0fft26YP69ttvU0JCAtnZ2ZEQgpYtW2awuHl5edLA0atXL/r111/pww8/lB4XFyS5ZWZm0vz586lDhw704MEDg8TQ5q+//pLeX0lJSYrFTU9Pp7Zt29KSJUsoOzvbIDFu3LghHVOZO3cubdu2jZ544gkCQB06dDDYoKhWq+mzzz4jANSwYUO6fPkyTZgwQdqwCAsLM0jc8owYMYJee+01unbtmqJxMzMzpffXyZMnq9RHjS4AJQ9sZWVl6fz/9fHCCy8QABJC0Pnz5xWN/eqrr2occC35eMqUKQb7AO/cuZMOHTokbalNnjyZnnnmGYMNVCWZ8paarpR4rRkZGTRo0CCaPXs2ERUOzpGRkbRnzx6DxMvLy6MJEyaUet8W3VlQWoKDgw0Suyx//vmnFHv8+PGKxk5OTpZiV/WECS4ARioA58+fl968zz//vKKx1Wo1RURE0NNPP13mGR3PPfcc5eTkGDyXQ4cOUUFBgcHjMMPIzc2lo0ePGjxOVlYWBQQElPl+DQwMLLVhoRR/f38CQNbW1orvAXABMOECQEQ0ZswYKf6ZM2cUj09EtH79+jI/VD4+PnT37l2j5MVYsVu3bknTS9qWn376ySh5HTp0SMph4sSJisdXugDwWUAymzNnDiwtLQEAn3zyidSenp6uSPzNmzdLpw1qc/ToUfj4+ODatWuK5MPY4y5fvoynnnoKp0+fLnOd8ePH43//+58i+RR/Non+/Va5ra0tZs6cqUh8Y+ICIDN3d3eMHVt4m+SQkBD8+eefuHr1KsaPH2/w2AUFBWjWrBmWL1+OqVOnYtiwYejUqRPs7OxKrRcfH48nn3xS+gYtY0o5cuQIevfujevXr5dqt7e3h6enJ4YPH44PPvgAy5Ytg7OzM9RqtcFzGjNmDG7evIn9+/dL93yYOHEiXF1dDR7b2PhSEAYwe/Zs/PDDD8jPz8fMmTORmZmJxMREEJHsFw4rycrKCn5+fvDz8yvVrlarkZycrHFNnilTpuCLL76Ar6+vwXJirFhERAQ+/fRT+Pn5SdcfKl6aNGli0M9GWQoKCnDgwAG8/PLL0hVlHRwcMGPGDMVzMQYuAAbQqlUrTJgwAWvWrMG+ffuk9hs3bqB58+aK52NhYQFXV1e4urqiX79+isdnDIDWjRNji4+PR25ubqkvoAUHB6Nx48bGS0pBZjMFVHjsQ/fn5KRWq/Hqq6/C3t4e69ev13g+Li5OkTwYY5Wj7TO5bNkyODg44K233lJs7DDW+GU2BWDChAkICQkp9ctKT0/H7NmzFTuYZGFhgfXr16Nfv34oKCjQeJ4LAGPVi7bPZEFBAYYOHYqVK1cqNi21adMmfP7558jMzJTa1Go1duzYUe5JHXqr7OlC2hYALwC4AEANwLuc9QYDiAeQAGBGZfvX5TTQ4q/KOzo6SqdRFf988+bNyp9DJYOcnBzq37+/xqltSn83gDFWPm3fmwkMDFT8KgKXL18uc/zS9VInUPA00PMARgA4VNYKQghLAKsBBADwADBKCOGhZ1wNQUFBAIDs7GypLTs7G927d0fTpk3lDlcue3t77Nq1C7179y7VznsAjFUvj38mBw0ahJ9//lnvGxXpqn379nB3d9cYv4B/xzZD0KsAENElIoqvYLWeABKI6CoR5QHYBkD2+9r16NFD62VrDfnLK0+tWrUQFhaG7t27S20JCQl48OCBUfJhjJV2+/Zt3LlzR3rs6+uL33//XeO0aaVoG6vc3NzQoUMHg8VU4hhAMwD/lHh8o6hNVhYWFhg6dKhGu7EKAADUrVsX4eHh8PT0BFA43Xbu3Dmj5cMY+1fJrf9evXohNDRUuuy2MWgbq4KCggx6HKLCAiCEiBRCnNeyGOTu1EKIN4UQ0UKIaG23EyxPYGBgqcdNmzZFt27d5ExPZ87OzoiMjJSqOE8DMVY9FH8Wu3btit27d6NOnTpGzcfHx0fjvgyPj2lyq7AAENEgIuqsZdlZyRg3AZQ8+d21qK2seOuIyJuIvBs0aFDJEIX8/PxKzd0FBgYa5cslj2vYsCEiIyPRpk0bxMbGGjsdxhiA2NhYeHh4YO/evdLd5IzJ2toaAQEB0uM6deoY/EuaSkwBnQLgLoRoLYSwAfASgF2GCFSrVq1SN8025vTP45o1a4Z9+/bh3r17xk6FMQYgJycHkZGR0HVD05BKjlmDBw82+MFovQqAEOJZIcQNAE8B+EMIEV7U3lQIEQYARFQAIBhAOIBLAP5HRBf0S7tsxb9Ae3t7DBw40FBhqqRVq1ZYvny5sdNgrMYjInz99ddo0qSJsVMpZfDgwdLFJJXYgBWk0DfdqsLb25t0vWBZUlISWrVqhcDAQISEhBgoM8YYM4y+ffviyJEjSElJgbOzs87/XwgRQ0TelVnX7K4F1LJlS3Tp0qVaTf8wxlhlBQUFQa1WV2nw15XZFQAAeOaZZwx+9JwxxgyhuAAoweymgAAgNTW1Wh3YYYwxXegzhukyBWQ2F4MriQd/xpgpU2oMM8sCwBhjrGJcABhjrIbiAsAYYzUUFwDGGKuhuAAwxlgNxQWAMR3Rv3e6Y8ykcQFgTEfR0dGIj6/oPkiMVX9m+U1gxgwpJCQEtWvXNuidmhhTAu8BMKajkJAQvtAgMwu8B8BqDCLS+wZB//zzD2JjY2FhYYH09HTUr1/f6DkxVlW8B8BqhH/++QeLFy/Wu5/Q0FAAgFqtxu7du/Xub8GCBbh9+7be/TBWFVwAmNnbs2cPunXrhlatWundV3EBACDLNFCjRo3QrVs3REVF6d0XY7riAsDMlkqlwieffIIhQ4YgNzdX70uEZ2dnY9++fdLjPXv2ID8/X68+n3vuOaSlpWHAgAH48ssvFbsMMGMAFwBmpu7cuQN/f3/MmzcPRIThw4fDwcFBrz4jIyORm5srPc7IyMDhw4f16tPJyQkBAQFQq9X46KOPMGzYMKSnp+vVJ2OVxQWAmZ3Dhw+jW7dupbbWR48erXe/2qZ8Sk4JVVXJ3EJDQ/HEE0/g1KlTevfLWEW4ADCzoVarsWjRIvTv3x+3bt2S2p2dnTFo0CC9+9Y22IeEhOj9reCgoCA4OjpKj5OSkuDj44NVq1bxN46ZQXEBYGbh3r17GD58OKZPnw6VSlXquRdffBHW1tZ69R8dHY07d+5otCckJOj9rWAHBwc8++yzpdry8/Px7rvvYtSoUcjKytKrf8bKwgWAmbyrV6/C29u7zLNy5Jr+6datG6ZOnSq1TZ8+HZ07d5blbKBRo0Zpbd++fTuefPJJ3LhxQ+8YjD2OCwAzeW3atMG5c+cQFRUFe3v7Us81b94cvXv31jvGyJEjERMTAz8/P6ktKCgIcXFxGDJkiN79+/n5wdnZuVRb3bp1ceLECcTExMDV1VXvGIw9Tq8CIIR4QQhxQQihFkKUeRNiIcR1IcQ5IUSsEEL3u7wzVgEHBwds3LgRDx8+BADY2NgAKNyytrDQfzunc+fOWr+xa2FhgU6dOundv7W1NV588UUA/+aekZGBH3/8EXZ2dnr3z5g2+l4K4jyAEQC+qcS6/Ynorp7xmAkhIuzfvx+ZmZnIy8srteTn52tt8/DwwLhx43S+PMKWLVvw/fffAwD69++PV199Fa+99pos0z9KGT16NNasWYMff/wRS5YswcmTJ7Fq1Sr4+fnhmWee0akvtVqN9evXIyEhATY2NrCxsYG1tbX0s7a2+vXro2/fvgZ6daxaKr62uT4LgIMAvMt5/joAF1377d69OzHTtmfPHrKxsSEA5S4WFhb00UcfUW5urs4xEhISqHbt2gSAnJ2d6caNG1RQUEDPP/88qdVqWV/P7t27pZyPHDkia98qlYqee+45UqvVlJiYqPGadJWTk0OTJ0+u8HcPgBwcHOjQoUOyvh5mHACiqbJjd2VXLLeTigvANQCnAcQAeLOCvt4EEA0gukWLFob7LTHFrFy5stzBp02bNlUeTHNzc6lHjx5SX7t27ZKey8nJkeslSAxZAIhK57xlyxYpVv/+/amgoKBKfUZGRpKrq2u5f4PNmzfL9RKYkclaAABEonCq5/FlWIl1KioAzYr+bQggDkCfyiTHewCm68KFC/TZZ5+Rl5dXuQPPG2+8QVlZWVWOM336dKmvd999V8ZXoJ2hC8Djxo4dK8WbP39+lfu5d+8evfzyy+X+Lby9vemLL76gK1euyPgKmNKq3R7AY+t+CuCDyqzLBcB0qNVqOnPmDM2aNYs6duxY4ZRDo0aNKCQkRK+Ye/fulfrr0qULPXz4UKZXUzalC0BmZia5u7sTALK0tKRjx47p1d///vc/ql+/foV/H09PT5o7dy6dP39e9mk0ZljVqgAAcARQu8TPxwAMrky/XAAqT61W09WrVytcR24pKSk0bdo0cnNz0zqQeHh40OzZs2nu3LlS27PPPkspKSl6xb1z5w41btyYAJC9vT1dvHhRpldUPqULABFRdHQ0WVtbEwBq2bIl3bt3T6/+bt68SYMHD5Zex6JFi+ijjz6idu3aaf0btm/fnj7++GO94z5OrVZX+J5MTEyUNWZNoFgBAPAsgBsAcgHcARBe1N4UQFjRz22Kpn3iAFwAMLOy/VelAGRnZ1NKSgqlpKQYZcslLCyMVqxYoXjcU6dOEQDy8vKiefPmaR0Qs7Ky6LXXXqOMjAzZ4t67d08anIqXbt260eeff06XLl2S1ps/fz7Vrl2bNm/eLMvf5dKlS9KAtW7dOr37q6zIyEhycXEhFxcXOnHihGJxly5dKm2ZyzEoqtVqWrNmDTk4ONDKlSultnPnztGcOXOoc+fOpf6mjo6Osu5hpaWl0SuvvEL5+fkaecXFxdEnn3xCnTp1IgAUHx8vW9zKWrRoEe3fv1/xuAUFBdL4VdVjWIrvARhqqUoBWLhwofSm1WduuSoWL14s2666rmbMmKGx5daxY0eaNWsWxcbGSoNus2bNqGXLlnTw4EHZYg8ZMoR69epFixYtKnNw2rhxI12/fl22mESFBe3rr7+uEVMUKpWKVq1aJfuB7StXrtD333+v9bn4+HhasGABde/enUaOHClbzN27d1OTJk2oY8eORFQ46J86dYpmzJghTXeVXBYsWCBb7MoonloUQtCaNWsUjZ2cnCy97tWrV1epDy4ARioAcu+q6+LChQs0e/Zs8vDw0Lob7+bmRtOmTaMOHTpIb+73339flq26qpy6yUyLHH/jBw8e0FtvvSW9J7t27UpTp06lli1ban3PdunShT777DNFD0o/PrV44cIFxWITcQEw6QJARLRkyRIp/siRI42ydXrp0iX6/PPPqWvXrhUe7OvcuTOdOXNG8RxZzXLixAmtW/ePLz169KAvv/yS/vrrL8VzVKlUFBAQIOWi5NRiMS4AJl4AVCoV+fv7Szls3LhR8RxKSkxMpEWLFlGvXr3K/NBZW1vTF198UeXzzBkrS15eHs2ePZssLS21vveEEOTj40PLly+npKQko+a6fPlyKS9DfImwMrgAmHgBICK6ffs2NWzYkIDCb1iWPBh6+/Zto+SUnp6ucWDv8cXHx4fPumCyuXjxInXv3r3c95y3t7dRP6fFYmJipOnbFi1aUHp6ulFyUroA8NVADaBRo0bSdWlycnLw0ksv4dGjR4iKisLcuXMVz+f+/fvw9/fH+fPny13v6NGj6NKlC9avX1+4dcBYFajVaqxYsQJPPPEEYmJiyl03OjoaQ4cORXZ2tkLZ/evDDz/E6dOn8eDBA4waNQr5+fmwsLDA1q1b4eTkpHg+xqDvxeBYGfz9/fHBBx9gyZIliIuLw/Tp0xETE4P4+Hh89dVXet+gRBdWVlbYtm1buRdje/xxQkIC3N3dFcuRmY/Lly+jVq1aWL9+fZkXnnu8TWk5OTn49ddfcffuXTRu3BhXrlwBAHz66afw8fFRPB9j4QJgQPPnz8eBAwcQExODFStWSO2RkZEICAhQLI9atWqhVq1aisVjNZuHhwc8PDyMnUa5QkNDkZ2djd27d0ttffr0wccff2zErJRnNlNAv//+O/Lz8zXak5KS8OeffyqeT2pqKo4fP6718ro//fST4vkwxv61detWjTZfX18cO3YMaWlpiudz9OhRJCcna7Tn5uZi165dBotrNgVgx44daN++PSIiIqS24OBguLu7Q61WK57P+fPnMXLkSCxbtkzjud9++w05OTmK58QYK7x/dFhYmEb7/PnzMXr0aGk6SEkPHjxAmzZtMH36dKktJCQEbdu2xZ49ewwWV1Tng31CiFQASQbo2gWAKd+cxtTzB0z/NXD+xmfqr8FQ+bckogaVWbFaFwBDEUJEE1GZt7Cs7kw9f8D0XwPnb3ym/hqqQ/5mMwXEGGNMN1wAGGOshqqpBWCdsRPQk6nnD5j+a+D8jc/UX4PR86+RxwAYY4zV3D0Axhir8WpkARBCzBNCnBVCxAoh9gohmho7J10JIRYLIS4XvY7fhBD1jJ2TLoQQLwghLggh1EIIkzmTQwgxWAgRL4RIEELMMHY+uhJCfCuESBFClH9hqGpKCNFcCHFACHGx6P0z2dg56UoIYSeEOCmEiCt6DcpfIKw4l5o4BSSEqENEmUU/TwLgQUQTjZyWToQQTwPYT0QFQoiFAEBE0yv4b9WGEKIjADWAbwB8QETRRk6pQkIISwBXAPih8FaopwCMIqKLRk1MB0KIPgAeAPieiDobOx9dCSGaAGhCRKeFELUBxAAYbmJ/AwHAkYgeCCGsARwBMJmITiidS43cAyge/Is4ovDyqyaFiPYSUUHRwxMAXI2Zj66I6BIRxRs7Dx31BJBARFeJKA/ANgDDjJyTTojoEIB0Y+dRVUR0i4hOF/2cBeASgGbGzUo3RVdtflD00LpoMcoYVCMLAAAIIeYLIf4BMAbAJ8bOR0+vA9hd4VpMX80A/FPi8Q2Y2OBjToQQrQB0A6D8xb70JISwFELEAkgBEEFERnkNZlsAhBCRQojzWpZhAEBEM4moOYAtAIKNm612Fb2GonVmAihA4euoViqTP2NVIYSoBeAXAFMe26M3CUSkIqKuKNxz7ymEMMp0nNleDpqIBlVy1S0AwgDMMWA6VUGVfboAAAFESURBVFLRaxBCjAMQCGAgVcODOTr8DUzFTQDNSzx2LWpjCiqaN/8FwBYi+tXY+eiDiO4LIQ4AGAxA8QPzZrsHUB4hRMk7nQwDcNlYuVSVEGIwgGkAniEivrSoMk4BcBdCtBZC2AB4CYDhrtXLNBQdQN0I4BIRaV5q1wQIIRoUn7UnhLBH4UkFRhmDaupZQL8AaI/Cs1CSAEwkIpPakhNCJACwBVB88fITpnQmkxDiWQArATQAcB9ALBH5GzerigkhhgD4PwCWAL4lovlGTkknQoifAPRD4ZUo7wCYQ0QbjZqUDoQQ/wFwGMA5FH5+AeBjItK8vnM1JYToAuA7FL6HLAD8j4g+M0ouNbEAMMYYq6FTQIwxxrgAMMZYjcUFgDHGaiguAIwxVkNxAWCMsRqKCwBjjNVQXAAYY6yG4gLAGGM11P8Dr5yt8tV4VCsAAAAASUVORK5CYII=\n",
- "text/plain": [
- "<Figure size 432x288 with 1 Axes>"
- ]
- },
- "metadata": {
- "needs_background": "light"
- },
- "output_type": "display_data"
- }
- ],
+ "outputs": [],
"source": [
"plt.quiver(vdir[..., __spherical_indeces__['azimuth']],\n",
" vdir[..., __spherical_indeces__['elevation']],\n",
@@ -374,7 +367,7 @@
},
{
"cell_type": "code",
- "execution_count": 179,
+ "execution_count": null,
"metadata": {},
"outputs": [],
"source": [
@@ -385,9 +378,9 @@
" names=['position',\n",
" 'orientation'])\n",
"velocity = pd.Series(index=index_posvel, data=0)\n",
- "#velocity.loc[(convention,'alpha_0')] = 2*np.pi*(np.random.rand()-0.5)\n",
- "#velocity.loc[(convention,'alpha_1')] = np.pi*(np.random.rand()-0.5)\n",
- "#velocity.loc[(convention,'alpha_2')] = 2*np.pi*(np.random.rand()-0.5)\n",
+ "velocity.loc[(convention,'alpha_0')] = 2*np.pi*(np.random.rand()-0.5)\n",
+ "velocity.loc[(convention,'alpha_1')] = np.pi*(np.random.rand()-0.5)\n",
+ "velocity.loc[(convention,'alpha_2')] = 2*np.pi*(np.random.rand()-0.5)\n",
"velocity.loc[(convention,'dalpha_2')]= np.random.rand()\n",
"\n",
"elevation = np.linspace(-np.pi/2,np.pi/2,5)\n",
@@ -412,32 +405,9 @@
},
{
"cell_type": "code",
- "execution_count": 180,
+ "execution_count": null,
"metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- "<matplotlib.quiver.Quiver at 0x7f00b173e358>"
- ]
- },
- "execution_count": 180,
- "metadata": {},
- "output_type": "execute_result"
- },
- {
- "data": {
- "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADx0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wcmMyLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvMCCy2AAAIABJREFUeJzt3XdYFFf3B/DvpQkKGIoVFYLY0MSo2FuMnYi+thh7jWmWWGJNYolGzc9INJq8lhjbK4gxalCMMRoNUQmCYqKIimLFgijSpCx7fn8AExYWBXZ2BtjzeZ55sjs7uefssN4zc6cJIgJjjDHTY6Z2AowxxtTBBYAxxkwUFwDGGDNRXAAYY8xEcQFgjDETxQWAMcZMFBcAxhgzUVwAGGPMRHEBYIwxE2WhdgLP4+zsTG5uborGTEtLw6VLl6DVanXmW1tbw9HREY6OjqhQoYKiOeV14cIFZGRkwNzcHBqNBjY2NvD09IRGo4EQAubm5ormo9FoYGZmBq1Wi/PnzwMArKyskJGRAWdnZ7i6uiqaT14pKSl4/Pgxnjx5gszMTJ3PrK2t0ahRI5iZ8TYQK1/Cw8MfEVGVIi1MRKV2atGiBalh165dBKDQqW3btrRnzx5VcmvRogUdOHCAFixYIOVz5coVVXLJ6/vvv5fy+e2332jBggX0/vvvK56HRqOh5cuXk4eHR6F/P1tbW4qKilI8N8aUACCMitjHqt7JP29SqwAQEU2aNElv5+Hs7ExfffUVpaamqpLXpUuXiIjo77//lnJatmyZKrnk5e3tTQDIycmJMjMziejfXJV27949mjJlCllZWen9G/r7+6uSF2NKKE4B4P3fQqxcuRItW7YsMN/W1hbt2rWDjY2NClkBDRs2BAA0adIE9evXBwDs2bNHlVxyPX36FEeOHAEA9OvXDxYW2SOLubkqrVq1amjevDksLS0LfPbhhx9iyJAhKmTFWOnDBaAQFSpUwO7du+Hg4KAz/8aNG+jQoQOWLVtW4DiBkoQQGDhwIAAgLCwMN2/eVC2XAwcOSGPsuTmpJTExESNGjMCYMWOQkpKi81nLli3x1VdfqZQZY6UPF4DncHV1xfbt2wEAL7/8MlauXAlLS0tkZWVh3rx56N69O2JjY1XLL29nq+ZeQG5se3t7dO3aVbU8/vrrLzRr1gw7d+4EAFSpUgXff/89AMDBwQEBAQGqHsBnrNQp6liRGlPeYwAajYaCg4Pp8uXL8g2WFdGcOXPIx8eHiIjOnDmjc4DRycmJAgMDFc+JiEir1ZKbmxsBoHbt2qmSQ3JyMtnY2BAAGjFihCo5ZGVl0fLly8nCwkL6u3Tr1o1iY2NJq9Wq9jf6888/6eLFi6TVahWPzUqfmJgYOnLkiNHjoLwcBH7ttdcoICCARo4cSU5OTuTs7ExJSUmyr7AXyczMpAMHDkjvExMTadSoUToHFqdMmULPnj1TPLfp06dLOdy9e1fx+Lt375bi7927V/H4sbGx1K1bNykHCwsLWrFiBWVlZUnL7Nu3T/G8iIjCwsIIALm7u9PUqVPpyJEjlJ6erkouTHkajYZOnjxJc+bMocaNGxMA+vHHH40et9wUACGETif7xRdfkEajKdKkhO3bt5Otra2UX9OmTSkyMlKR2LlOnTolxV+7dq2isYmI3n77bQJAlSpVUvzMqIMHD5Kzs7P0/d3d3emvv/4yelytVlvk3+Gbb76p8xu2t7enwYMH07Zt2yguLs7ouTJlPX36lHbv3k2jRo3S+W0CoMaNG1NGRkaRfjeG7DWWmwKQd+UVZ/Ly8irxyiuuq1evkpeXlxS7Z8+eisUmyh7+qFmzJgEgb29vRWNnZmaSo6MjAaDBgwcrGjs5OZmqVasmrfdhw4bR06dPFYmdnp5eot9l/snMzIzat29Py5cvpwsXLhTpH/2NGzdo4cKFlJGRocA31ZWWlqbYxlWulJQUReMRZRf4NWvWUGhoaJGWv379Oq1Zs4a6d+9OlpaWsvw2DDmFmguAggWAKLtDmDVrFjk5OdGtW7cUjU1EtGPHDgoMDFRlCCo+Pp5++OEH+uOPPxSPfejQIbKzs6OtW7cqOs4uVwHInV566SUaOnQonT59+rlxMzMzqX379gSA2rRpQ9evX1fk+16/fp2GDRtGdnZ2dOzYMUVi5sa1tramfv360dmzZxWJGRcXRz4+PgSA6tatS4mJic9d/tKlSzR+/HidjZGyVABE9vKGEUJsBtAHwEMiaqLn89cB7AcQkzPrJyJaXIR2dZIbOHAgmjQp0HwBNWvWxMSJE4uQubweP34MR0dHxeOaMjXWeVZWFj7//PMiLRsUFIQzZ84UmF+/fn34+PjAx8cH7du3l66deJ7Y2Fj06NEDFy9eBJB91tWGDRuMfl1DXFwcqlevDq1Wiw8//BBr1641arxcX331FWbOnAkACA8PR/PmzY0a79ixYxg5cqR0Zl/t2rURFBRUpD5Hq9UiPDwcgYGBCAwMRERERIFlPDw8MHz48CLlMmnSJDg7OxfvC+QQQoQTkVeRFi5qpXjeBKATgOYALhTy+esADhS33erVq9Orr74qVcVevXqVuCoyprTk5GSqUqUKASBzc3N6/fXXaeXKlQadyZaamkrvv/++ztbi+PHjKTk5WWe5F225FleXLl0IANWoUUPnALsxtW3blgCQm5ubrHt4iYmJOu1lZGTQvHnzKO8xx/79+1N8fHyJY9y6dYu+/fZb6t27N1WoUIEAkJWVFd2+fVuOr/BcUGMICICb3AUg9zTQGzdu0Nq1a6lnz5507tw5+dcYY0awadMmGjp0KO3cuZMeP34sa9s//fQTOTg4SB1WgwYNpH8bGo2G3njjDVmPE6xdu1aKdfLkSdnaLcydO3ekeDNmzJCtXa1WS4MGDZLOmLt+/Tq1adNGimVtbU3fffedrAUnOTmZ9u3bR+PHj6fPP/9ctnYLU1oLQDyA8wAOAWhclDb13QuIz6lmZYWxf6u3bt2ijh07Sp2XlZUVrV69mi5fvkwAaNq0abLFunv3rhRn+vTpsrVbmG+++UaK96JjI8WxevVqAkC//vor+fv7k729vRSncePG9M8//8gWSx8l+q/SWADsAdjmvPYGcPU57UwEEAYgrE6dOsZbS4yVAxqNhhYtWkRmZmZSR1avXj3ptZx3rW3Xrh0BIFdXV6N3ZJ07dyYA5OLiItuQU0hIiHSWTt51BIDef/991W7wKLdSVwD0LHsDgPOLllPzbqCMlSXBwcFUu3btAmeT2Nvb09WrV2WJsWrVKqndsLAwWdrU58GDB1JBmzx5sixtPnr0SO/6cXBwoJ9++kmWGKVFcQqAIvcCEkJUF0KInNetkH0PonglYjNmCjp06AB/f/8C9zpKTEzE4MGD8ezZM4NjDBgwQHqd995T2X2OYfK2sW/fPulGi3LcXFCr1WLkyJG4fft2gc/Wr1+P/v37GxyjrJKlAAgh/ACcBtBACHFHCDFeCPGeEOK9nEUGAbgghDgPYA2At0mOXw1jDGlpaZg4cSI6duyI9PT0Ap9HRERg6tSpBsdxdXWFl1f22YU//vgjiAgHDx5EUFCQwW1HRERg48aNAP4tLlWrVkWHDh0Mbnv58uU4dOiQ3s/eeust/Oc//8Hjx48NjlMmFXVXQY2Jh4AYK5rHjx/Tpk2bqGvXrpT/Fiq507Zt2wyOs2zZMqm9ZcuWkaWlJfn5+Rnc7unTp6U2c2/qN3HiRIPbPXbsmM7xkbxT69atafXq1XTv3j2D45QmUPpCMGPx8vKisLAwtdNgrEyJjY1FQEAA/Pz8EBoaKs2vWLEiQkND0bhx42K3uWnTJkRHR8PS0hJLlizR+Wzbtm0YOXKkQTkHBwejU6dOOvM++OADODo6wtbWFrNnzy52m/fu3UOzZs3w4MEDaV6jRo0wfPhwvP3226hbt65BOZdWil8IZqyJ9wAYM8zVq1dp8eLF1LBhQ+l6gZJcJHb//v0CNzfLnTZt2mRwnr/99pvetitUqEAXL14sdnuZmZnSmUR16tSh2bNnU0REhEmcRo7SdhCYMaYODw8PfPrpp4iMjMS5c+fQt29fLFmypNgHbqtVq4YNGzbo/SwjI8PgPAtrY9myZfD09Cx2e9988w0aN26M4OBgxMTEYPny5WjatClyzkVhOV58AxLGWJknhMBrr72G1157DVqtNnv8t5idYf/+/TFmzBhs2bJFZ76xCkCXLl1KfPB6ypQpMDc3NzStco/3ABgzMWZmZjAzK9k//dWrV8PV1VVnnjEKgL29PbZs2VLiPLnzLxouAIyxIrO3t8fWrVt19h4yMzMNbjd/G2vXrkWdOnUMbpc9HxcAxlixdO7cGdOnT5fey70HMHDgQIwYMcLgNtmLcQFgjBXbkiVLpPvky1kAqlevjv/+9798sFYhXAAYY8VmbW2NHTt2wNLSUtYCsHnz5hI/CIUVHxcAxliJNG3aFIsXL5atALz77rvo3bu3DJmxouICwBgrsY8//hitW7c2uB1XV1esXLlShoxYcfCtIBhjBinJNQXGaINlK86tIHgPgDFmEDk6bu781cEFgDHGTBQXAMYYM1FcABhjzERxAWCMMRPFBYAxxkwUF4AyKPeB2aVVac+PMZaNC0AZdPnyZfz+++9qp1Go9evXq50CY6wIuACUQU+ePMG0adOQlZWldioFxMbGYs6cObwXwFgZwAWgDEpISMD58+cLPJmpNJg/fz4SExORnJysdiqMsReQpQAIITYLIR4KIS4U8rkQQqwRQkQLIf4WQjSXI66pevLkCYB/O9vSIjw8XCpKuTkyxkovufYAtgDo9ZzPewOolzNNBPCdTHFNUm7n+uDBAyxbtkzlbLIRET766CPpPRcAxko/WQoAEf0B4PFzFukHYBtlCwHwkhCihhyxTVHeztXX1xcxMTEqZpPtxx9/xJ9//im95wLAWOmn1DEAFwC387y/kzOPlUDezjU9PR2zZ89WMRsgLS0Ns2bN0pnHBYCx0q/UHQQWQkwUQoQJIcLi4uJUzSUjIwOl8XbZCQkJOu93796N4OBglbIBvv76a9y4cUNnXv4cS4PS+vdkTC1KFYC7AGrneV8rZ14BRLSBiLyIyKtKlSqKJFeYrKwsDB8+HKGhoarmkZ++retp06apcurl/fv3sXTp0gLzS9MeABEhICAAc+fO5dsOM5aHUgXgZwCjcs4GagPgKRHdUyh2idnY2KBTp05o3bo1hg0bVmArVy36Otfw8HBs375d8Vw++eQTvad8lpYCcOrUKbRr1w5vv/02xowZo3Y6jJUqcp0G6gfgNIAGQog7QojxQoj3hBDv5SwSBOA6gGgAGwF8IEdcJYwbNw7u7u7w8/NDw4YNMXv2bNWHN/J3ri4uLmjRogU+/fRTRc+/j4iIgL+/P9544w1YWVk9N0elXbt2DYMHD0b79u0REhKCoUOH4pVXXlE1J8ZKHSIqtVOLFi2oNNi2bRsBkCYnJyf65ptvKCMjQ5V8ateuTb1796a2bdsSAKpfvz4RET19+pQSEhIUy+P+/fuUnp5ODx48kNbNyJEjqV69ejR06FDF8sgrPj6epk2bRpaWllJO5ubmdPXqVVXyYUxpAMKoiH1sqTsIXBoNGzYMnp6e0vv4+HhMnjwZjRs3xr59+xQ/sDhhwgQEBgaie/fuAIArV64gLi4O9vb2qFy5smJ5VKtWDVZWVjh16pQ0b9SoUfjrr7/QsmVLxfIAss+G8vX1hYeHB3x9fZGZmSl9NnbsWHh4eCiaD2NlQlErhRpTadkDICLas2ePzl5A3qlTp04UGhqqeE6HDx+Wcti3b5/i8XPNnDmTAJCZmRklJiYqGlur1dLu3bvJ3d1d79/GysqKbt26pWhOjKkJvAcgv/79+6NFixYF5ltZWcHW1hb//PMPNBqNojm1adMGZmbZf8KTJ08qGjuv3NhNmzaFnZ2dorGfPHmC69evo7Azxt5//33Url1b72eMmTouAEUkhNB7umPt2rWxceNGjBs3DhYWFormZG9vLx3YVKsApKWlITw8HADQvn17xeM7Ojpi1qxZ+L//+z9UrFhR57OKFSti7ty5iufEWFlRbgpA3jFfY+nRowc6duwIAKhQoQKA7LNNunbtigcPHhg9vj65nW5YWBjS0tIUjx8WFoaMjAydXJQWEhICb29vpKamQgghFeKPPvoI1apVM2psrVZbKm/LzcoWIlJ8BAEo4wXg3r172LhxI/r27YtPPvnE6PHy7gUsXboU77zzDgAgKioKXbt2hRpXLud2uhkZGdKWuJLy7nmoUQDCwsLQq1cv6fTXzZs3w8fHB5UrV8bMmTONHl8IAW9vb4wZMwZ79uxBUlKS0WOy8iEtLQ1BQUH44IMP0KZNG0U2Ygso6sECNab8B4G1Wi2Fh4fTwoULycvLS+dA3507d2Q5gFIU3t7edOfOHcrKyqLRo0dLeTRt2pTi4+MVy4OIKCYmRor/5ZdfKhqbiKhv374EgGrVqqV47HPnzpGDg4P0/devX09ERAEBAbRkyRLF8ti7d6+Ug6WlJXXv3p3WrFlDMTExiuXAyoZ79+7Rpk2bqF+/flSxYkXpd+Pr6ytbDBTjILAghU9hLA4vLy8KDg7GsWPHEBgYiAMHDuDu3YJ3kHBxcUGfPn2k966urkYd+42Li5MOOmZlZWHkyJHw8/MDALRs2RInT56EpaWl0eLnRUSoVasWYmNj0a9fP+zbt0+RuLmxq1atikePHmHIkCHw9/dXLPbNmzfh5eWFR48eAQC++eYbTJo0CQDw7NkzZGVlwdbW1iixNRqNFAvIXg+bN2/WuwvfpEkT+Pj4oE+fPmjdujXMzc0Njp+cnGy078ayJScno1KlSgbfOoSIcP78eQQGBiIwMBBnzpzRu9zIkSN1jmEtXLgQ1atXL1FMIUQ4EXkVOcHSOjk4OJCNjU2hp18WNnl5eclWTYsiMzOTBg4cqLMVqqR3332X+vTpQ+vWrVM0bmpqKs2dO5c6d+5MGzZsUDS2RqOhMWPGEABauXKlorHT09OL/ZsEQFWqVKHp06cbdLHe6tWrqU6dOoruXaSkpNBvv/1GixYtovv37ysWNyMjg+bPn08HDhygJ0+eKBY3Pj6emjZtStOmTSOtVluiNp49e0aLFy+m2rVrl+i3cunSpRLnj2LsAajeyT9vatasGe3du5fGjRtHVatWLXRlWVlZkbOzszR169atxCuvpNLT0+nw4cOKxzVlWVlZFBQUpHjc9PR0nd+bs7MzCSEK/X02aNCAZs6cSSdOnKDMzMwSxz1x4oTUppubm871DVlZWXJ8Nb1+//13Ke7u3buNFie/M2fOSHE3b95stDh5111CQoLO8PKWLVsMavfMmTP02WefUbNmzZ7b4Ts5Oen8nq5cuVLiuOWmAOQ9BpCVlUUhISE0f/58evXVV3VWnouLC6WlpZV4hTFmiMjISJ0CYG5uTl26dKFVq1YZ9A85P61WS9OmTZPieHh40N27d4mIyM/Pj86ePStbrLxSUlLIwsKCANBHH31klBj6rF69Wvquly9fNkqM2NhY+vrrr4mIKDExkdq0aSPFHDhwoEEFO7/bt2/Td999R97e3lShQgWdPuyXX36RLU65LAD53bhxg9atW0e9evUiKysrxYc/GMs1dOhQcnBwoGHDhpGfn59Rhyu0Wi19+OGHUsfRsGFDun//Pvn4+NCbb75ptLgtW7YkANSqVSujxcjvrbfeIgDk7Oxc4qGYF5k0aRI1adKEkpOTqWPHjtJ67du3r1Hv9ZWcnEz79++nCRMmUPXq1alNmzayfUeTKAB5JSUlUUhISPHWEmMy0Gg0FBwcLOuW4otkZWXRhAkTpM7K09NTuvndyZMnjRJz6tSpBIAsLCwoJSXFKDHy0mq15OLiInXGxnDjxg1pveUdoundu7eiIwpZWVkUGhoq221UilMAyvR1ALlsbW3RunVrtdNgJsjc3BwdOnRQ9CpwMzMzrF+/HqNGjQIAREZGSueQz58/P3vLTma513hoNJpCz2SR061bt6Qz/ox1fcmiRYuk9Xbu3DkAQLdu3bBnzx7pQk8lmJmZoWXLlorfRgUo4xeCMWaKiAj//PMP6tSpI90LKtfx48dx9OhR2WPm7YSVuO1I3jvMGqMAXL58GVu3bi0wv379+oiIiFDl6Xpq4ALAWBlz8+ZNzJo1C0uWLNHbUc2bN0/2vYCaNWvCzc0NgDIFIDeGlZWV3pswGuqzzz7Tu+6+/fZbTJ06FX///bfsMUsjLgCMlTFubm44fPgwfvnlFzRp0qTA52fOnMHPP/8se9zcLfFTp05Bq9UiLS1NGjqRQ1RUlPS0vdwC4OXlBWtra9liANlPsgsICCgw39XVFTt37kRISAhee+01WWOWVlwAGCujevbsiYiICGzcuLHAVaOffPKJ7Depa9euHQAgISEBx48fR5cuXXD69GnZ2o+JiUHr1q0RHh4ubYHnxpRT/vuG2dvbY8WKFYiKisLQoUMLDKuVZ6bzTRkrh8zNzTFhwgRcvXoVCxYskG4ncOHCBezatUuWGCtWrMDEiRMRGxsrzevZsydCQkLg4OAgSwwAcHBwwJUrV9CmTRtpeCYlJQVTp07F9OnTZYlx+vRpHDx4EABgYWGBSZMmITo6GrNmzZJ9T6NMKOrpQmpMpemJYIyVBXfu3KFx48aREII8PDxkOZc9NjaWKlWqpPcKVjmvxL506ZLeGEIIOnfunMHta7Vaev311wkA9evXj6KiomTIuvSBqZ0GyhjL5uLigu+//x7nzp2Dm5sbtmzZYnCbNWrUwLx58/R+JvcegD7jx4+XZUz+6NGjSEpKwvHjx7Fv3z40aNDA4DbLulJ/N9CwsDC102CsTCIiXLx4Ue+B4uJ69uwZGjVqhJs3b+rMj4qKkq0jzcjIKHD+vZ2dHa5evSrLg30uXryIRo0alfsx/uLcDbR8rwnGTJgQQpbOHwBsbGywYsWKAvPl3AOwsrIq8FjPefPmyfZUt8aNG5f7zr+4ZFkbQoheQojLQohoIcQcPZ+PEULECSEicqYJcsRljCnnrbfeKnBWjpwFIH97bm5u+Oijj2Rtn+kyuAAIIcwBrAPQG4AngKFCCE89i+4iotdypk2GxmWMKUsIga+//lp6X6lSJdkffJS3APzf//2faZ6ZoyA59gBaAYgmoutElAHAH0A/GdpljJUyLVu2lO5B9NJLL8nefm6bHTt2xMCBA2Vvn+mSowC4ALid5/2dnHn5DRRC/C2E+FEIUVuGuIwxFXzxxReoWLGi7MM/QPYegBACvr6+Bj+Okb2YUkdEAgG4EdGrAI4AKHgXphxCiIlCiDAhRFhcXJxC6THGisrFxQWzZ882WgEYPXq0Ue7/wwqSowDcBZB3i75WzjwJEcUTUXrO200ACv3rEtEGIvIiIq/cB68zxkqXmTNnomnTprK3W7t2bSxdulT2dpl+chSAMwDqCSFeFkJYAXgbgM6dqIQQNfK87QvgkgxxGWMqqVixIhYtWiR7u9OnT0fNmjVlb5fpZ/BTLIhII4SYBOAwAHMAm4noohBiMbIvSf4ZwBQhRF8AGgCPAYwxNC5jTF2Ojo5lok1WOL4SmDHGyhG+EpgxxtgLcQFgjDETxQWAMcZMFBcAxhgzUVwAGGPMRHEBKCcuXVLn0oqMjAxcu3ZNldiMMcNwASgn5s2bh8TERMXjnjhxArt371Y8LmPMcFwAyoGUlBQcOnQIhw8fVjx2YGAgAgMDFY/LGDMcF4By4OjRo0hPT1e8IyYiHDhwAKdPn8ajR48Ujc0YMxwXgHIgt+MPCgpCVlaWYnEjIyMRExMDIkJQUJBicRlj8uACUMZptVocOHAAABAfH4+QkBDFYufd4+BhIMbKHi4AZVx4eDju378vvVeyI84b6/Dhw8jIyFAsNmPMcFwAyrjcrf9cShWAR48e4fTp09L7pKQk/PHHH4rEZozJgwuAjIgI58+fVzRm/g4/MjIS169fN3rcoKAg5L+TrNLDQEqva8bKGy4AMrp8+TLWrFmjWLw7d+7g3LlzBeYr0RHrixEYGFigKBjT5MmTFT3ozVh5wwVARoGBgTh48CC0Wq0i8fIP/+TNw5gyMjL0XnMQExODyMhIo8bOdeXKFQQHByM0NFSReIyVR1wAZBQYGIgHDx7gzJkzisWrVKkSXn75ZQBAhQoVUK9ePZw4cQJPnz41WtwTJ04gKSkJr7zyijTP09MTQgjFhoFy4/DZR4yVHBcAmTx+/BgnT54EUPiWuZw0Gg1at26NGzduoGfPngAAW1tbREZGYvPmzbh8+bLRYj9+/Bi///67zi0g5s6di4sXL6JixYpGi5tX7jpWYl0zVl4Z/Exglu3QoUPS0E9gYCA+//xzo8azsLDAZ599pnf+yJEjjRp7yJAhAFCgyDRq1AiNGjUyamwAePLkCYKDgwEA//zzD27evAlXV1ejx2WsvOE9AJnkHYo4f/48bt26pWI25dsvv/yic/CXh4EYKxkuADLIzMzEL7/8ojOPhyaMJ3+HzwWAsZLhAiCDP//8s8BBVy4AxqHRaHDo0CGdecePH0dSUpJKGTFWdslSAIQQvYQQl4UQ0UKIOXo+ryCE2JXz+V9CCDc54pYW+rZAjx07hpSUFBWyKd9OnjyJhIQEnXkZGRk4cuSIShkxVnYZXACEEOYA1gHoDcATwFAhhGe+xcYDeEJEHgB8AawwNG5pQUR6C0B6ejp3SkZQ2HAPDwMxVnxy7AG0AhBNRNeJKAOAP4B++ZbpB2BrzusfAXQVQggZYuul5HDA5cuXER0djbZt20rzmjRpgkqVKnGnZASBgYFwcXFBrVq1AABVqlSBu7s7Dh48qNhVwUSk6BXPzLQQEZKTkxWJJUcBcAFwO8/7Oznz9C5DRBoATwE4yRC7gPPnz8PV1RU+Pj7GaL6A2NhYHD16FMeOHZPmDRkyBDdu3NC5UIoZLj4+HlOnTsW1a9fQuHFjAICrqyuioqLw5Zdf4vbt2y9owXDJyclwcXGBh4cH0tPTjR6PmZ4uXbrA3d0dV65cMX6w3K2Zkk4ABgHYlOf9SABr8y1zAUCtPO+vAXAupL2JAMIAhNWpU4eK688//yRzc3Nq3rx5sf9NSarQAAAgAElEQVRfQ2i1WoqLi6O4uDhKSUlRNHZSUhLFxcXRo0ePFI2bmZkpfednz54pGjshIYHi4uLoyZMnisZ98OAB2dnZkbOzMyUkJBg93rVr1yg8PJyIiFq3bk1OTk40duxY0mq1dOjQIUpKSjJ6Dj/88AM5OTmRk5MT3bt3z+jx8qpRowY5OTnRggULjB5Lq9XS3r17SaPR0MKFC8nJyYmqV69ORERXr16lc+fOGT0HIiIPDw+ytLQscTwAYVTE/luOC8HuAqid532tnHn6lrkjhLAAUBlAfCEFaQOADQDg5eVV7P3s9u3b4+LFi6hdu/aLF5aREALOzs6Kxsxla2sLW1tbxeNaWFio9p0rV66sStyqVaviwoULMDc3VySHKlWqoFGjRujVqxdiYmIQHx+PsLAwtG3bFgkJCYiKijJ6DmlpaYiPz/7nqtR9rnLFx8cjIyMDqampRo8lhMD69esxd+5cODs7Iz4+HmZmZhgzZgx27NiBixcvGj0HADh79iwePHgADw8Po8eSowCcAVBPCPEysjv6twEMy7fMzwBGAziN7D2GYzmVyigaNGhgrKYZQ506dRSLZWdnh9dffx0///yzNO+ff/4BAMycOVOxPEyFj48PPvzwQ+m9VqvF1q1bUa9ePcX6FTs7O9jZ2SkSy+BjAJQ9pj8JwGEAlwAEENFFIcRiIUTfnMW+B+AkhIgGMB1AgVNFGWP69enTR+98pY5zmZI333xT7/zyuq5luQ6AiIKIqD4R1SWipTnzPiOin3NepxHRYCLyIKJWRGT8J5YwVk7o63wcHBzQrl07FbIp31xdXfHqq68WmM8FgDGmCjc3NzRp0kRnXu/evWFhwfdyNIb8nX3lypXRvn17lbIxLi4AjJUB+Tul8rpFWhrkX7e9e/eGpaWlStkYFxcAxsqAvJ2ShYUFevXqpWI25VvLli1RtWpV6X15LrZcABgrA1q1aiWdctuxY0e89NJLRo/p7++v94rUa9eu4ffffzda3PDwcERERBSYn5GRge3btxv9KmwzMzPpYLC5uTl69+5t1Hhq4gLAWBlgbm4udUpKbZGeOnUKL7/8MgICAqR5kyZNQoMGDYx6PYCDgwOaNWuGAQMGSLf3OHHiBOrVq4eAgAAY8S4yktx13KFDBzg4OBg9nlq4ADBWRuR2SkoVAB8fHzx69Ehna3/v3r2oVKkSOnbsaLS47u7u8PT0xN69e6UCEBoailu3bin23bt37w4rK6tyPfwDcAFgrMzo0aMHmjZtqsgVogDQuXNnvVeY9+rVC1ZWVkaNXVjHW9h5+nKztbXFG2+8wQWAMVY62NnZYeXKlYrFs7KyQs+ePQvML+zCNDnpi9G8eXO4uOS/z6TxzJgxA/Xr11csnhq4ADBWhnTr1k3RePm3gM3MzODt7W30uG3btoWTk+4Ng5XeGld6XauBCwBjrFDe3t46B13btWtXoGM2BnNz8wKFprwPx6iBCwBjrFBVqlTRediRkp1w3lg1a9ZE8+bNFYttKrgAMMaeK29HrGQB6Nmzp3S7iz59+ihy+qep4QLAGHuu3AOy7u7uaNiwoWJx7e3t0blzZ50cmLy4ADDGnqtx48Zwc3ODj4+P4lvhPj4+sLa2RteuXRWNayq4ADDGnksIAR8fH1UOwvr4+KBbt26oWLGi4rFNAd9PljH2QsOHD0ezZs0Uj+vu7o7Zs2crHtdUcAFgjL1Q69atVYvdoUMH1WKXdzwExBhjJooLAGOMmSguAIwxZqK4ADDGmIniAsAYYyaKCwBjjJkogwqAEMJRCHFECHE15796n50mhMgSQkTkTD8bEpMVXVpaWplokzGmDkP3AOYAOEpE9QAczXmvzzMiei1n6mtgTFYEcXFxWL9+veztLl261OgP5WaMKcPQAtAPwNac11sB/MfA9phMFixYgJiYGNnbPXToEPz8/GRvlzGmPEMLQDUiupfz+j6AaoUsZy2ECBNChAghuEgY2YULF7B+/Xo8efJE9rafPHmCOXPmIDU1Vfa2GWPKemEBEEL8JoS4oGfql3c5yh4XKGxswJWIvAAMA/C1EKLuc+JNzCkWYXFxccX5LgwAEWHGjBnQarVGKwC3b9/GV199JXvbjDFlvbAAEFE3ImqiZ9oP4IEQogYA5Pz3YSFt3M3573UAxwEUelcpItpARF5E5FWlSpUSfCXTdujQIfz6668AIHsB0Gq1SEhIAAAsX74csbGxsrbPGFOWoUNAPwMYnfN6NID9+RcQQjgIISrkvHYG0B5ApIFxmR6ZmZmYPn269F7uApCYmCgdAE5NTcW8efNkbZ8xpixDC8ByAN2FEFcBdMt5DyGElxBiU84yjQCECSHOA/gdwHIi4gJgBN999x0uX74svc/dWpdL/va2bt2KsLAwWWMwxpQjSvMpfV5eXlTeO5iHDx+iatWqBrfz+PFjeHh46Gz1V6xYESkpKQa3nevcuXMFHszdoUMH/PHHH7I8KUqudcGYKRNChOccc30hvhJYJdHR0Rg0aBDOnDkjS3uLFi0qMOSTmpqKjIwMWdoH9A8p/fnnn/jxxx9laX/dunWYNWuW7HsujDH9uAAo7PHjx5g2bRo8PT0RGxsLb29vg9uMiorCunXr9H4m53GAwtqaNWuWLFcIT506FevXr4eHhwe++eYbWYsXY6wgLgAKSU9Px6pVq1C3bl18/fXXyMzMxNKlS2UZOlm8eDF69OiBPn36FPjM2AVg3LhxqFmzJr799luD23d0dMTMmTMRHx+PKVOmoEmTJti7dy9fecyYkXABMDIiQkBAABo1aoQZM2ZIwxvdunVDly5dZGl/+/btCAoKgoeHBwDA3Nwc/v7+sLW1lXU4JSEhAR4eHjp7G97e3jh58iSmTJkiS4yPPvoIzs7OAICrV69iwIAB6Ny5s2xDZYyxf3EBMKJTp06hXbt2GDJkSIHbMixdulSWGEIImJubAwBOnjwJAGjWrBmGDBmCkJAQVKxYUZY4APDKK68gNDQUo0ePLhDTwkKex0vb2dlh7ty5OvOCg4PRqlUrDBs2DDdu3JAlDmMM2VuQpXVq0aIFlUXR0dE0ePDg3CujC0x9+/aVPWZKSgpZWFgQAJo6dars7efXokULAkCtW7eWve3U1FSqWbOm3nVXoUIFmjVrFj158kT2uIyVBwDCqIh9LO8BGIGtrS26deuGbt26FfhMCIHPP/9c9pihoaHQaDQAgPbt28vefn65Mc6ePYtnz57J2raNjQ0+/fTTAvMtLS0xaNAgtGrVCpaWlrLGZMwUcQEwgmrVqmH8+PHSWHZeb7/9Nl599VXZY546dUp63a5dO9nbzy83RmZmplHG58eNGwd3d3edeZmZmWjVqhUGDhyISpUqyR6TMVNT7gtAUlISnj59qmhMrVaLcePGwd/fHwBQvXp1ANkHZxctWmSUmLlj8a6urnBxcTFKjLzy7mXkxpaTlZUVFi5cCCB7y9/R0RHAv6eKKi02NpbPRmKquHv3rtHaLvcFwMrKCs2bN0e3bt2wevVqXL9+3ajxtFot3n33XWzbtg1A9oHTc+fOwdnZGWPHjkW9evWMEjN3D0CJ4R8AqFWrFurUqQNAd+9DTsOGDYOnpyd69uyJI0eOoHLlygCA9957Dz/88INRYubSarUICwvDggUL0Lx5c0yePFmWU3YZKy5fX180bNgQH3/8Mf744w9pqFcWRT1YoMYk10Hgb7/9VudAoqenJ82aNYuCg4NJo9HIEoOISKvV0gcffCDFadiwIT148ICIiKZNm0Y3b96ULVZeFy5ckGKuW7fOKDH0GTp0KAEgR0dHysrKMkqMPXv20M6dO4mIKCQkhOzs7AgACSFo+/btssZKSUmh/fv304QJE6hGjRo6v5mIiAhZYzFWVPfv3ycbGxvpt+jg4EDDhw8nf39/vSdDoBgHgVXv5J83lbQAREZGUseOHaWpVatWhZ6R4+TkRCNHjqSAgAB6+vRpieLlWrdundSuh4cH3b17V/rs2bNnBrX9POvXr1elo1q7dq0UNzIy0igxtFotpaWlSe+Dg4OpUqVKBIDMzMzo7NmzBrV/584d+u9//0tvvvkmWVtbF3rmUd7fU5cuXQz9WiVy5coVVeKaMrXW+fDhw3V+c5UrV9b727SwsKAuXbrQqlWr6OrVq0TEBYBCQ0ML7fCfN1WqVIm+/PJLSk9PL1Hc5ORk6tq1K7388st069atErVRElFRUbRixQoaPHiwrHs0L3Lx4kUaMWIEfffddxQXF6dY3GPHjpG1tTXNmjWLtFptidp49OgRvfvuuySEKPbvxNLSUuZv9GLbt28nMzMzWrt2reKxd+7cSR988AHt2LFD0d8XEdH+/ftp2rRp5Ovrq2hcIqK9e/eSpaUlLV26VPHYDRo0KFEfNmDAAC4Aly5dok6dOklT69atC11hzs7ONGrUKNq9e7fBewBE2eew37592+B22PNdvXq1xJ1/Xnfu3KH169dTnz59nrsHkPf39MYbb8jwDQr37Nkzne+2a9cuMjMzkzZSYmNjjRo/v4EDBxIAqlatmizrvDhy/+02b95c0bipqalUq1Yt6TewcuVKnc+NuUdPRDRixAid39zz9gC6du1Kvr6+FB0dTUS8B1BA3qEZANS4cWOaM2cOnTx5UvEtGlZ6paSkUGBgIE2cOLHAhWjnzp1TLI+vvvqKTp06RUTZW6Hm5uYEgKytrenYsWOK5UGUPQRXvXp1aetSadOnTycAZG5uTklJSYrG/ueff8jZ2Vn6DaxZs4aIsjcapkyZolge9+7d0zkG4OjoSCNGjKBdu3ZRQkJCgeW5AOSRlpZG7u7u1L17d1qzZg1dv37d4DZZ+afVaiksLIwWLlxILVq0UKzzS0xMJGdnZ5o8eTIdOHCALC0tCQBZWVnR4cOHFckhr2vXrkkdz1dffaV4/D179kjxf/vtN8XjR0REkIODg5TD+vXradWqVWRpaUkxMTGK5DB9+nRq1KiRdOJKZmbmc5fnApBHYmKiLEM7zLTFxsYqMvyxePFiAkCVK1emChUqSMccDhw4YPTY+mzbtk3q/EJCQhSPf//+fSn+okWLFI9PRHTmzBmyt7eX8qhSpQoBoLFjxyoSP+/JJEVRnAJQ7q8DsLOzg729vdppsDKuRo0aRr8O4PHjx1i5ciUA4OnTp0hPT4e5uTl27dqFN99806ixC5N7kZ+1tTWaNWumePxq1aqhbt26OrkozcvLC4cPH4atrS0AIC4uDkD2I1GjoqKMHr9mzZpGa7vcFwDGyoovv/wSiYmJOvOsra1x8OBBBAUFIT09XfGccjvdVq1awcrKSvH4wL8XN54+fRpZWVmKxs7KysLx48fh5+cHa2trnc+0Wi0WLFigaD5y4wLAWClw7949rFmzpsD8lJQU/Pjjj7h48WL2mK2CEhIScPHiRQDK3F+qMLmxk5KSpHyUdO3aNQQEBODRo0cFPgsICMC5c+cUz0kuXAAYKwW++OKLAndVtbCwwOTJkxEdHY2PP/64wBaosezatQtZWVkICQmRio5StxjRJ/99p4gIfn5+isQ2NzfH+PHjcfXqVSxYsEDv8zX03bm2rBBKb1UUh5eXF4WFhamdBmNGdfPmTdSrVw+ZmZnSvP79+2P58uWoX7++4vl4e3vDzMwMdevWlfZK4uPj4ejoCI1GI9vDf14kN5ZWq4WTkxMSEhIwePBg2Nvb49SpU4iMjFQkj7xiY2Px2WefYfPmzTp7ZCdPnlR1LykvIUQ4EXkVaeGiHi3WNwEYDOAiAC0Ar+cs1wvAZQDRAOYUtf2y+kAYxopj7Nix0hkmLVu2pBMnTqiaT+49nnIne3t7mjhxIjVq1EjRM4Fu3LhBnp6eNGrUKHJ1ddXJqW3btorloc/58+epe/fuUj6vv/664hfJFQZKnQYKoBGABgCOF1YAAJgDuAbAHYAVgPMAPIvSPhcAVt5FRUWRmZkZubq60s6dO412U73iyHtDw7zTsGHDFM9l/vz5enPx9vZWPBd9Dh06RI0bNyYAdOTIEbXTISIFTwMloktEdPkFi7UCEE1E14koA4A/gH6GxGWsvFi1ahWWLVuGqKgoDB06FGZm6h+Wc3BwKDDPxsYGy5cvVzyXOXPmSM/TyEtfjmro1asXIiIisHHjRqxdu1bxA/WGUuLX5gLgdp73d3LmMWbStFotvvjiC8yaNUuxA7xFoa9z/fjjj1G7dm3Fc7G1tcUXX3xRYH5pKQBA9sH6CRMmYMeOHUhLS1M7nWJ5YQEQQvwmhLigZzLKVrwQYqIQIkwIEZZ7wQVj5ZGZmRmcnJzUTqOA/J1rzZo1MWvWLJWyAUaPHl3gIrTSVABy2drawsbGRu00iuWFh/OJqOCTzYvnLoC8mw61cuYVFm8DgA1A9llABsZmjBVT/s512bJlqj6D2czMDF9//TU6d+4szSuNBaAsUmII6AyAekKIl4UQVgDeBvCzAnEZYyXw0ksvSa+9vLwwYsQIFbPJ1qlTJwwcOFB6nzdHVnIGFQAhRH8hxB0AbQEcFEIczplfUwgRBABEpAEwCcBhAJcABBCR8pfzMcaKJO/Wta+vb6k4MA1k3yoj93YUvAcgD0PPAtpLRLWIqAIRVSOinjnzY4nIO89yQURUn4jqEtFSQ5NmjBlPbuf61ltvoUOHDipn8y93d3d89NFHALgAyKV0lHbGWKnh4OCAChUqYMWKFWqnUsD8+fNRtWpVLgAy4QLAGNNhZ2eHWbNmwc3NTe1UCrC3t8eSJUu4AMiE7wXEGCsgLS2tVF2bkFdWVha0Wi0sLS3VTqVUKs69gJS5qxNjrEwprZ0/kH2HTnNzc7XTKBd4CIgxxkwUFwDGGDNRXAAYY8xEcQFgjDETxQWAMcZMFBcAxhgzUVwAyrn09HQEBwcb3M5vv/0mQzaMsdKEC0A59+mnnyIqKsrgdrZv3479+/fLkBFjrLTgAlCO/fHHH1i5cqV0B0VDWFlZ4Z133sHDhw9lyIwxVhpwASinEhMTMXr0aBCRbAUgLi4OEydOLHPPPWWM6ccFoJyaNm0abty4AQCyFQAA2L9/P7Zs2WJwe4wx9XEBKIf279+PzZs3S+/lLAAAMGXKFMTExBjcJmNMXVwAypmHDx/inXfe0Zknx10T87aRnJyM0aNHIysry+B2GWPq4QJQjhAR3nnnHcTFxenMl3sPAACCg4OxatUqg9tljKmHC0ApcvfuXaxatQrff/99if7/H374AT///HOB+cYoAADwySef4O+//y52WxqNBp999hm2b9+OpKQkg3NjjJUMFwCVPX78GBs3bkSXLl1Qu3Zt/O9//8Pw4cOL3U5MTAymTp2q9zNjFYCMjAyMHDkS6enpxWrLwsICY8aMweTJk1G1alUMGTIE+/btK3Y7jDHDcAFQQUpKCvz9/dG3b19Ur14dEydOxPHjx2Fvb4/du3cX+2EcRARfX194e3ujT58+BT43VgGYMGEC3NzcsG3btmK35+7ujq1btyItLQ0BAQHo378/qlWrhvHjx+Po0aN8fIExJRBRqZ1atGhB5UVGRgYdOHCAhg0bRpUqVSIABaZ9+/YZHGfixIkEgCwsLGj48OEEgC5cuGBwu//9739JCEGjRo2S8t2+fbvB7c6cOVPvuqhevTpNnTqV/vrrL9JqtQbHYcxUAAijIvaxqnfyz5vKSwG4f/8+tWrVSm9HlzvNnDnT4DgajYaqVKlCAKhHjx6UlZVF7777Ll25csXgtrds2UI//PADpaenU+XKlQkA9evXz+B2MzIyqEOHDs9dN2PHjqWMjAyDYzFmChQrAAAGA7gIQAvA6znL3QDwD4CI4iRXXgoAUXZH9+WXX+rt4Dp06CBLB/f7779LbW7YsIGIiLRaLaWnpxvcdlpamvQ6dy/A2tqakpKSDG77zp07UuHKO1WqVIkOHjzIewCMFUNx+lhDjwFcADAAwB9FWLYLEb1GRXxafXlz9uxZfPfddwXmV6lSBf7+/rKcq79nzx4AgJmZGf7zn/8AAIQQshwDqFChgvR64MCBAIC0tDQEBQUZ3LaLiwt27twJIYTO/JSUFKxcuRKxsbEGx2CMFWRQASCiS0R0Wa5kyiOtVovly5ejQ4cO0tWzlStXBpDdOe/cuRMuLi6yxPnpp58AAJ06dUKVKlUMbrMwPXr0gK2tLYB/i46hunXrhoULF0rvc9fR77//jqZNmyIwMFCWOIyxfyl1FhAB+FUIES6EmPi8BYUQE4UQYUKIsPwXNJU19+7dQ48ePTB37lxoNBpYWFjgyy+/xLp16wAACxYsQLdu3WSJ9ddff0lbyrlb6MZibW2NN998EwBw8OBBPHv2TJZ2P/nkE/To0QMAcOrUKYwePRoAEB8fj759+2LKlClIS0uTJRZjDC8+BgDgN2QP9eSf+uVZ5jiefwzAJee/VQGcB9CpKONTxj4GsGLFCvL19aXo6GjZ2z5w4AA5OztL49l169al0NBQIiKKiIig7t27k0ajkS3ejBkzpFh3796Vrd3C7N69W4q3d+9e2dp9+PAheXp6SuP+O3bsIFtbWynWq6++SpGRkbLFIyLKysqi0NBQ+vTTT2X9LowZ08mTJ2nGjBl0/PhxyszMlOZD6bOAXlQA8i27EMDMoixr7AJw+vRpqWNp2LAhffzxx/THH3/orMziSktLo6lTp+oczBw+fDg9ffpUZ5mHDx/K8RWIKPtAr5ubGwGgdu3aydbu8yQnJ5ONjQ0BoBEjRsja9u3bt3XeR0dHU8uWLaX1aWNjQxs3bjTo4HBycjLt27ePxo8fT9WrVycA5OjoqPN3Yqw002g01KBBAwJADg4ONGzYMPLz8ytWARBEht/bXQhxPKdTD9PzWSUAZkSUlPP6CIDFRPTLi9r18vKisLACTb5QfHx8kR9hOGfOHOm2ybkcHR3h7e0NHx8f9OzZUxqPfpHo6GgMHjwYERERAABbW1t8++23GDlyZLHyL67w8HB4eWUfW//qq68wffp0o8bLNWDAAOzduxeVK1fGgwcPdA4Uyy0jIwOffvopvvzyS2ne4MGDsWnTJtjb2xepjTt37uDAgQMIDAzEsWPHCgwn9erVC2PGjHlhO2ZmZhg8eHCx8pdDVlYWzM3NFY9rytRa50FBQUW6TUpQUJC+CzHDqagn2xS1UuibAPQHcAdAOoAHAA7nzK8JICjntTuyh33OI/uU0flFbb+kewChoaHPPa+8OJOFhQV17dqVfH196datW8+Ne+3aNbKzsyMA1KJFC1nOvy+KQ4cOkbu7OwGgmJgYRWISZQ/PAKBWrVrR9evXFYn566+/UrVq1QgANW3alJ49e/bc5c+dO0effvopNWvWTLbfhKWlpSLfNa+jR49SvXr1KCoqSvHYERERFBYWptrpuJGRkXTy5EnF4z58+JAaNmxI/v7+isfO3bIv4WTaF4LJWQBefvllmjJlCh05cqRI59Pv3LmTZsyYIcu598Wh1Wrp0qVLisZMTk6mmzdvKhqTiOjBgwc0YMCAIh0LSEhIoF27dtGIESPI0dGxzBWAjIwMmjt3LgkhCAC99tpriv+2Bg0aRADolVdeUbwI+Pj4EADy8vJSNK5Wq6U+ffpIf/Nx48ZRcnKyYvGVKgCyDAEZS0mHgFJTU3H9+vUXLkdE6Nevn87DTczMzNC2bVv4+PigT58+8PT0LHB+OiubNBoNTp8+LQ0DXbp0SefzAQMGYNGiRS9sRwiBxo0bGytNSUxMDIYOHYq//voLQPbZV19//TUmTpyo2G8yNTUVVapUQWpqKt5++234+fkpEjfXvHnzsGzZMgDAjRs34Orqqljsffv2Ydy4cXjy5AkAoH79+vD390ezZs2MHvvKlSvIyMh44XIbNmzAN998k3+2MkNAxp6MfRA4MDCQAJCdnR0NGjSItm7dSnFxcUaNyUqP6Oho8vX1pTfeeIMsLCzI2tqaYmNjVcklLCxM572fnx/Z29tLW3VNmjSR5Z5OxfXTTz9JOezevVvx+GFhYVL8VatWKR7/1q1b1KlTJykHKysr8vX11dkTyv+3U0pqaipVr16dhBDUpk0bWrp0Kf399988BFRUW7ZsKfLQDivfcoeKfv31V8Vjnzhxgjw9PYmIKCkpicaOHauzS//BBx9Qamqq4nkRkXRDQRsbG0WHQHLlPcOtffv2iscnyj7bZvHixWRmZib9Tby9venBgwdERNSzZ0/63//+p3heZ86coR9++EHKIxcXAMbKiPv371ONGjXIzMyMTp06RfXr15c6GQcHB1WvS0hLS5P2Qvr3769aHtOnT5fWiRLXuBTmzz//pDp16ki5VK9enY4cOUK1atWiSpUqyX59SkkVpwDw8wAYU0lWVhaGDRuGe/fuQavVon379rhy5QqA7Nt5nD9/XrqnkxqOHj2KxMREAMCgQYNUyyNv7L1796qWR/v27RERESHlc//+ffTo0QN37txBSkoKBg0ahJSUFNXyKwkuAIypZOHChTh27Jj0nohgZmaGxYsX49ixY6hdu7aK2f17nycrKyu9DxpSSuvWrVGzZk2dnNTi4OCAgIAAbNiwATY2NtnDKDkiIyPx3nvv6cwr7bgAMKaCX375BUuWLCkw38nJCTVr1lStE3n06BEAIDMzE/v27QMAdO/evcgX2xmDmZkZBgwYAAA4ceIEcu8RlpurGpydnVGjRo0C83fs2IGNGzeqkFHJcAFgTGG3b9/GiBEj9H4WFxeHKVOmYOnSpaoUgdmzZ2P//v04ceIEHj9+DMD4NxcsitwctFot9u3bh5CQEOlmgUrbtGkTxo8fX+ip5lOmTMHZs2cVzqqEinqwQI2JDwKz8iY9PZ3atGlT4OIdCwsL8vHxIT8/P1XOtsn11ltvkYWFBXl5eREAMjc3p0ePHhERyXrzwqLKjZn3aXeNGzcmOzs7atmypeL55EpPT6eff/6Zhg4dShUrVizw93R3d6cnT9AnwD0AAAOrSURBVJ6okhvKy4VgQog4ADeN0LQzAPX2Hw1X1vMHyv534PzVV9a/g7HydyWiIj0QpFQXAGMRQoRRGX4yWVnPHyj734HzV19Z/w6lIX8+BsAYYyaKCwBjjJkoUy0AG9ROwEBlPX+g7H8Hzl99Zf07qJ6/SR4DYIwxZrp7AIwxZvJMsgAIIT4XQvwthIgQQvwqhKipdk7FJYT4PyFEVM732CuEeEntnIpDCDFYCHFRCKEVQpSZMzmEEL2EEJeFENFCiDlq51NcQojNQoiHQogLaudSEkKI2kKI34UQkTm/n6lq51RcQghrIUSoEOJ8znd48UMojJWLKQ4BCSHsiSgx5/UUAJ5E9J7KaRWLEKIHgGNEpBFCrAAAIpqtclpFJoRoBEALYD0KeZ50aSOEMAdwBUB3ZD8K9QyAoUQUqWpixSCE6AQgGcA2Imqidj7FJYSoAaAGEZ0VQtgBCAfwnzL2NxAAKhFRshDCEsCfAKYSUYjSuZjkHkBu55+jErKv3itTiOhXItLkvA0BUEvNfIqLiC4R0WW18yimVgCiieg6EWUA8AfQT+WcioWI/gDwWO08SoqI7hHR2ZzXSQAuAXBRN6viyblgNznnrWXOpEofZJIFAACEEEuFELcBDAfwmdr5GGgcgENqJ2ECXADczvP+DspY51OeCCHcADQD8Je6mRSfEMJcCBEB4CGAI0SkyncotwVACPGbEOKCnqkfABDRfCKqDeB/ACapm61+L/oOOcvMB6BB9vcoVYqSP2MlIYSwBbAHwEf59ujLBCLKIqLXkL3n3koIocpwnIUaQZVARN2KuOj/AAQBWGDEdErkRd9BCDEGQB8AXakUHswpxt+grLgLIO9N+mvlzGMKyhk33wPgf0T0k9r5GIKIEoQQvwPoBUDxA/Pldg/geYQQ9fK87QcgSq1cSkoI0QvALAB9iShV7XxMxBkA9YQQLwshrAC8DeBnlXMyKTkHUL8HcImIVqmdT0kIIarknrUnhLBB9kkFqvRBpnoW0B4ADZB9FspNAO8RUZnakhNCRAOoACA+Z1ZIWTqTSQjRH8A3AKoASAAQQUQ91c3qxYQQ3gC+BmAOYDMRLVU5pWIRQvgBeB3Zd6J8AGABEX2valLFIIToACAYwD/I/vcLAPOIKEi9rIpHCPEqgK3I/g2ZAQggosWq5GKKBYAxxpiJDgExxhjjAsAYYyaLCwBjjJkoLgCMMWaiuAAwxpiJ4gLAGGMmigsAY4yZKC4AjDFmov4fsafew+ynta4AAAAASUVORK5CYII=\n",
- "text/plain": [
- "<Figure size 432x288 with 1 Axes>"
- ]
- },
- "metadata": {
- "needs_background": "light"
- },
- "output_type": "display_data"
- }
- ],
+ "outputs": [],
"source": [
"plt.quiver(vdir[..., __spherical_indeces__['azimuth']],\n",
" vdir[..., __spherical_indeces__['elevation']],\n",