diff --git a/doc/images/external-fields-implicit-field.svg b/doc/images/external-fields-implicit-field.svg new file mode 100644 index 000000000..a56edc8dd --- /dev/null +++ b/doc/images/external-fields-implicit-field.svg @@ -0,0 +1,1982 @@ + + + + + + + + 2026-04-07T15:46:02.965607 + image/svg+xml + + + Matplotlib v3.10.8, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc/images/external-fields-morph.svg b/doc/images/external-fields-morph.svg new file mode 100644 index 000000000..7738b9ba7 --- /dev/null +++ b/doc/images/external-fields-morph.svg @@ -0,0 +1,6952 @@ + + + + + + + + 2026-04-07T13:20:59.401673 + image/svg+xml + + + Matplotlib v3.10.8, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc/images/external-fields-single-cell.svg b/doc/images/external-fields-single-cell.svg new file mode 100644 index 000000000..bd1e60ba1 --- /dev/null +++ b/doc/images/external-fields-single-cell.svg @@ -0,0 +1,1499 @@ + + + + + + + + 2026-04-07T12:04:23.578217 + image/svg+xml + + + Matplotlib v3.10.8, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc/python/morphology.rst b/doc/python/morphology.rst index 8ae095dc3..b6d459f47 100644 --- a/doc/python/morphology.rst +++ b/doc/python/morphology.rst @@ -606,19 +606,19 @@ region. .. py:class:: loaded_morphology - .. py:attr:: segment_tree + .. py:attribute:: segment_tree Raw segment tree, identical to morphology. - .. py:attr:: morphology + .. py:attribute:: morphology Morphology constructed from description. - .. py:attr:: labels + .. py:attribute:: labels Regions and locsets defined in the description as ``label_dict`` - .. py:attr:: metadata + .. py:attribute:: metadata Loader specific metadata, see below in the individual sections. diff --git a/doc/tutorial/cosim.rst b/doc/tutorial/cosim.rst index 2c3732dd1..b3a32fd02 100644 --- a/doc/tutorial/cosim.rst +++ b/doc/tutorial/cosim.rst @@ -12,7 +12,7 @@ tutorial, we are assuming you are comfortable with the basics of Arbor (cells, recipes, and networks), Python package management, and MPI, as wells as installing software on your system. All source code for all intermediate steps can be in the directory -`python/example/cosim `__ +`python/example/cosim `__ of the Arbor source tree. Setup diff --git a/doc/tutorial/index.rst b/doc/tutorial/index.rst index 761e1f85c..bb3c50c5b 100644 --- a/doc/tutorial/index.rst +++ b/doc/tutorial/index.rst @@ -1,10 +1,9 @@ - .. _tutorial: Tutorials ========= -Grouped loosely by primary (but not exclusive!) focus, we have a set of tutorials to help you learn by doing. +Grouped loosely by primary (but not exclusive!) focus, we have a set of tutorials to help you learn by doing. You can find some examples of full Arbor simulations in the ``python/examples`` directory of the `Arbor repository `_. @@ -33,7 +32,7 @@ Recipes .. toctree:: :maxdepth: 1 - + single_cell_recipe single_cell_detailed_recipe @@ -47,7 +46,7 @@ Networks network_two_cells_gap_junctions brunel - + Probes ------ @@ -84,6 +83,7 @@ Advanced connectivity full-feature-diffusion cosim + reading_external_fields Demonstrations -------------- diff --git a/doc/tutorial/reading_external_fields.rst b/doc/tutorial/reading_external_fields.rst new file mode 100644 index 000000000..1c875ef98 --- /dev/null +++ b/doc/tutorial/reading_external_fields.rst @@ -0,0 +1,418 @@ +.. _tutorial_external_fields: + +Interacting with External Fields +================================ + +In this tutorial we are going to discuss how to use Arbor and an external simulation +of smooth fields to facilitate interaction at the microscopic level, e.g. + +- stimulation by electric fields as induced by direct currents or magnetic fields [#f1]_ +- simulations of the extra-cellular medium [#f2]_ + +As a concrete example, we will model the influence of a time-varying electric +field on a single cell. + +This is an experimental technique leaning into the implementation details of +Arbor and ``modcc`` and as such might be changed sooner or later. For a valuable +outcome, readers should be familiar with modeling cable cells in Arbor, NMODL +and its use in Arbor, the basics of electro-dynamics, and some low-level (C++) +programming. + +All source code for all intermediate steps can be in the directory +`python/example/reading-external-fields +`__ +of the Arbor source tree. + + +Single Cell Model +----------------- + +As the first step, we create a single cell model with a sufficiently complex +geometry, following the standard method outlined in for example :ref:`tutorialsinglecellrecipe` + +.. literalinclude:: ../../python/example/reading-external-fields/single-cell.py + :language: python + +For convenience and the purposes of this tutorial, we store the actual cable +cell data as global variables outside the recipe + +.. literalinclude:: ../../python/example/reading-external-fields/single-cell.py + :linenos: + :lineno-match: + :lines: 9-23 + :language: python + +We add a spiking input to the soma centre + +.. literalinclude:: ../../python/example/reading-external-fields/single-cell.py + :linenos: + :lineno-match: + :lines: 9 + :language: python + +with a regular rate of :math:`\nu = 1/50ms` + +.. literalinclude:: ../../python/example/reading-external-fields/single-cell.py + :linenos: + :lineno-match: + :lines: 42-47 + :language: python + +For reference, this is the voltage trace after simulation: + +.. figure:: ../images/external-fields-single-cell.svg + :width: 600 + :align: center + +Theory: Injecting Currents by Electric Fields +--------------------------------------------- + +Now, we want to apply the effects of an electric field :math:`\vec E(\vec x, t)` +to each segment of the morphology. We follow the quasi-potential approach [#f3]_ + +.. math:: + + \phi = - \int \mathrm{d}\vec l\cdot \vec E + +and add an extra term to the cable equation + +.. math:: + + I_\mathrm{ext} = -\sigma \frac{\partial^2}{\partial x^2} \phi + +This term can then be absorbed into the cable equation by adding an extra +mechanism per segment, feeding it the segment's endpoints :math:`\vec x` and +:math:`\vec x'` in 3d space and a handle to the vector field :math:`\vec E(\vec +x, t)`. Then, we compute + +.. math:: + + I_\mathrm{ext} \simeq -\frac{1}{r_\mathrm{L}} \frac{\vec E \cdot \vec a}{||a||}\\ + \vec a = \vec x - \vec x' + +which converges to the proper result iff the segment lengths are sufficiently +small [#f4]_, where :math:`r_\mathrm{L}` is the familiar axial resistivity. + +Let's step back for a moment and take a look at the morphology + +.. figure:: ../images/external-fields-morph.svg + :width: 600 + :align: center + +which has been coloured according to the angle to a uniform vector field :math:`\vec E = (1, 0, 0)` + +Stepping Stone: Implicit Electric Fields +---------------------------------------- + +We will now show how to manipulate a standard NMODL mechanism into an injected +current for a homogeneous, varying electric field :math:`\vec E = (E_0 +\sin(\omega t), 0, 0)`. This allows us to demonstrate the technique without +delving into the implementation details of transferring data from and to C++. + +First, we will create a point mechanism in NMODL, then compile it to C++ using +``modcc``. In the next section, we use the compiled C++ code and change the +implementation to produce the current we want. Without further ado + +.. literalinclude:: ../../python/example/reading-external-fields/mod/efield.mod + :linenos: + :lineno-match: + +This code formalises what we discussed above, if you are unfamiliar with writing +mechanisms, please refer to :ref:`tutorial_nmodl_density`. We pass the desired +amplitude and frequency of the E-field, as well as, the distal and proximal +coordinates as parameters + +.. literalinclude:: ../../python/example/reading-external-fields/mod/efield.mod + :linenos: + :lineno-match: + :lines: 7-12 + +The time is modelled via an ODE + +.. literalinclude:: ../../python/example/reading-external-fields/mod/efield.mod + :linenos: + :lineno-match: + :lines: 39 + +which is solved to compute the induced current as outlined above + +.. literalinclude:: ../../python/example/reading-external-fields/mod/efield.mod + :linenos: + :lineno-match: + :lines: 30-37 + +Compile the catalogue + +.. code-block:: + + $ arbor-build-catalogue efields mod + Building catalogue 'efields' from mechanisms in [...]/mod + * NMODL + * efield + Catalogue has been built and copied to [...]/efields-catalogue.so + +and it is ready to use. We create a new simulation derived from the one above + +.. literalinclude:: ../../python/example/reading-external-fields/implicit-field.py + :linenos: + :lineno-match: + :language: python + :lines: -73 + +As promised, we add one point process per segment, passing in the endpoints and +field parameters, after loading our catalogue into the model + +.. literalinclude:: ../../python/example/reading-external-fields/implicit-field.py + :linenos: + :lineno-match: + :language: python + :lines: 23-44 + +We plot the membrane potential and induced current over time + +.. figure:: ../images/external-fields-implicit-field.svg + :width: 600 + :align: center + +and clearly observe the influence of the driving current and non-linearity of +the HH-mechanism. + +If you can -- analytically or otherwise -- determine the value of :math:`\vec +E(\vec r, t)` at the coordinates of the segment centres without reliance on +external tools, this approach is sufficient and there is no need for the +techniques demonstrated in the rest of this tutorial. The same holds for other +quantities interacting with the simulation. + +(Ab)using NMODL to Create Custom Mechanisms +------------------------------------------- + +Now, this gives us a current on every segment, but for interfacing with an +external simulation, we need to enable the mechanism to read from that source. +For this, we create a new point mechanism with a similar layout as above + +.. literalinclude:: ../../python/example/reading-external-fields/mod/template.mod + :linenos: + :lineno-match: + +This will not be used in a simulation, but rather we will compile it to C++ + +.. code-block:: + + $ arbor-build-catalogue --debug tmp efields mod + Building catalogue 'efields' from mechanisms in [...]/mod + * NMODL + * template + * efield + Building debug code in '[...]/tmp'. + Catalogue has been built and copied to [...]/efields-catalogue.so + +By using the ``debug`` mode, the intermediate C++ files will be preserved in the +``tmp`` directory, and we can edit and add them back to the catalogue. Open +``tmp/build/generated/efields/template_cpu.cpp`` and take a look at the +``compute_currents`` function + +.. code-block:: c++ + + static void compute_currents(arb_mechanism_ppack* pp) { + PPACK_IFACE_BLOCK; + for (arb_size_type i_ = 0; i_ < _pp_var_width; ++i_) { + auto node_indexi_ = _pp_var_node_index[i_]; + arb_value_type current_ = 0; + arb_value_type i = 0; + arb_value_type e; + e = 42.0; + i = e*(_pp_var_xd[i_]-_pp_var_xp[i_])*_pp_var_da[i_]; + current_ = i; + _pp_var_vec_i[node_indexi_] = fma(_pp_var_weight[i_], current_, _pp_var_vec_i[node_indexi_]); + } + } + +Note that if you have SIMD enabled, this might look differently. We recommend +disabling SIMD for this tutorial, since the C++ is simpler to follow. + +We are going to continue, as before, under the assumption that the electric +field is homogeneous, so we need only to concern ourselves with three ``double`` +numbers. Generalizing to complete fields is straightforward, if tedious [#f5]_. + +Copy the two files ``tmp/build/generated/efields/template_cpu.cpp`` and +``tmp/build/generated/efields/template.hpp`` into the ``mod`` directory and +rename them to ``mod/reader.hpp`` and ``mod/reader_cpu.hpp``. + +.. code-block:: + + cp tmp/build/generated/efields/template.hpp mod/reader.hpp + cp tmp/build/generated/efields/template_cpu.cpp mod/reader_cpu.cpp + +then, replace the word ``template`` with ``reader`` in both files, either using +a text editor or + +.. code-block:: + + sed -i '' 's/template/reader/g' mod/reader_cpu.cpp + sed -i '' 's/template/reader/g' mod/reader.hpp + +on MacOS while Linux users drop the ``''``. Now, you should be able to compile +the catalogue again + +.. code-block:: + + $ arbor-build-catalogue fields mod --raw reader + Building catalogue 'efields' from mechanisms in [...]/mod + * NMODL + * template + * efield + * Raw + * reader + Building debug code in '[...]/tmp'. + Catalogue has been built and copied to [...]/efields-catalogue.so + +This has built an extra point process ``reader`` from a raw C++ implementation. +Let's take it for a test drive using a new simulation + +.. literalinclude:: ../../python/example/reading-external-fields/external-field.py + :linenos: + :lineno-match: + :language: python + :lines: -79 + +We added the electric field vector + +.. literalinclude:: ../../python/example/reading-external-fields/external-field.py + :linenos: + :lineno-match: + :language: python + :lines: 10 + +and swapped the point process for ``reader`` passing nono-sense data to ``field`` + +.. literalinclude:: ../../python/example/reading-external-fields/external-field.py + :linenos: + :lineno-match: + :language: python + :lines: 32-46 + +and our simulation is now stepping in increments of 1ms after which we change the +x-component of the electric field + +.. literalinclude:: ../../python/example/reading-external-fields/external-field.py + :linenos: + :lineno-match: + :language: python + :lines: 76-79 + +Final Step: Getting Data from NumPy into C++ +-------------------------------------------- + +Now we need to wire our 'simulation' of the electric field into the custom mechanism. +To that end, we need to acquire the pointer to the backing data of the NumPy array. + +.. literalinclude:: ../../python/example/reading-external-fields/external-field-data.py + :linenos: + :lineno-match: + :language: python + :lines: 7-11 + +We use ``ctypes`` to fetch the pointer from the array. However, this is a +somewhat brittle process and the machinery relies on this pointer never +changing again. This means, you can write to the array like this + +.. code-block:: python + + E[0] = 31.0 + +which does not change the backing memory, just the values contained within. +This, however, does + +.. code-block:: python + + x = np.array(...) + y = np.array(...) + E = x + y + +and thus breaks the programm, as do all operations that assign a new array to +``E``. + +Next, we pass the pointer to the mechanism + +.. literalinclude:: ../../python/example/reading-external-fields/external-field-data.py + :linenos: + :lineno-match: + :language: python + :emphasize-lines: 6 + :lines: 34-48 + +Note that we pass an integral value as a floating point number, which is +questionable, but a fix requires tweaking the NMODL language and might happen in +the future. For demonstration purposes, we shift the sine wave by 50ms versus +the implicit field example + +.. literalinclude:: ../../python/example/reading-external-fields/external-field-data.py + :linenos: + :lineno-match: + :language: python + :emphasize-lines: 3 + :lines: 77-81 + +All that's left to do is to read the field value from the ``reader`` mechanism + +.. literalinclude:: ../../python/example/reading-external-fields/mod/reader_cpu.cpp + :linenos: + :lineno-match: + :language: c++ + :lines: 89-98 + +First, we cast the floating point value to a double-valued array + +.. literalinclude:: ../../python/example/reading-external-fields/mod/reader_cpu.cpp + :linenos: + :lineno-match: + :language: c++ + :lines: 91 + +then we read the first value :math:`E_x`, use it to compute the current ``i`` as :math:`I = \frac{E_x a_x}{|a|}` + +.. literalinclude:: ../../python/example/reading-external-fields/mod/reader_cpu.cpp + :linenos: + :lineno-match: + :language: c++ + :lines: 94-95 + +which is written back the data consumed by Arbor's cable equation solver + +.. literalinclude:: ../../python/example/reading-external-fields/mod/reader_cpu.cpp + :linenos: + :lineno-match: + :language: c++ + :lines: 96 + +Conclusions +----------- + +We have seen how custom C++ mechanisms can source data from external data and +they can be created, relatively simply, from NMODL templates. There are many +ways to use this, examples include driving two simulations exchanging data over +a shared memory segment, using a network connection to obtain values from +external sources, and many more. The largest hurdle here is the synchronisation +between Arbor and the external data source. + +It might be prudent to note that this approach will be slow, both in terms of +setup and execution. We are explicitly adding one point process to each segment +and setting its coordinates. Also, technically, setting both proximal and distal +coordinates is redundant as segments (mostly) share their endpoints. The way the +simulation is set up, the time is stepped at the granularity of the electric +field evolution. We would much rather use asynchronous simulations and a +lightweight signalling mechanism. + +References and Remarks +---------------------- + +.. [#f1] Yonezawa et al *Biophysical Simulation of the Effects of Transcranial Magnetic Stimulation (TMS) on a Cortical Column*, ICNIP 2025 +.. [#f2] McDougal et al *Reaction-diffusion in the NEURON simulator*, FNINF 2013 +.. [#f3] Wang et al, *Coupling Magnetically Induced Electric Fields to Neurons: Longitudinal and Transverse Activation*, Biophys 2018 +.. [#f4] Pashut et al, *Mechanisms of Magnetic Stimulation of Central Nervous System Neurons*, Plos CompBio 2011 +.. [#f5] A general outline might be + + 1. Add parameters for the extents both in terms of colocation points and physical dimensions. + 2. Shift the segment endpoints into the reference frame of the field layout above. + 3. Compute the indices into the array and retrieve the value at index. + 4. Continue as outlined in the tutorial. diff --git a/python/example/reading-external-fields/Acker2008.swc b/python/example/reading-external-fields/Acker2008.swc new file mode 100644 index 000000000..3a6ae208b --- /dev/null +++ b/python/example/reading-external-fields/Acker2008.swc @@ -0,0 +1,1185 @@ +# SWC neuron morphology file exported by +# Data History - +# +# + 1 1 -53.420 7.500 -5.950 7.500 -1 + 2 2 -54.060 -4.250 -5.960 1.250 1 + 3 2 -57.100 -6.710 1.340 0.655 2 + 4 2 -57.420 -14.480 1.540 0.613 3 + 5 2 -59.660 -18.370 2.460 0.480 4 + 6 2 -62.220 -20.960 4.940 0.480 5 + 7 2 -63.180 -25.180 7.200 0.480 6 + 8 2 -65.100 -29.060 10.100 0.480 7 + 9 3 -54.060 -1.660 -5.960 1.250 1 + 10 3 -46.010 -2.930 -0.020 1.120 9 + 11 3 -40.250 -4.540 -1.600 0.800 10 + 12 3 -36.368 -6.407 -1.016 0.800 11 + 13 3 -33.866 -7.550 0.690 0.800 12 + 14 3 -29.246 -9.330 -9.550 0.640 13 + 15 3 -25.336 -8.930 -9.030 0.480 14 + 16 3 -18.626 -8.930 -10.910 0.480 15 + 17 3 -13.826 -6.990 -10.910 0.480 16 + 18 3 -9.356 -6.020 -8.590 0.480 17 + 19 3 -4.556 -5.040 -8.590 0.480 18 + 20 3 -0.396 -4.070 -10.830 0.480 19 + 21 3 5.044 -1.800 -7.450 0.320 20 + 22 3 8.234 1.440 -4.490 0.320 21 + 23 3 11.434 3.050 -5.770 0.320 22 + 24 3 14.314 7.590 -3.250 0.320 23 + 25 3 8.874 -3.100 -7.570 0.320 21 + 26 3 12.394 -4.720 -7.570 0.320 25 + 27 3 16.874 -6.020 -3.270 0.320 26 + 28 3 21.024 -4.720 -3.270 0.320 27 + 29 3 27.104 -5.690 -3.270 0.320 28 + 30 3 34.144 -5.040 -3.270 0.320 29 + 31 3 42.134 -5.370 -5.950 0.320 30 + 32 3 -28.926 -14.190 -13.170 0.640 14 + 33 3 -25.836 -15.300 -16.870 0.320 32 + 34 3 -26.156 -17.890 -16.870 0.320 33 + 35 3 -24.024 -22.264 -18.550 0.320 34 + 36 3 -20.516 -26.376 -20.924 0.320 35 + 37 3 -19.246 -30.256 -23.784 0.320 36 + 38 3 -19.886 -33.496 -26.964 0.320 37 + 39 3 -17.326 -36.416 -28.544 0.247 38 + 40 3 -15.726 -39.336 -32.184 0.265 39 + 41 3 -15.086 -40.306 -42.104 0.306 40 + 42 3 -11.566 -38.036 -42.944 0.320 41 + 43 3 -8.366 -37.716 -42.944 0.320 42 + 44 3 -5.486 -39.006 -46.344 0.320 43 + 45 3 -3.896 -41.596 -46.344 0.320 44 + 46 3 -1.976 -44.196 -46.344 0.320 45 + 47 3 0.264 -46.786 -50.644 0.260 46 + 48 3 2.824 -48.726 -50.684 0.238 47 + 49 3 5.064 -50.026 -55.644 0.250 48 + 50 3 10.814 -52.616 -61.884 0.228 49 + 51 3 14.514 -55.326 -64.804 0.211 50 + 52 3 16.434 -57.266 -64.804 0.152 51 + 53 3 18.994 -60.186 -70.144 0.160 52 + 54 3 -14.126 -42.571 -44.114 0.252 41 + 55 3 -12.728 -46.195 -45.327 0.238 54 + 56 3 -11.732 -50.174 -47.162 0.224 55 + 57 3 -10.199 -52.853 -48.131 0.197 56 + 58 3 -5.399 -56.093 -50.231 0.194 57 + 59 3 -4.759 -58.363 -50.331 0.214 58 + 60 3 -2.839 -60.633 -52.211 0.185 59 + 61 3 -1.889 -63.873 -52.211 0.212 60 + 62 3 0.671 -65.813 -52.211 0.203 61 + 63 3 1.311 -71.643 -60.271 0.191 62 + 64 3 3.551 -74.883 -63.351 0.207 63 + 65 3 5.151 -78.123 -67.451 0.148 64 + 66 3 5.151 -81.693 -71.031 0.160 65 + 67 3 -30.846 -17.110 -16.150 0.541 32 + 68 3 -32.766 -18.730 -16.150 0.540 67 + 69 3 -32.126 -23.270 -17.470 0.551 68 + 70 3 -33.406 -25.210 -17.470 0.534 69 + 71 3 -35.646 -27.800 -18.050 0.549 70 + 72 3 -39.156 -30.720 -18.770 0.547 71 + 73 3 -40.116 -34.600 -20.110 0.566 72 + 74 3 -39.906 -40.240 -24.570 0.480 73 + 75 3 -39.906 -45.750 -24.570 0.480 74 + 76 3 -38.668 -51.661 -25.382 0.480 75 + 77 3 -38.572 -56.527 -25.946 0.480 76 + 78 3 -39.342 -61.021 -26.780 0.480 77 + 79 3 -40.302 -64.911 -26.780 0.504 78 + 80 3 -39.662 -71.391 -28.040 0.480 79 + 81 3 -39.982 -76.571 -29.980 0.320 80 + 82 3 -41.216 -80.587 -31.934 0.320 81 + 83 3 -42.920 -84.501 -32.904 0.320 82 + 84 3 -45.160 -87.741 -30.504 0.320 83 + 85 3 -46.440 -90.651 -31.344 0.320 84 + 86 3 -47.400 -94.541 -31.344 0.320 85 + 87 3 -47.400 -98.431 -32.864 0.320 86 + 88 3 -47.720 -102.641 -34.304 0.320 87 + 89 3 -49.000 -105.561 -34.304 0.320 88 + 90 3 -48.680 -109.441 -35.684 0.320 89 + 91 3 -47.400 -113.331 -35.684 0.320 90 + 92 3 -46.760 -116.251 -34.264 0.320 91 + 93 3 -46.760 -120.141 -31.004 0.223 92 + 94 3 -46.120 -123.381 -26.584 0.193 93 + 95 3 -44.200 -126.621 -26.584 0.199 94 + 96 3 -41.716 -35.410 -23.290 0.565 73 + 97 3 -45.096 -37.479 -24.752 0.562 96 + 98 3 -47.299 -41.094 -25.735 0.518 97 + 99 3 -48.612 -44.234 -27.368 0.546 98 + 100 3 -49.572 -48.454 -30.748 0.543 99 + 101 3 -51.812 -53.634 -30.748 0.495 100 + 102 3 -52.020 -58.504 -31.286 0.504 101 + 103 3 -51.485 -63.294 -31.630 0.529 102 + 104 3 -51.508 -68.690 -33.076 0.506 103 + 105 3 -53.108 -72.260 -35.296 0.479 104 + 106 3 -54.388 -76.470 -37.036 0.387 105 + 107 3 -56.938 -80.360 -39.256 0.386 106 + 108 3 -56.298 -84.890 -40.956 0.399 107 + 109 3 -57.698 -89.240 -45.596 0.369 108 + 110 3 -58.968 -92.800 -52.736 0.320 109 + 111 3 -57.698 -97.980 -51.616 0.320 110 + 112 3 -62.488 -94.740 -53.496 0.274 110 + 113 3 -66.008 -97.660 -55.516 0.226 112 + 114 3 -68.888 -100.250 -57.636 0.195 113 + 115 3 -72.408 -102.520 -59.716 0.160 114 + 116 3 -33.454 -9.753 0.736 0.690 13 + 117 3 -32.081 -13.934 0.337 0.614 116 + 118 3 -30.188 -16.751 0.108 0.621 117 + 119 3 -26.988 -17.081 2.608 0.560 118 + 120 3 -26.988 -20.321 3.708 0.562 119 + 121 3 -27.628 -23.231 8.008 0.486 120 + 122 3 -26.668 -26.471 10.388 0.516 121 + 123 3 -24.758 -30.361 14.628 0.500 122 + 124 3 -22.838 -32.951 14.628 0.481 123 + 125 3 -22.518 -36.191 14.208 0.508 124 + 126 3 -21.238 -38.461 14.208 0.461 125 + 127 3 -20.918 -41.701 14.208 0.456 126 + 128 3 -21.238 -44.291 15.168 0.445 127 + 129 3 -19.638 -47.861 15.168 0.436 128 + 130 3 -17.458 -53.202 15.019 0.410 129 + 131 3 -16.847 -57.954 13.855 0.433 130 + 132 3 -14.927 -62.494 13.775 0.394 131 + 133 3 -13.007 -63.794 12.175 0.339 132 + 134 3 -12.687 -67.354 12.175 0.336 133 + 135 3 -9.817 -69.624 12.175 0.270 134 + 136 3 -8.574 -74.473 11.450 0.350 135 + 137 3 -5.827 -78.113 11.213 0.349 136 + 138 3 -3.080 -81.753 10.977 0.346 137 + 139 3 -1.160 -85.963 10.977 0.330 138 + 140 3 2.680 -89.523 10.977 0.310 139 + 141 3 4.600 -91.793 10.977 0.325 140 + 142 3 4.600 -96.003 10.977 0.296 141 + 143 3 4.600 -100.543 9.397 0.274 142 + 144 3 7.504 -105.610 8.349 0.277 143 + 145 3 12.467 -109.418 8.221 0.263 144 + 146 3 14.387 -112.008 8.221 0.262 145 + 147 3 15.337 -116.548 7.741 0.295 146 + 148 3 16.297 -121.728 7.741 0.289 147 + 149 3 18.537 -125.618 9.161 0.246 148 + 150 3 20.457 -131.128 9.221 0.252 149 + 151 3 20.457 -135.018 7.521 0.206 150 + 152 3 19.497 -136.958 6.101 0.216 151 + 153 3 20.457 -146.678 7.981 0.160 152 + 154 3 20.607 -149.468 9.521 0.160 153 + 155 3 -54.060 -4.250 -5.860 1.250 1 + 156 3 -49.950 -7.220 3.600 1.078 155 + 157 3 -48.670 -12.080 3.780 0.735 156 + 158 3 -48.670 -17.270 4.700 0.603 157 + 159 3 -47.390 -19.860 5.600 0.559 158 + 160 3 -47.523 -24.137 7.008 0.544 159 + 161 3 -47.787 -28.401 7.217 0.538 160 + 162 3 -46.397 -32.726 7.605 0.548 161 + 163 3 -46.523 -36.824 8.233 0.555 162 + 164 3 -44.603 -41.034 6.453 0.518 163 + 165 3 -42.043 -46.864 5.833 0.410 164 + 166 3 -40.931 -51.092 1.762 0.360 165 + 167 3 -39.371 -55.006 0.727 0.320 166 + 168 3 -36.979 -57.042 0.736 0.320 167 + 169 3 -34.782 -59.944 0.399 0.320 168 + 170 3 -32.088 -62.488 0.371 0.320 169 + 171 3 -29.848 -65.078 0.371 0.320 170 + 172 3 -28.248 -68.968 -2.169 0.320 171 + 173 3 -25.368 -69.938 -2.169 0.320 172 + 174 3 -22.808 -71.558 -2.969 0.320 173 + 175 3 -22.168 -74.798 -1.929 0.320 174 + 176 3 -19.928 -78.368 -0.229 0.320 175 + 177 3 -17.698 -82.898 -0.229 0.320 176 + 178 3 -17.378 -88.088 -3.449 0.320 177 + 179 3 -17.058 -92.618 -3.449 0.320 178 + 180 3 -16.418 -96.508 0.131 0.320 179 + 181 3 -15.178 -99.754 -0.075 0.320 180 + 182 3 -13.399 -101.604 -0.737 0.320 181 + 183 3 -12.171 -104.963 -1.239 0.320 182 + 184 3 -10.571 -108.853 1.121 0.320 183 + 185 3 -10.891 -113.393 3.021 0.320 184 + 186 3 -11.519 -117.770 4.264 0.320 185 + 187 3 -12.120 -122.248 4.078 0.320 186 + 188 3 -13.423 -127.170 5.023 0.320 187 + 189 3 -13.743 -132.360 7.783 0.320 188 + 190 3 -16.133 -137.050 8.683 0.320 189 + 191 3 -16.105 -139.642 7.714 0.320 190 + 192 3 -15.710 -142.279 7.455 0.320 191 + 193 3 -16.350 -146.809 8.095 0.320 192 + 194 3 -19.230 -149.729 9.435 0.320 193 + 195 3 -19.870 -153.609 10.595 0.320 194 + 196 3 -21.790 -156.849 13.535 0.320 195 + 197 3 -20.830 -160.419 13.535 0.320 196 + 198 3 -21.470 -163.979 14.795 0.246 197 + 199 3 -22.110 -169.169 16.355 0.202 198 + 200 3 -21.790 -174.029 16.375 0.201 199 + 201 3 -23.710 -178.559 14.595 0.202 200 + 202 3 -24.990 -182.129 14.595 0.187 201 + 203 3 -25.950 -186.009 14.595 0.160 202 + 204 3 -27.230 -188.279 14.595 0.160 203 + 205 3 -26.590 -191.519 17.275 0.160 204 + 206 3 -29.610 -196.579 17.275 0.160 205 + 207 3 -33.130 -200.469 19.655 0.160 206 + 208 3 -54.060 -4.250 -5.950 1.250 1 + 209 3 -58.330 -5.390 -5.200 0.957 208 + 210 3 -60.570 -9.920 -19.840 0.800 209 + 211 3 -65.520 -9.130 -22.500 0.544 210 + 212 3 -69.680 -7.510 -24.560 0.480 211 + 213 3 -74.160 -4.600 -25.900 0.480 212 + 214 3 -77.990 -2.330 -26.920 0.480 213 + 215 3 -79.730 -0.240 -30.980 0.480 214 + 216 3 -83.098 2.124 -35.508 0.480 215 + 217 3 -85.955 4.956 -40.509 0.480 216 + 218 3 -87.875 8.836 -45.249 0.480 217 + 219 3 -89.795 12.076 -47.129 0.480 218 + 220 3 -90.435 15.646 -48.609 0.480 219 + 221 3 -92.675 16.616 -53.969 0.480 220 + 222 3 -94.275 19.526 -59.469 0.320 221 + 223 3 -95.875 22.126 -60.809 0.320 222 + 224 3 -98.755 24.386 -60.869 0.320 223 + 225 3 -102.265 27.626 -63.749 0.320 224 + 226 3 -104.825 31.196 -66.049 0.320 225 + 227 3 -107.385 33.786 -70.129 0.320 226 + 228 3 -82.790 -2.650 -28.060 0.370 214 + 229 3 -87.590 -2.980 -26.000 0.369 228 + 230 3 -90.763 -4.278 -25.829 0.402 229 + 231 3 -95.400 -5.559 -25.739 0.389 230 + 232 3 -99.767 -7.014 -25.498 0.450 231 + 233 3 -102.790 -7.459 -25.391 0.365 232 + 234 3 -106.310 -6.809 -24.611 0.320 233 + 235 3 -108.550 -7.459 -24.031 0.320 234 + 236 3 -111.110 -8.759 -23.331 0.320 235 + 237 3 -114.048 -12.051 -23.429 0.302 236 + 238 3 -116.649 -15.992 -23.457 0.324 237 + 239 3 -118.676 -18.162 -23.352 0.305 238 + 240 3 -122.506 -19.462 -22.932 0.307 239 + 241 3 -124.930 -20.824 -23.543 0.301 240 + 242 3 -128.255 -22.108 -24.483 0.256 241 + 243 3 -131.712 -24.009 -25.028 0.243 242 + 244 3 -133.632 -25.629 -27.428 0.234 243 + 245 3 -136.832 -28.869 -27.748 0.202 244 + 246 3 -140.672 -28.869 -28.448 0.200 245 + 247 3 -145.782 -29.519 -28.508 0.244 246 + 248 3 -149.622 -30.809 -28.948 0.246 247 + 249 3 -155.382 -31.779 -27.388 0.214 248 + 250 3 -159.852 -32.429 -26.628 0.224 249 + 251 3 -164.652 -31.139 -24.128 0.205 250 + 252 3 -168.492 -29.189 -22.088 0.198 251 + 253 3 -177.122 -26.919 -22.688 0.250 252 + 254 3 -182.102 -26.779 -22.988 0.213 253 + 255 3 -186.262 -27.749 -24.568 0.207 254 + 256 3 -191.702 -28.719 -26.068 0.200 255 + 257 3 -199.052 -30.989 -26.648 0.189 256 + 258 3 -210.562 -32.609 -27.068 0.195 257 + 259 3 -60.570 -14.460 -21.620 0.640 210 + 260 3 -59.297 -18.328 -22.197 0.640 259 + 261 3 -58.861 -22.600 -22.725 0.640 260 + 262 3 -59.203 -26.538 -23.530 0.640 261 + 263 3 -59.843 -29.778 -25.810 0.640 262 + 264 3 -63.233 -30.638 -25.590 0.480 263 + 265 3 -66.433 -32.578 -26.670 0.387 264 + 266 3 -68.033 -35.818 -27.070 0.411 265 + 267 3 -71.350 -38.118 -27.057 0.409 266 + 268 3 -75.262 -40.341 -27.546 0.408 267 + 269 3 -79.098 -43.335 -27.250 0.404 268 + 270 3 -82.288 -45.605 -27.250 0.391 269 + 271 3 -86.128 -47.875 -28.730 0.403 270 + 272 3 -90.405 -51.038 -28.876 0.362 271 + 273 3 -95.435 -53.990 -29.034 0.320 272 + 274 3 -101.195 -55.280 -31.314 0.320 273 + 275 3 -104.352 -56.462 -32.367 0.320 274 + 276 3 -106.847 -58.815 -33.579 0.292 275 + 277 3 -108.592 -62.898 -34.595 0.316 276 + 278 3 -109.728 -66.517 -36.074 0.352 277 + 279 3 -112.928 -68.457 -38.554 0.364 278 + 280 3 -116.128 -69.437 -38.554 0.341 279 + 281 3 -119.638 -71.377 -41.954 0.335 280 + 282 3 -123.848 -73.662 -41.956 0.302 281 + 283 3 -127.498 -76.496 -41.933 0.303 282 + 284 3 -131.338 -78.116 -41.453 0.273 283 + 285 3 -134.858 -78.766 -41.453 0.239 284 + 286 3 -138.698 -79.086 -41.453 0.240 285 + 287 3 -141.248 -77.146 -39.793 0.276 286 + 288 3 -145.728 -76.496 -39.793 0.239 287 + 289 3 -150.528 -76.496 -39.793 0.250 288 + 290 3 -153.408 -78.766 -40.533 0.276 289 + 291 3 -154.998 -82.326 -44.853 0.239 290 + 292 3 -157.878 -85.246 -44.853 0.260 291 + 293 3 -162.038 -85.246 -48.813 0.243 292 + 294 3 -165.878 -86.216 -48.813 0.253 293 + 295 3 -170.028 -85.896 -50.953 0.253 294 + 296 3 -174.628 -86.036 -52.873 0.217 295 + 297 3 -177.508 -85.716 -55.873 0.191 296 + 298 3 -183.268 -87.016 -58.953 0.213 297 + 299 3 -189.018 -87.986 -58.953 0.171 298 + 300 3 -194.778 -87.986 -58.953 0.204 299 + 301 3 -199.898 -88.306 -56.093 0.176 300 + 302 3 -59.523 -31.728 -24.950 0.640 263 + 303 3 -60.033 -34.848 -27.770 0.507 302 + 304 3 -61.194 -39.364 -27.360 0.464 303 + 305 3 -61.369 -43.906 -26.715 0.421 304 + 306 3 -60.729 -46.496 -24.115 0.452 305 + 307 3 -60.089 -51.036 -21.235 0.456 306 + 308 3 -59.139 -54.596 -20.835 0.454 307 + 309 3 -58.898 -58.816 -20.302 0.445 308 + 310 3 -61.025 -62.868 -19.689 0.436 309 + 311 3 -61.938 -66.912 -18.401 0.417 310 + 312 3 -61.618 -70.802 -16.441 0.347 311 + 313 3 -61.618 -74.692 -16.441 0.374 312 + 314 3 -63.563 -78.060 -15.268 0.368 313 + 315 3 -66.626 -80.751 -14.195 0.363 314 + 316 3 -68.546 -82.041 -13.355 0.361 315 + 317 3 -69.186 -86.581 -13.355 0.362 316 + 318 3 -71.746 -89.501 -13.035 0.391 317 + 319 3 -72.386 -93.061 -12.835 0.406 318 + 320 3 -73.986 -97.601 -12.995 0.330 319 + 321 3 -76.546 -101.161 -12.995 0.306 320 + 322 3 -77.186 -104.401 -12.995 0.334 321 + 323 3 -78.466 -109.581 -15.175 0.333 322 + 324 3 -77.826 -114.771 -15.175 0.195 323 + 325 3 -58.243 -33.018 -27.750 0.482 302 + 326 3 -57.851 -36.548 -27.926 0.407 325 + 327 3 -57.978 -41.481 -28.520 0.407 326 + 328 3 -58.529 -45.965 -29.761 0.406 327 + 329 3 -58.344 -50.639 -30.780 0.406 328 + 330 3 -57.064 -53.559 -31.080 0.362 329 + 331 3 -57.569 -56.735 -32.163 0.358 330 + 332 3 -57.314 -61.817 -33.036 0.353 331 + 333 3 -57.471 -66.283 -34.835 0.344 332 + 334 3 -56.991 -69.038 -35.295 0.369 333 + 335 3 -56.881 -73.860 -35.515 0.382 334 + 336 3 -57.276 -78.837 -35.240 0.394 335 + 337 3 -57.596 -82.717 -35.240 0.381 336 + 338 3 -58.556 -85.317 -36.200 0.389 337 + 339 3 -58.236 -89.847 -36.220 0.372 338 + 340 3 -57.916 -93.087 -37.120 0.394 339 + 341 3 -57.156 -97.827 -37.420 0.391 340 + 342 3 -56.836 -102.037 -38.700 0.343 341 + 343 3 -55.556 -105.927 -40.460 0.363 342 + 344 3 -52.676 -109.487 -40.460 0.392 343 + 345 3 -50.436 -113.057 -42.020 0.305 344 + 346 3 -48.196 -116.297 -43.040 0.272 345 + 347 3 -45.316 -118.887 -44.180 0.220 346 + 348 3 -42.766 -121.157 -40.480 0.253 347 + 349 3 -40.846 -124.397 -40.480 0.230 348 + 350 3 -38.286 -125.367 -40.480 0.239 349 + 351 3 -36.686 -128.927 -37.240 0.250 350 + 352 3 -33.806 -131.517 -37.240 0.205 351 + 353 3 -33.166 -135.407 -37.240 0.238 352 + 354 3 -30.286 -137.997 -36.240 0.160 353 + 355 3 -29.656 -141.887 -36.240 0.160 354 + 356 3 -53.420 3.520 -5.960 1.250 1 + 357 3 -60.300 3.990 0.280 0.640 356 + 358 3 -64.028 3.787 1.455 0.640 357 + 359 3 -68.616 2.577 1.405 0.640 358 + 360 3 -72.550 1.133 1.864 0.640 359 + 361 3 -77.030 0.483 3.784 0.640 360 + 362 3 -80.680 2.633 0.564 0.480 361 + 363 3 -84.200 3.613 -0.576 0.480 362 + 364 3 -88.771 4.452 -1.634 0.480 363 + 365 3 -93.902 6.048 -2.616 0.480 364 + 366 3 -96.462 6.688 -4.516 0.480 365 + 367 3 -100.622 5.718 -5.996 0.480 366 + 368 3 -102.852 5.718 -7.876 0.480 367 + 369 3 -102.852 3.128 -17.756 0.320 368 + 370 3 -101.262 2.478 -21.296 0.320 369 + 371 3 -100.622 -1.082 -24.456 0.320 370 + 372 3 -101.892 -1.732 -27.696 0.320 371 + 373 3 -103.172 -2.702 -35.536 0.320 372 + 374 3 -105.092 -3.032 -41.596 0.320 373 + 375 3 -105.412 -2.052 -46.196 0.320 374 + 376 3 -107.012 -1.412 -47.816 0.320 375 + 377 3 -108.932 -1.412 -50.276 0.320 376 + 378 3 -110.212 0.538 -52.336 0.248 377 + 379 3 -110.212 1.508 -56.156 0.238 378 + 380 3 -104.452 5.398 -7.876 0.480 368 + 381 3 -108.932 5.718 -9.356 0.402 380 + 382 3 -113.092 4.428 -12.676 0.320 381 + 383 3 -115.332 3.128 -12.676 0.391 382 + 384 3 -118.842 3.128 -14.996 0.369 383 + 385 3 -121.402 2.158 -17.416 0.365 384 + 386 3 -124.282 1.188 -17.416 0.320 385 + 387 3 -127.162 0.858 -20.256 0.320 386 + 388 3 -129.082 0.858 -20.256 0.320 387 + 389 3 -131.632 -1.732 -22.016 0.320 388 + 390 3 -134.192 -3.352 -24.736 0.320 389 + 391 3 -138.352 -4.322 -29.736 0.320 390 + 392 3 -140.592 -5.622 -32.256 0.320 391 + 393 3 -143.472 -6.912 -32.256 0.320 392 + 394 3 -146.342 -6.912 -35.356 0.320 393 + 395 3 -149.222 -7.562 -37.116 0.320 394 + 396 3 -152.102 -7.562 -37.936 0.320 395 + 397 3 -154.022 -7.562 -37.916 0.320 396 + 398 3 -157.222 -8.212 -39.196 0.320 397 + 399 3 -159.782 -10.152 -39.636 0.320 398 + 400 3 -162.332 -11.452 -43.116 0.320 399 + 401 3 -165.852 -13.392 -44.816 0.258 400 + 402 3 -168.092 -14.042 -47.096 0.243 401 + 403 3 -126.522 -2.052 -19.176 0.290 386 + 404 3 -130.042 -4.972 -19.176 0.294 403 + 405 3 -133.872 -5.942 -17.176 0.311 404 + 406 3 -137.072 -8.212 -17.176 0.332 405 + 407 3 -139.312 -9.832 -17.176 0.328 406 + 408 3 -142.512 -12.422 -17.596 0.309 407 + 409 3 -146.342 -12.752 -16.076 0.324 408 + 410 3 -149.862 -13.072 -17.236 0.257 409 + 411 3 -153.062 -12.752 -18.816 0.306 410 + 412 3 -155.942 -12.422 -20.536 0.249 411 + 413 3 -159.142 -11.772 -23.916 0.202 412 + 414 3 -114.692 7.338 -3.716 0.360 381 + 415 3 -119.482 8.638 1.824 0.379 414 + 416 3 -122.362 9.288 5.544 0.329 415 + 417 3 -125.882 10.578 7.024 0.338 416 + 418 3 -131.002 10.908 9.624 0.299 417 + 419 3 -133.552 13.498 10.304 0.326 418 + 420 3 -136.752 14.788 11.724 0.273 419 + 421 3 -138.992 16.088 13.144 0.307 420 + 422 3 -144.112 18.358 14.244 0.238 421 + 423 3 -80.230 -0.487 1.264 0.394 361 + 424 3 -82.470 -0.487 1.084 0.426 423 + 425 3 -85.350 -2.107 0.664 0.375 424 + 426 3 -88.084 -4.512 0.566 0.424 425 + 427 3 -90.408 -7.229 0.826 0.476 426 + 428 3 -93.590 -9.125 0.336 0.459 427 + 429 3 -96.470 -9.445 -0.304 0.366 428 + 430 3 -98.070 -11.395 -1.384 0.426 429 + 431 3 -99.350 -13.985 -1.924 0.403 430 + 432 3 -101.590 -16.895 -2.324 0.417 431 + 433 3 -104.470 -20.135 -3.244 0.436 432 + 434 3 -106.502 -23.151 -3.947 0.392 433 + 435 3 -108.055 -26.209 -4.515 0.399 434 + 436 3 -110.728 -29.518 -5.631 0.399 435 + 437 3 -113.278 -32.758 -9.111 0.378 436 + 438 3 -116.798 -35.348 -10.931 0.388 437 + 439 3 -119.038 -37.288 -14.551 0.393 438 + 440 3 -121.278 -39.558 -17.851 0.366 439 + 441 3 -124.798 -42.148 -20.151 0.394 440 + 442 3 -127.150 -46.524 -19.989 0.374 441 + 443 3 -130.234 -50.645 -19.025 0.356 442 + 444 3 -132.082 -54.729 -18.852 0.320 443 + 445 3 -134.952 -57.639 -18.852 0.320 444 + 446 3 -137.192 -60.229 -22.812 0.320 445 + 447 3 -140.072 -62.499 -24.332 0.320 446 + 448 3 -143.272 -64.449 -23.972 0.320 447 + 449 3 -145.192 -66.709 -23.972 0.320 448 + 450 3 -147.752 -69.309 -26.752 0.320 449 + 451 3 -149.982 -72.219 -25.852 0.320 450 + 452 3 -150.302 -77.079 -26.312 0.320 451 + 453 3 -152.222 -79.999 -27.112 0.320 452 + 454 3 -153.182 -83.889 -27.772 0.320 453 + 455 3 -156.062 -87.119 -28.132 0.320 454 + 456 3 -157.342 -91.009 -28.052 0.320 455 + 457 3 -158.822 -95.029 -26.292 0.320 456 + 458 3 -160.732 -99.569 -26.272 0.320 457 + 459 3 -162.012 -102.479 -25.372 0.320 458 + 460 3 -164.572 -106.049 -24.672 0.160 459 + 461 3 -167.132 -110.579 -24.672 0.160 460 + 462 3 -168.732 -116.409 -26.352 0.160 461 + 463 3 -53.420 3.520 -5.950 1.150 1 + 464 3 -60.270 4.330 -11.160 1.173 463 + 465 3 -63.460 5.300 -18.840 0.803 464 + 466 3 -65.700 6.920 -21.700 0.749 465 + 467 3 -68.260 7.570 -25.300 0.821 466 + 468 3 -69.860 9.190 -27.700 0.790 467 + 469 3 -71.700 10.878 -31.447 0.726 468 + 470 3 -73.919 12.850 -34.624 0.507 469 + 471 3 -75.519 14.470 -35.724 0.502 470 + 472 3 -76.799 14.790 -38.744 0.484 471 + 473 3 -76.799 16.410 -41.704 0.416 472 + 474 3 -77.759 18.680 -44.324 0.431 473 + 475 3 -77.119 19.650 -46.684 0.416 474 + 476 3 -78.079 23.210 -47.684 0.473 475 + 477 3 -79.989 23.860 -49.304 0.427 476 + 478 3 -80.629 25.160 -50.144 0.439 477 + 479 3 -81.269 26.780 -52.924 0.388 478 + 480 3 -82.229 27.750 -56.024 0.321 479 + 481 3 -83.509 30.670 -57.384 0.369 480 + 482 3 -86.389 31.960 -59.124 0.369 481 + 483 3 -87.669 33.260 -63.364 0.319 482 + 484 3 -88.309 35.530 -63.364 0.345 483 + 485 3 -89.589 35.850 -67.304 0.343 484 + 486 3 -90.229 36.820 -69.564 0.299 485 + 487 3 -90.869 36.170 -72.364 0.383 486 + 488 3 -90.549 35.530 -75.464 0.265 487 + 489 3 -53.100 14.220 -5.960 2.880 1 + 490 3 -53.420 16.810 -6.900 1.440 489 + 491 3 -53.230 19.230 -8.220 1.440 490 + 492 3 -52.270 27.010 -6.940 1.440 491 + 493 3 -51.950 34.130 -7.320 1.440 492 + 494 3 -49.710 42.560 -7.840 1.440 493 + 495 3 -48.110 43.210 -8.960 1.280 494 + 496 3 -50.670 45.470 -6.740 0.901 495 + 497 3 -50.990 48.390 -4.960 0.656 496 + 498 3 -53.550 49.680 -3.460 0.624 497 + 499 3 -54.830 51.300 -2.100 0.508 498 + 500 3 -56.110 51.950 0.580 0.585 499 + 501 3 -57.390 51.630 0.940 0.545 500 + 502 3 -58.990 51.950 2.100 0.469 501 + 503 3 -60.270 51.950 3.860 0.504 502 + 504 3 -61.550 52.600 5.760 0.500 503 + 505 3 -63.140 53.570 6.580 0.494 504 + 506 3 -65.560 55.340 7.980 0.492 505 + 507 3 -66.200 55.660 8.680 0.506 506 + 508 3 -67.480 58.580 8.840 0.483 507 + 509 3 -70.040 57.930 9.280 0.464 508 + 510 3 -71.000 59.550 10.100 0.482 509 + 511 3 -72.280 60.200 10.880 0.463 510 + 512 3 -73.560 60.200 12.080 0.415 511 + 513 3 -75.480 60.520 12.460 0.501 512 + 514 3 -75.800 62.140 12.920 0.453 513 + 515 3 -77.400 63.120 13.640 0.488 514 + 516 3 -79.000 64.090 13.740 0.429 515 + 517 3 -80.590 65.710 14.020 0.418 516 + 518 3 -82.510 65.710 14.320 0.418 517 + 519 3 -84.110 65.380 15.140 0.427 518 + 520 3 -84.430 64.740 15.940 0.406 519 + 521 3 -87.310 64.090 16.880 0.455 520 + 522 3 -88.270 62.790 17.580 0.446 521 + 523 3 -89.230 62.140 17.680 0.428 522 + 524 3 -91.790 62.140 16.400 0.450 523 + 525 3 -92.750 63.440 16.080 0.414 524 + 526 3 -94.350 63.120 17.600 0.458 525 + 527 3 -96.260 63.440 17.980 0.445 526 + 528 3 -98.180 64.740 20.160 0.459 527 + 529 3 -101.060 64.410 20.580 0.438 528 + 530 3 -102.340 65.710 20.960 0.444 529 + 531 3 -103.940 66.360 21.900 0.452 530 + 532 3 -106.180 67.000 23.000 0.424 531 + 533 3 -108.100 67.000 23.060 0.403 532 + 534 3 -109.700 68.950 23.780 0.341 533 + 535 3 -111.290 70.240 23.780 0.298 534 + 536 3 -113.210 71.210 24.380 0.323 535 + 537 3 -117.050 72.830 23.760 0.344 536 + 538 3 -121.210 75.750 24.220 0.265 537 + 539 3 -47.660 48.210 -11.680 0.960 495 + 540 3 -46.380 52.750 -13.220 0.960 539 + 541 3 -45.100 56.960 -15.660 0.960 540 + 542 3 -44.140 62.140 -17.120 0.960 541 + 543 3 -42.540 67.330 -18.720 0.960 542 + 544 3 -39.980 72.190 -17.060 0.960 543 + 545 3 -38.380 76.720 -16.440 0.960 544 + 546 3 -36.630 85.290 -18.760 0.800 545 + 547 3 -32.800 91.450 -21.020 0.800 546 + 548 3 -30.240 98.250 -20.720 0.800 547 + 549 3 -26.080 102.790 -20.760 0.800 548 + 550 3 -24.160 108.940 -22.280 0.800 549 + 551 3 -22.240 116.400 -25.640 0.800 550 + 552 3 -20.810 122.800 -27.140 0.800 551 + 553 3 -20.320 125.470 -27.140 0.800 552 + 554 3 -17.010 132.750 -29.880 0.800 553 + 555 3 -16.830 133.990 -29.880 0.800 554 + 556 3 -16.010 135.760 -29.880 0.649 555 + 557 3 -13.810 139.880 -30.160 0.800 556 + 558 3 -13.170 146.680 -30.300 0.800 557 + 559 3 -9.650 152.830 -30.580 0.800 558 + 560 3 -5.180 158.340 -32.160 0.800 559 + 561 3 -2.940 164.500 -33.280 0.800 560 + 562 3 1.220 169.360 -34.460 0.800 561 + 563 3 3.460 172.270 -34.800 0.800 562 + 564 3 7.810 176.000 -35.220 0.800 563 + 565 3 11.330 182.160 -35.240 0.800 564 + 566 3 14.200 186.050 -35.940 0.800 565 + 567 3 17.870 190.390 -37.400 0.800 566 + 568 3 19.790 191.690 -40.000 0.800 567 + 569 3 21.390 193.310 -40.100 0.800 568 + 570 3 22.020 198.530 -41.220 0.800 569 + 571 3 23.620 204.030 -42.740 0.800 570 + 572 3 25.540 211.490 -44.760 0.800 571 + 573 3 29.060 218.290 -46.280 0.800 572 + 574 3 32.900 223.800 -47.980 0.800 573 + 575 3 34.810 233.840 -49.140 0.800 574 + 576 3 37.370 241.290 -49.140 0.800 575 + 577 3 42.810 249.390 -47.980 0.800 576 + 578 3 45.690 258.140 -48.560 0.800 577 + 579 3 47.290 265.910 -47.180 0.800 578 + 580 3 47.930 270.450 -47.180 0.800 579 + 581 3 52.080 276.280 -47.180 0.800 580 + 582 3 51.440 280.490 -47.180 0.800 581 + 583 3 53.360 288.270 -46.660 0.800 582 + 584 3 56.240 294.750 -46.660 0.800 583 + 585 3 57.670 300.750 -46.660 0.800 584 + 586 3 60.550 304.320 -46.660 0.800 585 + 587 3 60.550 309.500 -46.660 0.800 586 + 588 3 61.190 315.650 -48.600 0.800 587 + 589 3 62.470 320.840 -48.620 0.800 588 + 590 3 65.030 327.640 -50.460 0.800 589 + 591 3 66.310 333.470 -51.560 0.800 590 + 592 3 66.310 339.310 -48.780 0.800 591 + 593 3 68.550 344.160 -53.580 0.800 592 + 594 3 71.250 348.210 -48.960 0.800 593 + 595 3 73.490 352.420 -52.120 0.800 594 + 596 3 77.650 356.630 -49.480 0.800 595 + 597 3 78.930 363.430 -46.460 0.800 596 + 598 3 82.120 370.890 -45.560 0.800 597 + 599 3 85.640 377.040 -44.400 0.800 598 + 600 3 88.520 383.850 -42.820 0.800 599 + 601 3 91.890 385.640 -39.600 0.640 600 + 602 3 94.450 388.560 -39.280 0.640 601 + 603 3 98.290 391.470 -39.280 0.640 602 + 604 3 101.490 393.740 -37.060 0.640 603 + 605 3 104.680 397.300 -37.060 0.640 604 + 606 3 107.350 400.340 -41.900 0.480 605 + 607 3 109.910 404.220 -42.880 0.480 606 + 608 3 112.470 406.490 -44.420 0.480 607 + 609 3 114.390 410.060 -46.560 0.480 608 + 610 3 115.670 413.290 -47.580 0.480 609 + 611 3 120.460 415.240 -49.380 0.480 610 + 612 3 122.700 417.510 -50.740 0.480 611 + 613 3 126.540 417.830 -52.320 0.480 612 + 614 3 130.060 419.770 -54.260 0.480 613 + 615 3 131.980 422.370 -55.180 0.480 614 + 616 3 135.170 423.340 -55.780 0.480 615 + 617 3 139.010 423.340 -56.820 0.480 616 + 618 3 143.490 423.660 -58.320 0.480 617 + 619 3 148.960 422.650 -58.480 0.320 618 + 620 3 150.240 421.680 -58.480 0.320 619 + 621 3 152.800 422.970 -57.840 0.320 620 + 622 3 155.040 419.410 -56.260 0.320 621 + 623 3 158.240 419.410 -56.260 0.320 622 + 624 3 160.150 420.380 -56.260 0.320 623 + 625 3 163.030 420.380 -59.260 0.320 624 + 626 3 168.470 422.330 -58.880 0.320 625 + 627 3 173.590 423.300 -57.940 0.320 626 + 628 3 175.500 425.240 -56.220 0.320 627 + 629 3 176.780 428.480 -55.580 0.320 628 + 630 3 178.700 428.160 -52.340 0.320 629 + 631 3 180.940 428.480 -50.920 0.320 630 + 632 3 182.860 429.780 -49.480 0.320 631 + 633 3 184.780 428.810 -46.620 0.320 632 + 634 3 186.380 428.480 -45.120 0.320 633 + 635 3 189.250 428.160 -45.120 0.320 634 + 636 3 191.810 428.810 -44.060 0.320 635 + 637 3 193.730 428.160 -42.600 0.320 636 + 638 3 195.970 430.100 -41.440 0.320 637 + 639 3 197.250 432.370 -41.440 0.320 638 + 640 3 202.370 433.020 -41.580 0.320 639 + 641 3 205.240 432.040 -40.520 0.320 640 + 642 3 207.480 434.310 -40.520 0.320 641 + 643 3 210.040 434.960 -40.420 0.320 642 + 644 3 212.920 435.280 -40.420 0.320 643 + 645 3 214.840 434.640 -40.420 0.320 644 + 646 3 218.670 436.260 -39.740 0.320 645 + 647 3 221.230 437.880 -39.740 0.320 646 + 648 3 224.750 437.880 -38.100 0.320 647 + 649 3 226.350 439.820 -38.100 0.320 648 + 650 3 229.870 439.820 -36.480 0.320 649 + 651 3 232.750 440.790 -36.480 0.320 650 + 652 3 234.980 444.680 -36.460 0.320 651 + 653 3 237.860 446.300 -36.440 0.320 652 + 654 3 240.420 447.270 -37.640 0.320 653 + 655 3 244.900 447.270 -37.640 0.320 654 + 656 3 246.500 449.860 -37.800 0.208 655 + 657 3 247.780 449.540 -37.800 0.199 656 + 658 3 250.010 452.780 -36.600 0.199 657 + 659 3 251.930 453.430 -36.600 0.233 658 + 660 3 254.490 453.750 -35.440 0.219 659 + 661 3 256.090 455.370 -35.440 0.237 660 + 662 3 257.690 457.320 -35.440 0.231 661 + 663 3 260.570 459.580 -35.900 0.210 662 + 664 3 263.760 460.880 -35.900 0.227 663 + 665 3 267.600 462.820 -36.980 0.204 664 + 666 3 270.160 464.120 -36.980 0.196 665 + 667 3 272.080 466.060 -36.980 0.150 666 + 668 3 275.280 469.300 -38.320 0.160 667 + 669 3 279.430 468.650 -38.320 0.160 668 + 670 3 280.870 471.410 -38.320 0.160 669 + 671 3 284.070 473.680 -38.340 0.160 670 + 672 3 245.220 445.650 -36.280 0.320 655 + 673 3 247.460 443.060 -35.380 0.160 672 + 674 3 147.010 424.310 -58.180 0.480 618 + 675 3 148.920 425.280 -66.420 0.480 674 + 676 3 152.120 426.900 -68.100 0.480 675 + 677 3 156.600 427.550 -70.460 0.480 676 + 678 3 159.480 427.230 -70.460 0.320 677 + 679 3 162.990 429.170 -72.920 0.320 678 + 680 3 165.750 430.290 -74.740 0.320 679 + 681 3 170.230 430.940 -76.100 0.320 680 + 682 3 171.510 433.860 -76.780 0.320 681 + 683 3 174.070 433.530 -79.040 0.320 682 + 684 3 177.580 433.210 -79.380 0.320 683 + 685 3 179.500 435.480 -80.700 0.320 684 + 686 3 183.020 434.830 -77.780 0.320 685 + 687 3 186.860 436.120 -79.180 0.320 686 + 688 3 191.010 435.480 -80.600 0.320 687 + 689 3 192.290 437.100 -81.220 0.320 688 + 690 3 199.970 439.360 -81.220 0.320 689 + 691 3 203.170 439.040 -81.220 0.320 690 + 692 3 205.090 442.280 -81.220 0.320 691 + 693 3 209.880 442.930 -83.320 0.320 692 + 694 3 212.440 445.840 -86.080 0.320 693 + 695 3 215.640 447.140 -86.080 0.320 694 + 696 3 220.440 447.790 -86.820 0.320 695 + 697 3 223.470 449.860 -87.220 0.320 696 + 698 3 227.310 451.810 -88.580 0.320 697 + 699 3 233.070 452.130 -92.020 0.320 698 + 700 3 236.260 454.400 -96.560 0.320 699 + 701 3 241.700 456.670 -97.960 0.320 700 + 702 3 244.580 453.750 -99.440 0.320 701 + 703 3 248.100 454.400 -101.440 0.320 702 + 704 3 250.650 454.080 -101.720 0.320 703 + 705 3 253.210 456.020 -102.740 0.320 704 + 706 3 256.730 457.640 -104.560 0.320 705 + 707 3 260.250 457.960 -106.300 0.320 706 + 708 3 263.440 457.960 -105.780 0.160 707 + 709 3 265.360 456.340 -105.780 0.160 708 + 710 3 268.880 456.340 -107.480 0.160 709 + 711 3 270.480 457.960 -107.480 0.160 710 + 712 3 262.810 460.230 -106.980 0.160 707 + 713 3 263.760 463.790 -107.160 0.160 712 + 714 3 105.320 400.220 -38.400 0.320 605 + 715 3 107.560 401.840 -35.900 0.320 714 + 716 3 109.160 402.160 -36.120 0.320 715 + 717 3 110.760 401.190 -34.580 0.320 716 + 718 3 111.080 402.160 -25.200 0.320 717 + 719 3 111.720 401.840 -22.760 0.320 718 + 720 3 111.080 399.250 -19.640 0.320 719 + 721 3 112.360 398.280 -19.640 0.320 720 + 722 3 113.640 398.280 -16.340 0.320 721 + 723 3 115.560 398.280 -13.060 0.320 722 + 724 3 117.790 397.630 -13.060 0.320 723 + 725 3 119.390 398.280 -13.060 0.320 724 + 726 3 123.230 397.950 -11.840 0.320 725 + 727 3 128.030 398.920 -10.860 0.320 726 + 728 3 130.910 397.950 -9.380 0.320 727 + 729 3 133.140 397.300 -8.260 0.320 728 + 730 3 137.300 396.980 -8.260 0.320 729 + 731 3 140.820 397.630 -8.900 0.320 730 + 732 3 144.020 399.570 -8.920 0.320 731 + 733 3 149.950 402.640 -7.540 0.320 732 + 734 3 153.150 406.200 -5.060 0.320 733 + 735 3 158.260 408.470 -3.240 0.320 734 + 736 3 161.140 410.740 -1.440 0.320 735 + 737 3 164.340 412.680 0.000 0.320 736 + 738 3 166.260 413.000 0.000 0.320 737 + 739 3 170.420 414.620 -2.720 0.320 738 + 740 3 173.610 416.570 -4.680 0.320 739 + 741 3 177.450 418.830 -5.660 0.342 740 + 742 3 180.970 420.780 -6.100 0.340 741 + 743 3 183.530 423.690 -6.680 0.281 742 + 744 3 186.730 426.930 -8.420 0.311 743 + 745 3 188.320 432.120 -10.040 0.305 744 + 746 3 191.520 434.710 -11.160 0.280 745 + 747 3 193.440 437.300 -12.040 0.255 746 + 748 3 195.040 438.920 -13.060 0.247 747 + 749 3 196.960 441.840 -14.020 0.248 748 + 750 3 199.840 442.810 -15.360 0.254 749 + 751 3 201.910 443.990 -15.360 0.252 750 + 752 3 204.150 448.200 -15.940 0.256 751 + 753 3 205.750 449.820 -16.460 0.237 752 + 754 3 208.630 451.120 -16.820 0.222 753 + 755 3 88.520 390.000 -42.740 0.640 600 + 756 3 89.480 395.510 -40.160 0.640 755 + 757 3 89.160 401.340 -39.280 0.640 756 + 758 3 87.060 405.410 -37.120 0.480 757 + 759 3 83.860 409.940 -36.580 0.480 758 + 760 3 82.900 412.540 -36.580 0.480 759 + 761 3 81.300 415.780 -35.340 0.480 760 + 762 3 80.660 419.010 -35.340 0.480 761 + 763 3 79.060 420.310 -35.340 0.480 762 + 764 3 77.780 423.230 -35.140 0.480 763 + 765 3 74.900 424.200 -35.140 0.480 764 + 766 3 72.340 425.490 -34.460 0.480 765 + 767 3 69.150 427.440 -33.300 0.480 766 + 768 3 67.870 429.710 -32.800 0.480 767 + 769 3 66.270 432.300 -32.800 0.320 768 + 770 3 65.950 436.830 -33.960 0.320 769 + 771 3 64.350 441.690 -33.960 0.314 770 + 772 3 62.750 444.930 -34.880 0.272 771 + 773 3 61.790 449.470 -34.900 0.310 772 + 774 3 60.510 453.360 -34.920 0.280 773 + 775 3 59.550 456.600 -36.320 0.299 774 + 776 3 57.770 460.610 -36.320 0.283 775 + 777 3 57.770 464.500 -37.160 0.288 776 + 778 3 57.140 468.060 -37.160 0.335 777 + 779 3 59.050 470.980 -37.160 0.297 778 + 780 3 58.090 474.210 -36.200 0.293 779 + 781 3 58.410 478.750 -36.200 0.245 780 + 782 3 72.380 346.430 -56.460 0.640 593 + 783 3 74.620 351.290 -57.900 0.640 782 + 784 3 78.630 353.360 -57.900 0.640 783 + 785 3 81.830 358.870 -60.800 0.640 784 + 786 3 83.430 363.080 -63.000 0.640 785 + 787 3 83.590 367.110 -59.640 0.640 786 + 788 3 81.240 370.140 -65.340 0.480 787 + 789 3 81.240 374.020 -65.340 0.480 788 + 790 3 79.960 376.620 -65.340 0.480 789 + 791 3 80.600 380.500 -65.340 0.480 790 + 792 3 78.680 381.800 -65.340 0.480 791 + 793 3 78.680 386.660 -63.980 0.480 792 + 794 3 77.720 389.900 -63.760 0.480 793 + 795 3 76.440 391.520 -63.620 0.480 794 + 796 3 77.080 393.460 -63.420 0.480 795 + 797 3 75.480 397.030 -63.200 0.480 796 + 798 3 76.760 400.270 -64.320 0.480 797 + 799 3 75.160 402.860 -67.320 0.480 798 + 800 3 73.560 407.070 -69.660 0.480 799 + 801 3 72.280 409.660 -69.680 0.480 800 + 802 3 72.920 414.200 -73.540 0.480 801 + 803 3 72.130 417.260 -72.220 0.320 802 + 804 3 71.170 419.850 -73.740 0.320 803 + 805 3 72.130 424.060 -74.580 0.320 804 + 806 3 70.850 426.980 -75.700 0.320 805 + 807 3 70.210 431.520 -75.860 0.320 806 + 808 3 71.170 435.400 -77.020 0.320 807 + 809 3 71.810 437.350 -77.020 0.320 808 + 810 3 71.490 440.590 -76.180 0.320 809 + 811 3 70.210 442.210 -76.180 0.320 810 + 812 3 70.850 446.420 -76.180 0.320 811 + 813 3 69.260 449.010 -76.180 0.320 812 + 814 3 69.570 450.950 -76.180 0.320 813 + 815 3 68.940 452.900 -76.180 0.320 814 + 816 3 68.940 456.460 -77.780 0.320 815 + 817 3 71.170 459.700 -75.780 0.320 816 + 818 3 69.890 461.000 -75.780 0.320 817 + 819 3 70.530 464.560 -75.940 0.320 818 + 820 3 68.940 469.100 -77.600 0.320 819 + 821 3 69.570 471.360 -76.640 0.320 820 + 822 3 67.020 475.580 -76.640 0.320 821 + 823 3 66.700 479.790 -79.160 0.320 822 + 824 3 86.470 368.080 -61.360 0.640 787 + 825 3 86.790 371.320 -68.280 0.480 824 + 826 3 89.350 374.560 -70.640 0.480 825 + 827 3 90.630 377.800 -72.320 0.480 826 + 828 3 91.590 380.710 -72.320 0.480 827 + 829 3 94.150 379.740 -75.240 0.480 828 + 830 3 95.110 382.980 -77.200 0.480 829 + 831 3 95.430 384.930 -79.460 0.480 830 + 832 3 97.660 387.190 -80.800 0.480 831 + 833 3 97.980 389.790 -80.820 0.480 832 + 834 3 96.700 391.730 -82.600 0.480 833 + 835 3 97.020 394.650 -83.800 0.480 834 + 836 3 98.620 396.910 -85.560 0.480 835 + 837 3 100.540 398.860 -86.380 0.480 836 + 838 3 100.220 402.420 -88.560 0.480 837 + 839 3 101.820 405.660 -90.300 0.480 838 + 840 3 105.020 407.610 -92.000 0.480 839 + 841 3 105.660 410.850 -93.420 0.480 840 + 842 3 105.340 415.060 -95.120 0.480 841 + 843 3 106.300 417.320 -98.040 0.480 842 + 844 3 109.180 420.560 -100.260 0.480 843 + 845 3 111.210 424.280 -101.600 0.480 844 + 846 3 112.490 425.580 -107.860 0.480 845 + 847 3 115.050 426.550 -108.100 0.480 846 + 848 3 118.890 428.170 -109.040 0.480 847 + 849 3 122.720 430.120 -109.900 0.480 848 + 850 3 123.360 431.740 -112.580 0.480 849 + 851 3 123.360 433.360 -112.580 0.480 850 + 852 3 123.490 436.760 -112.540 0.480 851 + 853 3 124.450 441.300 -114.120 0.480 852 + 854 3 125.410 445.180 -114.120 0.480 853 + 855 3 124.770 450.370 -116.140 0.480 854 + 856 3 124.130 453.610 -120.420 0.320 855 + 857 3 123.170 456.200 -126.760 0.320 856 + 858 3 120.300 459.760 -126.760 0.320 857 + 859 3 120.620 464.300 -130.660 0.320 858 + 860 3 119.980 468.190 -130.660 0.320 859 + 861 3 119.980 472.070 -132.560 0.320 860 + 862 3 127.010 452.960 -116.800 0.295 855 + 863 3 129.570 453.930 -116.800 0.283 862 + 864 3 129.890 457.490 -115.700 0.278 863 + 865 3 132.130 461.710 -113.840 0.256 864 + 866 3 133.730 463.330 -113.840 0.278 865 + 867 3 135.010 466.890 -112.420 0.266 866 + 868 3 136.600 469.160 -112.420 0.244 867 + 869 3 141.080 470.450 -113.320 0.242 868 + 870 3 142.040 473.690 -114.360 0.278 869 + 871 3 142.680 478.230 -113.420 0.289 870 + 872 3 142.360 480.820 -113.420 0.278 871 + 873 3 126.240 432.710 -114.120 0.480 851 + 874 3 128.800 431.410 -119.440 0.462 873 + 875 3 131.360 431.410 -122.220 0.418 874 + 876 3 133.920 430.440 -122.220 0.444 875 + 877 3 136.800 432.380 -122.220 0.448 876 + 878 3 139.030 434.000 -124.580 0.456 877 + 879 3 141.590 433.680 -125.360 0.437 878 + 880 3 144.790 436.920 -127.060 0.471 879 + 881 3 147.030 438.860 -128.480 0.494 880 + 882 3 149.910 437.570 -129.500 0.476 881 + 883 3 152.460 437.570 -131.100 0.430 882 + 884 3 155.020 438.540 -132.700 0.452 883 + 885 3 157.900 439.190 -135.540 0.512 884 + 886 3 160.780 438.860 -137.220 0.452 885 + 887 3 163.930 440.200 -132.860 0.441 886 + 888 3 168.730 441.170 -133.200 0.350 887 + 889 3 172.570 442.790 -133.960 0.301 888 + 890 3 176.720 445.380 -135.340 0.287 889 + 891 3 178.960 450.240 -136.860 0.286 890 + 892 3 182.480 451.220 -138.060 0.269 891 + 893 3 184.080 453.480 -138.900 0.243 892 + 894 3 186.960 455.750 -141.740 0.253 893 + 895 3 187.920 458.020 -142.560 0.211 894 + 896 3 190.800 459.640 -145.680 0.247 895 + 897 3 192.710 465.150 -145.700 0.223 896 + 898 3 195.590 469.360 -144.220 0.262 897 + 899 3 198.470 473.250 -145.400 0.277 898 + 900 3 201.670 473.250 -146.880 0.263 899 + 901 3 204.550 476.810 -146.660 0.228 900 + 902 3 209.020 480.370 -146.660 0.277 901 + 903 3 161.740 436.920 -143.520 0.504 886 + 904 3 163.660 437.890 -144.800 0.422 903 + 905 3 165.260 438.540 -146.640 0.370 904 + 906 3 167.310 438.400 -149.760 0.334 905 + 907 3 168.910 439.370 -149.760 0.339 906 + 908 3 170.510 437.750 -151.020 0.335 907 + 909 3 173.710 437.750 -154.660 0.351 908 + 910 3 177.230 440.990 -157.460 0.374 909 + 911 3 179.140 442.280 -158.620 0.378 910 + 912 3 183.300 443.900 -159.300 0.391 911 + 913 3 187.140 445.520 -160.800 0.319 912 + 914 3 190.020 446.500 -156.260 0.342 913 + 915 3 193.850 448.120 -160.220 0.339 914 + 916 3 198.650 448.120 -163.180 0.366 915 + 917 3 203.130 448.440 -162.660 0.276 916 + 918 3 206.650 449.740 -166.040 0.296 917 + 919 3 210.160 451.030 -171.540 0.316 918 + 920 3 213.360 451.360 -168.860 0.331 919 + 921 3 219.120 451.680 -173.520 0.293 920 + 922 3 222.960 449.740 -176.200 0.283 921 + 923 3 225.010 448.300 -180.860 0.282 922 + 924 3 226.610 450.240 -181.740 0.301 923 + 925 3 230.770 448.950 -185.260 0.219 924 + 926 3 234.610 447.980 -189.360 0.187 925 + 927 3 88.220 364.700 -61.000 0.452 786 + 928 3 92.060 366.970 -61.540 0.470 927 + 929 3 93.660 369.240 -62.220 0.479 928 + 930 3 96.220 371.830 -63.180 0.439 929 + 931 3 100.060 372.150 -67.120 0.420 930 + 932 3 104.530 374.750 -70.740 0.451 931 + 933 3 108.050 376.040 -73.280 0.406 932 + 934 3 110.610 378.960 -73.280 0.476 933 + 935 3 116.050 381.870 -73.280 0.436 934 + 936 3 119.880 384.140 -74.080 0.443 935 + 937 3 125.320 386.090 -75.640 0.417 936 + 938 3 130.120 389.330 -78.500 0.352 937 + 939 3 134.130 391.790 -78.580 0.362 938 + 940 3 139.890 393.410 -82.180 0.391 939 + 941 3 146.600 396.320 -89.900 0.424 940 + 942 3 151.400 398.270 -92.560 0.429 941 + 943 3 156.200 400.210 -99.540 0.413 942 + 944 3 161.310 401.510 -104.000 0.440 943 + 945 3 165.470 399.890 -106.560 0.429 944 + 946 3 168.670 403.130 -107.620 0.445 945 + 947 3 173.460 401.510 -110.540 0.441 946 + 948 3 178.260 405.720 -114.120 0.412 947 + 949 3 181.780 407.660 -118.460 0.462 948 + 950 3 185.300 411.230 -120.880 0.456 949 + 951 3 187.860 414.790 -124.580 0.477 950 + 952 3 193.110 417.840 -125.280 0.416 951 + 953 3 196.630 418.480 -125.280 0.429 952 + 954 3 198.870 420.750 -126.660 0.440 953 + 955 3 198.030 423.170 -124.340 0.463 954 + 956 3 199.620 426.740 -124.340 0.405 955 + 957 3 202.820 429.650 -127.700 0.442 956 + 958 3 204.420 432.240 -128.800 0.484 957 + 959 3 205.060 435.480 -128.800 0.473 958 + 960 3 207.940 436.130 -129.960 0.432 959 + 961 3 207.940 439.370 -129.960 0.444 960 + 962 3 212.740 440.670 -132.140 0.455 961 + 963 3 214.010 444.230 -133.720 0.419 962 + 964 3 214.330 448.440 -130.060 0.469 963 + 965 3 215.290 452.330 -128.840 0.425 964 + 966 3 219.130 453.950 -127.300 0.435 965 + 967 3 219.770 458.160 -128.140 0.404 966 + 968 3 223.610 459.780 -128.500 0.420 967 + 969 3 226.170 461.730 -129.820 0.425 968 + 970 3 229.360 463.670 -128.000 0.402 969 + 971 3 229.680 467.880 -128.000 0.330 970 + 972 3 233.200 470.470 -128.240 0.363 971 + 973 3 237.040 475.330 -130.380 0.236 972 + 974 3 204.620 423.020 -129.440 0.320 954 + 975 3 210.700 425.610 -130.440 0.320 974 + 976 3 215.500 425.290 -131.180 0.320 975 + 977 3 220.930 427.230 -132.260 0.320 976 + 978 3 225.730 427.560 -134.020 0.320 977 + 979 3 229.890 425.290 -135.760 0.320 978 + 980 3 233.720 426.580 -137.100 0.320 979 + 981 3 237.240 429.180 -134.880 0.320 980 + 982 3 242.040 430.150 -135.940 0.320 981 + 983 3 245.870 431.770 -137.240 0.320 982 + 984 3 250.470 435.480 -138.860 0.320 983 + 985 3 255.590 435.480 -139.740 0.320 984 + 986 3 259.100 437.430 -140.580 0.320 985 + 987 3 262.940 438.720 -141.140 0.320 986 + 988 3 268.060 440.020 -141.140 0.320 987 + 989 3 271.580 440.990 -144.060 0.320 988 + 990 3 276.370 440.340 -145.480 0.320 989 + 991 3 23.310 191.040 -43.660 0.540 569 + 992 3 25.870 190.390 -44.460 0.493 991 + 993 3 29.070 189.740 -44.460 0.551 992 + 994 3 32.580 188.450 -44.880 0.462 993 + 995 3 36.420 189.420 -45.340 0.477 994 + 996 3 39.940 190.390 -45.340 0.491 995 + 997 3 41.100 193.810 -41.680 0.480 996 + 998 3 41.730 197.690 -40.100 0.480 997 + 999 3 44.610 201.910 -39.840 0.480 998 + 1000 3 45.570 207.410 -37.080 0.480 999 + 1001 3 46.530 213.570 -36.720 0.480 1000 + 1002 3 48.130 218.110 -36.720 0.480 1001 + 1003 3 52.290 222.640 -37.600 0.480 1002 + 1004 3 53.890 227.180 -36.740 0.480 1003 + 1005 3 57.720 230.090 -36.380 0.480 1004 + 1006 3 60.600 231.060 -36.280 0.480 1005 + 1007 3 62.520 235.920 -36.220 0.464 1006 + 1008 3 65.080 240.140 -35.020 0.448 1007 + 1009 3 67.110 244.860 -38.280 0.508 1008 + 1010 3 70.310 249.720 -37.260 0.411 1009 + 1011 3 73.510 254.900 -37.260 0.386 1010 + 1012 3 76.070 259.430 -36.780 0.334 1011 + 1013 3 77.990 263.970 -34.260 0.368 1012 + 1014 3 81.500 264.940 -33.340 0.309 1013 + 1015 3 85.660 265.910 -31.440 0.320 1014 + 1016 3 90.140 266.890 -29.460 0.276 1015 + 1017 3 42.820 188.770 -44.240 0.456 996 + 1018 3 46.970 187.480 -46.440 0.485 1017 + 1019 3 49.210 184.880 -46.720 0.422 1018 + 1020 3 52.410 184.560 -47.700 0.453 1019 + 1021 3 55.610 182.940 -48.660 0.403 1020 + 1022 3 59.130 184.560 -49.940 0.365 1021 + 1023 3 61.360 182.940 -51.140 0.415 1022 + 1024 3 65.200 181.640 -52.160 0.392 1023 + 1025 3 70.000 180.670 -53.640 0.451 1024 + 1026 3 73.200 181.640 -54.800 0.385 1025 + 1027 3 75.950 181.170 -55.820 0.428 1026 + 1028 3 79.470 179.550 -60.080 0.418 1027 + 1029 3 81.390 177.280 -60.600 0.358 1028 + 1030 3 82.670 174.370 -61.620 0.355 1029 + 1031 3 86.500 174.040 -63.660 0.272 1030 + 1032 3 88.740 173.070 -64.940 0.233 1031 + 1033 3 1.860 175.840 -37.060 0.480 563 + 1034 3 -0.700 178.750 -37.520 0.480 1033 + 1035 3 0.130 181.840 -36.740 0.480 1034 + 1036 3 -1.470 184.430 -36.660 0.320 1035 + 1037 3 -3.380 185.400 -36.080 0.320 1036 + 1038 3 -3.380 188.320 -36.700 0.320 1037 + 1039 3 -4.340 190.260 -37.520 0.320 1038 + 1040 3 -5.940 190.910 -38.620 0.320 1039 + 1041 3 -6.580 193.820 -38.720 0.320 1040 + 1042 3 -8.500 196.090 -38.720 0.320 1041 + 1043 3 -9.460 199.010 -39.300 0.320 1042 + 1044 3 -10.100 202.890 -40.000 0.320 1043 + 1045 3 -10.740 204.840 -40.500 0.320 1044 + 1046 3 -11.700 208.080 -41.040 0.320 1045 + 1047 3 -12.660 210.670 -41.400 0.320 1046 + 1048 3 -13.620 212.290 -42.260 0.320 1047 + 1049 3 -14.580 213.910 -43.420 0.320 1048 + 1050 3 -17.140 215.210 -44.860 0.320 1049 + 1051 3 -18.090 218.120 -46.200 0.320 1050 + 1052 3 -18.090 219.420 -48.840 0.320 1051 + 1053 3 -20.330 220.390 -48.840 0.320 1052 + 1054 3 -20.330 222.980 -52.500 0.320 1053 + 1055 3 -22.890 224.280 -52.500 0.320 1054 + 1056 3 -4.220 180.370 -38.100 0.320 1034 + 1057 3 -6.660 181.170 -38.100 0.320 1056 + 1058 3 -7.620 183.110 -37.400 0.320 1057 + 1059 3 -10.820 183.760 -37.100 0.320 1058 + 1060 3 -13.060 183.760 -36.140 0.320 1059 + 1061 3 -15.620 187.320 -34.400 0.320 1060 + 1062 3 -17.540 191.210 -34.060 0.320 1061 + 1063 3 -19.460 195.420 -33.400 0.320 1062 + 1064 3 -20.740 198.340 -31.680 0.320 1063 + 1065 3 -19.460 200.610 -31.680 0.320 1064 + 1066 3 -21.370 202.230 -29.920 0.320 1065 + 1067 3 -22.330 203.850 -26.440 0.320 1066 + 1068 3 -24.250 205.790 -23.340 0.320 1067 + 1069 3 -24.250 207.730 -23.340 0.320 1068 + 1070 3 -25.530 212.590 -21.620 0.320 1069 + 1071 3 -23.610 214.860 -20.240 0.320 1070 + 1072 3 -24.250 217.450 -19.200 0.320 1071 + 1073 3 -24.250 220.050 -19.200 0.320 1072 + 1074 3 -25.530 221.670 -18.440 0.320 1073 + 1075 3 -24.890 225.550 -17.360 0.320 1074 + 1076 3 -26.810 228.150 -16.080 0.320 1075 + 1077 3 -27.770 230.740 -14.200 0.320 1076 + 1078 3 -29.050 232.030 -12.020 0.320 1077 + 1079 3 -30.010 231.380 -12.020 0.320 1078 + 1080 3 -31.610 230.740 -10.680 0.320 1079 + 1081 3 -32.570 232.030 -7.840 0.320 1080 + 1082 3 -32.890 234.950 -6.180 0.320 1081 + 1083 3 -34.630 236.420 -4.420 0.320 1082 + 1084 3 -34.630 238.040 -1.700 0.320 1083 + 1085 3 -36.220 239.980 -0.280 0.320 1084 + 1086 3 -37.180 242.250 0.360 0.320 1085 + 1087 3 -38.460 244.190 2.240 0.320 1086 + 1088 3 -38.140 246.780 2.260 0.320 1087 + 1089 3 -41.020 247.110 3.840 0.320 1088 + 1090 3 -42.300 247.110 7.540 0.320 1089 + 1091 3 -43.260 250.350 7.540 0.320 1090 + 1092 3 -46.460 250.670 9.920 0.320 1091 + 1093 3 -46.780 252.610 12.600 0.320 1092 + 1094 3 -12.810 134.470 -39.180 0.559 556 + 1095 3 -9.300 134.140 -40.200 0.546 1094 + 1096 3 -6.740 131.880 -40.840 0.604 1095 + 1097 3 -5.140 129.930 -42.480 0.635 1096 + 1098 3 -2.900 129.280 -46.680 0.631 1097 + 1099 3 -1.300 127.020 -48.520 0.540 1098 + 1100 3 -1.620 124.420 -50.620 0.628 1099 + 1101 3 0.620 123.130 -54.380 0.480 1100 + 1102 3 1.900 124.420 -56.680 0.480 1101 + 1103 3 2.860 123.450 -59.260 0.480 1102 + 1104 3 4.910 122.320 -60.980 0.480 1103 + 1105 3 6.510 120.700 -62.880 0.503 1104 + 1106 3 9.070 119.410 -66.060 0.449 1105 + 1107 3 11.310 119.410 -70.600 0.336 1106 + 1108 3 12.910 115.840 -74.860 0.376 1107 + 1109 3 14.190 114.220 -76.820 0.372 1108 + 1110 3 17.700 113.580 -76.820 0.307 1109 + 1111 3 20.900 116.820 -80.660 0.315 1110 + 1112 3 21.220 120.060 -76.500 0.322 1111 + 1113 3 25.380 121.030 -75.080 0.335 1112 + 1114 3 -18.110 132.040 -39.000 0.585 555 + 1115 3 -20.350 130.750 -39.920 0.506 1114 + 1116 3 -20.990 129.130 -41.260 0.463 1115 + 1117 3 -22.270 127.510 -42.440 0.467 1116 + 1118 3 -22.590 125.890 -43.220 0.467 1117 + 1119 3 -23.870 123.620 -44.580 0.436 1118 + 1120 3 -25.790 123.300 -47.880 0.459 1119 + 1121 3 -28.030 122.650 -50.140 0.414 1120 + 1122 3 -28.990 121.030 -52.700 0.393 1121 + 1123 3 -29.620 118.760 -55.280 0.421 1122 + 1124 3 -32.180 118.440 -57.700 0.428 1123 + 1125 3 -32.500 113.900 -59.660 0.444 1124 + 1126 3 -33.460 110.010 -59.700 0.456 1125 + 1127 3 -35.060 109.360 -59.920 0.419 1126 + 1128 3 -35.700 110.340 -64.700 0.374 1127 + 1129 3 -36.660 108.390 -68.000 0.444 1128 + 1130 3 -37.300 106.120 -68.000 0.412 1129 + 1131 3 -38.260 103.210 -74.760 0.320 1130 + 1132 3 -18.250 122.480 -21.480 0.483 552 + 1133 3 -18.250 125.400 -21.480 0.496 1132 + 1134 3 -20.490 127.340 -20.300 0.437 1133 + 1135 3 -21.450 128.960 -19.020 0.401 1134 + 1136 3 -24.010 130.260 -17.520 0.385 1135 + 1137 3 -24.970 131.880 -14.940 0.378 1136 + 1138 3 -26.240 135.440 -13.340 0.390 1137 + 1139 3 -27.520 137.380 -11.780 0.360 1138 + 1140 3 -29.760 139.650 -9.320 0.375 1139 + 1141 3 -31.040 140.300 -8.540 0.383 1140 + 1142 3 -31.360 143.860 -7.540 0.320 1141 + 1143 3 -33.920 145.810 -6.680 0.297 1142 + 1144 3 -36.160 146.780 -7.200 0.287 1143 + 1145 3 -38.080 151.310 -5.160 0.314 1144 + 1146 3 -40.000 152.290 -4.780 0.283 1145 + 1147 3 -40.630 155.530 -4.260 0.286 1146 + 1148 3 -43.190 157.150 -3.920 0.261 1147 + 1149 3 -45.430 158.440 -1.060 0.242 1148 + 1150 3 -47.670 160.060 0.700 0.224 1149 + 1151 3 -41.580 78.990 -19.000 0.480 545 + 1152 3 -44.460 81.260 -20.620 0.480 1151 + 1153 3 -46.380 81.910 -20.000 0.480 1152 + 1154 3 -48.620 85.470 -19.940 0.480 1153 + 1155 3 -51.490 87.410 -19.580 0.480 1154 + 1156 3 -53.730 88.060 -20.700 0.480 1155 + 1157 3 -56.140 88.860 -21.420 0.320 1156 + 1158 3 -58.380 89.180 -21.360 0.320 1157 + 1159 3 -60.620 88.860 -23.580 0.320 1158 + 1160 3 -63.810 88.860 -24.000 0.320 1159 + 1161 3 -66.050 87.560 -26.240 0.320 1160 + 1162 3 -66.690 86.910 -30.900 0.320 1161 + 1163 3 -68.610 86.590 -33.200 0.308 1162 + 1164 3 -70.850 85.620 -34.320 0.310 1163 + 1165 3 -72.450 85.940 -38.100 0.275 1164 + 1166 3 -72.450 87.560 -40.120 0.242 1165 + 1167 3 -74.690 87.560 -41.780 0.240 1166 + 1168 3 -76.290 86.910 -43.120 0.225 1167 + 1169 3 -78.530 84.970 -44.540 0.249 1168 + 1170 3 -80.440 83.350 -47.780 0.217 1169 + 1171 3 -55.330 91.300 -16.860 0.480 1156 + 1172 3 -57.250 92.270 -13.840 0.432 1171 + 1173 3 -60.450 93.570 -13.820 0.320 1172 + 1174 3 -63.330 95.510 -14.320 0.322 1173 + 1175 3 -65.240 97.460 -14.420 0.350 1174 + 1176 3 -69.080 98.750 -14.420 0.330 1175 + 1177 3 -70.680 101.020 -15.680 0.300 1176 + 1178 3 -71.320 102.960 -16.420 0.295 1177 + 1179 3 -72.920 104.580 -16.760 0.242 1178 + 1180 3 -75.480 105.560 -16.760 0.261 1179 + 1181 3 -76.760 108.470 -16.780 0.211 1180 diff --git a/python/example/reading-external-fields/external-field-data.py b/python/example/reading-external-fields/external-field-data.py new file mode 100644 index 000000000..8875fe5e9 --- /dev/null +++ b/python/example/reading-external-fields/external-field-data.py @@ -0,0 +1,98 @@ +import arbor as A +from arbor import units as U + +import matplotlib.pyplot as plt +from pathlib import Path +import numpy as np +import ctypes + +E = np.zeros(shape=(3,), dtype=np.float64) +buf = ctypes.c_double * 3 +addr = ctypes.addressof(buf.from_buffer(E)) + +ctr = "(location 0 0.5)" + +here = Path(__file__).parent +mrf = A.load_swc_neuron(here / "Acker2008.swc") + +dec = ( + A.decor() + .set_property(Vm=-55 * U.mV) + .paint('"soma"', A.density("hh")) + .paint('"dend"', A.density("pas/e=-70")) +) + +cvp = A.cv_policy("(every-segment)") + +prp = A.neuron_cable_properties() +cat = A.load_catalogue(here / "efields-catalogue.so") +prp.catalogue.extend(cat, "") + +e0 = 30.0e-3 +omega = 0.05 # 1/ms + +for id, seg in enumerate(mrf.segment_tree.segments): + loc = f"(support (on-components 0.5 (segment {id})))" + dec.place( + loc, + A.synapse( + f"reader/field={addr}", + xp=seg.prox.x, + xd=seg.dist.x, + yp=seg.prox.y, + yd=seg.dist.y, + zp=seg.prox.z, + zd=seg.dist.z, + ), + label=f"ef{id}", + ) + + +class recipe(A.recipe): + def __init__(self): + A.recipe.__init__(self) + + def num_cells(self): + return 1 + + def cell_kind(self, _): + return A.cell_kind.cable + + def global_properties(self, kind): + return prp + + def cell_description(self, _): + return A.cable_cell(mrf.morphology, dec, mrf.labels, cvp) + + def probes(self, gid): + return [A.cable_probe_membrane_voltage(ctr, tag="Um")] + + +if __name__ == "__main__": + T = 1000 + dt = 0.01 + rec = recipe() + sim = A.simulation(rec) + hdl = sim.sample((0, "Um"), A.regular_schedule(10 * dt * U.ms)) + t = 0 + while t < T: + E[0] = e0 * np.sin(omega * (t + 50)) + sim.run((t + 1) * U.ms, dt * U.ms) + t += 1 + fg, ax = plt.subplots() + for data, meta in sim.samples(hdl): + ax.plot(data[:, 0], data[:, 1], label=meta) + ax.legend() + ax.set_xlabel("Time $(ms)$") + ax.set_ylabel("Membrane potential $(mV)$") + ax.set_xlim(0, T) + ax.set_ylim(-120, -20) + + bx = ax.twinx() + bx.plot(data[:, 0], e0 * np.sin(omega * (50 + data[:, 0])), label="E(t)", color="r") + bx.set_xlim(0, T) + bx.set_ylim(-0.05, 0.05) + bx.set_ylabel("Induced Current $(nA)$") + + fg.savefig("external-fields-external-field-data.svg", bbox_inches="tight") + fg.savefig("external-fields-external-field-data.pdf", bbox_inches="tight") diff --git a/python/example/reading-external-fields/external-field.py b/python/example/reading-external-fields/external-field.py new file mode 100644 index 000000000..5b07cef2b --- /dev/null +++ b/python/example/reading-external-fields/external-field.py @@ -0,0 +1,95 @@ +import arbor as A +from arbor import units as U + +import matplotlib.pyplot as plt +from pathlib import Path +import numpy as np + +here = Path(__file__).parent + +E = np.zeros(shape=(3,), dtype=np.float64) + +ctr = "(location 0 0.5)" + +mrf = A.load_swc_neuron(here / "Acker2008.swc") + +dec = ( + A.decor() + .set_property(Vm=-55 * U.mV) + .paint('"soma"', A.density("hh")) + .paint('"dend"', A.density("pas/e=-70")) +) + +cvp = A.cv_policy("(every-segment)") + +prp = A.neuron_cable_properties() +cat = A.load_catalogue(here / "efields-catalogue.so") +prp.catalogue.extend(cat, "") + +e0 = 30.0e-3 +omega = 0.05 # 1/ms + +for id, seg in enumerate(mrf.segment_tree.segments): + loc = f"(support (on-components 0.5 (segment {id})))" + dec.place( + loc, + A.synapse( + f"reader/field={0xC0FFEE}", + xp=seg.prox.x, + xd=seg.dist.x, + yp=seg.prox.y, + yd=seg.dist.y, + zp=seg.prox.z, + zd=seg.dist.z, + ), + label=f"ef{id}", + ) + + +class recipe(A.recipe): + def __init__(self): + A.recipe.__init__(self) + + def num_cells(self): + return 1 + + def cell_kind(self, _): + return A.cell_kind.cable + + def global_properties(self, kind): + return prp + + def cell_description(self, _): + return A.cable_cell(mrf.morphology, dec, mrf.labels, cvp) + + def probes(self, gid): + return [A.cable_probe_membrane_voltage(ctr, tag="Um")] + + +if __name__ == "__main__": + T = 1000 + dt = 0.01 + rec = recipe() + sim = A.simulation(rec) + hdl = sim.sample((0, "Um"), A.regular_schedule(10 * dt * U.ms)) + t = 0 + while t < T: + E[0] = e0 * np.sin(omega * t) + sim.run((t + T) * U.ms, dt * U.ms) + t += 1 + fg, ax = plt.subplots() + for data, meta in sim.samples(hdl): + ax.plot(data[:, 0], data[:, 1], label=meta) + ax.legend() + ax.set_xlabel("Time $(ms)$") + ax.set_ylabel("Membrane potential $(mV)$") + ax.set_xlim(0, T) + ax.set_ylim(-120, -20) + + bx = ax.twinx() + bx.plot(data[:, 0], e0 * np.sin(omega * data[:, 0]), label="E(t)", color="r") + bx.set_xlim(0, T) + bx.set_ylim(-0.05, 0.05) + bx.set_ylabel("Induced Current $(nA)$") + + fg.savefig("external-fields-external-field.svg", bbox_inches="tight") diff --git a/python/example/reading-external-fields/implicit-field.py b/python/example/reading-external-fields/implicit-field.py new file mode 100644 index 000000000..59e0cb7db --- /dev/null +++ b/python/example/reading-external-fields/implicit-field.py @@ -0,0 +1,89 @@ +import arbor as A +from arbor import units as U + +import matplotlib.pyplot as plt +from pathlib import Path +import numpy as np + +here = Path(__file__).parent + +ctr = "(location 0 0.5)" + +mrf = A.load_swc_neuron(here / "Acker2008.swc") + +dec = ( + A.decor() + .set_property(Vm=-55 * U.mV) + .paint('"soma"', A.density("hh")) + .paint('"dend"', A.density("pas/e=-70")) +) + +cvp = A.cv_policy("(every-segment)") + +prp = A.neuron_cable_properties() +cat = A.load_catalogue(here / "efields-catalogue.so") +prp.catalogue.extend(cat, "") + +e0 = 30.0e-3 +omega = 0.05 # 1/ms + +for id, seg in enumerate(mrf.segment_tree.segments): + loc = f"(support (on-components 0.5 (segment {id})))" + dec.place( + loc, + A.synapse( + f"efield/e0={e0},omega={omega}", + xp=seg.prox.x, + xd=seg.dist.x, + yp=seg.prox.y, + yd=seg.dist.y, + zp=seg.prox.z, + zd=seg.dist.z, + ), + label=f"ef{id}", + ) + + +class recipe(A.recipe): + def __init__(self): + A.recipe.__init__(self) + + def num_cells(self): + return 1 + + def cell_kind(self, _): + return A.cell_kind.cable + + def global_properties(self, kind): + return prp + + def cell_description(self, _): + return A.cable_cell(mrf.morphology, dec, mrf.labels, cvp) + + def probes(self, gid): + return [A.cable_probe_membrane_voltage(ctr, tag="Um")] + + +if __name__ == "__main__": + T = 1000 + dt = 0.01 + rec = recipe() + sim = A.simulation(rec) + hdl = sim.sample((0, "Um"), A.regular_schedule(10 * dt * U.ms)) + sim.run(T * U.ms, dt * U.ms) + fg, ax = plt.subplots() + for data, meta in sim.samples(hdl): + ax.plot(data[:, 0], data[:, 1], label=meta) + ax.legend() + ax.set_xlabel("Time $(ms)$") + ax.set_ylabel("Membrane potential $(mV)$") + ax.set_xlim(0, T) + ax.set_ylim(-120, -20) + + bx = ax.twinx() + bx.plot(data[:, 0], e0 * np.sin(omega * data[:, 0]), label="E(t)", color="r") + bx.set_xlim(0, T) + bx.set_ylim(-0.05, 0.05) + bx.set_ylabel("Induced Current $(nA)$") + + fg.savefig("external-fields-implicit-field.svg", bbox_inches="tight") diff --git a/python/example/reading-external-fields/main.py b/python/example/reading-external-fields/main.py new file mode 100644 index 000000000..63f77b6be --- /dev/null +++ b/python/example/reading-external-fields/main.py @@ -0,0 +1,2 @@ +#!/usr/bin/env python3 + diff --git a/python/example/reading-external-fields/mod/efield.mod b/python/example/reading-external-fields/mod/efield.mod new file mode 100644 index 000000000..a895dd405 --- /dev/null +++ b/python/example/reading-external-fields/mod/efield.mod @@ -0,0 +1,41 @@ + NEURON { + POINT_PROCESS efield + RANGE xp, yp, zp, xd, yd, zd + NONSPECIFIC_CURRENT i + } + + PARAMETER { + xp yp zp : proximal coordinates + xd yd zd : distal coordinates + e0 : amplitude + omega : frequency + } + + ASSIGNED { da } + + STATE { t } + + INITIAL { + LOCAL dx, dy, dz + : vector spanned by segment endopints + dx = xd - xp + dy = yd - yp + dz = zd - zp + : unit length of the segment + da = 1.0/sqrt(dx*dx + dy*dy + dz*dz) + : time + t = 0 + } + + BREAKPOINT { + SOLVE state METHOD cnexp + LOCAL e + : electrical field along x + e = e0 * sin(omega*t) + : induced current + i = e*(xd - xp)*da + } + + DERIVATIVE state { t' = 1 } + + NET_RECEIVE(weight) {} diff --git a/python/example/reading-external-fields/mod/reader.hpp b/python/example/reading-external-fields/mod/reader.hpp new file mode 100644 index 000000000..a855ad5a0 --- /dev/null +++ b/python/example/reading-external-fields/mod/reader.hpp @@ -0,0 +1,60 @@ +#pragma once + +#include +#include +#include + +extern "C" { + arb_mechanism_type make_arb_efields_catalogue_reader_type() { + // Tables + static arb_field_info globals[] = {{ "field", "", NAN, -179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000000, 179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000000 } }; + static arb_size_type n_globals = 1; + static arb_field_info state_vars[] = {{ "da", "", NAN, -179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000000, 179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000000 } }; + static arb_size_type n_state_vars = 1; + static arb_field_info parameters[] = {{ "xp", "", NAN, -179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000000, 179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000000 }, + { "yp", "", NAN, -179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000000, 179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000000 }, + { "zp", "", NAN, -179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000000, 179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000000 }, + { "xd", "", NAN, -179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000000, 179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000000 }, + { "yd", "", NAN, -179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000000, 179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000000 }, + { "zd", "", NAN, -179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000000, 179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.000000 } }; + static arb_size_type n_parameters = 6; + static arb_ion_info* ions = NULL; + static arb_size_type n_ions = 0; + static arb_random_variable_info* random_variables = NULL; + static arb_size_type n_random_variables = 0; + + arb_mechanism_type result; + result.abi_version=ARB_MECH_ABI_VERSION; + result.fingerprint=""; + result.name="reader"; + result.kind=arb_mechanism_kind_point; + result.is_linear=true; + result.has_post_events=false; + result.globals=globals; + result.n_globals=n_globals; + result.ions=ions; + result.n_ions=n_ions; + result.state_vars=state_vars; + result.n_state_vars=n_state_vars; + result.parameters=parameters; + result.n_parameters=n_parameters; + result.random_variables=random_variables; + result.n_random_variables=n_random_variables; + return result; + } + + arb_mechanism_interface* make_arb_efields_catalogue_reader_interface_multicore(); + arb_mechanism_interface* make_arb_efields_catalogue_reader_interface_gpu(); + + #ifndef ARB_GPU_ENABLED + arb_mechanism_interface* make_arb_efields_catalogue_reader_interface_gpu() { return nullptr; } + #endif + + arb_mechanism make_arb_efields_catalogue_reader() { + static arb_mechanism result = {}; + result.type = make_arb_efields_catalogue_reader_type; + result.i_cpu = make_arb_efields_catalogue_reader_interface_multicore; + result.i_gpu = make_arb_efields_catalogue_reader_interface_gpu; + return result; + } +} // extern C diff --git a/python/example/reading-external-fields/mod/reader_cpu.cpp b/python/example/reading-external-fields/mod/reader_cpu.cpp new file mode 100644 index 000000000..ae4a8a8d7 --- /dev/null +++ b/python/example/reading-external-fields/mod/reader_cpu.cpp @@ -0,0 +1,132 @@ +#include +#include +#include +#include +#include +#include + +#include + +namespace arb { +namespace efields_catalogue { +namespace kernel_reader { + +using ::arb::math::exprelr; +using ::arb::math::safeinv; +using ::std::abs; +using ::std::cos; +using ::std::exp; +using ::std::max; +using ::std::min; +using ::std::pow; +using ::std::sin; +using ::std::sqrt; +using ::std::tanh; + +static constexpr unsigned simd_width_ = 1; +static constexpr unsigned min_align_ = std::max(alignof(arb_value_type), alignof(arb_index_type)); +using ::std::log; + +#define PPACK_IFACE_BLOCK \ +[[maybe_unused]] auto _pp_var_width = pp->width;\ +[[maybe_unused]] auto _pp_var_n_detectors = pp->n_detectors;\ +[[maybe_unused]] auto _pp_var_dt = pp->dt;\ +[[maybe_unused]] arb_index_type * __restrict__ _pp_var_vec_ci = pp->vec_ci;\ +[[maybe_unused]] arb_value_type * __restrict__ _pp_var_vec_v = pp->vec_v;\ +[[maybe_unused]] arb_value_type * __restrict__ _pp_var_vec_i = pp->vec_i;\ +[[maybe_unused]] arb_value_type * __restrict__ _pp_var_vec_g = pp->vec_g;\ +[[maybe_unused]] arb_value_type * __restrict__ _pp_var_temperature_degC = pp->temperature_degC;\ +[[maybe_unused]] arb_value_type * __restrict__ _pp_var_diam_um = pp->diam_um;\ +[[maybe_unused]] arb_value_type * __restrict__ _pp_var_area_um2 = pp->area_um2;\ +[[maybe_unused]] arb_value_type * __restrict__ _pp_var_time_since_spike = pp->time_since_spike;\ +[[maybe_unused]] arb_index_type * __restrict__ _pp_var_node_index = pp->node_index;\ +[[maybe_unused]] arb_index_type * __restrict__ _pp_var_peer_index = pp->peer_index;\ +[[maybe_unused]] arb_index_type * __restrict__ _pp_var_multiplicity = pp->multiplicity;\ +[[maybe_unused]] arb_value_type * __restrict__ _pp_var_weight = pp->weight;\ +[[maybe_unused]] auto& _pp_var_events = pp->events;\ +[[maybe_unused]] auto _pp_var_mechanism_id = pp->mechanism_id;\ +[[maybe_unused]] arb_size_type _pp_var_index_constraints_n_contiguous = pp->index_constraints.n_contiguous;\ +[[maybe_unused]] arb_size_type _pp_var_index_constraints_n_constant = pp->index_constraints.n_constant;\ +[[maybe_unused]] arb_size_type _pp_var_index_constraints_n_independent = pp->index_constraints.n_independent;\ +[[maybe_unused]] arb_size_type _pp_var_index_constraints_n_none = pp->index_constraints.n_none;\ +[[maybe_unused]] arb_index_type* __restrict__ _pp_var_index_constraints_contiguous = pp->index_constraints.contiguous;\ +[[maybe_unused]] arb_index_type* __restrict__ _pp_var_index_constraints_constant = pp->index_constraints.constant;\ +[[maybe_unused]] arb_index_type* __restrict__ _pp_var_index_constraints_independent = pp->index_constraints.independent;\ +[[maybe_unused]] arb_index_type* __restrict__ _pp_var_index_constraints_none = pp->index_constraints.none;\ +[[maybe_unused]] auto _pp_var_field = pp->globals[0];\ +[[maybe_unused]] auto const * const * _pp_var_random_numbers = pp->random_numbers;\ +[[maybe_unused]] arb_value_type* __restrict__ _pp_var_da = pp->state_vars[0];\ +[[maybe_unused]] arb_value_type* __restrict__ _pp_var_xp = pp->parameters[0];\ +[[maybe_unused]] arb_value_type* __restrict__ _pp_var_yp = pp->parameters[1];\ +[[maybe_unused]] arb_value_type* __restrict__ _pp_var_zp = pp->parameters[2];\ +[[maybe_unused]] arb_value_type* __restrict__ _pp_var_xd = pp->parameters[3];\ +[[maybe_unused]] arb_value_type* __restrict__ _pp_var_yd = pp->parameters[4];\ +[[maybe_unused]] arb_value_type* __restrict__ _pp_var_zd = pp->parameters[5];\ +//End of IFACEBLOCK + + +// interface methods +static void init(arb_mechanism_ppack* pp) { + PPACK_IFACE_BLOCK; + for (arb_size_type i_ = 0; i_ < _pp_var_width; ++i_) { + arb_value_type dy, dz, dx; + dx = _pp_var_xd[i_]-_pp_var_xp[i_]; + dy = _pp_var_yd[i_]-_pp_var_yp[i_]; + dz = _pp_var_zd[i_]-_pp_var_zp[i_]; + _pp_var_da[i_] = 1.0/sqrt(dx*dx+dy*dy+dz*dz); + } + if (!_pp_var_multiplicity) return; + for (arb_size_type ix = 0; ix < 0; ++ix) { + for (arb_size_type iy = 0; iy < _pp_var_width; ++iy) { + pp->state_vars[ix][iy] *= _pp_var_multiplicity[iy]; + } + } +} + +static void advance_state(arb_mechanism_ppack* pp) { +} + +static void compute_currents(arb_mechanism_ppack* pp) { + PPACK_IFACE_BLOCK; + auto e_field_ptr = (double*)((uint64_t) _pp_var_field); + for (arb_size_type i_ = 0; i_ < _pp_var_width; ++i_) { + auto node_indexi_ = _pp_var_node_index[i_]; + arb_value_type e = e_field_ptr[0]; + arb_value_type i = e*(_pp_var_xd[i_]-_pp_var_xp[i_])*_pp_var_da[i_]; + _pp_var_vec_i[node_indexi_] = fma(_pp_var_weight[i_], i, _pp_var_vec_i[node_indexi_]); + } +} + +static void write_ions(arb_mechanism_ppack* pp) { +} + +static void apply_events(arb_mechanism_ppack* pp, arb_deliverable_event_stream* stream_ptr) { + PPACK_IFACE_BLOCK; + auto [begin_, end_] = *stream_ptr; + for (; begin_