diff --git a/.gitignore b/.gitignore index a93d24249..e51e9d4c9 100644 --- a/.gitignore +++ b/.gitignore @@ -50,6 +50,7 @@ build/ # Tex files *.aux *.log +*.synctex.gz *.out *.toc diff --git a/Notes_Upgrades/SolidAngle/Solid_Angle_Sphere_Approx.ipynb b/Notebooks/Solid_Angle_Sphere_Approx.ipynb similarity index 100% rename from Notes_Upgrades/SolidAngle/Solid_Angle_Sphere_Approx.ipynb rename to Notebooks/Solid_Angle_Sphere_Approx.ipynb diff --git a/Notebooks/solid_angle_triangle.ipynb b/Notebooks/solid_angle_triangle.ipynb new file mode 100644 index 000000000..b9bb2272e --- /dev/null +++ b/Notebooks/solid_angle_triangle.ipynb @@ -0,0 +1,386 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "guided-triple", + "metadata": {}, + "source": [ + "# Benchmarking computation of solid angle for triangular aperture" + ] + }, + { + "cell_type": "markdown", + "id": "suffering-buddy", + "metadata": {}, + "source": [ + "From [wikipedia](https://en.wikipedia.org/wiki/Solid_angle#Tetrahedron):\n", + "\n", + "Let OABC be the vertices of a tetrahedron with an origin at $O$ subtended by the triangular face $ABC$ where $\\vec a\\ ,\\, \\vec b\\ ,\\, \\vec c$ are the vector positions of the vertices $A$, $B$ and $C$. Define the vertex angle $\\theta_a$ to be the angle $BOC$ and define $\\theta_b$, $\\theta_c$ correspondingly. Let $\\phi_{ab}$ be the dihedral angle between the planes that contain the tetrahedral faces $OAC$ and $OBC$ and define $\\phi_{ac}$, $\\phi_{bc}$ correspondingly. The solid angle $\\Omega$ subtended by the triangular surface $ABC$ is given by\n", + "\n", + "$$\n", + "\\Omega = \\left(\\phi_{ab} + \\phi_{bc} + \\phi_{ac}\\right)\\ - \\pi\\\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "biological-cedar", + "metadata": {}, + "source": [ + "In ``tofu`` we have the coordinates of the points $O$, $A$, $B$, and $C$. If we want to compute the solid angle this ways, we would have to compute the dihedral angles $\\phi_{ab}$, $\\phi_{bc}$, and $\\phi_{ac}$ first. To do so, the computation of the latter is described as such:" + ] + }, + { + "cell_type": "markdown", + "id": "public-blend", + "metadata": {}, + "source": [ + "\\[...\\] if $\\mathbf{n}_A$ and $\\mathbf{n}_B$ are normal vector to the planes (defining the angle), one has:\n", + "\n", + "$$\\cos \\varphi = \\frac{ \\left\\vert\\mathbf{n}_\\mathrm{A} \\cdot \\mathbf{n}_\\mathrm{B}\\right\\vert}{|\\mathbf{n}_\\mathrm{A} | |\\mathbf{n}_\\mathrm{B}|}$$\n", + "\n", + "where $\\mathbf{n}_A \\cdot \\mathbf{n}_B$ is the dot product of the vectors and $|\\mathbf{n}_\\mathrm{A} | |\\mathbf{n}_\\mathrm{B}|$ is the product of their lengths.\n", + "\n", + "The absolute value is required in above formulas, as the planes are not changed when changing all coefficient signs in one equation, or replacing one normal vector by its opposite.\n" + ] + }, + { + "cell_type": "markdown", + "id": "cross-costs", + "metadata": {}, + "source": [ + "## Example of computation with original formula" + ] + }, + { + "cell_type": "markdown", + "id": "alternative-baker", + "metadata": {}, + "source": [ + "Let's generate four arbitrary points in space and compute the solid angle on the triangle using this formula" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "id": "chicken-fleet", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy\n", + "from mpl_toolkits import mplot3d" + ] + }, + { + "cell_type": "code", + "execution_count": 76, + "id": "demonstrated-aruba", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 0., 1., -5., 1.],\n", + " [ 0., -5., 4., -1.],\n", + " [ 0., -6., 6., 3.]])" + ] + }, + "execution_count": 76, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# We create randomly the coordinates of A, B, C\n", + "coord = numpy.random.randint(-10, 10, size=(3, 4)) + 0.0\n", + "# We set O to orign\n", + "coord[:, 0] = 0\n", + "coord" + ] + }, + { + "cell_type": "code", + "execution_count": 77, + "id": "typical-separation", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 78, + "id": "executive-memorabilia", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 78, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPUAAADyCAYAAACVtIPdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABkX0lEQVR4nO19eXhkZZn9ubWlqrJXKns6WyfdSXf2dELTLeCwCzTdCMiICoogLoCIuDIqMy4og4PoT8FREMQGh+5m3wQRFdl6zb519j2VWlL7euv+/gjf5VallnurbiWddJ3nmWewq+5SqTr3e7/3Pe95KYZhkEQSSWwcSNb6BpJIIglxkSR1EklsMCRJnUQSGwxJUieRxAZDktRJJLHBkCR1EklsMMiivJ6sdyWRROJBiXmy5EqdRBIbDElSJ5HEBkOS1EkkscGQJHUSSWwwJEmdRBIbDElSJ5HEBkOS1EkkscGQJHUSSWwwJEmdRBIbDElSJ5HEBkOS1EkkscGQJHUSSWwwJEmdRBIbDElSJ5HEBkOS1EkkscGQJPUagGEYeDwe+Hw+JC2akxAb0UwSkhAZfr8fHo8HLpeL/TepVAq5XA6ZTAapVAqKErVnPonTDFSUlSK5jIgEhmHg8/ng8/lAURS8Xi/77wzDwO/3s2R2u91IT0+HQqFIkvz0gKhfcHKlXgWQcJtLXAKKokBRFCQSCfvekZERlJeXQ61WA0iu5EkIQ5LUCYbP58P09DRomkZxcTEoimJX51DkJCSXSqWQSqXsKu50Otn3y2Qy9v+SJE8iGElSJwjccNvv97Nht1CEWslpmobP52PfI5PJ2JVcIpEkSX6aI0nqBMDv98Pr9bLhNlmdCRwOB06ePInU1FRkZ2cjLS2NJS2AFe/ngpyPIJjkFEUFrORJkp9+SJJaRBCCkSQYISqXpHNzcxgdHUVlZSU8Hg+mp6dhs9mgVCqRlZWF7OxsQWWuUCT3+Xw4evQoGhsbkyQ/DZEktUhgGAZerxc0Ta8gGkVRoGkaPT098Pl8aG9vZ/fVhYWFYBgGLpcLJpMJk5OTMBqN8Hq9yM3NRXZ2NtRqNW8ikmvTNM3uyX0+H/ugSZJ84yNJahFAas8k+RVMEpfLhenpaVRVVaGkpAQURcHj8cBiAd5/XwqTiUJRkRRtbSoUFRWht7cXeXl5cLvdGB0dhcPhQFpaGruSq1QqwSQnIA+fYJLL5XJIpdIkyTcAkqSOA8G1Z+6+mLw+PT2NyclJ5ObmYtOmTexrHg/wl7/I4HZTkMudGBpSwW6ncNFFNCQSCZRKJXJzc1FSUgKGYWC322EymTA8PAyXy4W0tDRkZ2cjOzsbSqWS9z2TzDr3HoNJTpJuMpks5EMqiVMbSVLHiODac/AP3+v1ore3FzKZDFu3boXZbA543WymYLFQmJh4B//7v/+L//7vezEzkw+3m2bPT0BRFNLS0pCWloZNmzaBYRhYrVaYTCYMDAzA4/EgIyODXclTUlJ4f45QJPd4PHC73ZidnUV+fj7UajW7kidJfuojSeoYQJJh4cJts9mMnp4eVFZWorCwEHq9fkXySyplwDDA5s3VgKYcTz73V1x96acglSIqaSiKQkZGBjIyMlBWVga/3w+r1Qqj0Yi5uTn4fD64XC7odDpkZ2dDLpfz/mxckptMJuTl5bEkJ9FIcLiexKmFJKkFgE+4PT4+joWFBTQ3N7OKsFAlquxsoKaGRl9fHraVNaHHVYHP509CLi+MWNIKBYlEgszMTGRmZgJYfugcPnwYNpuNFb6QVTwrKwsyGf+vXSKRsCQn9+TxeODxeNjXkyQ/tZAkNU/4/X7Y7XYMDQ1h+/btK1ZTj8eD7u5uqNVqtLe3R607UxSwe7cf5eU+bNnSipufeh/3HkvB2Wdr4r5XqVQKmUyGyspKAMskX1pagslkwvj4OAAEkJwbfnOx8p4p9vzc15MkP7WQJHUUcGvPwXJNAqPRiP7+flRXVyMvL2/FOcKtvBQFbNrEYNOmQlS8bsCstBqf++MJ3H12FvJFbMmUSqXIyclBTk4OgGXp6tLSEoxGI8bGxkBRFJt0y8jICCB5pK1AKJKTPTkhOSmtpaamJkm+SkiSOgKCa8+k7st9fXh4GCaTCa2trWGz0HzC6c9evBM/ecuAySUPfvhPI359VS6yxPwwHMhkMmi1Wmi1WgDLST2TyQSdTofh4WHIZDJkZWWxDzK+CFU+W1xchNPpRGlpKYDllZzbnJIkufhIkjoMQtWeJRIJ+yN3uVzo6upCdnY22traoq5o0Uh9+Rlbcc8/3wJDSTG25MM3XhjFI9dlQa0IHRqLCblcjry8PDbKcLvdWFpagtvtRkdHB1JSUtiVPD09XVCNHEBAcwrDMHC73XC73WAYJiBUJyW0JOJDktRB4IbbwckwQk6dToeTJ0+ipqaGDWkjgQ+pU1NkqM1ToW9xOWztXXDilqd68Jtr6qCUJ57YXKSkpCA/Px8zMzNoaGiAz+eDyWQKkLQSkqempkYkIrcbLdRK7vf7k4YRIiNJag6i1Z4ZhoHT6cTU1BTa2tqgUCh4nTeY1F6vF0NDQ5DL5dBoNEhPT4dEIsF524rQ949x9n1HJsz42qE+/PLq7ZBL1yZMpSgKSqUShYWFrKTV6XRiaWkJExMTsNlsUKvVLMmDJa3hWkzJuZMkFx9JUn+AaFJPh8OBrq4uUBSFlpYWQT8uLqmXlpbQ29uL0tJSMAyD2dlZWK1WpKSkYHOqesWx/xox4VvPDuC/r6iFVLK6P+jQyT0KarUaarUaRUVFYBgGDocDJpMpQNJKSB7KGCIcwpHc6XRifHwcmzZtQkpKSpLkUXDakzpa7Rn4sLOqrq4Ovb29gn9IFEXB7/djfHwc8/PzaG5uhlwuh9/vR0FBAQDA6XQix2hEuoKC1RNIptcH9Pjei4P48Z6tq/4j5iOESU1NRWpqKitptdlsMJlMOHnyJKxWK5RKJSQSSUySVvKdWK1W9u/IrUAkV/KVOK1JHarvmQuaptHf3892VglRZgWfZ2lpCSkpKWwNm2itCVQqFUqKi7G7yoJX+xZXnOOFbh0Yrxt3X1YjSAYaD2JxOqUoCunp6UhPT0dpaSkmJyfh8/ng8XgCJK1kJee7hfH7/Wy2nGsYkXSFWYnTktTBfc+hCG21WtHT04OSkhK2syoWmM1mdHd3IyUlBdu2bYv6/l2V2SFJDQAvDphBu4/j6moZMjMzWWIIUYgJhRikUKvVKCgoYCWtFosFJpMJMzMzoGma/SxZWVlhH5yh9ubhXGFIkhM4PV1hTjtSMwwDs9kMh8MBjUYTMhk2PT2Nqakp1NfXIz09PebrTE5OYm5uDtu3b8fo6Civ43ZVZEd8/ZUxL0qLCnDt5kyYTCZMTEyAoihWIZaZmRlWISYUYniS+/3+gC2NRCJBVlYWsrKyUFFRAZqmWZJPTk6CYZgAkpMHVqSEG0Ek6ydyvNvtZs+7UUl+WpGahNukwym4HOXz+dDT0wOZTIYzzjgjZnKQDi2FQoG2tjZBpv35GSnYlCnDlNkX9j2/fXsKaUoZPrtzM3u9paUl6PV6jIyMQCaTwe12w2KxCKorh0K8P/poZJRKpWzEASx/B2azmZW0kgcW8XoT8p2ESrz19PSgtbWVfX0jGkacFqQODstkMtkKpZTZbEZvby8qKipQWFgY87UsFgt6enoCzkPTtKBVr6VQhSmzNeJ7fv7GGFRyKa5pLYJcLkdubi5yc3MBLItHjh07hpmZGTZRxbeuLDb4rLBcyGSyAEmr1+uF2WzG7Owsjh8/DqlUGhCVCFGkhZK1bkRXmA1P6lA2QySLSl4nnVWNjY1ITU2N+TpTU1OYmZlZcR6hXVfNhUo8NxCZ1ADw41eHoZJLcXlDfsC/k7JPbW0tW1cmK5/dbg8oOalUKv4fMgYIJXUw5HI5tFotlEol2tra4PF4sLS0FCBp5ardhJJ8I7rCbGhSh6s9SyQSVmgSrrNKCEjYLpfL0d7eviJEFErq+nwl5FIKXjryMQyA7784CJVCggtqckO+h1tXLi4uDig5DQ0Nwe12x5SN5ot4SR0MhUKxQtJqMpkC6v1kPy506yHEFeZUJvmGJHW02rNEIoHT6cSRI0fCdlZFOz/5Mkm4XV5ejqKiopDvF0rqFJkEdfkqnJh1RH0vzQDffGYAD1wtxdlV0ds2g0tOobLRsbiahoPYpA5GSkoKCgoK2Ho/MXAMJ2kVAr6uMKea9dOGIzUfqefU1BSsVit27dolSAwBBBJ0enoa09PTaGhoQFpaGq9jCLjNIaHev6NYzYvUAODzM7jjUB9+c00d2suz+H0Qzn0EZ6NJ77XD4cCxY8dYUgjdwwKJJ3UwQklauVsPt9uNmZkZZGVlCXJpBQJJbjQakZ+fH+AKQ1byte4l31B9b8TGJxyhXS4Xjhw5AoqikJmZKZjQwDIJPB4Purq6YDab0d7eHpHQQGhSR1sFW4uFrSpunx+3HuhF54xF0HHBIL3XVVVVSE1NZR9YOp0OR48eRUdHByYmJmCxWHit5KtNai7I1qO4uBh1dXVoaWmBUqkEwzAYHR3F4cOH0dvbi9nZWTidTkHnJn3ipLuMhOIejwd2ux3XXXcdBgYGEvTJImNDrNTB4XaoH9Hi4iKGhoZQU1ODtLQ0dHd3x3Qtv9+PY8eOoby8HMXFxTGdw+PxoKenBwzDQKPRQKPRBCSsKIpCRbYM2lQF9HYP7/M6PDS+9OcePPKphpjuKxSC2zK54a3Vao3YzAGIQ2oxtgHkPHK5nBUUhcovpKens58nknIvuLwWnFk3Go0JT0KGw7ondTSpp9/vx9DQEOx2O9tZJbT5H1j+QczMzMBut6OlpQUaDX/bIe49kYaOzZs3Q6FQsBppl8uFjIwMaDQa0DQNmUyGMyuz8EK3TtB9Wl0+3PxkN+5oSkwQFhzeBjdzEFJoNBqkpKSIQupgAUs85wkmYnB+wWazsU420SStkT4XqTKsBdYtqcONuOGCdFbl5+dj69YPmyEi7WdDwefzoa+vjxVCEENBofc7MTGBubk5NDc3Q6FQwOfzsba/3ITV/Pw8GIZBpTq2H4XR4cV9x4CWJic2ZSdutQjVzEGEPX19fezDMzU1FUqlMmbtvFghfLSHg0QiYV1ayfvDSVqj/X5It9paYF2SmpQauru7Q5oAAh92Vm3fvh1ZWVkBrwkhtdVqRXd3N0pLS1FSUoITJ04IXuV9Ph+cTidsNhva2toglUoDplaSeyIJK7lcDoZhoGXk+OX7S4gl+FxyAzft78Kj1zWhIGN1GkAoaqV1cWdnJ/twZRiGl+FhMMRcqYWcJ1QSkajdnE4njh49yr4e7NLq9XpFKQ9SFJUF4PcA6rBcxbyBYZh3Ix2z7kjNrT1bLJaYOqv4PvVnZmYwMTERoAEXWp6y2Wzo6uqCTCbD9u3beR1DthFbiguwNX8GAwt23tcLuH+zGzft78IfPtMIbZq49Wc+IH5kpaWlUKvVAYaHo6OjrESUaxQRCgzDrAmpgyGVStkciNFoRFNTU0hJ68TEhJiZ7wcAvMowzFUURSkARA0T1w2p+fQ9i9VZRdM0+vr6wDAM2tvbA57AQlb5ubk5jI2Nob6+PubE3JmV2TGTGgDGjU584clu/OHTDchUxRb+xgNu6BxseOjxeGAymTA3N4fBwcEAL7S0tDT2OCFGC5FAMtbxgjzUQ0laTSYTDh48iMnJSfzbv/0bLrzwQnz3u9+N6ToURWUAOBvAZz+4rgdA1MzpuiA1n9qzGJ1VwPLK2t3djU2bNqG4uDhku1+0ldrv97O9w8EPBT7gylh3VWTjD+9OC/sQQTips+OLT/bg95+qR2rK6n7lkVZZhUKB/Px85Ocvy1xJTXlychI2m42d361UKldlT80X4R4OpFLw61//GmeffTYOHDiAzs7OeC5VCWARwB8oimoEcAzAVxmGifiUP+Xr1DRNw+12hyW0z+dDZ2cnzGYzzjjjjLgIPTs7i66uLtTV1YVd6aOt1ESpplar0djYGFOvM/e6LZsyoZLH/zX1zFlx8xOdcHjCd38lAkKSXCrV8tTP7du3o729HeXl5fD7/ZiYmGAz0vPz83C73THdS6JJTeDz+SCVSpGbm4vzzz8/nkvJALQAeJBhmGYAdgDf5nPQKQk+4TZN03j//ffZmVWxguzDaZqOurJGWqn1ej0GBwexbds2tpUw0nkigVxDIZOgtTQT/xoxRfkU0dE5a8cND7+Lb+/OQp42Z1UaOmINnSnqw6GAmZmZmJmZQVFREZtZ9/l8go0ixCR1pPPY7faYG4OCMA1gmmGY9z/43wexXkkdrfZMykMulwtnnnlmXKUDu92Orq4u3vtwbmjMvZ+RkRGYTCbs2LEjbruh4AfHrspsUUgNAL1GBr/r8+POnV5WcMElR6xlp3AQS3wilUrZeWHl5eUBmejJyUkAH44SCmcUEY2MfBGtr1usGjXDMPMURU1RFLWVYZhBAOcB6It23ClFaj42Q9zOqoyMjLgIRBJZdXV1bG0yGkiHF/d+urq6kJmZiR07diREErm7UgOAn3MKH/xzxIw0pQL37G0APnCCMZlMmJqaAsMwbEb6VGnoCLXaczPRQGijiOCWTKEmC+EQLfwWcaUGgFsB7P8g8z0K4HPRDjhlSB2q7zkYwTOrjh49KrhmDCx/KS6XC/Pz84ITWdxVlKjDYun04nsNAKjUqlGQkYJ5S2x7yVB4uXcRKrkUP7ikeoXziMlkwuLiIhwOBzo6OqDRaFZkpPlitRRloYwiuC2ZROefkZER9z1FI7XNZhNNeMIwTAeAHUKOOSVIHc1zm4S3RqMxYGaVVCoFTdOCrkXCbZlMhrq6OsGJLIlEApqmA9RhsSjMaJoWFA6eWZGFZzoXBF8nEg51zEMll+JbF25m/00mk7HkMJvNqKmpWZGRJiskn4YYscJvoefgtmSSbq3h4WHo9XosLCzEZRQR7Xsjf6e1wpqSOthmKNQfyuVyobu7G1lZWStmVgmVe87Pz7Mqs6GhoZjCS+KUkp6ezqrDhMJqtaKrqwvAMokISchKGCoZt6tSIzqpAeBPR2agUkhw20crQr4erPW22+0wGo0YHBxk9+MajSaiE+haa79Jt1ZaWhoyMjKQk5MTl1FEtJV6LSWiwBqSmrh6Li4uYtOmTVE7q0LNrOK7UpO6sdvtRltbG+RyueAHArD8BJ6ZmUFubi5vdRgXDMNgbm4O4+PjqKurQ0pKCrxeL4xGI7sSpqWlsc4aXJxZkQUJBfjFm3DL4ndvT0GtkOLGXaUR38fNSJMGiFDJKo1GE1PvdTiILRONZhTh9/tDOppyz7OKe2rBWBNSk3Cb7N/ImFPu68GdVaHAh5hEd1xYWIja2tqYmzpIUq2wsDCmpzDDMOjv74fb7UZ7ezsoioLX60VKSkrASmiz2TA1NQWDwQCTycSSJCsrC9sK0tEzF927LBY88OY41HIprm3j305Kpm6Q/ThJVpEBggqFAh6Phw1HY12xxWroCLfCRjKKIPJPrlEE6aILh7Xs0AJWmdTBtWeZTLZipQ3XWRUK0Vbq+fl5jIyMoK6uDpmZmQGvhSpNhQJ3lW9vb8fs7KzgFd7tdsPhcKCgoIB9sIS6b7KC5ObmQqVSobS0NEArXZlKo0fQlYXhp6+NQKWQ4orGgpiOD05WuVwuHD9+nB2kR/axfPfjBKvd0EGMIoLln+RhRdM0MjIykJqaGtIHzW63C2rNFRurRupQtedgUkbqrAqFcKut3+/H4OAgnE5n2KaO4NJUKDidTvYBQ8gotKGDiCWUSiXKy8t5rTjkGsE/LlqziOeH+3lfWygYAHe/NASVXIKLt8WfzSftltu3bw8wJCASWj77cWD5+xRjCkmsD4dgo4jBwUFQFBXWKMJut2PTpk1x3y8AUBQlBXAUwAzDMJfxOSbhpI5UeyakjnVmVShSk5W+oKAANTU1YUkULfwm+/lgdVioOVihwJ3Q0dLSwrYexoMdFVqkKqSwe4Rl/IXAzwDfeW4QKTIJxAwgQ+1jzWYzjEYjO2UknBfaavVT84VEIoFWq2XNGblGEU8++SQOHz4Ms9mMs846CyUlJfFe7qsA+gHwE1IgwaSOZjNECHL48OGYOquCV/qFhQUMDw/zWunDkTqaOozPSk3TNHp7eyGRSNgMefBxkT5nuGvIJBTay7Pw5pAh4vXjhc/P4M6n+3FbkxxtCbpGqP04N8QlHVvECeZUaL0k4O7Ng40iamtrceutt8LpdOKGG27Af/3Xf2Hnzp0xXYeiqBIAlwL4MYA7+B6XMFLzkXpOT0/D6XTizDPPjKkRg9SMueE232HwoUjNRx0WbYV3OBzo7OzEpk2bAp7SQsL2SITfVZGdcFIDgIdm8KsODxq2m9G8KTP6AXEiOMTluoCaTCakpqaCpmnB43C5EKv1MtJ5SOXixhtvRFtb3I/EXwD4JgBB5BCd1Hxqz9yZVSTZEAukUikcDgeOHDmCvLy8iOF2MILJyVcdFomcJGQPl5gTEn6He++uysiNImLCTQNf+b8e/O5TDdheGHv3WyxQqVRs19bw8DCUSiW8Xi+7Hyc6byFTP1eroUMM8QlFUZcB0DEMc4yiqI8KOVZUUvORepKZVcT8/p133on5elarFXNzc2hqaoraFRUMQmru3pePOiwUObkhe7hIQehKHe69pRoVSrKUmF5y8TpXvLC6aXzxyW784TONqMpdm9orwzBsb3VpaWlAMwd3P67RaJCRkRHRQUWsvXm0OnU8LcAfYDeAyymKugSAEkAGRVF/Yhjm09EOFI3UDMPA7XZHlHpOTExgfn4+rplVwId1bKPRiMLCQsGEBj7czxPJKF91WPAKT86Rnp4esaFD6EodCbsqs/HU8TlRzsUHS04fbtrfjceua0SpZvVtb4ONFkI1c5hMJiwsLGBoaAgpKSmsXj24Pr4aDirEVTUeMAzzHQDfAYAPVuo7+RAaEJHUhMih/miRZlaRejHfsIiUmXJzc1FdXQ2j0RjT/Xq9XkxNTaGqqkqQfzeXnMSUcPPmzax7B5/jhFwjFHZVrC6pAUBv9+DG/V147LpGFGbGtqeNFdF6siPtx4kQRKPRxNT8E+5+VqmfOiaIGn6HSiIFd1YFg2Sw+ZA6uMxkMpkEN3QAy/Xwubk5FBcXCzbkJ59xdnYW4+PjUUfuEAQT1WAwwGAwICcnR/C0xvbyLMgkFHyJ0IxGwJzFjc88cgT/fXERNhfnRh1AJ1ZkInQvzN2Pk/q40WhkJ7Rw6+Ox1L+jhfHRFGcxXO/vAP7O9/0Jy36H66wKBiF1NPHByZMnYbVaA/asQqWeXHVYRUXoBoZoYBgGRqNRsP8YITVpCNHpdCgoKMDs7CwsFgvUajVycnKg0WiirtTpShnqi9JxYjq+ETuxYMHB4Pt/W8R/7HICHnvEri2x9rDxnIdbH9fpdGhpaVnhAMpnPx58zkj3utZICKm5nVU7duyI+IeKJvV0uVzo6upCTk4OWltbA/6gQlovg9Vhc3Nzgr2uXC4XBgYGIJFI0NjYKOiHRlEUfD4furq6IJfL0draCp/Pt6L7aWBgAC6XCxRFwWAwhPXHPrMye01IDQDjSx78/LgHv7+2BRLazd631+sNWAXDbceEQqysNRB+Pz4/Px91P84XYn3uWCEqqSmKitpZFYxIxCTnqq2tDaml5btSh1KHCV3lidyztLQUJpNJ8JdG0zQ78rakpCRgnxjc/USyukTzLZfL2R8i+aHtqszGb/45IegexET/vA23PNWLhz5Zj9LS0oCstNFoZFdBt9sNi8UieFY0F2L4fodbQUPtx8n9kyw2d4xQpHNFu9ZqQVRSj42NYXFxkbcABAhNar/fj+HhYVgsloieX9FWaoZhMDw8jKWlpRXn4UtqkrVfWFhgV1eDQZj4Q6/Xw2g0YuvWrbxkg1KpFEqlEtXV1QCWIwSj0YixsTE4HA5kZGQgNysbGUopLK7ESUaj4cS0Bbcd6MVvrqmDQiZZsQra7XZ0d3ezGmmhBgsEYvh+813tVSoVm2sJHiPk8/mQkZEBmqbh8/lCbr08Ho/oPm9CISqpi4qKwvZGh0MwMUm4rdFoVoTbwYhETKIOy8jICFlq4kNqn8+H3t5etuQlkUjgcDh4P4nJ/nlxcRG5ubm82/GC99RKpRJFRUUoKiqC3++H1WqFwWDAlkzg6OqUq8Pi/fEl3HGoD/dftQ1yaSBp5HI5lEoltm3btmKLERyqR8pNiBF+x3IOigocI0TTNIxGI/R6PTo6OgKkrmQ/Llbmm6KoTQD+CKAAgB/A/zIM8wCfY0UltVKpXDEjKhq4pCYWu/GG7nzUYdFITWyPguWefLq7gA/Dbblcjh07dmBgYECUsEwikbCumpeaVTj68sm4zxkv/jFsxHeeG8S9V9RAwnl4chNcwVuM4FBdIpGwq3hwqC6Wz1m8ElGpVIqMjAykpaWhsbGRnTBC9uM0TeOVV16BQqEQ4559AL7OMMxxiqLSARyjKOp1hmFOfTdRMizu5MmTIcPkSAhe0YSowyKRkzQVxNqHHUr/HU6JFuqL51vT3r2KktFo+Ev/IlRyCf7rsi3sZ4r0ww4O1T0eD4xGY8hQXYyVWqymEK7wJHjCiE6nA0VRGBoaQmNjI66//np8/etfj+k6DMPMAZj74L+tFEX1AyjGerAIZhgGo6OjKCgoEGyxy31vcKgc7akcaqUme3Cz2Ryz3NNgMGBgYGBFp1io4+JdfQozlSjPUWHc4IzrPGLh2a4FqBRSfPeiKgDCVliFQhFgFMj1QrPb7Th58mRcteXV0H3n5eXh8ssvx9LSEn73u98Jzr2EA0VR5QCaAbwf5a0AEpD9FgKDwYDJyUlWHRYryGTJsrIy3mKS4BWXqN7S09Mj7uUjtWyS/TOfls1YWi9DYVdF9ilDagB48ugsVHIpvnZuRcwhaHCofvjwYWi1Wl6hejgkanB9MIiCTSqVimIbTVFUGoBDAG5nGIZXDXNNVmquMGXz5s3weKIO8gsLor0WOhiPS06LxYKenp6Y5Z5k/yyTycLW5cVqvQzG7spsPHF0lvf7VwOPvDsFtUKCTzdrRanXUhTFK1SPZPe7WnO0xLQHpihKjmVC72cY5mm+x606qd1ud0DPssFggNMpfKUh6jCv14vdu3fH5N/t9/sxMzODyclJ3k0mwSt1uP7pYIjVehmMHWVZkEspeOm1VzJx8f/+MQEpQ2NXjvgijOBQ3eFwwGg0sna/XLNG8rtYjV5qQDzTQWr5afgwgH6GYf5HyLGrSmqy39yyZQtrTheLIT9XHaZSqWLW2VosFtadRIjckyDc/jnccVyi+v1+uFyukGNahTwA1AopmkoycGTCzOv9q4kH/jkNV0s6tm1L3DW4ziObNm0KsEkioXp2djZrDRwvVovUWG69/AyAboqiOj74t+8yDPNytANXZU9NkmEGg2GFDlwoqYPVYXNzc4JDK5fLhc7OzpjknuTzjI+PY2FhgXe2nktUr9eLjo4OVsSQmZmJnJzlKZSxrCa7KrNPSVIDwO+OW1GxSYdL68QbSxQJwTZJpOw0MzMDh8PBOn1qNJqYJn7ymXgpxl6aYZh/AYjpKZTwlZqIQEi/cfAfhC+pw6nDSDjMl9Ska2zLli0YGRmJSe7pcrlYT3K+1yWkJkm9zZs3Iysrix1qYDAYMDY2BplMhoyMDHi9Xt6Jpl0V2XjgzXFBn2O14AfwHy8MQimX4Lyt2lW/Pik7EXutnJycgFCdCGCEjMON5vm9lm2XQIJJzSUQCbeDwYfUkdRh5PhoX0iw3DMlJQUnTwoTbjidTnR0dEAqlQqe0EFRFMxmM8bGxlBfX4+0tDR4PJ4VK4vL5cLCwgIcDgcOHz4csIqH+4y1BWnQqOUwOqK7nK4FfH4G33imH7+6ejt2bxbmhy1m+6ZUKg0bqk9MTLDfBenYCvVApWk6YmS21kb+QILCb4ZhWB14S0tLxDAnGqmjqcP4yj17enqgUCgEra5ccPfPfX1R6/8BYBgGBoMBLpcL7e3trOIoFJRKJQoKCrC0tIT6+npYLBYYDAb2R0faM7kTKCmKws6KLLzcuyj4c60WvDSD2w/24cFP1mFHaRbv48S0IAp+KIYL1WdnZzE4OAiVSrUiVD/V52gBCVipSb03LS2NF4HCkZqvOow4ioYDkXuWlpYKNkQg98Fd4YU6WdI0je7ubtA0jbKyMt6NLkDgOBhguXJAVhW73c4Oe8vOzsauyuxTmtQA4PL58eUnu3H/5RXYuaWQV/5AzJE70XIfXIVYqKx6ZmYmXC7XCpUhFxsu/GYYBidOnEBFRQXvZEGolVaIOkwqlYZdqYncs76+nvdQeS6If7dUKo1phSfhOgn1+CJc9jt47hZZxaemppDmEseqJ9Fw+hjc+eIYvtE6g0rNcu9yTk4O1Gp1SPKu9sgdgnBZ9aGhIYyMjGBycpIN1bnONSKWtC4G8AAAKYDfMwzzU77Hih5+k+FvQo7hQqg6LNRKz0fuGQ1OpxOdnZ0oLi6OaYQKadcjWfqpqSlR+2wpimIbO4DlCKmi8zjGTLELeVYLNg+DB7r8+O0nyiCDE6Ojo2xLKYk8SPuiGL3UQPwPBxKqp6Wloby8HHK5PGCovUqlwsjIiCimg9TyqJ1fA7gAwDSAIxRFPc+nmQNIQPjNt4spFIjvlxB1WPBKzzXkj9a6GQ5k/xw8cocvpqamMDMzExCu82kEIYjFeVShUODsLXkYe39a8P2uBYx2L245OIjHrmtCfVBL6dTUFACwuQMxILaiLFSo/tZbb2FwcBAf+9jHcNZZZ+G+++6LqWwGoB3AMMMwowBAUdSfAewFj2YOABDHI4aDWEjEMAz6+vqwsLCA9vb2uOSeR44cQWlpKaqrq2OuPw8PD6O1tTUioUORzu/3o6+vD0ajEW1tbQH7bzEtgsNhV2VWQs8vNhasyw6lOqubbSmtrKxEa2srGhoaoFarMT8/j6WlJfT09GB2dhYuV2wN5IlUlJFQ/ctf/jKys7Px3nvv4ZOf/GTMk0Sw3I01xfnf0x/8Gy+seZeW0+mEw+FAcXFxwPxoviDhN5F7NjU1xZSoELJ/JgTl3qvH40FnZydycnJCfg4xLYLDobU0C0qZBC7f+thfA8D0kgs3PdGNP3y6AZrUD7dJcrkc+fn5bJa/vLw8wGCBKwXlQ9bVms7h9/uhUqnwkY98JJ7LhCIB7x/EmpKaqMOUSiVKS0tj7uiZmppiy1VCJKMkJHa73YL2z8Gks1qt6OrqiliPjyfPwBcpMglaSjPxzqgppuPXCqN6B25+shsPf7oRGcrA74/Ul4MNFrhzu4mH21on3ESMxKYBcH+IJQB4d+2ITmo+P8hgdVhnZ2dMTewulwvT09PIyMhAQ0ODYDJIJBIYDIaQI2ujHUe+wIWFBYyMjKCxsTHi/k/InhqI/QeyqyJ73ZEaAAYW7PjSn7vxu2sboFZ8uPKGSpQFz+12uVwwGAwRE25ikTqaU6hITqJHAFRTFFUBYAbAvwO4lu/Bq75Sh1KH8fH+DgZJZuXn54dsiogGMiaI7J+F7H8oigJN0xgfH8fS0hLa2tqi3nvw/c3MzGBkZATp6ensD5Rk6eP5UeyqzAbeiPnwNUXXjBW3PNWD31xTB6V8mdh8TAeVSiVrFhgu4ebxeBJu20t83UU4j4+iqFsA/AXLJa1HGIbp5Xv8qpI6nDpMSFMH1wy/tbUVJpNJcOsm2T/7/X40NjYKTmhQFMXLUCH4GDKQb2hoCA6HA21tbXA6nTAYDOjp6QHDMKwOOdYRMdV5qchLU0BnO/VLW6FwZMKMrx3qwy+v3g65VJiuHwj0cAOWm2fIdI4TJ07E7GhKEIm0Ho+HtxUXj+u8DCBqR1YorEr2m6iyBgYG0NzcvEKYwpfUPp8PnZ2d7BxqpVIp2L/b6XTiyJEjyM7ORmZmpuCnt8PhwNLSErRaLbZu3cr7eELqEydOsN1hUqkU6enpKC8vR0tLC5vxnZ2dhd1uZysCXq8wTffOU8i7LBb8a8SEbz07ANrPxF2nJgk3pVKJtrY2lJeXw+fzYWBgAEeOHMHJkydhMBhiGt8UjFNBTQaswkpNdNdyuTysOowPqYkohYzAFXIsAWkwIftng8Eg6IFAQv6MjAzB7XVutxszMzOora1lFWHBID/AvLw82Gw2lJSUwGAwYHp6ufZMwnSu7jsUdlVk4/muBUH3d6rh9QE9vvfiIG5sz8IBtxuyuTl8NCMDrXGQRiKRxJVwixY1iOl6Eg8SSmq+6rBoxFxYWMDw8HBIuSeflZroyOfn5wP2z0IM/aempjA3N4fW1lYMDg4Knu4xMjICrVaLwsLCqO8nyRbiOV1RUcFa+ExOTsJms7HJII1GsyLjf2ZFFigIqIGEQKlpDpvMC1hI02BYWxrHmWLHc3067M+1w6MAfE4n7l9YwINlZdgTgyAoFCIl3JxOJ5vvIAm3VTRIiAsJy34LUYdFauogg/Ha29tDJqOiPRDI/pk4nHCftHxITQQlDMNgx44dkEqlgurI09PTmJ6expYtW7C0tBTwGvFIJz+USKtAsIUP0X1PTk6y3Vs5OTnLe8ZUBUrTKUxYY6P1hUPv4Pa3n4QfFCSMH/ubPob9LZfGdK54YC+UYYliIPH6IaEYeOVyfHt6WjRSByNawi0jI4PNi4SKlDZs+E3TNPr6+uB2u3lPhQxFTK7cs6WlRbC7J/ChfruoqAilpStXm2ikJvXrvLw8lJWVsffANzoYHByEy+VCW1sbLBZLwGtkbIvf7wdN0+wPhY/tDlf3XVlZyXZvjY+PsyWdWk1spE51O3Db23+GVaGCVyqH1E/jUx2v4M2qNsxmrI57CYFfBjASwM9IQHtdkMvlsIuw9+WDUAk3MlTx8OHDbMItJyeHTY6JvVJTFPXfAPYA8AAYAfA5hmGWoh0nOqntdjvUarUgdVgwqc1mM3p6eiJO2OAeG4pgwfvnUIhETovFgu7ubmzduhVabaBjR7SVmiT0MjMzA+ySCIFJmYY88MjT3+/3s95lfr8fPp8PEokkaqKI273l9/thsVhQm7WIV2OYn5fudoAC4JUuR0W0RApaIkW207rqpFYa/aD8ACUBVGo1JBSFCyO0PSYScrkc2dnZsNvtqKmpYX3J+/v74fV64XQ68e6778YjDQ2F1wF854MS188AfAfAt6IdJDqpMzMzBYvYpVIpm+Gdnp7G1NRU1AkbBMH91OH2z+GODUXqubk5jI2NhZWcRnoYEHfRiooKFBQUBLxGSBssUCCklUqlsNls6Ovrw+bNmwEsRz7k85FVPBLJSQ92bW4K1AovHB5hK5s+NQtLyjRkumwwK9Og9jjhk0gxlRnZOjkRSLH4oe1yw9eohFQqxQUZGbi/rEzwecRSepE9dbAvOU3T6OjowIkTJ9DV1YWOjg7cddddOOuss+K6HsMwr3H+53sAruJz3Jprv4HlH7PD4UBPTw9omkZ7eztv8T2XYCT0pyiKV/9zMDnJHt5ms0XcOoRbqcmTO3hcD8MwUCqVsNlsOH78OLRaLbRa7Yr5xwaDgR33Q8I47ipOyE3TNEvucJ9RJqGwozQT/xw2RvwbBMMnleGui76C/3z9t8i3GWBWpuG/zrsJFuXaJIDUizTucfmRCxfS5XJYFxch12gECZUS7fktlUrR2tqKj33sY7jkkkuwd+9e0eZpc3ADgP/j88ZTgtSkIaOiokKwBpyE7i6XCx0dHYImb3JJTQbCp6Wlobm5OaoUMHilDtVuCSwTmniotbW1wePxwGAwsL23WVlZ0Gq1cDgc0Ol0aGlpCej/5q7icrk8YA9O3EjJ68Gr+K6KLMGkBoBxTTGu/8R/QulzwyVLARKsxIqE0kwZzthWDq1WyyauSImP7GmjTelYLaMFm80WNn8TDueffz7eeOONnhAv3cUwzHMAQFHUXVgemLefzznXnNQGgwHDw8PIzMxEWQyhlUQigcfjwbFjxwT3PxNSk5C5vLycV8mJq/32+/0YHByE2+1eUYcnCTFyDLC8/+WOpTWZTDh58iScTicyMzOxsLAArVYbdgvDXZ3JKk726cCHqzjDMDgzHhEKRcElF3V/GBPO3KQKWeIjSjEypSMtLY0t8QWbYqzWdA6HwyE4+/3Xv/4VAOrCvU5R1PUALgNwHsNzH7EmDR1A4Oypbdu2YWFBuFiC7J/dbjfOOusswUkKiUQCi8WC8fHxkBMuIx3HMAw78iczMxM1NTUBxovh9s9c+P1+dpZYZWUlnE4n9Ho9m3zRaDTQarXIzMwM+aPkruLkfDRNw2azgaZplGTIUZiRgjmLW9Df5VTCGUXKkJ+dCHWISYHNZmPltn6/n5WCkiHxqzGdw2aziZ39vhjLibFzGIZx8D0uISs1n+xwd3c3lEolduzYAYfDIVimR/bPAKBWqwUTmmEYGI1GWCwWnHHGGYI0uxRFsXLTysrKgIQYX0K7XC52/jWJDtRqNUpLS9nki8FgwNzcHAYGBpCamsruxcPZM0kkEpjNZgwMDKChoQEKhQJnVmTh6c71qS7blK1EeZaMV4kvPT2dldz6fD4YjUbMzc1hcHAQKSkp8Pl8cLvdcWmzozUdJcBJ9P8BSAHw+gd/g/cYhvlitINWPfwOJfcUOqUjeP/87rvvCroHv9+P3t5eeDwelJSUCP6inU4nFhcX0dzcvCIhxofQZrMZfX19qK2tDTuuh0xNzMvLY1civV6Pzs5OAMuSUa1WG7CfnJ+fx+TkJJqbm9mH3FnV2nVL6vNrtPD7vYJDZ5lMFvC30+l0mJqaQl9fH3w+H7Kzs5GTkxM2AgoHmqYjLh5ii08YhqmK5bhVJfX8/DxGR0dXqMyEkDrY0E8o3G43Ojo6UFBQwCaohGBqagp6vR5lZWUrCM0VkYQj9MLCAsbHx9HU1MS79MddiYhklHiBE8mo3++Hx+NBS0tLQNZ+Z0U2pBQFOsFWSonABVu18Ftm4toPUxSFlJQUZGRkYMuWLfD5fDCZTKzTLPH2zsnJiRrt8RljG6/poBhYlfDb7/ezpaJQvcd8SE3017OzszH5bwMfilpqamqQk5MDnU7HW8PNnbJZVlYWNSEW6v7Hx8dhMpnQ0tIiqCQTDIVCwYpNyFbG5XJBIpGgq6uLDdPVajUyVXLUFaWjc4bXaONTBgUZKagrSke/OX7fb26iTCaTITc3F7m5uQHe3lybpJycHGRlZa34Lld7Tx0rVmWWVmdnJ7Kzs8PKPaPJLrn753CdXtFM32dnZzExMREgauHb0OH1etnPUFtbi5mZmYCaMZ+EWH9/PyQSCZqamkSrYRJCZ2Vloby8HBRFweVyQa/Xsxn57OxsNBcq1x2pz9+qZUuH8f69wp0j2NubdG3p9XoMDw8jJSWF1dSrVCpe2e8Nu1ITkJUxkncXEDljTiZUFhYWhq0/R5qnFWxKwH0PH1Lb7XZ0dnYGDKQnPzY+hCYTS3Jzc3nXz/kgVKINWG5KKCkpQUlJCWiahslkQvXiqTWUng/Or1mW5orh+803+x3ctRU8oYOmaWRkZECtVoc8n9PpjNUSWFQkjNRTU1OYnp7mLfcMBbJ/rq2thUYTfrBaOHJyS05NTU0rCBXNo5z0Twe3fHJJHakBw263o7u7G5s3b474UBMKm82G7u5u1NTURMwrSKVSaLVa7NmlwU/f/hes7tVphogXOalyNG9a/nvzsTOKhlhXe7VaDbVazT4gT5w4AYvFgtnZWcjl8oBVnGw5E6AkA0VRdwL4bwC5DMPoo70/IaQeHh5mpZax1Ae5++doA/aA0KQmK2xwySnacQRkjlfw/GmGYaBWqzE2NgaTyYTc3FxotdoVe3yj0YjBwUHU1dWJGpKRlYNMzuQDmUSCMyqy8deBqL+HUwINGgZTk5PIycmJyZAyGGKE8FKpFFKpFFVVVZDL5XA6nTAajRgeHobL5cJTTz3FavfF/L4pitqE5Ukdk3yPSQipy8rKeLUQhgK3fznaHC2C4EQbsR5uaGiI+AcORWqSEPP5fGz/NAHJcKempuKMM85gS1u9vb2gaZotM1mtVszNzaGlpUU0zypgOS9Aoh+h591VqVk3pL5qZxUUimVjR4vFguHhYWi1Wt4zpIMRbaY0X3DDeJVKFdB7rdfr8frrr+O8887D1q1b8fjjj8d9vQ9wP4BvAniO7wEJIbVCoRAsJqEoCg6HA93d3SgoKBCkASfkJBlmvV7Pa4ZWMKm9Xi86OjqQk5ODioqKgOtz98/kqa9Wq1FWVoaysjJ4vV4YDAb09vbC5XIhPz8fZrMZOTk5cauZGIbB6OgorFYrWltbYzrfRwTOhV4rZKpk2FWVC5mEQmFhIY4cOYL8/HyYTCZMTEwE7HvDeXwHg8/ESz4Il4yVSCS4+OKLcc899+Dw4cOw2WxxXwsAKIq6HMAMwzCdQhbIhJW0hMLv9+P48ePYtm1bxP1zKJDWze7ubshkMrS2tvIKt7iNGUQUw02IAfwFJRKJBAsLC8jLy0NlZSXMZjP0ej1GR0eRkpICrVaL3NxcwaU4kjmXSqUBvdlCUZKtQqlGhUmjMOfV1cZHq3MgkwR+RjJdElhOEBJPMaKXJ5ZD4R52Ymm/gfC/bZfLxW4ThZS1IjV0APgugAuF3uOaN3SQ/bPT6cSOHTvCKqwigYTsZWVlgiZUkpWalICCw3Whks+SkhJWJUeGmVdXV8PhcECv16O3txc+nw85OTnIzc1FRkZGRJKSB1VOTk7ME0y42F2pwaRxJq5zJBoX1GhX/Bv3cyuVyoCGGLPZDIPBgLGxsYDkFTc5KyapwyFWNVm4hg6KouoBVAAgq3QJgOMURbUzDDMf6ZxrSmpCRr/fj5ycnJgEGUtLS1hcXER5ebngkbMSiQQOhwMjIyMhE2J8FGJWq5UVtITLRHM13T6fj/W9slqtyMjIQG5u7goDQVLKKy8vD4gc4sHuzRo8efTUJXWqQoodJWm8HV/IeFnydyce6iR5ReSgYiTbokFsJ1GGYboBsFYzFEWNA9ixZtlvPisK+dGS/TMxSBACMhSvsLBQ8B+UtEx6vV585CMfCfjS+RJ6cXERIyMjaGho4H19mUwW0F1kNpuxuLjIrjS5ublQqVQ4efJkRG14LGgqVENKAfQpqhg9u0oDVYo8oFeclA6jOb4Ay8krbo2eCEkWFhZgs9lQUFAAjUYTUy05WtfjqWI6CKzRSh2q/ixE/83tYW5vb8fk5KSgBwJ3QiUZCEAQKiEWDNLyqdfr0draGrPkk6IoZGVlscR1Op0YHx/HyZMnoVQqodcvP5RjGToQDLvdjqG+btQVpqJz1h7XuRKFC7flsX9LmqYxMDAQMK2Ej+MLATeh5vP5kJeXB6fTyT7IhTZ1RKuXJ9oemGGYcr7vXXVST05Ohqw/8yU1V7JJepiFTOkgCbGqqirk5eVhfn55eyKkB3pgYAAMw6C5uVnUsM5gMMBms7GRg8FgwMzMDPr7+5GRkQGtVoucnBzB5RnSFVZXV4dzvAZ0zo6Jds9iQSmT4KyqZSUXSQ4qlUps3rw5QOxDVnHy36EcX4JBypBarZaVgwY3dZAHQLgsOZ9mjlNB9w2sYvjN3T+Hqj8LmdJBCCnkWCB8/ZovoUniSqPRBFgGxwuGYTAyMgK73Y6Wlhb2b8MN0y0WCxYXFzE+Pg65XM5m06OFkouLixgdHWW7wnZXMvjlm6ceqdtL06CULZO3u7sbmZmZKC8vZ1+P5vjC3YcHEzw4UUaUdlqtlm3qMBgMbGsm6drKyMhgj+Nj5H9ahd9k/5yfnx+WDNGIqdPpMDw8jIaGhhVPRIlEEnHeFJnlpdPpQtav+eyfSQ29oqJC8MidSCC93SkpKWHH8XJ9vquqqgIcUjweD5tNDw7TZ2ZmWBEMCWu3F6UjSyXHklPYfK5EozlXgvfffx8ejwcajSairVQ4xxfy/4PdVyNlv7lNHSSRaTKZMD8/j8HBQaSmprJS0PXQoQWsAqn56rcjTekYHR2FyWTCjh07QgpKwnl/AysnbAQnxPx+P8bHx5GXlxf2Sbu0tIT+/n5s3759xdifeEC06Xl5eYIy9yqVCps2bcKmTZtYlw8Spqenp0Or1cJms8Fms6G5uTngxyihKJxZmY1XenWifY54IZdS+MTuWgz3LwuPKGp5qiixJYpW/gtexbn/R2ynyHcdbbsU3Jppt9vZaoXb7cbo6Ci7inPvJxHhN0VRtwK4Bcumgy8xDPNNPscllNTEYZOPfpvr/U1A0zS6u7uRkpKClpaWsF9IsPc3AUmIabVatjWRgPRANzc3s9a8Lpdrxao3NzfH+pCLadROpodUVlbGtfIHu3wQOyOXy4X09HTMzMywvdUEuyo1pxSpzyjPwnB/N0pLS9nyXXl5OWsuSMp/5IEVqfzJJThN0xgbG2OtmElGnaKWZ6JHIzjX3zszMxPz8/NIS0vD7OwsBgYGAswOHQ6HYNFUlGv/G4C9ABoYhnFTFMX7R5IQUjMMwxrACdFvu1wu9n+THz0pUURCKFKT/Xd1dXVAh1Tw/pk7P4mm6YBVj5y7qalJVA03SVxt27aNt9khH/j9fkxMTCAvLw8VFRVwu90BvdVEm35mRZZo1xQD1Uo7ystXtucGmwtaLBbo9Xp2flg4/3Rg+XseGxuD2+1GfX09G4ZzQ3QhQxKIP1koi6mf/exneO6559DS0oL6+nqxeua/BOCnDMO4P/g8vJ/CVJT6W0wVTYZhMD09jdzcXN7JpMXFRZhMJmzZskWwZdHS0hJmZmawfft29lzxJMRomkZPTw8kEgmUSiUMBgNSUlLYsCwegut0OoyOjqKxsVHU3ltSFSgoKAj5ECQPrMXFRZjNZnz/PS9mrGvfiimhgOdv2I7KYuGjgfV6PfR6fYB/ukajgUQiwfDwMLxeb9jxT6GGJAAIWzLT6XRwOBwByTsu7rzzTmRnZ2NmZgZ33nknGhsbhXycFTdIUVQHlps4LgbgAnAnwzBH+JwsYdlv8nTlC7KnDmeKH+1Y8iWFS4jxJbTb7UZXVxcKCwtZchCp5+LiIrvXIwQPtUqEw9TUFHQ6XVy17VAgicjKysqwfdtSqTRgr3iWvg9/PrH2IXhzcbpgQgPL/uncLimiLBwZGYHX64VSqURdXV3EfTgQfkhCcMksWvabpmns2bMHu3fv5v0Zzj//fMzPz6O3tzdY+30XlrmZDWAngDYAT1EUVcnH+3vNtd8EFEVBr9fD4/HwDtkJJBIJfD4fent7ASBkQoxPhttms7GD+Yj7BQG3Iyt4ygZJ5oTytSLXP3nyJNxut+i1bWKYIER9RlEUPlpTcEqQ+pL66MMTokEikUCj0SA7OxuDg4Pw+XxIT0/n7Z9OzhGpZObxeFg9RKhzxFLS+kD3DYTWfn8JwNMfkPgwRVF+AFoAi9HOmzBSR/P+5sLj8WBwcBASiSRsWScSSGi5efPmFSUzPqaAAFhfqvr6+qhfDtf4z+/3w2g0YmFhAYODg0hLS2ONE2QyGRvKp6amRlw5YoHJZGIbUYT+oNrKs6CQSuCh+Yl2EgEKH9oWxQuGYTAwMACpVIrt27eDoiiUlZUFeIAL8U8HPiyZWSwWzM/PY9u2bWH34glwEn0WwLkA/k5R1BYACgC8GuLXfKW2Wq3o7l7OfOp0OsE/eqvViq6uLiiVyoD9Dt9wG1gOixcWFlbMseIDbsKGYRhYrVYsLi6yyRyXy4Xi4mJUVFQIOm80LCwsYGJiAk1NTTFl5VVyKVpKM/HemEnU+xKC5k2ZyE0Xp8+5r68PCoUCVVVVAd91cHWAj386F3a7Hb29vaw+IpzwxWg0il3SegTAIxRF9WB5PvX1azZ2RwgWFhbYhgiFQoG5uTlBxxNBSn19PQYGBth/FyL5HBoags/ni1gy4wuK+nDeU2FhIWu4YDKZsLi4yJbLog10iwayN29ubo5rb767MntNSX1Bbfy+bQzDoLe3N0BSGg5UFP/0zMxMNtkmk8lgt9vR1dUVYB0VSvjyzjvvYGRkRNRtFcMwHgCfjuXYhIbf4UBkkWazmfUB54Y10cB1ONmxYwd7PHmND6GJvW5mZia2bt0qalhMxCrcoQXEGYX8gLKystiWS74/BvJ3czgcouzNd23W4OdvjMZ1jnhwfk18pCZqvNTUVFRWVgo+PngbRYwtxsbG2Lbcbdu2RVyBjx8/jm9+85t47733RDWXjAervlITMqnV6gAfcL5NGeSLlEgkKxxO+CbEnE4nurq6UFZWFtaUMFaQCRzBYhW5XI6CggIUFBQEZGtPnjyJ1NRUdh8ebuUlTQ4ymYytu8aLmvw05KQqYLB74j6XUGwvTEdxVuxiHr/fj56eHnbVjRfc3myHw4ETJ06guLgY09PTGBkZQXZ2NuuTRn5zHR0duPXWW/H000/HNLE1UVhVUpORsWVlZaxDCAGfH6nH48GJEyfCepjxSYglSvhB2jENBkPUCRwkW6vRaFgp4uLiIk6cOAGJRMKWnogKjKZpdHV1ITs7W9RGEoqi0F6ajlf6DaKcTwjiCb3DNX2IAYfDga6urgCNA+nqIvqHEydOQKfT4YUXXsBzzz2HzZs3i3oP8WLVSE08tIWMjOWCJMS2bt0KrTYwY0qM/Ds7O5GXl4fc3NyQCa/5+Xk2uSSm8INhGAwODoKmacFqIq4UkajAFhcXWRVYVlYWTCZTyAdhvFhaWkKxdG0md/xbVVZMx/n9fnR1dUGj0Qga7s4HJILbtm3billv3GTo0tISnnjiCWRmZuL666/H/v37Ywr/E4WEKMoAsDpb4kE2NzeHxsbGiJnad955B7t27Vrx7+E6tIL3z0Qgsri4CIqikJubi7y8PCiVSoyNjcFsNqO+vl4Uu1ju5+zu7kZ6ejoqKytF3ZtbrVZ0dHQgNTUVbrcbmZmZ7D48XodSUsLbVL0NF/76WOxfdAwoy5Ljv85cznqTFlI+Ih4SsZC+aDFBZMm1tbURF52hoSFcd9112L9/P+rr62E2m6FWq+MVE4n3o0GCV+pgD2+hiR2i3zUYDCs6tEIlxEgLXXl5Obvi9ff3w2KxQKVSoaamRpTh4wREfVZcXCz6KmqxWNDb24vGxkZkZGQEWB+NjIxAqVSyYbrQMtz8/DympqbYEt6W/DQMLohja8sHlzYUo61tOfus1+vZ5F92djZyc3MD9q0ENE2zkVi0XgChIIq8aIQeGxvDddddh8ceewz19fUAIOoWTiwkjNRutxvHjh1DXl4e730gt/eVJEJCWf7ySYilpKSwziZlZWVQqVSYmJiA3W6PqgDjg0jqs3hBusYaGxvZfTXX+qi6uprdh3d1dYFhGN4rHrccRiKW3ZXZq0rqCz/YTysUigBnUO6+Va1Ws8lDqVTKjh8uLi4W9V7IrPOampqIBJ2cnMS1116Lhx9+GM3NzaLeg9hIWPhtNBrhdDpX7H8j4fDhw2hubobf70dHRwcKCwtX7Ju4hI5ESDLHqqqqKuAeiAJscXERS0tLSE9PR15eniDTfTJSR8joG76Ym5vD9PQ0Ghsbea/AZMVbXFyE0+lkH1pcWSSJeqxWK+rr6wP+du+OGvH5P3WK+jnCoUyjwiu37Iz4HpI81Ov10Ol07IO4srISaWlpom1xuISOJLGdmZnB1Vdfjd/85jcht4ciYH2E31lZWYKli1KpFGazGYODg2ETYtFMAYEPV7q6urqQLincpIfFYmE7p5RKJfLy8iJKCBPVXw2AnV/NXUX5gLviEckskUWSHmSTyQSGYULKcFtLs6CSS+D0Jl4yegGP2jRJHhLzxerqakilUoyNjcFut7M1/kgG/tHgdrvR0dGBrVu3RiT0/Pw8rrnmGjzwwAOJIrToSNhK7ff7I1oMhcJ7773HGhdwHwhCJJ/T09OYm5tDQ0OD4BZJu90OnU4HvV4PiqLYTLpKpWIFL0tLS6In28i4Xa/Xi23btommTCL7cGJ7lJ6ezu7Dgx9IX9jfiX+NGEW5biT83+dbUV8c3T2GjEDiGicACKjxm0wmNreg1Wp5f99utxsnTpzA1q1bI7b26nQ6XHnllbj33ntx3nnn8Tp3jBB1pU4YqRmGgcfDT9RALIsmJyfR2NgY4CDB7XmNNHSPdEK5XC5s37497oSYy+ViM+nEDkelUqGurk7UZBsR06hUqqgyR6Hg1rfLy8vZgX6Li4ugaZrdh6elpeHx96fx09eGRbt2KBRmpuCNr0Zf7bxeL06cOIHy8vKIrjDENHBxcRF6vR5+vz/gM4X6WxKtQ3V1dUSnEr1ejyuvvBI//OEPcfHFF/P7gLFjY5Gapmn09vZCJpOBYRgUFRWxT0++CjFSVkpLSxOdGD6fD52dnWzJQqxEGzk3KdGIXXMlpgmFhYUhk0ter5fdh9vtdtgkafjqXxI7FfO6M0rw7YuqI77H4/Ggo6MDFRUVgmWX5DPp9XpWiku03FKplCV0VVVVxOSmyWTCxz/+cfzHf/wH9uzZI+geYsT62FPzAdnXkIQYEXAA/Akdao6VWCDnLi0tZeWkJNFG3CbJ2Byh0y3dbjerrhNrrE7wuSOtdHK5PED3bDKZkPNPIwzOxO2ro6nICOk2b94sKMFKEPyZyISOkZERKBQK2O32qNUKs9mMq6++Gt/61rdWi9CiI2ErNbD84woHi8WC7u5u1NTUsH/k4eFhNrHDZ/9MarmR5ljFCj4zsri1Y4PBAJVKxe7vImWuSWZ+y5YtoprVAR+KKGI59388P4CnO4R1yvGFNk2Bv39tFyQRHGc6OjqirqKxwOPx4NixY8jMzITT6WRniQd3zFmtVlx11VW45ZZbcM0114h6D1Gw/lfq+fl51mCemxAjjqJ8CM31+uI6ZYoBkj2PZj7ArR1XVVWxteOOjo4A+yCuJJU7LUPkpnr2QRSrlfHuzdkJI/X5NdqohI62z40FJOFWXV3Nrv7BHXMmkwl6vR5PPvkkvvCFL6w2oUVHQkkd7H7C9fAmLZfc1xQKBSYmJlgPsFCrHfEhMxqNont9AcDs7CxmZmbQ3NwsKHserOEmibb+/n74fD7WCWV2dlZ07TmwrOMeGBiIyQWF4MwKDSQU4E+AZjRcmyWpFUfLRMcCknCrqKgICOeDO+YOHz6M//mf/8HMzAyeeuopbN26Fe3t7aLey2oioeG3x+NhSU1sfRQKBbZu3RpWIcbVb0ulUraspFQq2fZDiqJQU1MjdlM6RkdHWXGGmBlur9eLkydPQqfTQaFQICcnB3l5ecjKyhIlqUf2jdG09Xxwze+PonvWGvc9cZGpkuGtr++GLOj7IluFaOKPWEBW6PLy8ogJN5fLhWuvvRb79u3DzTffjJGREahUKtGVazRNY8eOHSguLsaLL74Y/PL6C7+JtraoqGiFED84IRa82ul0OnZYu9frRX5+/grLmnhBNOoymQyNjY2inps0tHg8Hpx11lkAECAOycjIQF5eXsxNGkSB1tzcLFgDHgq7KjWik/rcrdqwhI6mt44FPp8PHR0dKCsri0hoj8eD66+/HpdccgluvvlmUBSFqqoqUe+F4IEHHkBtbS0slsR3xSV0EjdFUTCbzTh27Biqq6tXENrv97OdXKFq0EqlEqWlpaitrWWTGzabDYcPH8bIyAisVqsgG+JQICFaeno6O0VTLBAzPLfbjcbGRkilUnavvW3bNuzcuRPFxcUwmUw4cuQIOjs7MTs7y1u0Q7rfxCI0sDyYXmzUZXixtLTEflcOhwMdHR2i97QDy4Q+ceIESktLI9a4vV4vbrjhBpxzzjm49dZbRf3egzE9PY2XXnoJN954Y8KuwUVCV+q5uTmMjIygubk5IJklRCFmMpnYPmyS/CH1yLGxMTgcDjacjTRvKRRI/2x5ebnoZSWy3SDOHOEG35FEG9E763Q6nDhxgiU/aR3lgmwVbDabWNMgWDSWZEAtl8AhkmQ0LUWKc2ryMTs7i/7+fqjValgsFtTX14s6lwz4cIUOVqGFet9NN92E1tZWfP3rX08ooQHg9ttvx7333gurVdwIKBwS3nrZ3t4eIKkUQujZ2Vk2tAy2BiL1SJqm2QFmVqsV2dnZ7H410o+dlMOE+GXzRbRpGaHA3XpUVlayibbe3l5W/ZWXlwe1Wo2hoSH4/f6Y7JSjwbCow9ZsCmJZgp9TrUVRQT6KCvJhtVrR2dkJjUaDgYEB0aaeAB8SuqSkJCKhaZrGl7/8ZdTW1uK73/1uwgn94osvIi8vD62trfj73/+e0GsRJDRR5vP5AswE+RKaO69ZiCyTiCh0Oh2WlpbC7ldJT3JDQ4Po5TA+0zKEgkQmOp0ORqMRarUaW7ZsES3RRjAzM4P5+Xn0e3Pxk9dGRDnnA1fX4YLaXHboALezjZsUjXXqCbBM1I6ODhQVFUUcgUvTNG677Tbk5+fjnnvuSTihAeA73/kOHn/8cchkMrhcLlgsFnz84x/Hn/70J+7b1odMFAgktRDJJ9FCx5MQI8IQnU4Hg8GA1NRU5Ofnw+12Q6fTsbbEYiKWaRl8QXTcmZmZSE9PZ2diieWGQvzVGhoaMLXkxiW/fj/ue1bJJfjXnR+Bz+VAT09PxHIbseslM6v4SnH5Etrv9+OOO+5AWloa7rvvPlG3LHzx97//Hffdd9/GyH7znZJBnESKioriLikE71etVisGBgbY1j29Xo/c3FzR6tzxTMuIBhLOk/ZKAOxMLNKxNDw8DLVazbaOCvlcY2NjsFgsaGxshEQiQXmOGiVZSkwvuaIfHAEf2ZwDr9POOrhEioqETD0hIG4o5Lhw8Pv9+Pa3vw2FQrFmhF5NJHyl9ng8vPbPVqsVvb29CZFOEhcVpVLJDrvT6XQBtfC8vLyY93VkWkZDQ4PoPdZ8dNwAAhJter1+RY0/3DHDw8Nwu90rWj7vfmkQTx2bjeve776oHKXMQlzbHO7UE71eD7lczkYmg4ODyM/Pj7gA+P1+/OAHP4DVasVDDz10qhJ6/YTfjz32GCorK9HU1BQxNCR7XD5zrISCrHL5+fkhzepIO6JOpwPDMGzGme+PkNgDNTQ0iK5uI3a1sTzoQrVZ5uXlsftV4oDKMEzIUt7r/Yv46oHgYYz8IZdQ+MVHU7CzVVz1nNPphE6nw9jYGKRSKYqKisJOPWEYBj/60Y8wNzeHhx9+WFRBkchYP6R+5pln8MQTT2BwcBDnnnsu9u7dG2BASIQZi4uLqK+vF32PS0ixefNmXkkrj8cDnU4HnU4Hr9fLEjzcUHNimFdXVyf6ChCvjpsLr9fLEpzYHZEpjdXV1SEjKKvLh933/Qu+GDWjDVoJHr3hDNEjF7/fj87OTuTm5iI/Px8GgwGLi4srpp5QFIV7770Xw8PDeOyxx0Q1tUgA1g+pCZxOJ1599VUcPHgQnZ2dOOecc3DppZfixRdfxL//+7+LMscqGPGa9nMzzk6nM6AWzjAMOy1jy5YtomdRiY47UZFLR0cHaw3FJULwSvapR47hxHRsCqgffGwzrmkTt0ecEDqURTDXEeX+++9HX18fJBIJXnjhBdE1CAnA+kuUqVQqXHHFFbjiiivgdrvx7LPP4uabb0ZeXh48Hg9sNht2794tWvhKwrN4GifC1cItFgurbgu3ysUDouOOdZplJBBBDHF4JYk24qsePP5n12ZNTKSWSShcVBf/3GkuiIl/OM9v7ozq+vp6LC4uYvfu3bjqqqvYbaBYmJqawnXXXYf5+XlIJBJ84QtfwFe/+lXRzh8vVmWlDsb3v/99NDY2Ys+ePXjzzTdx6NAhvP3222hvb8e+fftwzjnnxByKT05OYnFxMSF7XOLKkZWVBZ/PB7PZzNbCc3Jy4o42YnES5Qvi4JKfnx9SEEPGvJJEm1wuh55Jx+0vTQu+1pkV2Xj4M00i3PUyyJid7OzsiA4xDMPg4YcfxmuvvYZDhw7FLWgJh7m5OczNzaGlpQVWqxWtra149tlnsW3btlhPuf7Cbz7w+Xx46623cODAAfzjH/9Ac3Mz9u3bh3PPPZfXikXM+zweD7Zv3y56OE/251znDO5KR+YTE4IL3cOR3EJDQ4Po+z8ScpeUlEQs/XDhdDoxv6DDJ54chV2YfyS+f8kW/PsOcbqcCKGzsrKiDqF79NFH8dxzz+G5554TPcqJhL179+KWW27BBRdcEOspNiapuaBpGu+88w4OHjyIv/3tb9i2bRv27duHCy64IGRWmoSVqamponuUAR9KSiMlrUjphax0xG44Wi2c6LiJek7shxGJLqKVxMLh9gM9eK1/kff7JRTw5td2ITct/lWSlCIzMjKiDsLbv38//vznP+OFF14QXSUYCePj4zj77LPZ+4wRG5/UXPj9fhw5cgQHDhzA66+/jqqqKlx++eW4+OKLkZ6ezopKxBCshALXBUXIj8Vms63oCw+uhZOykt/vR21tregPIyJZjcci6MDxWfzgxUHe72/ZlIk/fa4lpmtxIWRU7YEDB/DII4/gpZdeEn24QiTYbDacc845uOuuu/Dxj388nlOdXqTmgkzuOHjwIF555RXk5ORgbGwMf/zjHxMyCkWsPS6prS4uLoJhGFb1NTo6mhBrYHLNzs7OuB1FZpZcuOCX7/J+/5fPzMfNH62OK5/BMAx6enrYvvpIePbZZ/Hggw/ixRdfXNW5Vl6vF5dddhkuuugi3HHHHfGe7vQlNRdvv/02Pv/5z+Pcc8/F4cOHodVqsW/fPlx66aWiGNeRaRliG/e73W4sLCxgdHQUUqkUxcXFYWvhscJut6Orq0uUGjcAXPLr9zBucPJ678P7iiBxLrHKr0iKtlBgGAa9vb1Qq9VRM9YvvfQS7r//frz00kuiWyFFAsMwuP7666HRaPCLX/xCjFOuv5JWovDXv/4VJSUlbBh78OBBXHXVVcjIyMDll1+OPXv2IDc3VxBZuNMyiBZaTEgkEuh0OmzZsgW5ublsCSu4Fh4rwYloRcw5X7s3azBumIn6vvqidJzZsBXAh9FJT08PO8CPPLzCgWEY9PX1QaVSRSX0a6+9hvvuuw8vv/zyqhIaWF5QHn/8cdTX16OpqQkA8JOf/ASXXHLJqt5HOKzblTociNLr0KFDeO6556BQKHD55Zdj7969KCgoiEgWso9Tq9UJCYmJa2ZFRcWKpBVN06zYxWazQaPRCPYx45oPipks+vuQHl/+c3fU991xXiVu3L0yQ00G+Ol0OrhcrpAPL0LolJSUqH/7N998E3fffTdeeumlmJJ/pyCS4TdfMAyDyclJHDp0CM8++yz8fj/27NmDffv2oaSkJOCHQ+q4eXl5og80B4TpuEmXkk6nY9srSV94uMjBaDRiaGgIjY2NojuVOjw0zrz3n4hmhvLyV85AeU7khwkR8uh0OlitVlbRtrCwAIVCEbXd9q233sJ3v/tdvPTSS+yAhQ2AJKljAcMwmJubw6FDh/DMM8/A6XTi0ksvxd69eyGVSnHkyBGcffbZCZEUxqPjDlcLJ3ObgeWGGOKjngjBxeTkJO54YQJ9Bl/Y92zJS8WzXxRmq0tMLYaGhuB2u9noJLjFkuDdd9/FnXfeiRdffDEhlY41RJLUYkCn0+GZZ57B448/jqGhIVx11VW46aabRNdyLy0tob+/X5Q+6+BauEqlQkpKCsxms6jmg1xMTExgaWkJ71ky8cCbY2Hf95VzyvGVcyJnqoNBjBmlUimqqqrYMiBRtJE6f0pKCo4ePYrbbrsNzz//vOhzx04BJEktFmZmZnDppZfiV7/6FU6ePIlDhw5hfn4eF110Ea644grU1tbGlSgjK6gYftyhMDY2hunpaSgUihUkEOv8VqsVdXV16J+34+rfHw373ue+2IbqPP6JOZLcpCgq5IOUWB298cYbePDBB2Gz2fDII48keqTsWiFJarHAMAxMJlPAHndpaQkvvPACnn76aYyNjeGCCy7Avn37BGfCE6njBj6UlRLrYZJt1ul0oCiKbRuNdX9NVG5EcsswDM76+dswOlZqRstzVHj5Kzt5n5tUGBiGwdatWyNGRj09PfjSl76EK6+8Eu+//z7OO+883HbbbTF9pkh49dVX8dWvfhU0TePGG2/Et7/9bdGvEQFJUq8WrFYrXnrpJRw6dAiDg4M477zzsHfvXuzYsSMiwScnJ6HX6xOi4waWa+hmsxn19fUh78PtdrPGDz6fL6AvPBqIbNXpdGL79u0BhPvG0714qWelzehNu0vxtfM287p3Mkfc7/dHJXR/fz8+97nP4c9//nM8zRJRQdM0tmzZgtdffx0lJSVoa2vDk08+mdBrBiFJ6rUAtye8q6sL55xzDvbu3YudO3eyCatE67hJuY4Qjs/5iUECKSeRenE4p5Dh4WF4PB5s27ZtxevPdMzhrucHVlzjwI07sL0o+rA/cn6fzxd1cMLQ0BCuu+467N+/H/X19VHPHQ/effdd3H333fjLX/4CALjnnnsALDuBrhJEJXXCDZvuvvtuFBcXo6mpCU1NTXj55ZcTfcmEgPSE79+/H0eOHMFFF12EP/3pTzjzzDNx++23429/+xvuvPNOGAyGsCtoPOB2oQl5YMjlchQVFaGpqQltbW1IT0/HxMQE3nvvPQwODsJkMoFhGHYF9Xq9IQkNhJ7eUZSpFERor9cbldBjY2O4/vrr8dhjjyWc0MByboVbxiwpKcHMTHSxzamKVVGUfe1rX8Odd965GpdaFSiVSuzZswd79uyBx+PB66+/jttvvx1qtZr1dj777LNF20sTpxWpVBpX44dUKkV+fj7y8/PZWjiZ6QUsP7jq6+vDnj8vPQVlWXJMLH24r76gJvpweBJhhIsAuJicnMS1116L3//+9wnR84e7v2Cshid4onBKWiuuJygUCgwPD+MrX/kKjh07hs985jN45ZVX8JGPfAQ333wzXnnlFbhcsVvtEpVbSkqKqOU2iUQCrVbLDqhTq9VQKpU4fPgwenp6oNPpAgYxEOyqCJRkXlAb3fttdHSUdSyNdP8zMzP45Cc/iQcffBBtbW3CP1SMKCkpwdTUFPu/p6enWSvm9YiE76nvvvtuPProo8jIyMCOHTvw85//fNW1uokGGVDABU3TePvtt3Ho0CH87W9/w/bt27Fv3z6cf/75vCWcNE2zBgHR+oljvW/itUasmRiGgcViYYcgqFSqAC/xfw0b8IUnugAAuWkK/P1ruyISdXR0FA6HY0XSLRjz8/O46qqr8Itf/AJnn3226J81Enw+H7Zs2YI33ngDxcXFaGtrwxNPPIHt27ev1i2ceomy888/H/Pz8yv+/cc//jF27twJrVYLiqLwve99D3Nzc3jkkUdivN31CTLY/ODBg2xP+L59+3DRRReFbbogRvW5ubkJka3y0VpzvcQXFxchl8uRpdHi0kf6QckU+OSOYnzvki1hrzE2NgabzYa6urqIhNbpdLjyyitx7733rlkd+uWXX8btt98OmqZxww034K677lrNy596pOaL8fFxXHbZZejpid1Per2D9IQfOHAAr776KkpLS3H55ZfjkksuYfuBiVd5cXExb/shISDtjaQbim9IT4YgXHLvy0BBLX56USEubCwPKazhClciJfX0ej2uvPJK/PCHP8TFF18c82da51hf2e+5uTn2v5955hnU1dXFfK777rsPFEVBr9eLcWtrAolEgpaWFtxzzz04fvw4fvSjH2FiYgJ79uzBlVdeiQcffBAf//jHo86GihXE8yuWTjS1Wo1ubzdmi3+IqZRP4tDUPTjRfQKHDx/G2NgY7HY7gOWHNx9Cm0wmXH311fj+979/OhNadCR8pf7MZz6Djo4OUBSF8vJy/Pa3v43pxzo1NYUbb7wRAwMDOHbsGLTa6FnX9QSGYfDWW2/h2muvRWlpKdLS0nD55ZfjsssuE9wTHg5CLIJC4fDsYVzyf5fA6Vs2TFDJVLhm2zX4xbm/YGvhFosFUqmUnSce7r7NZjOuvPJK3HnnnfFaAW0ErC+ThMcff1yU83zta1/Dvffei71794pyvlMNFEXhvffewxNPPIGzzjqL7Qm/9tprkZKSgj179vDqCQ8HIa6c4fDyyMtw+T7M5Dt9Tjw/9Dx+fdGvUVRUBJ/PB4ZhUFBQgImJCdjtdmg0GuTn5yMzM5O9b6vVik984hP46le/miR0ArAunE+ef/55FBcXo7Gxca1vJaH45je/yf53VVUVvvWtb+Gb3/wm2xP+2c9+FgBw2WWXhewJDwcy2SInJyeuDqcMRQbkEjk8fg/7b2r5ciZ/cnISRqOR1ciTIQhGoxEzMzPo7++HTqeD2+3Go48+ii984Qu45pprYr6XJMLjlJGJRsqg/+QnP8Frr72GzMxMlJeX4+jRoxsu/OYDbk/4008/DZfLhcsuuwx79+5FRUVFSIKTudbhJlsIgcFpQPsf2mF0GeGhPVDJVPj9Jb9Hq7oVer0+YtOL3+/HW2+9hR/84AeYm5vDueeei9tvvz2hApNvfOMbeOGFF6BQKLB582b84Q9/EH1uuEhYv9nvWNDd3Y3zzjuPre0SYcDhw4c3kvOFYDAMw/aEP/3001haWsIll1yCvXv3siIVUhbLy8sLOZUjFhicBjza+SgsHgs+tvljKGFK2G6xSEkxl8uFa6+9Fvv27cNNN92EY8eOISMjAzU1NaLcVyi89tprOPfccyGTyfCtb30LAPCzn/0sYdeLA6cXqYMRz0r9ve99D8899xwkEgny8vLw6KOPrmvlEBcGgwHPPvssnn76aSwsLODcc8/FW2+9hV/84hcJ009PT09Dp9Ox7Z/h4PF48OlPfxoXXnghbr311jWRYD7zzDM4ePAg9u/fv+rX5oEkqWMltcViYe2EfvnLX6Kvrw8PPfSQ2Le45piamsJFF12E3NxcWCwWtie8oaFBtEaTmZkZLCwsRCW01+vFZz/7WezevRtf//rX10xTvWfPHlxzzTX49Kc/vSbXj4L1lf0WG+Pj4zEfy/UHs9vt61q0HwkdHR34wQ9+gGuuuYbtCb///vsxNDTE9oS3trbGTPDZ2VlehPb5fLjxxhvR1taWMEJHysWQSsmPf/xjyGQyfOpTnxL9+qci1t1KHS/uuusu/PGPf0RmZibefPNNXsPoNwocDgdeeeUVHDp0CD09PWxP+BlnnBGRnFzMzc1hdnYWTU1NEY+haRpf/OIXUVVVhbvvvnvNHqCPPfYYHnroIbzxxhurOmNLIE7v8Dsa+Dy5geVGeJfLhf/8z/9czds7ZeByufD666/j4MGDOH78OHbt2oUrrrgCu3btCuvWIoTQt912GwoKCvCTn/xkzQj96quv4o477sA//vGPU/3hnSS1GJiYmMCll156WuvQCTweD/72t7/h0KFDePfdd3HGGWdg3759OOuss9ie8Pn5eUxPT6OpqSmiRZPf78cdd9yBtLQ03HfffaKbRQhBVVUV3G43O4Zp586dp2oOJUnqWHHy5ElUV1cDAH71q1/hH//4Bw4ePCjoHOuo9hkTfD4f/vnPf+LAgQN466230NLSgvz8fFitVtx7771RCU0M+375y1+uKaHXGZKkjhVXXnklBgcHIZFIUFZWhoceekiwKfw6qn3GDZqmcc899+C3v/0ttFotampqsHfv3pA94X6/H9///vdhs9nw0EMPJQktDKd39jseHDp0KO5zXHjhhex/79y5U/BKv55AjBRJE8jhw4dx4MAB3HPPPaiursa+fftw4YUXIjU1FT/60Y9gNBrx8MMPJwm9xjitVmqxcYrXPhMGv9+PEydO4MCBA/jLX/4Cj8eDLVu24ODBg7yz6EkEIBl+Jxp8a59Hjx7F008/vWHr3Xzg9/vx4osv4txzzxVtdO5piCSp1xrrpPaZxPrB+nI+2Wh49dVX8bOf/QzPP/98XIQ+cOAAa8h/9Gj4GVVJBGIjuN8kGklSC8Qtt9wCq9WKCy64AE1NTfjiF78Y03nq6urw9NNPr7pz5nrG1NQUXn/99Y049VJUnFbZbzEwPDwsynlqa2tFOc/phI3ufiMWkit1EusCp4v7jRhIrtQJBF8dehLL4ON+k0R0JEmdQPz1r39d61tYVwj39+ru7sbY2Bi7Sk9PT6OlpeW0d78Jh2T4vU7x6quvYuvWraiqqsJPf/rTtb6dhKK+vh46nQ7j4+MYHx9HSUkJjh8/niR0GCRJvUZ45plnUFJSgnfffReXXnopLrroIt7H0jSNr3zlK3jllVfQ19eHJ598En19fQm82yTWE5Lik3WIU2BIehLiIik+Od2x0YakJyEukqReh9hoQ9KTEBdJUq9DbLQh6UmIi2h76iROQVAUJQMwBOA8ADMAjgC4lmGY3hjO9QiAywDoGIaJfSRpEqcMkiv1OgTDMD4AtwD4C4B+AE/FQugP8CiA5BzZDYTkSp0EKIoqB/BicqXeGEiu1EkkscGQJHUSSWwwJEmdRBIbDElSJ5HEBkOS1Kc5KIp6EsC7ALZSFDVNUdTn1/qekogPyex3EklsMCRX6iSS2GBIkjqJJDYYkqROIokNhiSpk0higyFJ6iSS2GBIkjqJJDYYkqROIokNhiSpk0hig+H/A9j7sjgteszUAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig = plt.figure()\n", + "ax = plt.axes(projection='3d')\n", + "coord[:, 0] = np.sum(coord[:, 1:], axis=1) / 3.\n", + "# we plot the origin in red, A in green, B in blue and C in cyan.\n", + "ax.plot_trisurf(coord[0, 1:], coord[1, 1:], coord[2, 1:], shade=False)\n", + "ax.scatter3D(coord[0], coord[1], coord[2], c=[\"r\", \"g\", \"b\", \"c\"])\n", + "ax.plot([coord[0,0], coord[0,1]], [coord[1,0], coord[1,1]], [coord[2,0], coord[2,1]], c=\"black\")\n", + "ax.plot([coord[0,0], coord[0,2]], [coord[1,0], coord[1,2]], [coord[2,0], coord[2,2]], c=\"black\")\n", + "ax.plot([coord[0,0], coord[0,3]], [coord[1,0], coord[1,3]], [coord[2,0], coord[2,3]], c=\"black\")" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "bizarre-production", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 1, 3, 8],\n", + " [-1, -6, -6],\n", + " [-1, -5, 10]])" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "vectors = coord[:, 0, None] - coord[:, 1:] # [OA; OB; OC]\n", + "vectors" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "infinite-desire", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 3, 8, 1],\n", + " [-6, -6, -1],\n", + " [-5, 10, -1]])" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "others = np.roll(vectors, -1, axis=1) # [OB; OC; OA]\n", + "others" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "embedded-agent", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ -1, -90, 16],\n", + " [ 2, -70, 18],\n", + " [ -3, 30, -2]])" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "nA = np.cross(vectors, others, axis=0) # normal of [OAB; OBC; OCA]\n", + "nA" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "present-vaccine", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 16, -1, -90],\n", + " [ 18, 2, -70],\n", + " [ -2, -3, 30]])" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "nB = np.roll(nA, 1, axis=1) # normal of OCA, OAB, OBC\n", + "nB" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "driven-farmer", + "metadata": {}, + "outputs": [], + "source": [ + "norms = np.linalg.norm(nA, axis=0)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "generic-profile", + "metadata": {}, + "outputs": [], + "source": [ + "phiAB = np.arccos(np.abs(np.dot(nA[:, 2], nA[:,1]))/ (norms[2] * norms[1])) # OAC OBC" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "raised-yahoo", + "metadata": {}, + "outputs": [], + "source": [ + "phiBC = np.arccos(np.abs(np.dot(nA[:, 0], nA[:,2]))/ (norms[0] * norms[2])) # OAB OAC" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "latin-lounge", + "metadata": {}, + "outputs": [], + "source": [ + "phiAC = np.arccos(np.abs(np.dot(nA[:, 0], nA[:,1]))/ (norms[0] * norms[1])) # OAB OBC" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "naked-occurrence", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(0.2508036166329773, 1.2791357466336846, 1.247848641316634)" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "phiAB, phiBC, phiAC" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "inner-bibliography", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "-0.36380464900649745" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "res_f1 = phiAB + phiBC + phiAC - np.pi; res_f1" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "norwegian-scheme", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "78.5 µs ± 3.61 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n" + ] + } + ], + "source": [ + "%%timeit\n", + "vectors = coord[:, 0, None] - coord[:, 1:]\n", + "nA = np.cross(vectors, others, axis=0) # normal of [OAB; OBC; OCA]\n", + "norms = np.linalg.norm(nA, axis=0)\n", + "phiAB = np.arccos(np.abs(np.dot(nA[:, 2], nA[:,1]))/ (norms[2] * norms[1])) # OAC OBC\n", + "phiBC = np.arccos(np.abs(np.dot(nA[:, 0], nA[:,2]))/ (norms[0] * norms[2])) # OAB OAC\n", + "phiAC = np.arccos(np.abs(np.dot(nA[:, 0], nA[:,1]))/ (norms[0] * norms[1])) # OAB OBC\n", + "res_f1 = phiAB + phiBC + phiAC - np.pi; res_f1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "continuing-ferry", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Notes_Upgrades/SA_tetra/LM_version.pdf b/Notes_Upgrades/SA_tetra/LM_version.pdf new file mode 100644 index 000000000..563a57360 Binary files /dev/null and b/Notes_Upgrades/SA_tetra/LM_version.pdf differ diff --git a/Notes_Upgrades/SA_tetra/LM_version.tex b/Notes_Upgrades/SA_tetra/LM_version.tex new file mode 100644 index 000000000..52d3d080d --- /dev/null +++ b/Notes_Upgrades/SA_tetra/LM_version.tex @@ -0,0 +1,216 @@ +\documentclass[10pt,a4paper]{article} +\usepackage[margin=0.5in]{geometry} +\usepackage[utf8]{inputenc} +\usepackage[english]{babel} +\usepackage{amsmath} +\usepackage{amsfonts} +\usepackage{amssymb} + \usepackage{cancel} +\usepackage{graphicx} + +\title{Solid angle subtended by a tetrahedron computation} +\begin{document} + +\maketitle + + +\includegraphics[scale=0.4]{tetra.png} + + +\section{Notations} + +Let + +\begin{itemize} + \item $A$ be the vector $\vec{OA}$ + \item $B$ be the vector $\vec{OB}$ + \item $C$ be the vector $\vec{OC}$ + \item $a$ be the vector $\vec{GA}$ + \item $b$ be the vector $\vec{GB}$ + \item $c$ be the vector $c$ + \item $\mathbf{A}$ be the magnitude of the vector $\vec{OA}$ + \item $\mathbf{B}$ be the magnitude of the vector $\vec{OB}$ + \item $\mathbf{C}$ be the magnitude of the vector $\vec{OC}$ +\end{itemize} + +\section{Computation} + +The solid angle $\Omega$ subtended by the triangular surface $ABC$ is: + +$$ +{\displaystyle \tan \left({\frac {1}{2}}\Omega \right) + = {\frac {\left|{A}\ {B}\ {C}\right|} + {\mathbf{A}\mathbf{B}\mathbf{C} + \left({A}\cdot {B}\right)\mathbf{C} + + \left({A}\cdot {C}\right)\mathbf{B} + + \left({B}\cdot {C}\right)\mathbf{A}}}} +$$ + +where $ \left|{A}\ {B}\ {C}\right|={A}\cdot ({B}\times {C}) $ + + +\subsection{Numerator} + + +Given that $A = \vec{OA} = G + a$ (and resp. +with $B$ and $C$), we get + + +\begin{align*} +\left|{A}\ {B}\ {C}\right| + & = {A} \cdot ({B}\times {C}) \\ + & = {A} \cdot ({(G + b)} \times {(G + c)}) \\ + & = {(G + a)} \cdot (\cancel{G \times G} + + G \times c + + b \times G + + b \times c) \\ + & = \cancel{G \cdot (G \times c)} + + \cancel{G \cdot (b \times G)}\\ + & + G \cdot (b \times c) + + a \cdot (G \times c) + + a \cdot (b \times G) + + \cancel{a \cdot (b \times c)} \\ + & = G \cdot (b \times c) + + G \cdot (c \times a) + + G \cdot (a \times b) \\ + & = G \cdot \left(b \times c + +c \times a + + a \times b \right) +\end{align*} + +since $G$ is the centroid of $ABC$, + +\begin{equation} +a + b + c= 0 +\label{eq1} +\end{equation}. We obtain + +\begin{align} +\left|{A}\ {B}\ {C}\right| +& = G \cdot \left(b \times c + +c \times a + + a \times b \right) \nonumber \\ +& = G \cdot \left(b \times c + +c \times (-b - c) + + (-b - c) \times b \right) \nonumber \\ +& = G \cdot \left(b \times c + +c \times (-b) + \cancel{c \times (- c)} + + \cancel{(-b) \times b} + (- c) \times b \right) \nonumber \\ +& = G \cdot \left(3 b \times c \right) \nonumber \\ +& = 3 G \cdot \left( b \times c \right) +\end{align} + +\subsection{Denominator} + + +First let's see the term in $\mathbf{C}$ +\begin{align*} + \left({A}\cdot {B}\right)\mathbf{C} + =& \left( G + a \right) \cdot \left( G + b \right) \mathbf{C} \\ + =& (\mathbf{G}^2 + G \cdot b + a \cdot G + a \cdot b ) \mathbf{C} +\end{align*} + +Using \eqref{eq1} + +\begin{align} + \left({A}\cdot {B}\right)\mathbf{C} + =& (\mathbf{G}^2 + G \cdot b + (-b - c) \cdot G + (-b - c) \cdot b ) \mathbf{C} \nonumber \\ + =& (\mathbf{G}^2 - G \cdot c - c \cdot b -\mathbf{b}^2)\mathbf{C} +\end{align} + + + +Now, the term in $\mathbf{B}$ + + +\begin{align*} + \left({A}\cdot {C}\right)\mathbf{B} + =& \left( G + a \right) \cdot \left( G + c \right) \mathbf{B} \\ + =& (\mathbf{G}^2 + G \cdot c + a \cdot G + a \cdot c ) \mathbf{B} +\end{align*} + +Using \eqref{eq1} + +\begin{align} + \left({A}\cdot {C}\right)\mathbf{B} + =& (\mathbf{G}^2 + G \cdot c + (-b - c) \cdot G + (-b - c) \cdot c ) \mathbf{B} \nonumber \\ + =& (\mathbf{G}^2 - G \cdot b - c \cdot b - \mathbf{c}^2)\mathbf{B} +\end{align} + + +For the third term there is no simplification + +\begin{align} + \left({B}\cdot {C}\right)\mathbf{A} + =& \left( G + b \right) \cdot \left( G + c \right) \mathbf{A}\\ + =& (\mathbf{G}^2 + G \cdot b + G \cdot c + c \cdot b)\mathbf{A} +\end{align} + + +For the computation of the norms, we can use: + +\begin{align*} +\mathbf{A}^2 &= ||OG||^2 + \mathbf{a}^2 + 2 G \cdot a \\ + &= (G+a) \cdot (G+a) \\ + &= (G-b-c) \cdot (G-b-c) \\ +\end{align*} + +and respectively with B and C, we obtain + + +\begin{align} +(\mathbf{A}\mathbf{B}\mathbf{C}) &= ||G-b-c||^2\,||G+b||^2\,||G+c||^2 +\end{align} + +and + + +\begin{align} + \left({A}\cdot {B}\right)\mathbf{C} + +& \left({A}\cdot {C}\right)\mathbf{B} + + \left({B}\cdot {C}\right)\mathbf{A} \\ \nonumber + =& \mathbf{G}^2 (\mathbf{A} + \mathbf{B} + \mathbf{C}) + + c \cdot b (\mathbf{A} - \mathbf{B} - \mathbf{C}) + + G \cdot b (\mathbf{A} - \mathbf{B}) + + G \cdot c (\mathbf{A} - \mathbf{C}) + - ||b||^2\, \mathbf{C} - ||c||^2\, \mathbf{B} +\end{align} + + +\section{Final formula} + + +$$ +{\displaystyle \tan \left({\frac {1}{2}}\Omega \right) + = {\frac {3 G \cdot \left( b \times c \right)} + {\mathbf{A}\mathbf{B}\mathbf{C} + \left({A}\cdot {B}\right)\mathbf{C} + + \left({A}\cdot {C}\right)\mathbf{B} + + \left({B}\cdot {C}\right)\mathbf{A}}}} +$$ + +we get + +$$ +{\displaystyle \tan \left({\frac {1}{2}}\Omega \right) + = {\frac {3 G \cdot \left( b \times c \right)} + {\mathbf{A}\mathbf{B}\mathbf{C} + + (\mathbf{G}^2 - G \cdot c - c \cdot b)\mathbf{C} + + (\mathbf{G}^2 - G \cdot c - c \cdot b)\mathbf{B} + + \left( G + b \right) \cdot \left( G + c \right) \mathbf{A} + }}} +$$ + +or + +$$ +{\displaystyle \tan \left({\frac {1}{2}}\Omega \right) + = {\frac {3 G \cdot \left( b \times c \right)} + {\mathbf{A}\mathbf{B}\mathbf{C} + + \mathbf{G}^2 (\mathbf{A} + \mathbf{B} + \mathbf{C}) + + c \cdot b (\mathbf{A} - \mathbf{B} - \mathbf{C}) + + G \cdot b (\mathbf{A} - \mathbf{B}) + + G \cdot c (\mathbf{A} - \mathbf{C}) + }}} +$$ + + +\end{document} diff --git a/Notes_Upgrades/SA_tetra/SA_tetra.pdf b/Notes_Upgrades/SA_tetra/SA_tetra.pdf new file mode 100644 index 000000000..5c10085e3 Binary files /dev/null and b/Notes_Upgrades/SA_tetra/SA_tetra.pdf differ diff --git a/Notes_Upgrades/SA_tetra/SA_tetra.tex b/Notes_Upgrades/SA_tetra/SA_tetra.tex new file mode 100644 index 000000000..ff1d6ba63 --- /dev/null +++ b/Notes_Upgrades/SA_tetra/SA_tetra.tex @@ -0,0 +1,283 @@ +\documentclass[10pt,a4paper]{article} +\usepackage[margin=0.5in]{geometry} +\usepackage[utf8]{inputenc} +\usepackage[english]{babel} +\usepackage{amsmath} +\usepackage{amsfonts} +\usepackage{amssymb} +\usepackage{cancel} +\usepackage{xcolor} +\usepackage{graphicx} + +%\renewcommand\CancelColor{} +\newcommand\cc[2][black]{\renewcommand\CancelColor{\color{#1}}\cancel{#2}} +\newcommand{\ud}[1]{\underline{#1}} +\DeclareMathOperator{\tm}{\times} +\DeclareMathOperator{\A}{\ud{A}} +\DeclareMathOperator{\B}{\ud{B}} +\DeclareMathOperator{\C}{\ud{C}} +\DeclareMathOperator{\G}{\ud{G}} +\DeclareMathOperator{\av}{\ud{a}} +\DeclareMathOperator{\bv}{\ud{b}} +\DeclareMathOperator{\cv}{\ud{c}} + + +\title{Solid angle subtended by a tetrahedron computation} +\begin{document} + +\maketitle + + +\includegraphics[scale=0.4]{tetra.png} + + +\section{Notations} + +Let + +\begin{itemize} + \item $A$ be the vector $\vec{OA}$ + \item $B$ be the vector $\vec{OB}$ + \item $C$ be the vector $\vec{OC}$ + \item $a$ be the vector $\vec{GA}$ + \item $b$ be the vector $\vec{GB}$ + \item $c$ be the vector $\vec{GC}$ + \item $\mathbf{A}$ be the magnitude of the vector $\vec{OA}$ + \item $\mathbf{B}$ be the magnitude of the vector $\vec{OB}$ + \item $\mathbf{C}$ be the magnitude of the vector $\vec{OC}$ +\end{itemize} + +\section{Computation} + +The solid angle $\Omega$ subtended by the triangular surface $ABC$ is: + +$$ +{\displaystyle \tan \left({\frac {1}{2}}\Omega \right) + = {\frac {\left|{A}\ {B}\ {C}\right|} + {\mathbf{A}\mathbf{B}\mathbf{C} + \left({A}\cdot {B}\right)\mathbf{C} + + \left({A}\cdot {C}\right)\mathbf{B} + + \left({B}\cdot {C}\right)\mathbf{A}}}} +$$ + +where $ \left|{A}\ {B}\ {C}\right|={A}\cdot ({B}\times {C}) $ + + +\subsection{Numerator} + + +Given that $A = \vec{OA} = G + a$ (and resp. +with $B$ and $C$), we get + + +\begin{align*} +\left|{A}\ {B}\ {C}\right| + & = {A} \cdot ({B}\times {C}) \\ + & = {A} \cdot ({(G + b)} \times {(G + c)}) \\ + & = {(G + a)} \cdot (\cancel{G \times G} + + G \times c + + b \times G + + b \times c) \\ + & = \cancel{G \cdot (G \times c)} + + \cancel{G \cdot (b \times G)}\\ + & + G \cdot (b \times c) + + a \cdot (G \times c) + + a \cdot (b \times G) + + \cancel{a \cdot (b \times c)} \\ + & = G \cdot (b \times c) + + G \cdot (c \times a) + + G \cdot (a \times b) \\ + & = G \cdot \left(b \times c + +c \times a + + a \times b \right) +\end{align*} + +since $G$ is the centroid of $ABC$, + +\begin{equation} +a + b + c= 0 +\label{eq1} +\end{equation}. We obtain + + + +\begin{align} +\left|{\A}\ {\B}\ {\C}\right| +& = \G \cdot \left(\bv \tm \cv + \cv \tm \av + \av \tm \bv \right) \nonumber \\ +& = \G \cdot \left(\bv \tm \cv + \cv \tm (-\bv-\cv) + (-\bv-\cv) \tm \bv \right) \nonumber \\ +& = \G \cdot \left(\bv \times \cv + \cv \times (-\bv) + \cancel{\cv \times (-\cv)} + + \cancel{(-\bv) \times \bv} + (- \cv) \times \bv \right) \nonumber \\ +& = \G \cdot \left(3 \bv \times \cv \right) \nonumber \\ +& = 3 \G \cdot \left( \bv \times \cv \right) +\end{align} + + +\subsection{norms} + +We can write: +$$ +\begin{array}{ll} + A^2 + & = \A \cdot \A\\ + & = (\G+\av) \cdot (\G + \av)\\ + & = G^2 + 2\G \cdot \av + a^2\\ +\end{array} +$$ + +Hence: +$$ +\begin{array}{ll} + (ABC)^2 + & = (G^2 + 2\G \cdot \av + a^2) + (G^2 + 2\G \cdot \bv + b^2) + (G^2 + 2\G \cdot \cv + c^2)\\ + & = (G^2 + 2\G \cdot \av + a^2) + (G^4 + 2G^2\G \cdot \cv + G^2c^2 + + 2G^2\G \cdot \bv + 4(\G \cdot \bv)(\G \cdot \cv) + 2c^2\G \cdot \bv + + b^2G^2 + 2b^2\G \cdot \cv + b^2c^2)\\ +\end{array} +$$ + +Given that: $\av = -(\bv+\cv) \Rightarrow a^2=b^2 + 2\bv \cdot \cv + c^2$ + +We can write: +$$ +\begin{array}{lll} + (ABC)^2 + & = & (G^2 + 2\G \cdot \av + a^2) + (G^4 + 2G^2\G \cdot \cv + G^2c^2 + + 2G^2\G \cdot \bv + 4(\G \cdot \bv)(\G \cdot \cv) + 2c^2\G \cdot \bv + + b^2G^2 + 2b^2\G \cdot \cv + b^2c^2)\\ + & = & (G^2 - 2\G \cdot \bv - 2\G \cdot \cv + b^2 + 2\bv \cdot \cv + c^2)\\ + && (G^4 + 2G^2\G \cdot \cv + G^2c^2 + + 2G^2\G \cdot \bv + 4(\G \cdot \bv)(\G \cdot \cv) + 2c^2\G \cdot \bv + + b^2G^2 + 2b^2\G \cdot \cv + b^2c^2)\\ + & = & G^2 + \left(G^4 + \cc[black]{2G^2\G \cdot \cv} + G^2c^2 + + \cc[red]{2G^2\G \cdot \bv} + \cc[blue]{4(\G \cdot \bv)(\G \cdot \cv)} + + \cc[green]{2c^2\G \cdot \bv} + + b^2G^2 + \cc[yellow]{2b^2\G \cdot \cv} + b^2c^2\right)\\ + && - 2\G \cdot \bv + \left(\cc[red]{G^4} + \cc[blue]{2G^2\G \cdot \cv} + \cc[green]{G^2c^2} + + 2G^2\G \cdot \bv + 4(\G \cdot \bv)(\G \cdot \cv) + 2c^2\G \cdot \bv + + \cc[red]{b^2G^2} + \cc[cyan]{2b^2\G \cdot \cv} + \cc[blue]{b^2c^2}\right)\\ + && - 2\G \cdot \cv + \left(\cc[black]{G^4} + 2G^2\G \cdot \cv + \cc[magenta]{G^2c^2} + + 2G^2\G \cdot \bv + 4(\G \cdot \bv)(\G \cdot \cv) + + \cc[black]{2c^2\G \cdot \bv} + + \cc[yellow]{b^2G^2} + 2b^2\G \cdot \cv + \cc[green]{b^2c^2}\right)\\ + &&+ 2\bv \cdot \cv + \left(G^4 + 2G^2\G \cdot \cv + G^2c^2 + + 2G^2\G \cdot \bv + 4(\G \cdot \bv)(\G \cdot \cv) + 2c^2\G \cdot \bv + + b^2G^2 + 2b^2\G \cdot \cv + b^2c^2\right)\\ + && + b^2 + \left(G^4 + 2G^2\G \cdot \cv + G^2c^2 + + \cc[red]{2G^2\G \cdot \bv} + \cc[cyan]{4(\G \cdot \bv)(\G \cdot \cv)} + + \cc[blue]{2c^2\G \cdot \bv} + + b^2G^2 + 2b^2\G \cdot \cv + b^2c^2\right)\\ + && + c^2 + \left(G^4 + \cc[magenta]{2G^2\G} \cdot \cv + G^2c^2 + + 2G^2\G \cdot \bv + \cc[black]{4(\G \cdot \bv)(\G \cdot \cv)} + 2c^2\G \cdot \bv + + b^2G^2 + \cc[green]{2b^2\G \cdot \cv} + b^2c^2\right)\\ +\end{array} +$$ + + +\subsection{Denominator} + + +First let's see the term in $\mathbf{C}$ +\begin{align*} + \left({A}\cdot {B}\right)\mathbf{C} + =& \left( G + a \right) \cdot \left( G + b \right) \mathbf{C} \\ + =& (\mathbf{G}^2 + G \cdot b + a \cdot G + a \cdot b ) \mathbf{C} +\end{align*} + +Using \eqref{eq1} + +\begin{align} + \left({A}\cdot {B}\right)\mathbf{C} + =& (\mathbf{G}^2 + G \cdot b + (-b - c) \cdot G + (-b - c) \cdot b ) c \nonumber \\ + =& (\mathbf{G}^2 - G \cdot c - c \cdot b)\mathbf{C} +\end{align} + + + +Now, the term in $\mathbf{B}$ + + +\begin{align*} + \left({A}\cdot {C}\right)\mathbf{B} + =& \left( G + a \right) \cdot \left( G + c \right) \mathbf{B} \\ + =& (\mathbf{G}^2 + G \cdot c + a \cdot G + a \cdot c ) \mathbf{B} +\end{align*} + +Using \eqref{eq1} + +\begin{align} + \left({A}\cdot {C}\right)\mathbf{B} + =& (\mathbf{G}^2 + G \cdot c + (-b - c) \cdot G + (-b - c) \cdot c ) b \nonumber \\ + =& (\mathbf{G}^2 - G \cdot c - c \cdot b)\mathbf{B} +\end{align} + + +For the third term there is no simplification + +\begin{align} + \left({B}\cdot {C}\right)\mathbf{A} + =& \left( G + b \right) \cdot \left( G + c \right) \mathbf{A}\\ + =& (\mathbf{G}^2 + G \cdot b + G \cdot c + c \cdot b)\mathbf{A} +\end{align} + + +For the computation of the norms, we can use: + +\begin{align*} +\mathbf{A}^2 &= ||OG||^2 + \mathbf{a}^2 + 2 G \cdot a \\ + &= (G+a) \cdot (G+a) \\ + &= (G-b-c) \cdot (G-b-c) \\ +\end{align*} + +and respectively with B and C, we obtain + + +\begin{align} +(\mathbf{A}\mathbf{B}\mathbf{C}) &= ||G-b-c||\,||G+b||\,||G+c|| +\end{align} + +and + + +\begin{align} + \left({A}\cdot {B}\right)\mathbf{C} + + \left({A}\cdot {C}\right)\mathbf{B} + + \left({B}\cdot {C}\right)\mathbf{A} + =& \mathbf{G}^2 (\mathbf{A} + \mathbf{B} + \mathbf{C}) + + c \cdot b (\mathbf{A} - \mathbf{B} - \mathbf{C}) + + G \cdot b (\mathbf{A} - \mathbf{B}) + + G \cdot c (\mathbf{A} - \mathbf{C}) +\end{align} + + +\section{Final formula} + + +$$ +{\displaystyle \tan \left({\frac {1}{2}}\Omega \right) + = {\frac {3 G \cdot \left( b \times c \right)} + {\mathbf{A}\mathbf{B}\mathbf{C} + \left({A}\cdot {B}\right)\mathbf{C} + + \left({A}\cdot {C}\right)\mathbf{B} + + \left({B}\cdot {C}\right)\mathbf{A}}}} +$$ + + +$$ +{\displaystyle \tan \left({\frac {1}{2}}\Omega \right) + = {\frac {3 G \cdot \left( b \times c \right)} + {\mathbf{A}\mathbf{B}\mathbf{C} + + (\mathbf{G}^2 - G \cdot c - c \cdot b)\mathbf{C} + + (\mathbf{G}^2 - G \cdot c - c \cdot b)\mathbf{B} + + \left( G + b \right) \cdot \left( G + c \right) \mathbf{A} + }}} +$$ + +\end{document} diff --git a/Notes_Upgrades/SA_tetra/drawing.svg b/Notes_Upgrades/SA_tetra/drawing.svg new file mode 100644 index 000000000..87d03cb7c --- /dev/null +++ b/Notes_Upgrades/SA_tetra/drawing.svg @@ -0,0 +1,153 @@ + + + + + + + + image/svg+xml + + + + + + + + + A + B + C + + + + G + + O + + + diff --git a/Notes_Upgrades/SA_tetra/tetra.png b/Notes_Upgrades/SA_tetra/tetra.png new file mode 100644 index 000000000..c9f56ea4c Binary files /dev/null and b/Notes_Upgrades/SA_tetra/tetra.png differ diff --git a/Notes_Upgrades/SA_tetra/tetra.svg b/Notes_Upgrades/SA_tetra/tetra.svg new file mode 100644 index 000000000..fced958bd --- /dev/null +++ b/Notes_Upgrades/SA_tetra/tetra.svg @@ -0,0 +1,150 @@ + + + + + + + + image/svg+xml + + + + + + + + + A + B + C + + + + G + + O + + + diff --git a/setup.py b/setup.py index 22edbb15b..97370b51d 100644 --- a/setup.py +++ b/setup.py @@ -158,13 +158,13 @@ def get_version_tofu(path=_HERE): _README = [ ff for ff in os.listdir(_HERE) - if len(ff) <= 10 and ff[:7] == "README." + if len(ff) <= 10 and ff.startswith("README.") ] assert len(_README) == 1 _README = _README[0] with open(os.path.join(_HERE, _README), encoding="utf-8") as f: long_description = f.read() -if _README[-3:] == ".md": +if _README.endswith(".md"): long_description_content_type = "text/markdown" else: long_description_content_type = "text/x-rst" @@ -175,7 +175,7 @@ def get_version_tofu(path=_HERE): # Compiling files openmp_installed, openmp_flag = is_openmp_installed() -extra_compile_args = ["-O3", "-Wall", "-fno-wrapv"] + openmp_flag +extra_compile_args = ["-O3", "-Wall", "-fno-wrapv", "-ffast-math"] + openmp_flag extra_link_args = [] + openmp_flag extensions = [ diff --git a/tofu/geom/_GG.pyx b/tofu/geom/_GG.pyx index cd7a8af7a..1304f4eab 100644 --- a/tofu/geom/_GG.pyx +++ b/tofu/geom/_GG.pyx @@ -1621,7 +1621,7 @@ def _Ves_Smesh_Tor_SubFromInd_cython(double dL, double dRPhi, # Number of Phi per R dRPhirRef, dPhir, dRPhir = np.empty((Ln,)), np.empty((Ln,)), -np.ones((Ln,)) dLr, Rref = -np.ones((Ln,)), -np.ones((Ln,)) - NRPhi, NRPhi0 = np.empty((Ln,),dtype=int), np.empty((Ln,),dtype=int) + NRPhi, NRPhi0 = np.empty((Ln,),dtype=int), np.empty((Ln+1,),dtype=int) radius_ratio = int(c_ceil(np.max(RrefRef)/np.min(RrefRef))) for ii in range(0,Ln): NRPhi[ii] = (c_ceil(DPhiMinMax*RrefRef[ii]/dRPhi)) @@ -1637,8 +1637,8 @@ def _Ves_Smesh_Tor_SubFromInd_cython(double dL, double dRPhi, nRPhi0 = NRPhi0[Ln-1]+NRPhi[Ln-1] if Out.lower()=='(x,y,z)': - for ii in range(0,NP): - for jj in range(0,Ln+1): + for ii in range(NP): + for jj in range(Ln + 1): if ind[ii]-NRPhi0[jj]<0.: break iiL = jj-1 @@ -1653,8 +1653,8 @@ def _Ves_Smesh_Tor_SubFromInd_cython(double dL, double dRPhi, Rref[iiL] = RrefRef[iiL] else: - for ii in range(0,NP): - for jj in range(0,Ln+1): + for ii in range(NP): + for jj in range(Ln + 1): if ind[ii]-NRPhi0[jj]<0.: break iiL = jj-1 @@ -1690,7 +1690,7 @@ def _Ves_Smesh_TorStruct_SubFromD_cython(double[::1] PhiMinMax, double dL, for the desired resolution (dR,dZ,dRphi) """ cdef double Dphi, dR0r=0., dZ0r=0. - cdef int NR0=0, NZ0=0, R0n, Z0n, NRPhi0 + cdef int NR0=0, NZ0=0, R0n, Z0n cdef double[::1] phiMinMax = np.array([c_atan2(c_sin(PhiMinMax[0]), c_cos(PhiMinMax[0])), c_atan2(c_sin(PhiMinMax[1]), @@ -1835,7 +1835,7 @@ def _Ves_Smesh_TorStruct_SubFromInd_cython(double[::1] PhiMinMax, double dL, """ Return the desired surfacic submesh indicated by the limits (DR,DZ,DPhi) for the desired resolution (dR,dZ,dRphi) """ cdef double Dphi, dR0r, dZ0r - cdef int NR0, NZ0, R0n, Z0n, NRPhi0 + cdef int NR0, NZ0, R0n, Z0n cdef double[::1] phiMinMax = np.array([c_atan2(c_sin(PhiMinMax[0]), c_cos(PhiMinMax[0])), c_atan2(c_sin(PhiMinMax[1]), @@ -2743,6 +2743,17 @@ def LOS_isVis_PtFromPts_VesStruct(double pt0, double pt1, double pt2, # # ============================================================================== def triangulate_by_earclipping(np.ndarray[double,ndim=2] poly): + """ + Triangulates a polygon by the earclipping method. + Params + ===== + poly : (3, nvert) double array + Contains 3D coordinates of open polygon in counter clockwise order + Returns + ======= + ltri : (3*(nvert-2)) int array + Indices of triangles + """ cdef int nvert = poly.shape[1] cdef np.ndarray[long,ndim=1] ltri = np.empty((nvert-2)*3, dtype=int) cdef double* diff = NULL @@ -2750,10 +2761,10 @@ def triangulate_by_earclipping(np.ndarray[double,ndim=2] poly): # Initialization ........................................................... diff = malloc(3*nvert*sizeof(double)) lref = malloc(nvert*sizeof(bint)) - _vt.compute_diff3d(&poly[0,0], nvert, diff) - _vt.are_points_reflex(nvert, diff, lref) + _vt.compute_diff3d(&poly[0,0], nvert, &diff[0]) + _vt.are_points_reflex(nvert, diff, &lref[0]) # Calling core function..................................................... - _vt.earclipping_poly(&poly[0,0], <ri[0], diff, lref, nvert) + _vt.earclipping_poly(&poly[0,0], <ri[0], &diff[0], &lref[0], nvert) free(diff) free(lref) return ltri @@ -2768,7 +2779,7 @@ def vignetting(double[:, ::1] ray_orig, double[:, ::1] ray_vdir, LOS normalized direction vector vignett_poly : (num_vign, 3, num_vertex) double list of arrays Coordinates of the vertices of the Polygon defining the 3D vignett. - POLY CLOSED + POLY CLOSED in counter clockwise lnvert : (num_vign) long array Number of vertices for each vignett (without counting the rebound) Returns @@ -4885,7 +4896,7 @@ def compute_solid_angle_map(double[:,::1] part_coords, double[::1] part_r, num_threads = _ompt.get_effective_num_threads(num_threads) lstruct_lims_np = flatten_lstruct_lims(lstruct_lims) # .............. - _st.sa_assemble_arrays(block, + _st.sa_sphere_assemble(block, approx, part_coords, part_r, @@ -4932,3 +4943,481 @@ def compute_solid_angle_map(double[:,::1] part_coords, double[::1] part_r, free(ncells_rphi) return pts, sa_map, ind, reso_r_z + +# ============================================================================== +# +# Solid Angle Computation +# subtended by a polygon +# +# ============================================================================== +def compute_solid_angle_poly_map(list poly_coords, + double[:, ::1] poly_lnorms, + long[::1] poly_lnvert, + double rstep, double zstep, double phistep, + double[::1] RMinMax, double[::1] ZMinMax, + bint approx=True, + list DR=None, list DZ=None, DPhi=None, + double[:,::1] limit_vpoly=None, + bint block=False, + double[:, ::1] ves_poly=None, + double[:, ::1] ves_norm=None, + double[::1] ves_lims=None, + long[::1] lstruct_nlim=None, + double[::1] lstruct_polyx=None, + double[::1] lstruct_polyy=None, + list lstruct_lims=None, + double[::1] lstruct_normx=None, + double[::1] lstruct_normy=None, + long[::1] lnvert=None, + int nstruct_tot=0, + int nstruct_lim=0, + double rmin=-1, bint forbid=True, + double eps_uz=_SMALL, double eps_a=_VSMALL, + double eps_vz=_VSMALL, double eps_b=_VSMALL, + double eps_plane=_VSMALL, str ves_type='Tor', + double margin=_VSMALL, int num_threads=48, + bint test=True): + """ + Computes the 2D map of the integrated solid angles subtended by a list of + npoly polygons of coordinates poly_coords[npoly], not necessarily flat, + one sided and side defined by poly_lnorms[npoly], nomal vector of polygon. + If you want to see the details of the computation see: + $TOFU_DIR/Notes_Upgrades/SA_tetra/SA_tetra.tex + Parameters + ---------- + poly_coords: double (npoly, 3, npts) list + coordinates of the points defining the polygons. + not necessarily flat. + poly_coords[np, i] being the i-th coordinate of the np polygon + The coordinates should be given in anti clockwise order and the polygon + should be open + poly_lnorms: double array + normal vector that defines the visible face of the polygon + poly_lnorms[np] defines the norm of the np-th polygon + poly_lnvert: int array + number of vertices for each polygon + rstep: double + refinement along radius `r` + zstep: double + refinement along height `z` + phistep: double + refinement along toroidal direction `phi` + approx: bool + do you want to use approximation (8th order) or exact formula ? + default: True + RMinMax: double memory-view + limits min and max in `r` + ZMinMax: double memory-view + limits min and max in `z` + DR: double memory-view, optional + actual sub-volume limits to get in `r` + DZ: double memory-view, optional + actual sub-volume limits to get in `z` + DPhi: double memory-view, optional + actual sub-volume limits to get in `phi` + limit_vpoly: (3, npts) double memory-view, optional + if we only want to discretize the volume inside a certain flux surface. + Defines the `(R,Z)` coords of the poloidal cut of the limiting flux + surface. + block: bool, optional + check if ppolygon is viewable from viewing points or if there is a + structural element blocking visibility (False) + ves_poly : (2, num_vertex) double array + Coordinates of the vertices of the Polygon defining the 2D poloidal + cut of the Vessel + ves_norm : (2, num_vertex-1) double array + Normal vectors going "inwards" of the edges of the Polygon defined + by ves_poly + nstruct_tot : int + Total number of structures (counting each limited structure as one) + ves_lims : array + Contains the limits min and max of vessel + lstruct_polyx : array + List of x coordinates of the vertices of all structures on poloidal plane + If no structures : None + lstruct_polyy : array + List of y coordinates of the vertices of all structures on poloidal plane + If no structures : None + lstruct_lims : array + List of limits of all structures + If no structures : None + lstruct_nlim : array of ints + List of number of limits for all structures + If no structures : None + lstruct_normx : double memory-view, optional + List of x-coordinates of "inwards" normal vectors of the polygon of all + the structures + If no structures : None + lstruct_normy : double memory-view, optional + List of y-coordinates of "inwards" normal vectors of the polygon of all + the structures + If no structures : None + rmin : double, optional + Minimal radius of vessel to take into consideration + forbid : bool, optional + Should we forbid values behind visible radius ? (see rmin) + eps_ : double, optional + Small value, acceptance of error + margin: double, optional + tolerance error. Defaults to |_VSMALL| + num_threads : int + The num_threads argument indicates how many threads the team should + consist of. If not given, OpenMP will decide how many threads to use. + Typically this is the number of cores available on the machine. + test : bool, optional + Should we run tests? Default True + + Returns + ------- + pts: (2, npts) array of (R, Z) coordinates of viewing points in + vignette where solid angle is integrated + sa_map: (npts, sz_p) array approx solid angle integrated along phi + integral (sa * dphi * r) + ind: (npts) indices to reconstruct (R,Z) map from sa_map + rdrdz: (npts) volume unit: dr*dz + """ + cdef int jj + cdef int sz_r + cdef int sz_z + cdef int npoly + cdef int r_ratio + cdef int npts_pol + cdef int ind_loc_r0 + cdef int npts_disc = 0 + cdef int tot_num_tri + cdef int[1] max_sz_phi + cdef double min_phi, max_phi + cdef double min_phi_pi + cdef double max_phi_pi + cdef double abs0, abs1 + cdef double reso_r_z + cdef double twopi_over_dphi + cdef long[1] ncells_r0, ncells_r, ncells_z + cdef long[::1] ind_mv + cdef long[::1] first_ind_mv + cdef double[2] limits_dl + cdef double[1] reso_r0, reso_r, reso_z + cdef double[::1] reso_rdrdz_mv + cdef double[::1] lstruct_lims_np + cdef double[:, ::1] poly_mv + cdef double[:, ::1] pts_mv + cdef long[:, ::1] indi_mv + cdef long[:, ::1] ind_rz2pol + cdef long[:, ::1] is_in_vignette + cdef long** ltri = NULL + cdef long* lindex = NULL + cdef long* sz_phi = NULL + cdef long* lindex_z = NULL + cdef long* ncells_rphi = NULL + cdef double* disc_r0 = NULL + cdef double* disc_r = NULL + cdef double* disc_z = NULL + cdef double* step_rphi = NULL + cdef double** data = NULL + cdef np.ndarray[long, ndim=1] ind + cdef np.ndarray[double, ndim=1] reso_rdrdz + cdef np.ndarray[double, ndim=1] dot_GBGC + cdef np.ndarray[double, ndim=2] pts + cdef np.ndarray[double, ndim=2] sa_map + cdef np.ndarray[double, ndim=2] vec_GB + cdef np.ndarray[double, ndim=2] vec_GC + cdef np.ndarray[double, ndim=2] centroids + cdef np.ndarray[double, ndim=2] poly_lnorms_tot + cdef np.ndarray[double, ndim=2] cross_GBGC + cdef double[:, ::1] temp + print("begining.....................................") + # + # == Testing inputs ======================================================== + if test: + if block: + msg = "ves_poly and ves_norm are not optional arguments" + assert ves_poly is not None and ves_norm is not None, msg + bool1 = (ves_poly.shape[0]==2 and ves_norm.shape[0]==2 + and ves_norm.shape[1]==ves_poly.shape[1]-1) + msg = "Args ves_poly & ves_norm must be of the same shape (2, NS)!" + assert bool1, msg + bool1 = (lstruct_lims is None + or len(lstruct_normy) == len(lstruct_normx)) + bool2 = (lstruct_normx is None + or len(lstruct_polyx) == len(lstruct_polyy)) + msg = "Args lstruct_polyx, lstruct_polyy, lstruct_lims,"\ + + " lstruct_normx, lstruct_normy, must be None or"\ + + " lists of same len()!" + assert bool1 and bool2, msg + msg = "[eps_uz,eps_vz,eps_a,eps_b] must be floats < 1.e-4!" + assert all([ee < 1.e-4 for ee in [eps_uz, eps_a, + eps_vz, eps_b, + eps_plane]]), msg + msg = "ves_type must be a str in ['Tor','Lin']!" + assert ves_type.lower() in ['tor', 'lin'], msg + # ... + print("after tests.................................") + # .. Formatting polygon coordinates ........................................ + npoly = len(poly_coords) + for ii in range(npoly): + print(poly_coords[ii]) + poly_coords[ii] = format_poly(poly_coords[ii], Clock=False, + close=False, Test=True) + print(poly_coords[ii]) + print("......................... after format") + # .. Dividing polys in triangles ........................................... + ltri = malloc(sizeof(long*) * npoly) + # re writting_polygons coordinates to C type: + data = malloc(npoly * sizeof(double *)) + for ii in range(npoly): + temp = poly_coords[ii] + data[ii] = &temp[0, 0] + _vt.triangulate_polys( + &data[0], + &poly_lnvert[0], + npoly, + ltri, + num_threads + ) + print("...............after triangulate") + # cpumputing total number of triangles + tot_num_tri = np.sum(poly_lnvert) - 2 * npoly + # .. Getting centroids of triangles ....................................... + dot_GBGC = np.zeros(tot_num_tri) + vec_GB = np.zeros((3, tot_num_tri)) + vec_GC = np.zeros((3, tot_num_tri)) + centroids = np.zeros((3, tot_num_tri)) + cross_GBGC = np.zeros((3, tot_num_tri)) + _bgt.find_centroids_GB_GC_ltri( + data, + ltri, + &poly_lnvert[0], + npoly, + num_threads, + centroids, + vec_GB, + vec_GC, + ) + + print("............ got here 1") + print(centroids) + + poly_lnorms_tot = np.repeat(poly_lnorms, + np.asarray(poly_lnvert) - 2, + axis = 0) + assert np.shape(poly_lnorms_tot)[0] == tot_num_tri, ( + f"Total number of triangles = {tot_num_tri} " + + " and was expecting: " + + str(np.shape(poly_lnorms_tot))) + _bgt.compute_dot_cross_vec(vec_GB, + vec_GC, + cross_GBGC, + dot_GBGC, + tot_num_tri, + num_threads, + ) + # .. Check if points are visible ........................................... + # Get the actual R and Z resolutions and mesh elements + # .. First we discretize R without limits .................................. + _st.cythonize_subdomain_dl(None, limits_dl) # no limits + _ = _st.discretize_line1d_core(&RMinMax[0], rstep, limits_dl, + True, 0, # discretize in absolute mode + margin, &disc_r0, reso_r0, &lindex, + ncells_r0) + free(lindex) # getting rid of things we dont need + lindex = NULL + # .. Now the actual R limited ............................................. + _st.cythonize_subdomain_dl(DR, limits_dl) + sz_r = _st.discretize_line1d_core(&RMinMax[0], rstep, limits_dl, + True, 0, # discretize in absolute mode + margin, &disc_r, reso_r, &lindex, + ncells_r) + free(lindex) # getting rid of things we dont need + # .. Now Z ................................................................. + _st.cythonize_subdomain_dl(DZ, limits_dl) + sz_z = _st.discretize_line1d_core(&ZMinMax[0], zstep, limits_dl, + True, 0, # discretize in absolute mode + margin, &disc_z, reso_z, &lindex_z, + ncells_z) + # .. Preparing for phi: get the limits if any and make sure to replace them + # .. in the proper quadrants ............................................... + print("got here 2...........................") + if DPhi is None: + min_phi = -c_pi + max_phi = c_pi + else: + min_phi = DPhi[0] # to avoid conversions + min_phi = c_atan2(c_sin(min_phi), c_cos(min_phi)) + max_phi = DPhi[1] # to avoid conversions + max_phi = c_atan2(c_sin(max_phi), c_cos(max_phi)) + # .. Initialization ........................................................ + sz_phi = malloc(sz_r*sizeof(long)) + ncells_rphi = malloc(sz_r*sizeof(long)) + step_rphi = malloc(sz_r*sizeof(double)) + r_ratio = (c_ceil(disc_r[sz_r - 1] / disc_r[0])) + twopi_over_dphi = _TWOPI / phistep + ind_loc_r0 = 0 + min_phi_pi = min_phi + c_pi + max_phi_pi = max_phi + c_pi + abs0 = c_abs(min_phi_pi) + abs1 = c_abs(max_phi_pi) + # ... doing 0 loop before .................................................. + if min_phi < max_phi: + # Get the actual RPhi resolution and Phi mesh elements (! depends on R!) + ncells_rphi[0] = c_ceil(twopi_over_dphi * disc_r[0]) + loc_nc_rphi = ncells_rphi[0] + step_rphi[0] = _TWOPI / ncells_rphi[0] + inv_drphi = 1. / step_rphi[0] + # Get index and cumulated indices from background + for jj in range(ind_loc_r0, ncells_r0[0]): + if disc_r0[jj]==disc_r[0]: + ind_loc_r0 = jj + break + # Get indices of phi + # Get the extreme indices of the mesh elements that really need to + # be created within those limits + if abs0 - step_rphi[0]*c_floor(abs0 * inv_drphi) < margin*step_rphi[0]: + nphi0 = int(c_round((min_phi + c_pi) * inv_drphi)) + else: + nphi0 = int(c_floor((min_phi +c_pi) * inv_drphi)) + if abs1-step_rphi[0]*c_floor(abs1 * inv_drphi) < margin*step_rphi[0]: + nphi1 = int(c_round((max_phi+c_pi) * inv_drphi)-1) + else: + nphi1 = int(c_floor((max_phi+c_pi) * inv_drphi)) + sz_phi[0] = nphi1 + 1 - nphi0 + max_sz_phi[0] = sz_phi[0] + ind_i = -np.ones((sz_r, sz_phi[0] * r_ratio + 1), dtype=int) + indi_mv = ind_i + for jj in range(sz_phi[0]): + indi_mv[0, jj] = nphi0 + jj + npts_disc += sz_z * sz_phi[0] + else: + # Get the actual RPhi resolution and Phi mesh elements (! depends on R!) + ncells_rphi[0] = c_ceil(twopi_over_dphi * disc_r[0]) + loc_nc_rphi = ncells_rphi[0] + step_rphi[0] = _TWOPI / ncells_rphi[0] + inv_drphi = 1. / step_rphi[0] + # Get index and cumulated indices from background + for jj in range(ind_loc_r0, ncells_r0[0]): + if disc_r0[jj]==disc_r[0]: + ind_loc_r0 = jj + break + # Get indices of phi + # Get the extreme indices of the mesh elements that really need to + # be created within those limits + if abs0 - step_rphi[0]*c_floor(abs0 * inv_drphi) < margin*step_rphi[0]: + nphi0 = int(c_round((min_phi + c_pi) * inv_drphi)) + else: + nphi0 = int(c_floor((min_phi + c_pi) * inv_drphi)) + if abs1-step_rphi[0]*c_floor(abs1 * inv_drphi) < margin*step_rphi[0]: + nphi1 = int(c_round((max_phi+c_pi) * inv_drphi)-1) + else: + nphi1 = int(c_floor((max_phi+c_pi) * inv_drphi)) + sz_phi[0] = nphi1+1+loc_nc_rphi-nphi0 + max_sz_phi[0] = sz_phi[0] + ind_i = -np.ones((sz_r, sz_phi[0] * r_ratio + 1), dtype=int) + indi_mv = ind_i + for jj in range(loc_nc_rphi - nphi0): + indi_mv[0, jj] = nphi0 + jj + for jj in range(loc_nc_rphi - nphi0, sz_phi[0]): + indi_mv[0, jj] = jj - (loc_nc_rphi - nphi0) + npts_disc += sz_z * sz_phi[0] + print("Got heeeeeeere 3...........................") + # ... doing the others ..................................................... + npts_disc += _st.sa_disc_phi(sz_r, sz_z, ncells_rphi, phistep, + disc_r, disc_r0, step_rphi, + ind_loc_r0, + ncells_r0[0], ncells_z[0], &max_sz_phi[0], + min_phi, max_phi, sz_phi, indi_mv, + margin, num_threads) + print("Got heeeeeeere 4...........................") + # ... vignetting ........................................................... + is_in_vignette = np.ones((sz_r, sz_z), dtype=int) # by default yes + if limit_vpoly is not None: + npts_vpoly = limit_vpoly.shape[1] - 1 + # we make sure it is closed + if not(abs(limit_vpoly[0, 0] - limit_vpoly[0, npts_vpoly]) < _VSMALL + and abs(limit_vpoly[1, 0] + - limit_vpoly[1, npts_vpoly]) < _VSMALL): + poly_mv = np.concatenate((limit_vpoly, limit_vpoly[:,0:1]), axis=1) + else: + poly_mv = limit_vpoly + _ = _vt.are_in_vignette(sz_r, sz_z, + poly_mv, npts_vpoly, + disc_r, disc_z, + is_in_vignette) + print("Got heeeeeeere 5...........................") + # .. preparing for actual discretization ................................... + ind_rz2pol = np.empty((sz_r, sz_z), dtype=int) + npts_pol = _st.sa_get_index_arrays(ind_rz2pol, + is_in_vignette, + sz_r, sz_z) + # initializing arrays + reso_rdrdz = np.empty((npts_pol, )) + sa_map = np.zeros((npts_pol, npoly)) + pts = np.empty((2, npts_pol)) + ind = -np.ones((npts_pol, ), dtype=int) + pts_mv = pts + ind_mv = ind + reso_rdrdz_mv = reso_rdrdz + reso_r_z = reso_r[0]*reso_z[0] + ind_i = np.sort(ind_i, axis=1) + indi_mv = ind_i + first_ind_mv = np.argmax(ind_i > -1, axis=1).astype(int) + # initializing utilitary arrays + num_threads = _ompt.get_effective_num_threads(num_threads) + lstruct_lims_np = flatten_lstruct_lims(lstruct_lims) + print("Got heeeeeeere 6...........................") + # .............. + _st.sa_tri_assemble( + block, + approx, + data, + npoly, + poly_lnvert, + ltri, poly_lnorms_tot, + centroids, + vec_GB, + vec_GC, + cross_GBGC, + dot_GBGC, + is_in_vignette, + sa_map, + ves_poly, + ves_norm, + ves_lims, + lstruct_nlim, + lstruct_polyx, + lstruct_polyy, + lstruct_lims_np, + lstruct_normx, + lstruct_normy, + lnvert, + nstruct_tot, + nstruct_lim, + rmin, + eps_uz, eps_a, + eps_vz, eps_b, + eps_plane, + forbid, + first_ind_mv, + indi_mv, + tot_num_tri, sz_r, sz_z, + ncells_rphi, + reso_r_z, + disc_r, + step_rphi, + disc_z, + ind_rz2pol, + sz_phi, + reso_rdrdz_mv, + pts_mv, + ind_mv, + num_threads + ) + print("Got heeeeeeere 7...........................") + # ... freeing up memory .................................................... + free(lindex_z) + free(disc_r) + free(disc_z) + free(disc_r0) + free(sz_phi) + free(step_rphi) + free(ncells_rphi) + free(data) + return pts, sa_map, ind, reso_r_z diff --git a/tofu/geom/_basic_geom_tools.pxd b/tofu/geom/_basic_geom_tools.pxd index 4047fd300..4c7604f7b 100644 --- a/tofu/geom/_basic_geom_tools.pxd +++ b/tofu/geom/_basic_geom_tools.pxd @@ -66,6 +66,49 @@ cdef void compute_cross_prod(const double[3] vec_a, const double[3] vec_b, double[3] res) nogil +cdef void compute_dot_cross_vec(const double[:, ::1] lvec_a, + const double[:, ::1] lvec_b, + double[:, ::1] cross_p, + double[::1] dot_p, + const int npts, + const int num_threads, + ) nogil + +cdef void find_centroids_ltri(const double[:, :, ::1] poly_coords, + const long** ltri, + const long* lnvert, + const int npoly, + const int num_threads, + double[:, ::1] centroid) nogil + +cdef void find_centroids_GB_GC_ltri(const double** poly_coords, + const long** ltri, + const long* lnvert, + const int npoly, + const int num_threads, + double[:, ::1] centroid, + double[:, ::1] vec_GB, + double[:, ::1] vec_GC, + ) nogil + +cdef void compute_vec_ass_tri(const double pt0, const double pt1, + const double pt2, int npts, + const double[:, ::1] ptG, + const double[:, ::1] poly_norm, + const double[:, ::1] cross_bc, + const double[:, ::1] vecb, + const double[:, ::1] vecc, + double* side_of_poly, + double* num, + double* dot_Gb, + double* dot_Gc, + double* normG2) nogil + +cdef void compute_dist_pt_arr(const double pt0, const double pt1, + const double pt2, int npts, + const double* vec, + double* dist) nogil + cdef void compute_dist_pt_vec(const double pt0, const double pt1, const double pt2, const int npts, const double[:, ::1] vec, diff --git a/tofu/geom/_basic_geom_tools.pyx b/tofu/geom/_basic_geom_tools.pyx index 2fd91e399..17fe5719f 100644 --- a/tofu/geom/_basic_geom_tools.pyx +++ b/tofu/geom/_basic_geom_tools.pyx @@ -4,12 +4,15 @@ # cython: cdivision=True # cython: initializedcheck=False # -from cython.parallel import prange -from cpython.array cimport array, clone from libc.math cimport sqrt as c_sqrt from libc.math cimport NAN as CNAN from libc.math cimport pi as c_pi -from libc.stdlib cimport malloc, free +from libc.stdlib cimport free +from libc.stdlib cimport malloc +from cpython.array cimport clone +from cpython.array cimport array +from cython.parallel import prange +from cython.parallel cimport parallel # cdef double _VSMALL = 1.e-9 cdef double _SMALL = 1.e-6 @@ -54,6 +57,7 @@ cdef inline bint is_point_in_path(const int nvert, c = not c return c + cdef inline int is_point_in_path_vec(const int nvert, const double* vertx, const double* verty, @@ -184,6 +188,33 @@ cdef inline double compute_dot_prod(const double[3] vec_a, return vec_a[0] * vec_b[0] + vec_a[1] * vec_b[1] + vec_a[2] * vec_b[2] +cdef inline void compute_dot_cross_vec(const double[:, ::1] lvec_a, + const double[:, ::1] lvec_b, + double[:, ::1] cross_p, + double[::1] dot_p, + const int npts, + const int num_threads, + ) nogil: + cdef int ii + cdef double vec1[3] + cdef double vec2[3] + cdef double vec3[3] + + with nogil, parallel(num_threads=num_threads): + for ii in prange(npts): + vec1[0] = lvec_a[0, ii] + vec1[1] = lvec_a[1, ii] + vec1[2] = lvec_a[2, ii] + vec2[0] = lvec_b[0, ii] + vec2[1] = lvec_b[1, ii] + vec2[2] = lvec_b[2, ii] + compute_cross_prod(vec1, vec2, &vec3[0]) + cross_p[0, ii] = vec3[0] + cross_p[1, ii] = vec3[1] + cross_p[2, ii] = vec3[2] + dot_p[ii] = compute_dot_prod(vec1, vec2) + return + cdef inline double compute_norm(const double[3] vec) nogil: return c_sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]) @@ -287,9 +318,136 @@ cdef inline int find_ind_lowerright_corner(const double[::1] xpts, return res +cdef inline void find_centroid_tri(const double[3] xpts, + const double[3] ypts, + const double[3] zpts, + double[:] centroid) nogil: + centroid[0] = (xpts[0] + xpts[1] + xpts[2]) / 3. + centroid[1] = (ypts[0] + ypts[1] + ypts[2]) / 3. + centroid[2] = (zpts[0] + zpts[1] + zpts[2]) / 3. + return + + +cdef inline void find_centroids_ltri(const double[:, :, ::1] poly_coords, + const long** ltri, + const long* lnvert, + const int npoly, + const int num_threads, + double[:, ::1] centroid) nogil: + cdef int ipol + cdef int itri + cdef int wim1 + cdef int wi + cdef int wip1 + cdef double[3] xtri + cdef double[3] ytri + cdef double[3] ztri + + with nogil, parallel(num_threads=num_threads): + for ipol in prange(npoly): + for itri in range(lnvert[ipol] - 2): + wim1 = ltri[ipol][itri*3] + wi = ltri[ipol][itri*3+1] + wip1 = ltri[ipol][itri*3+2] + xtri[0] = poly_coords[ipol, 0, wim1] + ytri[0] = poly_coords[ipol, 1, wim1] + ztri[0] = poly_coords[ipol, 2, wim1] + xtri[1] = poly_coords[ipol, 0, wi] + ytri[1] = poly_coords[ipol, 1, wi] + ztri[1] = poly_coords[ipol, 2, wi] + xtri[2] = poly_coords[ipol, 0, wip1] + ytri[2] = poly_coords[ipol, 1, wip1] + ztri[2] = poly_coords[ipol, 2, wip1] + find_centroid_tri(xtri, ytri, ztri, + centroid[:, ipol + itri * npoly]) + return + + + +cdef inline void find_centroids_GB_GC_ltri(const double** poly_coords, + const long** ltri, + const long* lnvert, + const int npoly, + const int num_threads, + double[:, ::1] centroid, + double[:, ::1] vec_GB, + double[:, ::1] vec_GC, + ) nogil: + cdef int wi + cdef int wip1 + cdef int wim1 + cdef int ipol + cdef int itri + cdef int iglob + cdef double* xtri + cdef double* ytri + cdef double* ztri + + with nogil, parallel(num_threads=num_threads): + xtri = malloc(3*sizeof(double)) + ytri = malloc(3*sizeof(double)) + ztri = malloc(3*sizeof(double)) + for ipol in prange(npoly): + with gil: + print(f"{ipol} / ({npoly})") + for itri in range(lnvert[ipol] - 2): + with gil: + print(f"{itri}/({lnvert[ipol]}-2)") + print(ltri == NULL) + print(ltri[0] == NULL) + wim1 = ltri[ipol][itri*3] + wi = ltri[ipol][itri*3+1] + wip1 = ltri[ipol][itri*3+2] + with gil: + print("xxxxxxxxxxxxx 1", wim1, wi, wip1) + print(lnvert == NULL) + print(lnvert[ipol]) + print(poly_coords == NULL) + print(poly_coords[0] == NULL) + xtri[0] = poly_coords[ipol][0 * lnvert[ipol] + wim1] + ytri[0] = poly_coords[ipol][1 * lnvert[ipol] + wim1] + ztri[0] = poly_coords[ipol][2 * lnvert[ipol] + wim1] + xtri[1] = poly_coords[ipol][0 * lnvert[ipol] + wi] + ytri[1] = poly_coords[ipol][1 * lnvert[ipol] + wi] + ztri[1] = poly_coords[ipol][2 * lnvert[ipol] + wi] + xtri[2] = poly_coords[ipol][0 * lnvert[ipol] + wip1] + ytri[2] = poly_coords[ipol][1 * lnvert[ipol] + wip1] + ztri[2] = poly_coords[ipol][2 * lnvert[ipol] + wip1] + iglob = ipol + itri * npoly + with gil: + print("...........before find centroid") + find_centroid_tri(xtri, ytri, ztri, + centroid[:, iglob]) + with gil: + print("...........after find centroid") + vec_GB[0, iglob] = xtri[1] - centroid[0, iglob] + vec_GB[1, iglob] = ytri[1] - centroid[1, iglob] + vec_GB[2, iglob] = ztri[1] - centroid[2, iglob] + vec_GC[0, iglob] = xtri[2] - centroid[0, iglob] + vec_GC[1, iglob] = ytri[2] - centroid[1, iglob] + vec_GC[2, iglob] = ztri[2] - centroid[2, iglob] + return + # ============================================================================== # = Distance # ============================================================================== +cdef inline void compute_dist_pt_arr(const double pt0, const double pt1, + const double pt2, int npts, + const double* vec, + double* dist) nogil: + """ + Compute the distance between the point P = [pt0, pt1, pt2] and each point + Q_i, where vec = {Q_0, Q_1, ..., Q_npts-1} + """ + cdef int ii + for ii in range(0, npts): + dist[ii] = c_sqrt( + (pt0 - vec[0 * npts + ii]) * (pt0 - vec[0 * npts + ii]) + + (pt1 - vec[1 * npts + ii]) * (pt1 - vec[1 * npts + ii]) + + (pt2 - vec[2 * npts + ii]) * (pt2 - vec[2 * npts + ii])) + return + + cdef inline void compute_dist_pt_vec(const double pt0, const double pt1, const double pt2, int npts, const double[:, ::1] vec, @@ -306,6 +464,55 @@ cdef inline void compute_dist_pt_vec(const double pt0, const double pt1, return +cdef inline void compute_vec_ass_tri(const double pt0, const double pt1, + const double pt2, int npts, + const double[:, ::1] ptG, + const double[:, ::1] poly_norm, + const double[:, ::1] cross_bc, + const double[:, ::1] vecb, + const double[:, ::1] vecc, + double* side_of_poly, + double* num, + double* dot_Gb, + double* dot_Gc, + double* normG2) nogil: + """ + Computes: + - numerator of triangle SA comp : 3 G \dot (b \cross c) + - G \dot b = OG \dot Gb + - G \dot c = OG \dot Gc + - (norm G)**2 = G \dot G = OG \dot \OG + ptG are the coordinates of all G, centroids of all triangles + pt0, pt1, pt2 are the coordinates of O + cross_bc are b x c for all triangles + vecb = \vec Gb + vecc = \vec Gc + """ + cdef int ii + cdef double[3] vec_OG + + for ii in range(npts): + vec_OG[0] = ptG[0, ii] - pt0 + vec_OG[1] = ptG[1, ii] - pt1 + vec_OG[2] = ptG[2, ii] - pt2 + normG2[ii] = (vec_OG[0] * vec_OG[0] + + vec_OG[1] * vec_OG[1] + + vec_OG[2] * vec_OG[2]) + num[ii] = 3.0 * ( vec_OG[0] * cross_bc[0, ii] + + vec_OG[1] * cross_bc[1, ii] + + vec_OG[2] * cross_bc[2, ii]) + dot_Gb[ii] = (vec_OG[0] * vecb[0, ii] + + vec_OG[1] * vecb[1, ii] + + vec_OG[2] * vecb[2, ii]) + dot_Gc[ii] = (vec_OG[0] * vecc[0, ii] + + vec_OG[1] * vecc[1, ii] + + vec_OG[2] * vecc[2, ii]) + side_of_poly[ii] = (vec_OG[0] * poly_norm[ii, 0] + + vec_OG[1] * poly_norm[ii, 1] + + vec_OG[2] * poly_norm[ii, 2]) + return + + cdef inline void compute_dist_vec_vec(const int npts1, const int npts2, const double[:, ::1] vec1, const double[:, ::1] vec2, diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index 872c45a04..7dfa211d6 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -1970,7 +1970,8 @@ def from_txt( + "\t- 'object': return {} instance\n".format(cls.__name__) + "\t- 'dict' : return a dict with polygon, pos and extent") raise Exception(msg) - if pfe[-4:] not in ['.txt', '.csv']: + + if not (pfe.endswith('.txt') or pfe.endswith('.csv')): msg = ("Only accepts .txt and .csv files (fed to np.loadtxt) !\n" + "\t file: {}".format(pfe)) raise Exception(msg) diff --git a/tofu/geom/_sampling_tools.pxd b/tofu/geom/_sampling_tools.pxd index 25794a699..a5c25c211 100644 --- a/tofu/geom/_sampling_tools.pxd +++ b/tofu/geom/_sampling_tools.pxd @@ -253,7 +253,7 @@ cdef int sa_get_index_arrays(long[:, ::1] lnp, int sz_r, int sz_z) nogil -cdef void sa_assemble_arrays(int block, +cdef void sa_sphere_assemble(int block, int use_approx, double[:, ::1] part_coords, double[::1] part_rad, @@ -291,3 +291,52 @@ cdef void sa_assemble_arrays(int block, double[:, ::1] pts_mv, long[::1] ind_mv, int num_threads) + +cdef void sa_tri_assemble( + int block, + int use_approx, + double** poly_coords, + int npoly, + long[::1] lnvert_poly, + long** ltri, + double[:, ::1] poly_norm, + double[:, ::1] centroids, + double[:, ::1] vec_GB, + double[:, ::1] vec_GC, + double[:, ::1] cross_GBGC, + double[::1] dot_GBGC, + long[:, ::1] is_in_vignette, + double[:, ::1] sa_map, + double[:, ::1] ves_poly, + double[:, ::1] ves_norm, + double[::1] ves_lims, + long[::1] lstruct_nlim, + double[::1] lstruct_polyx, + double[::1] lstruct_polyy, + double[::1] lstruct_lims, + double[::1] lstruct_normx, + double[::1] lstruct_normy, + long[::1] lnvert, + int nstruct_tot, + int nstruct_lim, + double rmin, + double eps_uz, double eps_a, + double eps_vz, double eps_b, + double eps_plane, + bint forbid, + long[::1] first_ind, + long[:, ::1] indi_mv, + int num_tot_tri, + int sz_r, int sz_z, + long* ncells_rphi, + double reso_r_z, + double* disc_r, + double* step_rphi, + double* disc_z, + long[:, ::1] ind_rz2pol, + long* sz_phi, + double[::1] reso_rdrdz, + double[:, ::1] pts_mv, + long[::1] ind_mv, + int num_threads +) diff --git a/tofu/geom/_sampling_tools.pyx b/tofu/geom/_sampling_tools.pyx index 642f7f471..2c35d382c 100644 --- a/tofu/geom/_sampling_tools.pyx +++ b/tofu/geom/_sampling_tools.pyx @@ -1,8 +1,8 @@ # cython: language_level=3 -# cython: boundscheck=False +# cython: boundscheck=True # cython: wraparound=False # cython: cdivision=True -# cython: initializedcheck=False +# cython: initializedcheck=True # ################################################################################ # Utility functions for sampling and discretizating @@ -11,6 +11,7 @@ from libc.math cimport ceil as c_ceil, fabs as c_abs from libc.math cimport floor as c_floor, round as c_round from libc.math cimport sqrt as c_sqrt from libc.math cimport pi as c_pi, cos as c_cos, sin as c_sin +from libc.math cimport atan2 as c_atan2 from libc.math cimport isnan as c_isnan from libc.math cimport NAN as C_NAN from libc.math cimport log2 as c_log2 @@ -19,6 +20,7 @@ from cython.parallel import prange from cython.parallel cimport parallel from cython.parallel cimport threadid from cpython.array cimport array, clone +from libc.stdio cimport printf # for utility functions: import numpy as np cimport numpy as cnp @@ -239,8 +241,8 @@ cdef inline void discretize_vpoly_core(double[:, ::1] ves_poly, double dstep, #.. Filling arrays.......................................................... if c_abs(din) < _VSMALL: for ii in range(np-1): - v0 = ves_poly[0, ii+1]-ves_poly[0, ii] - v1 = ves_poly[1, ii+1]-ves_poly[1, ii] + v0 = ves_poly[0, ii+1] - ves_poly[0, ii] + v1 = ves_poly[1, ii+1] - ves_poly[1, ii] lminmax[1] = c_sqrt(v0 * v0 + v1 * v1) inv_norm = 1. / lminmax[1] discretize_line1d_core(lminmax, dstep, dl_array, True, @@ -2040,7 +2042,10 @@ cdef inline int sa_disc_phi(int sz_r, int sz_z, return npts_disc -cdef inline void sa_assemble_arrays(int block, +# ------------------------------------------------------------------------------ +# -- Solid Angle Computation subtended by a SPHERE +# ------------------------------------------------------------------------------ +cdef inline void sa_sphere_assemble(int block, int use_approx, double[:, ::1] part_coords, double[::1] part_rad, @@ -2124,124 +2129,124 @@ cdef inline void sa_assemble_arrays(int block, if use_approx: # if block and use_approx - assemble_block_approx(part_coords, part_rad, - sa_map, - ves_poly, ves_norm, - ves_lims, - lstruct_nlim, lstruct_polyx, lstruct_polyy, - lstruct_lims, lstruct_normx, lstruct_normy, - lnvert, vperp_out, - coeff_inter_in, coeff_inter_out, - ind_inter_out, sz_ves_lims, - ray_orig, ray_vdir, npts_poly, - nstruct_tot, nstruct_lim, - rmin, - eps_uz, eps_a, - eps_vz, eps_b, eps_plane, - forbid, - first_ind, indi_mv, - sz_p, sz_pol, - ncells_rphi, - disc_r, step_rphi, - disc_z, ind_pol2r, ind_pol2z, sz_phi, - num_threads) - - else: # if block and (not use_approx) - - assemble_block_exact(part_coords, part_rad, - sa_map, - ves_poly, ves_norm, - ves_lims, - lstruct_nlim, - lstruct_polyx, - lstruct_polyy, - lstruct_lims, - lstruct_normx, - lstruct_normy, - lnvert, vperp_out, - coeff_inter_in, coeff_inter_out, - ind_inter_out, sz_ves_lims, - ray_orig, ray_vdir, npts_poly, - nstruct_tot, nstruct_lim, - rmin, - eps_uz, eps_a, - eps_vz, eps_b, eps_plane, - forbid, - first_ind, indi_mv, - sz_p, sz_pol, - ncells_rphi, - disc_r, step_rphi, - disc_z, ind_pol2r, ind_pol2z, sz_phi, - num_threads) - - else: # if not block - - if use_approx: # if (not block) and use_approx - - assemble_unblock_approx(part_coords, part_rad, + sphr_asmbl_block_approx(part_coords, part_rad, sa_map, + ves_poly, ves_norm, + ves_lims, + lstruct_nlim, lstruct_polyx, lstruct_polyy, + lstruct_lims, lstruct_normx, lstruct_normy, + lnvert, vperp_out, + coeff_inter_in, coeff_inter_out, + ind_inter_out, sz_ves_lims, + ray_orig, ray_vdir, npts_poly, + nstruct_tot, nstruct_lim, + rmin, + eps_uz, eps_a, + eps_vz, eps_b, eps_plane, + forbid, first_ind, indi_mv, sz_p, sz_pol, ncells_rphi, disc_r, step_rphi, - disc_z, ind_pol2r, ind_pol2z, - sz_phi, + disc_z, ind_pol2r, ind_pol2z, sz_phi, num_threads) - else: # if (not block) and (not use_approx) + else: # if block and (not use_approx) - assemble_unblock_exact(part_coords, part_rad, + sphr_asmbl_block_exact(part_coords, part_rad, sa_map, + ves_poly, ves_norm, + ves_lims, + lstruct_nlim, + lstruct_polyx, + lstruct_polyy, + lstruct_lims, + lstruct_normx, + lstruct_normy, + lnvert, vperp_out, + coeff_inter_in, coeff_inter_out, + ind_inter_out, sz_ves_lims, + ray_orig, ray_vdir, npts_poly, + nstruct_tot, nstruct_lim, + rmin, + eps_uz, eps_a, + eps_vz, eps_b, eps_plane, + forbid, first_ind, indi_mv, sz_p, sz_pol, ncells_rphi, disc_r, step_rphi, - disc_z, ind_pol2r, ind_pol2z, - sz_phi, + disc_z, ind_pol2r, ind_pol2z, sz_phi, num_threads) + + else: # if not block + + if use_approx: # if (not block) and use_approx + + sphr_asmbl_unblock_approx(part_coords, part_rad, + sa_map, + first_ind, indi_mv, + sz_p, sz_pol, + ncells_rphi, + disc_r, step_rphi, + disc_z, ind_pol2r, ind_pol2z, + sz_phi, + num_threads) + + else: # if (not block) and (not use_approx) + + sphr_asmbl_unblock_exact(part_coords, part_rad, + sa_map, + first_ind, indi_mv, + sz_p, sz_pol, + ncells_rphi, + disc_r, step_rphi, + disc_z, ind_pol2r, ind_pol2z, + sz_phi, + num_threads) return -cdef inline void assemble_block_approx(double[:, ::1] part_coords, - double[::1] part_rad, - double[:, ::1] sa_map, - double[:, ::1] ves_poly, - double[:, ::1] ves_norm, - double[::1] ves_lims, - long[::1] lstruct_nlim, - double[::1] lstruct_polyx, - double[::1] lstruct_polyy, - double[::1] lstruct_lims, - double[::1] lstruct_normx, - double[::1] lstruct_normy, - long[::1] lnvert, - double[:, ::1] vperp_out, - double[:, ::1] coeff_inter_in, - double[:, ::1] coeff_inter_out, - int[:, ::1] ind_inter_out, - int sz_ves_lims, - double[:, :, ::1] ray_orig, - double[:, :, ::1] ray_vdir, - int npts_poly, - int nstruct_tot, - int nstruct_lim, - double rmin, - double eps_uz, double eps_a, - double eps_vz, double eps_b, - double eps_plane, - bint forbid, - long[::1] first_ind_mv, - long[:, ::1] indi_mv, - int sz_p, - int sz_pol, - long* ncells_rphi, - double* disc_r, - double* step_rphi, - double* disc_z, - int[::1] ind_pol2r, - int[::1] ind_pol2z, - long* sz_phi, - int num_threads) nogil: +cdef inline void sphr_asmbl_block_approx(double[:, ::1] part_coords, + double[::1] part_rad, + double[:, ::1] sa_map, + double[:, ::1] ves_poly, + double[:, ::1] ves_norm, + double[::1] ves_lims, + long[::1] lstruct_nlim, + double[::1] lstruct_polyx, + double[::1] lstruct_polyy, + double[::1] lstruct_lims, + double[::1] lstruct_normx, + double[::1] lstruct_normy, + long[::1] lnvert, + double[:, ::1] vperp_out, + double[:, ::1] coeff_inter_in, + double[:, ::1] coeff_inter_out, + int[:, ::1] ind_inter_out, + int sz_ves_lims, + double[:, :, ::1] ray_orig, + double[:, :, ::1] ray_vdir, + int npts_poly, + int nstruct_tot, + int nstruct_lim, + double rmin, + double eps_uz, double eps_a, + double eps_vz, double eps_b, + double eps_plane, + bint forbid, + long[::1] first_ind_mv, + long[:, ::1] indi_mv, + int sz_p, + int sz_pol, + long* ncells_rphi, + double* disc_r, + double* step_rphi, + double* disc_z, + int[::1] ind_pol2r, + int[::1] ind_pol2z, + long* sz_phi, + int num_threads) nogil: cdef int rr cdef int zz cdef int jj @@ -2299,7 +2304,7 @@ cdef inline void assemble_block_approx(double[:, ::1] part_coords, lstruct_lims, lstruct_normx, lstruct_normy, - lnvert, + lnvert, vperp_out[thid], coeff_inter_in[thid], coeff_inter_out[thid], @@ -2317,7 +2322,7 @@ cdef inline void assemble_block_approx(double[:, ::1] part_coords, for pp in range(sz_p): if is_vis[pp] and dist[pp] > part_rad[pp]: sa_map[ind_pol, - pp] += sa_approx_formula(part_rad[pp], + pp] += comp_sa_sphr_appx(part_rad[pp], dist[pp], vol_pi) free(dist) @@ -2326,21 +2331,21 @@ cdef inline void assemble_block_approx(double[:, ::1] part_coords, return -cdef inline void assemble_unblock_approx(double[:, ::1] part_coords, - double[::1] part_rad, - double[:, ::1] sa_map, - long[::1] first_ind_mv, - long[:, ::1] indi_mv, - int sz_p, - int sz_pol, - long* ncells_rphi, - double* disc_r, - double* step_rphi, - double* disc_z, - int[::1] ind_pol2r, - int[::1] ind_pol2z, - long* sz_phi, - int num_threads) nogil: +cdef inline void sphr_asmbl_unblock_approx(double[:, ::1] part_coords, + double[::1] part_rad, + double[:, ::1] sa_map, + long[::1] first_ind_mv, + long[:, ::1] indi_mv, + int sz_p, + int sz_pol, + long* ncells_rphi, + double* disc_r, + double* step_rphi, + double* disc_z, + int[::1] ind_pol2r, + int[::1] ind_pol2z, + long* sz_phi, + int num_threads) nogil: cdef int rr cdef int zz cdef int jj @@ -2383,14 +2388,14 @@ cdef inline void assemble_unblock_approx(double[:, ::1] part_coords, for pp in range(sz_p): if dist[pp] > part_rad[pp]: sa_map[ind_pol, - pp] += sa_approx_formula(part_rad[pp], + pp] += comp_sa_sphr_appx(part_rad[pp], dist[pp], vol_pi) free(dist) return -cdef inline double sa_approx_formula(double radius, +cdef inline double comp_sa_sphr_appx(double radius, double distance, double volpi, int debug=0) nogil: @@ -2423,46 +2428,46 @@ cdef inline double sa_approx_formula(double radius, # ----------------------------------------------------------------------------- # Exact formula computation # ----------------------------------------------------------------------------- -cdef inline void assemble_block_exact(double[:, ::1] part_coords, - double[::1] part_rad, - double[:, ::1] sa_map, - double[:, ::1] ves_poly, - double[:, ::1] ves_norm, - double[::1] ves_lims, - long[::1] lstruct_nlim, - double[::1] lstruct_polyx, - double[::1] lstruct_polyy, - double[::1] lstruct_lims, - double[::1] lstruct_normx, - double[::1] lstruct_normy, - long[::1] lnvert, - double[:, ::1] vperp_out, - double[:, ::1] coeff_inter_in, - double[:, ::1] coeff_inter_out, - int[:, ::1] ind_inter_out, - int sz_ves_lims, - double[:, :, ::1] ray_orig, - double[:, :, ::1] ray_vdir, - int npts_poly, - int nstruct_tot, - int nstruct_lim, - double rmin, - double eps_uz, double eps_a, - double eps_vz, double eps_b, - double eps_plane, - bint forbid, - long[::1] first_ind_mv, - long[:, ::1] indi_mv, - int sz_p, - int sz_pol, - long* ncells_rphi, - double* disc_r, - double* step_rphi, - double* disc_z, - int[::1] ind_pol2r, - int[::1] ind_pol2z, - long* sz_phi, - int num_threads) nogil: +cdef inline void sphr_asmbl_block_exact(double[:, ::1] part_coords, + double[::1] part_rad, + double[:, ::1] sa_map, + double[:, ::1] ves_poly, + double[:, ::1] ves_norm, + double[::1] ves_lims, + long[::1] lstruct_nlim, + double[::1] lstruct_polyx, + double[::1] lstruct_polyy, + double[::1] lstruct_lims, + double[::1] lstruct_normx, + double[::1] lstruct_normy, + long[::1] lnvert, + double[:, ::1] vperp_out, + double[:, ::1] coeff_inter_in, + double[:, ::1] coeff_inter_out, + int[:, ::1] ind_inter_out, + int sz_ves_lims, + double[:, :, ::1] ray_orig, + double[:, :, ::1] ray_vdir, + int npts_poly, + int nstruct_tot, + int nstruct_lim, + double rmin, + double eps_uz, double eps_a, + double eps_vz, double eps_b, + double eps_plane, + bint forbid, + long[::1] first_ind_mv, + long[:, ::1] indi_mv, + int sz_p, + int sz_pol, + long* ncells_rphi, + double* disc_r, + double* step_rphi, + double* disc_z, + int[::1] ind_pol2r, + int[::1] ind_pol2z, + long* sz_phi, + int num_threads) nogil: cdef int rr cdef int zz cdef int jj @@ -2539,7 +2544,7 @@ cdef inline void assemble_block_exact(double[:, ::1] part_coords, for pp in range(sz_p): if is_vis[pp] and dist[pp] > part_rad[pp]: sa_map[ind_pol, - pp] += sa_exact_formula(part_rad[pp], + pp] += comp_sa_sphr_ext(part_rad[pp], dist[pp], vol_pi) free(dist) @@ -2547,21 +2552,21 @@ cdef inline void assemble_block_exact(double[:, ::1] part_coords, return -cdef inline void assemble_unblock_exact(double[:, ::1] part_coords, - double[::1] part_rad, - double[:, ::1] sa_map, - long[::1] first_ind_mv, - long[:, ::1] indi_mv, - int sz_p, - int sz_pol, - long* ncells_rphi, - double* disc_r, - double* step_rphi, - double* disc_z, - int[::1] ind_pol2r, - int[::1] ind_pol2z, - long* sz_phi, - int num_threads) nogil: +cdef inline void sphr_asmbl_unblock_exact(double[:, ::1] part_coords, + double[::1] part_rad, + double[:, ::1] sa_map, + long[::1] first_ind_mv, + long[:, ::1] indi_mv, + int sz_p, + int sz_pol, + long* ncells_rphi, + double* disc_r, + double* step_rphi, + double* disc_z, + int[::1] ind_pol2r, + int[::1] ind_pol2z, + long* sz_phi, + int num_threads) nogil: cdef int rr cdef int zz cdef int jj @@ -2604,14 +2609,14 @@ cdef inline void assemble_unblock_exact(double[:, ::1] part_coords, for pp in range(sz_p): if dist[pp] > part_rad[pp]: sa_map[ind_pol, - pp] += sa_exact_formula(part_rad[pp], + pp] += comp_sa_sphr_ext(part_rad[pp], dist[pp], vol_pi) free(dist) return -cdef inline double sa_exact_formula(double radius, +cdef inline double comp_sa_sphr_ext(double radius, double distance, double volpi) nogil: """ @@ -2635,3 +2640,478 @@ cdef inline double sa_exact_formula(double radius, cdef double r_over_d = radius / distance return 2 * volpi * (1. - c_sqrt(1. - r_over_d**2)) + + + +# ------------------------------------------------------------------------------ +# -- Solid Angle Computation subtended by a POLYGON +# ------------------------------------------------------------------------------ +cdef inline void sa_tri_assemble( + int block, + int use_approx, + double** poly_coords, + int npoly, + long[::1] lnvert_poly, + long** ltri, + double[:, ::1] poly_norm, + double[:, ::1] centroids, + double[:, ::1] vec_GB, + double[:, ::1] vec_GC, + double[:, ::1] cross_GBGC, + double[::1] dot_GBGC, + long[:, ::1] is_in_vignette, + double[:, ::1] sa_map, + double[:, ::1] ves_poly, + double[:, ::1] ves_norm, + double[::1] ves_lims, + long[::1] lstruct_nlim, + double[::1] lstruct_polyx, + double[::1] lstruct_polyy, + double[::1] lstruct_lims, + double[::1] lstruct_normx, + double[::1] lstruct_normy, + long[::1] lnvert, + int nstruct_tot, + int nstruct_lim, + double rmin, + double eps_uz, double eps_a, + double eps_vz, double eps_b, + double eps_plane, + bint forbid, + long[::1] first_ind, + long[:, ::1] indi_mv, + int num_tot_tri, + int sz_r, int sz_z, + long* ncells_rphi, + double reso_r_z, + double* disc_r, + double* step_rphi, + double* disc_z, + long[:, ::1] ind_rz2pol, + long* sz_phi, + double[::1] reso_rdrdz, + double[:, ::1] pts_mv, + long[::1] ind_mv, + int num_threads +): + if block and use_approx: + # .. useless tabs ..................................................... + # declared here so that cython can run without gil + if ves_lims is not None: + sz_ves_lims = np.size(ves_lims) + else: + sz_ves_lims = 0 + npts_poly = ves_norm.shape[1] + ray_orig = np.zeros((3, num_tot_tri)) + ray_vdir = np.zeros((3, num_tot_tri)) + vperp_out = clone(array('d'), num_tot_tri * 3, True) + coeff_inter_in = clone(array('d'), num_tot_tri, True) + coeff_inter_out = clone(array('d'), num_tot_tri, True) + ind_inter_out = clone(array('i'), num_tot_tri * 3, True) + + tri_asmbl_block_approx( + poly_coords, + npoly, + lnvert_poly, + ltri, + poly_norm, + centroids, + vec_GB, + vec_GC, + cross_GBGC, + dot_GBGC, + is_in_vignette, + sa_map, + ves_poly, ves_norm, + ves_lims, + lstruct_nlim, + lstruct_polyx, + lstruct_polyy, + lstruct_lims, + lstruct_normx, + lstruct_normy, + lnvert, vperp_out, + coeff_inter_in, coeff_inter_out, + ind_inter_out, sz_ves_lims, + ray_orig, ray_vdir, npts_poly, + nstruct_tot, nstruct_lim, + rmin, eps_uz, eps_a, + eps_vz, eps_b, eps_plane, + forbid, + first_ind, indi_mv, + num_tot_tri, sz_r, sz_z, + ncells_rphi, + reso_r_z, disc_r, step_rphi, + disc_z, ind_rz2pol, sz_phi, + reso_rdrdz, pts_mv, ind_mv, + num_threads) + elif not block and use_approx: + tri_asmbl_unblock_approx( + poly_coords, + npoly, + lnvert_poly, + ltri, + poly_norm, + centroids, + vec_GB, + vec_GC, + cross_GBGC, + dot_GBGC, + is_in_vignette, + sa_map, + first_ind, + indi_mv, + num_tot_tri, + sz_r, sz_z, + ncells_rphi, + reso_r_z, + disc_r, + step_rphi, + disc_z, + ind_rz2pol, + sz_phi, + reso_rdrdz, + pts_mv, + ind_mv, + num_threads + ) + return + + +cdef inline void tri_asmbl_block_approx( + double** poly_coords, + int npoly, + long[::1] lnvert_poly, + long** ltri, + double[:, ::1] poly_norm, + double[:, ::1] centroids, + double[:, ::1] vec_GB, + double[:, ::1] vec_GC, + double[:, ::1] cross_GBGC, + double[::1] dot_GBGC, + long[:, ::1] is_in_vignette, + double[:, ::1] sa_map, + double[:, ::1] ves_poly, + double[:, ::1] ves_norm, + double[::1] ves_lims, + long[::1] lstruct_nlim, + double[::1] lstruct_polyx, + double[::1] lstruct_polyy, + double[::1] lstruct_lims, + double[::1] lstruct_normx, + double[::1] lstruct_normy, + long[::1] lnvert, + double[::1] vperp_out, + double[::1] coeff_inter_in, + double[::1] coeff_inter_out, + int[::1] ind_inter_out, + int sz_ves_lims, + double[:, ::1] ray_orig, + double[:, ::1] ray_vdir, + int npts_poly, + int nstruct_tot, + int nstruct_lim, + double rmin, + double eps_uz, double eps_a, + double eps_vz, double eps_b, + double eps_plane, + bint forbid, + long[::1] first_ind_mv, + long[:, ::1] indi_mv, + int num_tot_tri, + int sz_r, int sz_z, + long* ncells_rphi, + double reso_r_z, + double* disc_r, + double* step_rphi, + double* disc_z, + long[:, ::1] ind_rz2pol, + long* sz_phi, + double[::1] reso_rdrdz, + double[:, ::1] pts_mv, + long[::1] ind_mv, + int num_threads +) nogil: + cdef int rr + cdef int zz + cdef int jj + cdef int itri + cdef int ipoly + cdef int ind_pol + cdef int loc_first_ind + cdef int loc_size_phi + cdef long indiijj + cdef double loc_x + cdef double loc_y + cdef double loc_r + cdef double loc_z + cdef double loc_phi + cdef double loc_step_rphi + cdef long* is_vis + cdef double* dot_Gb = NULL + cdef double* dot_Gc = NULL + cdef double* norm_G2 = NULL + cdef double* dist_opts = NULL + cdef double* numerator = NULL + cdef double* side_of_poly = NULL + + is_vis = malloc(num_tot_tri * sizeof(long)) + dot_Gb = malloc(num_tot_tri * sizeof(double)) + dot_Gc = malloc(num_tot_tri * sizeof(double)) + norm_G2 = malloc(num_tot_tri * sizeof(double)) + numerator = malloc(num_tot_tri * sizeof(double)) + side_of_poly = malloc(num_tot_tri * sizeof(double)) + + for rr in range(sz_r): + loc_r = disc_r[rr] + loc_size_phi = sz_phi[rr] + loc_step_rphi = step_rphi[rr] + loc_first_ind = first_ind_mv[rr] + for zz in range(sz_z): + loc_z = disc_z[zz] + if is_in_vignette[rr, zz]: + ind_pol = ind_rz2pol[rr, zz] + ind_mv[ind_pol] = rr * sz_z + zz + reso_rdrdz[ind_pol] = loc_r * reso_r_z + pts_mv[0, ind_pol] = loc_r + pts_mv[1, ind_pol] = loc_z + for jj in range(loc_size_phi): + indiijj = indi_mv[rr, loc_first_ind + jj] + loc_phi = - c_pi + (0.5 + indiijj) * loc_step_rphi + loc_x = loc_r * c_cos(loc_phi) + loc_y = loc_r * c_sin(loc_phi) + # OG norm and other associated values .... + _bgt.compute_vec_ass_tri(loc_x, loc_y, loc_z, + num_tot_tri, centroids, + poly_norm, + cross_GBGC, + vec_GB, vec_GC, + &side_of_poly[0], + &numerator[0], + &dot_Gb[0], + &dot_Gc[0], + &norm_G2[0]) + # checking if visible ..... + _rt.is_visible_pt_vec_core(loc_x, loc_y, loc_z, + centroids, + num_tot_tri, + ves_poly, ves_norm, + &is_vis[0], norm_G2, + ves_lims, + lstruct_nlim, + lstruct_polyx, + lstruct_polyy, + lstruct_lims, + lstruct_normx, + lstruct_normy, + lnvert, vperp_out, + coeff_inter_in, + coeff_inter_out, + ind_inter_out, sz_ves_lims, + ray_orig, ray_vdir, + npts_poly, + nstruct_tot, nstruct_lim, + rmin, + eps_uz, eps_a, + eps_vz, eps_b, eps_plane, + 1, # is toroidal + forbid, 1) + + for ipoly in range(npoly): + if side_of_poly[ipoly] < 0.: + continue # point is not on right side of poly + dist_opts = malloc(lnvert_poly[ipoly] + * sizeof(double)) + # computing distances to vertices of poly (OA, OB, OC) + _bgt.compute_dist_pt_arr(loc_x, loc_y, loc_z, + lnvert_poly[ipoly], + poly_coords[ipoly], + &dist_opts[0]) + for itri in range(lnvert_poly[ipoly] - 2): + iglob = ipoly + itri * npoly + if rr < 2 and zz < 2 and jj < 2: + printf("r, z, p = %f, %f, %f\n", loc_r, loc_z, loc_phi) + printf("x, y, z = %f, %f, %f\n", loc_x, loc_y, loc_z) + printf("Poly %d at the %d itri, iglob = %d => is vis = %d\n", + ipoly, itri, iglob, is_vis[iglob]) + if is_vis[iglob]: + sa_map[ind_pol, + ipoly] += comp_sa_tri_appx( + itri, + ltri[ipoly], + numerator[iglob], + dot_Gb[iglob], + dot_Gc[iglob], + norm_G2[iglob], + dot_GBGC[iglob], + dist_opts, + ) + free(dist_opts) + free(dot_Gb) + free(dot_Gc) + free(norm_G2) + free(is_vis) + free(numerator) + free(side_of_poly) + return + +cdef inline void tri_asmbl_unblock_approx( + double** poly_coords, + int npoly, + long[::1] lnvert_poly, + long** ltri, + double[:, ::1] poly_norm, + double[:, ::1] centroids, + double[:, ::1] vec_GB, + double[:, ::1] vec_GC, + double[:, ::1] cross_GBGC, + double[::1] dot_GBGC, + long[:, ::1] is_in_vignette, + double[:, ::1] sa_map, + long[::1] first_ind_mv, + long[:, ::1] indi_mv, + int num_tot_tri, + int sz_r, int sz_z, + long* ncells_rphi, + double reso_r_z, + double* disc_r, + double* step_rphi, + double* disc_z, + long[:, ::1] ind_rz2pol, + long* sz_phi, + double[::1] reso_rdrdz, + double[:, ::1] pts_mv, + long[::1] ind_mv, + int num_threads +) nogil: + cdef int rr + cdef int zz + cdef int jj + cdef int itri + cdef int ipoly + cdef int ind_pol + cdef int loc_first_ind + cdef int loc_size_phi + cdef long indiijj + cdef double loc_x + cdef double loc_y + cdef double loc_r + cdef double loc_z + cdef double loc_phi + cdef double loc_step_rphi + cdef double* dot_Gb = NULL + cdef double* dot_Gc = NULL + cdef double* norm_G2 = NULL + cdef double* dist_opts = NULL + cdef double* numerator = NULL + cdef double* side_of_poly = NULL + + dot_Gb = malloc(num_tot_tri * sizeof(double)) + dot_Gc = malloc(num_tot_tri * sizeof(double)) + norm_G2 = malloc(num_tot_tri * sizeof(double)) + numerator = malloc(num_tot_tri * sizeof(double)) + side_of_poly = malloc(num_tot_tri * sizeof(double)) + + for rr in range(sz_r): + loc_r = disc_r[rr] + loc_size_phi = sz_phi[rr] + loc_step_rphi = step_rphi[rr] + loc_first_ind = first_ind_mv[rr] + for zz in range(sz_z): + loc_z = disc_z[zz] + if is_in_vignette[rr, zz]: + ind_pol = ind_rz2pol[rr, zz] + ind_mv[ind_pol] = rr * sz_z + zz + reso_rdrdz[ind_pol] = loc_r * reso_r_z + pts_mv[0, ind_pol] = loc_r + pts_mv[1, ind_pol] = loc_z + for jj in range(loc_size_phi): + indiijj = indi_mv[rr, loc_first_ind + jj] + loc_phi = - c_pi + (0.5 + indiijj) * loc_step_rphi + loc_x = loc_r * c_cos(loc_phi) + loc_y = loc_r * c_sin(loc_phi) + # OG norm and other associated values .... + _bgt.compute_vec_ass_tri(loc_x, loc_y, loc_z, + num_tot_tri, centroids, + poly_norm, + cross_GBGC, + vec_GB, vec_GC, + &side_of_poly[0], + &numerator[0], + &dot_Gb[0], + &dot_Gc[0], + &norm_G2[0]) + for ipoly in range(npoly): + if side_of_poly[ipoly] < 0.: + continue # point is not on right side of poly + dist_opts = malloc(lnvert_poly[ipoly] + * sizeof(double)) + # computing distances to vertices of poly (OA, OB, OC) + _bgt.compute_dist_pt_arr(loc_x, loc_y, loc_z, + lnvert_poly[ipoly], + poly_coords[ipoly], + &dist_opts[0]) + for itri in range(lnvert_poly[ipoly] - 2): + iglob = ipoly + itri * npoly + sa_map[ind_pol, + ipoly] += comp_sa_tri_appx( + itri, + ltri[ipoly], + numerator[iglob], + dot_Gb[iglob], + dot_Gc[iglob], + norm_G2[iglob], + dot_GBGC[iglob], + dist_opts, + ) + free(dist_opts) + free(dot_Gb) + free(dot_Gc) + free(norm_G2) + free(numerator) + free(side_of_poly) + return + +cdef inline double comp_sa_tri_appx( + int itri, + long* ltri, + double numerator, + double dot_Gb, + double dot_Gc, + double norm_G2, + double dot_bc, + double* dist_opts, +) nogil: + """ + Given by: + numerator = 3 G \dot (b \cross c) + denominator = G \dot G (norm A + norm B + norm C) + + (c \dot b) (norm A - norm B - norm C) + + (c \dot G) (norm A - norm C) + + (G \dot b) (norm A - norm B) + with G centroid of triangle + """ + cdef double normA + cdef double normB + cdef double normC + cdef double denumerator + + normA = dist_opts[ltri[itri*3]] + normB = dist_opts[ltri[itri*3 + 1]] + normC = dist_opts[ltri[itri*3 + 2]] + + denumerator = ( + normA * normB * normC + + norm_G2 * (normA + normB + normC) + + dot_bc * (normA - normB - normC) + + dot_Gb * (normA - normB) + + dot_Gc * (normA - normC) + ) + + if denumerator < 0. : + denumerator += c_pi + + if denumerator < 0. : + printf(">>>>>> denumerator is still negative: %f", denumerator) + + return c_atan2(c_abs(numerator), denumerator) diff --git a/tofu/geom/_vignetting_tools.pxd b/tofu/geom/_vignetting_tools.pxd index bc7cd77dd..4e6d9e050 100644 --- a/tofu/geom/_vignetting_tools.pxd +++ b/tofu/geom/_vignetting_tools.pxd @@ -29,6 +29,10 @@ cdef void earclipping_poly(double* vignett, bint* lref, int nvert) nogil +cdef void triangulate_poly(double* vignett_poly, + long nvert, + long** ltri) nogil + cdef int triangulate_polys(double** vignett_poly, long* lnvert, int nvign, diff --git a/tofu/geom/_vignetting_tools.pyx b/tofu/geom/_vignetting_tools.pyx index 322852081..922ab428b 100644 --- a/tofu/geom/_vignetting_tools.pyx +++ b/tofu/geom/_vignetting_tools.pyx @@ -26,24 +26,33 @@ from . cimport _sorted_set as _ss cdef inline bint is_reflex(const double[3] u, const double[3] v) nogil: """ + /\ + \ + v\ + \______> + P u Determines if the angle between U and -V is reflex (angle > pi) or not. + We suppose the angle is defined anticlockwise from u to v Warning: note the MINUS in front of V, this was done as this is the only form how we will need this function. but it is NOT general. + Is reflex if + u x -v has a negative 3rd coordinate + https://dai.fmph.uniba.sk/upload/8/89/Gm17_lesson05.pdf """ cdef int ii cdef double sumc cdef double[3] ucrossv # ... - _bgt.compute_cross_prod(u, v, ucrossv) - sumc = 0.0 - for ii in range(3): - # normally it should be a sum, but it is a minus cause is we have (U,-V) - sumc = ucrossv[ii] - return sumc >= 0. + _bgt.compute_cross_prod(u, v, &ucrossv[0]) + # reflexive if u x v 3rd coordinate is negative, since we computed u x -v + # the angle in P is reflexive if the 3rd coordinate is positive + return ucrossv[2] > 0. + cdef inline void compute_diff3d(double* orig, int nvert, double* diff) nogil: + # poly can be open ! cdef int ivert for ivert in range(nvert-1): diff[ivert*3 + 0] = orig[0*nvert+(ivert+1)] - orig[0*nvert+ivert] @@ -55,6 +64,7 @@ cdef inline void compute_diff3d(double* orig, diff[3*(nvert-1) + 2] = orig[2*nvert] - orig[2*nvert+(nvert-1)] return + cdef inline void are_points_reflex(int nvert, double* diff, bint* are_reflex) nogil: @@ -63,10 +73,8 @@ cdef inline void are_points_reflex(int nvert, (angle > pi) or not. """ cdef int ivert - cdef int icoord - cdef double[3] u1, v1, un, vn # .. Computing if reflex or not ........................................... - for ivert in range(1,nvert): + for ivert in range(1, nvert): are_reflex[ivert] = is_reflex(&diff[ivert*3], &diff[(ivert-1)*3]) # doing first point: are_reflex[0] = is_reflex(&diff[0], &diff[(nvert-1)*3]) @@ -113,7 +121,7 @@ cdef inline int get_one_ear(double* polygon, double* diff, bint* lref, _cl.ChainedList* working_index, - int nvert, int orig_nvert) nogil: + int nvert, int orig_nvert) nogil except-1: """ A polygon's "ear" is defined as a triangle of vert_i-1, vert_i, vert_i+1, points on the polygon, where the point vert_i has the following properties: @@ -155,12 +163,12 @@ cdef inline int get_one_ear(double* polygon, # edges of the triangle if (lref[wj] and wj != wim1 and wj != wip1 and wj != wi): if is_pt_in_tri(&diff[wi*3], &diff[wim1*3], - polygon[0*orig_nvert+wi], - polygon[1*orig_nvert+wi], - polygon[2*orig_nvert+wi], - polygon[0*orig_nvert+wj], - polygon[1*orig_nvert+wj], - polygon[2*orig_nvert+wj]): + polygon[0 * orig_nvert + wi], + polygon[1 * orig_nvert + wi], + polygon[2 * orig_nvert + wi], + polygon[0 * orig_nvert + wj], + polygon[1 * orig_nvert + wj], + polygon[2 * orig_nvert + wj]): # We found a point in the triangle, thus is not ear # no need to keep testing.... a_pt_in_tri = True @@ -168,10 +176,19 @@ cdef inline int get_one_ear(double* polygon, # Let's check if there was a point in the triangle.... if not a_pt_in_tri: return i # if not, we found an ear - # if we havent returned, either, there was an error somerwhere + # if we havent returned, there was an error somerwhere with gil: - assert False, "Got here but shouldnt have " - return -1 + for j in range(nvert): + wj = _cl.get_at_pos(working_index, j) + print("coordinates at i =", j, + " ", polygon[0 * orig_nvert + wj], + " ", polygon[1 * orig_nvert + wj], + " ", polygon[2 * orig_nvert + wj], + ) + print("... and btw, original polygon had : ", orig_nvert, + " vertices, now we are working on :", nvert) + raise ValueError("Didn't find a non reflex angle in polygon") + cdef inline void earclipping_poly(double* vignett, long* ltri, @@ -188,9 +205,16 @@ cdef inline void earclipping_poly(double* vignett, # init... cdef int loc_nv = nvert cdef int itri = 0 + cdef int ii cdef int wi, wim1, wip1 cdef int iear cdef _cl.ChainedList* working_index + # .. Checking if poly is a triangle ....................................... + if nvert == 3: + ltri[0] = 0 + ltri[1] = 1 + ltri[2] = 2 + return # .. First computing the edges coodinates ................................. # .. and checking if the angles defined by the edges are reflex or not..... # initialization of working index tab: @@ -237,11 +261,38 @@ cdef inline void earclipping_poly(double* vignett, # ============================================================================== # = Polygons triangulation and Intersection Ray-Poly # ============================================================================== +cdef inline void triangulate_poly(double* vignett_poly, + long nvert, + long** ltri) nogil: + """ + Triangulates a single 3d polygon using the earclipping techinque + https://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf + Returns + ltri: 3*(nvert-2) : + {tri_0_0, tri_0_1, ... tri_0_nvert0} + where tri_i_j are the 3 indices of the vertex forming a sub-triangle + on each vertex (-2) and for each vignett + """ + cdef int ii + cdef double* diff = NULL + cdef bint* lref = NULL + # ... + # -- Defining parallel part ------------------------------------------------ + diff = malloc(3*nvert*sizeof(double)) + lref = malloc(nvert*sizeof(bint)) + ltri[0] = malloc((nvert-2)*3*sizeof(long)) + compute_diff3d(vignett_poly, nvert, &diff[0]) + are_points_reflex(nvert, diff, &lref[0]) + earclipping_poly(vignett_poly, <ri[0][0], + &diff[0], &lref[0], nvert) + return + + cdef inline int triangulate_polys(double** vignett_poly, - long* lnvert, - int nvign, - long** ltri, - int num_threads) nogil except -1: + long* lnvert, + int nvign, + long** ltri, + int num_threads) nogil except -1: """ Triangulates a list 3d polygon using the earclipping techinque https://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf @@ -251,7 +302,7 @@ cdef inline int triangulate_polys(double** vignett_poly, where tri_i_j are the 3 indices of the vertex forming a sub-triangle on each vertex (-2) and for each vignett """ - cdef int ivign, ii + cdef int ivign cdef int nvert cdef double* diff = NULL cdef bint* lref = NULL @@ -277,6 +328,7 @@ cdef inline int triangulate_polys(double** vignett_poly, return 0 + cdef inline bint inter_ray_poly(const double[3] ray_orig, const double[3] ray_vdir, double* vignett, diff --git a/tofu/geom/openmp_enabled.py b/tofu/geom/openmp_enabled.py new file mode 100644 index 000000000..37fa7dbe3 --- /dev/null +++ b/tofu/geom/openmp_enabled.py @@ -0,0 +1,6 @@ +# Autogenerated by Tofu's setup.py on 2021-12-09 14:54:42 +def is_openmp_enabled(): + """ + Determine whether this package was built with OpenMP support. + """ + return True diff --git a/tofu/tests/tests01_geom/test_01_GG.py b/tofu/tests/tests01_geom/test_01_GG.py index c4de572d8..02d00ca7e 100644 --- a/tofu/tests/tests01_geom/test_01_GG.py +++ b/tofu/tests/tests01_geom/test_01_GG.py @@ -1571,6 +1571,12 @@ def test21_which_los_closer_vpoly_vec(): # ============================================================================== def test22_earclipping(): # .. First test ............................................................ + # + # (4,5) +------+(5,5) + # | | + # | | + # (4,4) +------+(5,4) + # ves_poly0 = np.zeros((3, 4)) ves_poly00 = [4, 5, 5, 4] ves_poly01 = [4, 4, 5, 5] @@ -1587,8 +1593,6 @@ def test22_earclipping(): ves_poly1[1] = y1 # ...computing out = GG.triangulate_by_earclipping(ves_poly1) - # out = out.reshape((7, 3)) - # print(out) assert np.allclose(out, [1, 2, 3, 1, 3, 4, 1, 4, 5, 0, 1, 5, 0, 5, 6, 0, 6, 7, 0, 7, 8]) @@ -1603,10 +1607,9 @@ def test22_earclipping(): ves_poly2[2] = z2 # ...computing out = GG.triangulate_by_earclipping(ves_poly2) + assert np.allclose(out, [0, 1, 2, 2, 3, 4, 2, 4, 5, 2, 5, 6, 2, 6, 7, 0, 2, 7, 0, 7, 8, 0, 8, 9]) - # out = out.reshape((npts-2, 3)) - # print(out) def test23_vignetting(): # .. First configuration .................................................... diff --git a/tofu/tests/tests01_geom/test_04_sampling.py b/tofu/tests/tests01_geom/test_04_sampling.py index 9e7358304..d217e0de1 100644 --- a/tofu/tests/tests01_geom/test_04_sampling.py +++ b/tofu/tests/tests01_geom/test_04_sampling.py @@ -1,6 +1,7 @@ """Tests of the functions in `sampling_tools.pxd` or their wrappers found in `tofu.geom`. """ +import itertools import numpy as np import tofu as tf import tofu.geom._GG as GG @@ -185,9 +186,7 @@ def test04_ves_vmesh_lin(): # ============================================================================= # Ves - Solid angles # ============================================================================= - - -def test05_sa_integ_map(ves_poly=VPoly, debug=0): +def test05_sa_integ_map(ves_poly=VPoly, debug=3): import matplotlib.pyplot as plt print() @@ -363,6 +362,202 @@ def test05_sa_integ_map(ves_poly=VPoly, debug=0): equal_nan=True) assert np.allclose(sa_map_py, sa_map_py_ex, atol=1e-14, rtol=1e-16, equal_nan=True) + # ... + return + + +def test06_sa_integ_poly_map(ves_poly=VPoly, debug=3): + import matplotlib.pyplot as plt + print() + + config = tf.load_config("A1") + + # lightening up config by removing some PFCs + for name_pfc in ['BaffleV0', "DivUpV1"]: + config.remove_Struct("PFC", name_pfc) + + kwdargs = config.get_kwdargs_LOS_isVis() + ves_poly = kwdargs["ves_poly"] + + if debug > 2: + config.plot() + fig = plt.gcf() + fig.savefig("configuration") + + # coordonnées en x,y,z: + # Tester : enlever des points au rectangle et puis un rectangle avec plus de points + poly_coords = [ + # np.array([ + # [2.7, 0, -0.4], + # [2.75, 0, -0.4], + # [2.75, 0, -0.1], + # [2.70, 0, -0.1], + # [2.60, 0, -0.1], + # [2.50, 0, -0.1], + # [2.50, 0, -0.4], + # ]).T, # 1st polygon + # np.array([ + # [2.50, 0, -0.1], + # [2.50, 0, -0.4], + # [2.7, 0, -0.4], + # [2.75, 0, -0.4], + # [2.75, 0, -0.1], + # [2.70, 0, -0.1], + # [2.60, 0, -0.1], + # ]).T, # 2nd polygon + # np.array([ + # [2.60, 0, -0.1], + # [2.50, 0, -0.1], + # [2.50, 0, -0.4], + # [2.7, 0, -0.4], + # [2.75, 0, -0.4], + # [2.75, 0, -0.1], + # [2.70, 0, -0.1], + # ]).T, # 3rd polygon + np.array([ + [2.5, 0, -0.35], + [2.65, 0, -0.4], + [2.70, 0, -0.15], + ]).T, # 4th polygon + # np.array([ + # [2.75, 0, -0.1], + # [2.70, 0, -0.1], + # [2.60, 0, -0.1], + # [2.50, 0, -0.1], + # [2.55, 0, -0.35], + # [2.65, 0, -0.35], + # ]).T, # 5th polygon + # np.array([ + # [2.2, 0., 0.25], + # [2.6, 0., 0.25], + # [3.0, 0., 0.25], + # [3.0, 0., 0.50], + # [2.7, 0., 0.50], + # [2.5, 0., 0.50], + # [2.2, 0., 0.50], + # ]).T, # 6th polygon + ] + # print(poly_coords) + # poly_coords = [GG.coord_shift(np.ascontiguousarray(poly), + # in_format="(r,z,phi)", + # out_format="(x,y,z)") for poly in poly_coords] + poly_coords = [np.ascontiguousarray(poly) for poly in poly_coords] + # print(poly_coords) + poly_lnorms = np.array([ + [0, 1., 0], + # [0, 1., 0], + # [0, 1., 0], + # [0, 1., 0], + # [0, 1., 0], + # [0, 1., 0], + ]) + poly_lnvert = np.array([poly.shape[1] for poly in poly_coords]) + limits_r, limits_z = compute_min_max_r_and_z(ves_poly) + + if debug > 0: + for pp in range(len(poly_coords)): + plt.clf() + fig = plt.figure() + ax = plt.subplot(111) + poly = poly_coords[pp] + poly = np.concatenate((poly, poly[:, 0].reshape(-1, 1)), axis=1) + xpoly = poly[0] + zpoly = poly[2] + ax.plot( + xpoly, zpoly, + "r-", marker='o', + linewidth=2, + ) + for xzi, (x,z) in enumerate(zip(xpoly, zpoly)): + #ax.annotate(f"({x},{z})", (x, z)) + ax.annotate(f"{xzi}", (x,z)) + ax.plot() + plt.savefig("poly" + str(pp)) + print("...saved!\n") + + lblock = [False] + lstep_rz = [ + 0.005, + 0.01, + ] + + test_cases = list(itertools.product(lblock, lstep_rz)) + + for (block, step_rz) in test_cases: + + rstep = zstep = step_rz + zstep = step_rz + phistep = 0.01 + DR = [2.45, 2.8] + DZ = [-0.5, 0.] + DPhi = [-0.1, 0.1] + + # -- Getting cython APPROX computation -------------------------------- + res = GG.compute_solid_angle_poly_map( + poly_coords, + poly_lnorms, + poly_lnvert, + rstep, zstep, phistep, + limits_r, limits_z, + DR=DR, DZ=DZ, + DPhi=DPhi, + block=block, + limit_vpoly=ves_poly, + **kwdargs, + ) + pts, sa_map_cy, ind, reso_r_z = res + + # check sizes + npts_ind = np.size(ind) + dim, npts = np.shape(pts) + npts_sa, npoly = np.shape(sa_map_cy) + + if debug > 2: + print(f"sa_map_cy is of size {npts_sa},{npoly}") + + # Checking shapes, sizes, types + assert dim == 2 + assert npoly == len(poly_coords), str(len(poly_coords)) + assert npts_ind == npts + assert npts == npts_sa + assert isinstance(reso_r_z, float) + + if debug > 0: + print(f">>>>>>>>>> THERE IS {npoly=}") + for pp in range(npoly): + print(f"...... {pp} / {npoly}") + plt.clf() + fig = plt.figure() + ax = plt.subplot(111) + print("begining") + im = ax.scatter(pts[0, :], pts[1, :], + marker="s", edgecolors="None", + s=40, c=sa_map_cy[:, pp], + vmin=sa_map_cy[:, pp].min(), + vmax=sa_map_cy[:, pp].max()) + poly = poly_coords[pp] + poly = np.concatenate((poly, poly[:, 0].reshape(-1, 1)), axis=1) + xpoly = poly[0] + print(f"{xpoly=}") + zpoly = poly[2] + ax.plot( + xpoly, zpoly, + "r-", marker='o', + linewidth=2, + ) + print("so far so good") + for xzi, (x, z) in enumerate(zip(xpoly, zpoly)): + #ax.annotate(f"({x},{z})", (x, z)) + ax.annotate(f"{xzi}", (x,z)) + ax.plot() + print("sfsg 2") + ax.set_title("cython function") + fig.colorbar(im, ax=ax) + print("sfsg 3") + plt.savefig("sa_map_poly" + str(pp) + + "_block" + str(block) + + "_steprz" + str(step_rz).replace(".", "_")) + print("...saved!\n") # ... return