pvlib_gh656

gh656.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Deep dive into pvlib python GH656\n",
    "This [issue](https://github.com/pvlib/pvlib-python/issues/656) was about `NaN` returned when the sun is still above the horizon. The [patch](https://github.com/pvlib/pvlib-python/pull/697) was a change to this line:\n",
    "\n",
    "```diff\n",
    "- temp = np.minimum(axes_distance*cosd(wid), 1\n",
    "+ temp = np.clip(axes_distance*cosd(wid), -1, 1)\n",
    "```\n",
    "\n",
    "The test case was a low sun angle:\n",
    "\n",
    "| solar zenith | solar azimuth | axis tilt | axis azimuth | max angle | backtrack |  gcr |\n",
    "|--------------|---------------|-----------|--------------|-----------|-----------|------|\n",
    "| 80           | 338           | 30        | 180          | 60        | True      | 0.35 |\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# thinking about gh656\n",
    "%matplotlib inline\n",
    "import copy\n",
    "import pvlib\n",
    "import shapely\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "from shapely.geometry.polygon import LinearRing\n",
    "from shapely import affinity\n",
    "from shapely.geometry import LineString\n",
    "import matplotlib as mpl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'tracker_theta': array([-50.31051385]),\n",
       " 'aoi': array([61.35300178]),\n",
       " 'surface_azimuth': array([112.53615425]),\n",
       " 'surface_tilt': array([56.42233095])}"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check the test case\n",
    "low_sun = dict(\n",
    "    apparent_zenith=80, apparent_azimuth=338, axis_tilt=30,\n",
    "    axis_azimuth=180, max_angle=60, backtrack=True, gcr=0.35)\n",
    "result_back60 = pvlib.tracking.singleaxis(**low_sun)\n",
    "result_back60"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test Case\n",
    "With the patch, the test case now returns -50.3[deg] rotation and an AOI of 61.4[deg].\n",
    "\n",
    "* This means that the trackers backtracked from facing west past zero, and are now facing east.\n",
    "* This AOI is actually for the back side of the PV surface which is still facing west.\n",
    "\n",
    "## New Test Case\n",
    "So what would happen if we removed backtracking and the rotation limits. Lets' do it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'tracker_theta': array([129.68948615]),\n",
       " 'aoi': array([61.35300178]),\n",
       " 'surface_azimuth': array([292.53615425]),\n",
       " 'surface_tilt': array([56.42233095])}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# no backtracking and no rotation limit, max_angle=180\n",
    "low_sun_noback180 = copy.copy(low_sun)\n",
    "low_sun_noback180.update(max_angle=180, backtrack=False)\n",
    "result_noback180 = pvlib.tracking.singleaxis(**low_sun_noback180)\n",
    "result_noback180"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### No Backtracking or Rotation Limits\n",
    "So the AOI is also 61.4[deg]. That's how I knew it couldn't be the AOI for both test cases. And look at the tracker rotation, that's insane. I never ever thought that a tracker would turn past 90[deg]! What does this even mean? Why would the trackers turn so far they're practically facing down?\n",
    "\n",
    "## Tilted Trackers\n",
    "Remember that the trackers are tilted 30[deg], and we are looking at the trackers in their refernce frame, not the global. Let's check the solar vector to make sure this really does make sense. A quick sanity check on the solar angles should help"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-0.3689154775254824, 0.9130978484451157, 0.17364817766693041]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# solar vector\n",
    "x_sun = (np.sin(np.radians(low_sun['apparent_zenith']))\n",
    "     * np.sin(np.radians(low_sun['apparent_azimuth'])))\n",
    "y_sun = (np.sin(np.radians(low_sun['apparent_zenith']))\n",
    "     * np.cos(np.radians(low_sun['apparent_azimuth'])))\n",
    "z_sun = np.cos(np.radians(low_sun['apparent_zenith']))\n",
    "sv = [x_sun, y_sun, z_sun]\n",
    "sv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "338.0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# azimuth? CHECK!\n",
    "np.degrees(np.arctan2(x_sun, y_sun)) % 360"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.767631730218922"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# is sun higher than tracker tilt? NO!\n",
    "np.degrees(np.arctan2(z_sun, y_sun))\n",
    "# the track is tilted 30-degress\n",
    "# so the sun is below the tracker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "tracker_theta = np.radians(result_noback180['tracker_theta'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-50.31051385])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# they are exactly PI apart - interesting\n",
    "result_noback180['tracker_theta']-180"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.63862662])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# should it backtrack?\n",
    "lrot_noback180 = np.cos(np.radians(result_noback180['tracker_theta']))\n",
    "lrot_noback180"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "180-R = [50.31051385]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([0.63862662])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# should it backtrack?\n",
    "print(f'180-R = {180-result_noback180[\"tracker_theta\"]}')\n",
    "lrot_noback180 = np.cos(np.radians(180-result_noback180['tracker_theta']))\n",
    "lrot_noback180"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.82464749])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lrot_noback180/low_sun_noback180['gcr']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.141592653589793"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.arccos(-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pitch = 4.571428571428572\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[(0.5109012967999508, -0.6156134054161985),\n",
       " (4.571428571428572, 3.793856281918045)]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD7CAYAAAB37B+tAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3xUVf7/8deZZCZt0icEUuklBASiiGIBK4qKvWIDZRVde93V3d/uWnZdu6vurmIv2BuyYgFUFESakBBAWiqQXiZtMjPn98fEr6hA2p3cmeTzfDzmQWAm575JZj5z595zP0dprRFCCBG8LGYHEEII0T1SyIUQIshJIRdCiCAnhVwIIYKcFHIhhAhyUsiFECLIGVbIlVIhSqm1SqkFRo0phBCifUbukV8P5Bs4nhBCiA4INWIQpVQaMB24F7ipvcc7HA49cOBAIzZ9QA0NDURFRfl9O0aQrP4hWf1DsvpHe1lXr15dobVO+vW/G1LIgUeB24Do/T1AKTUHmAOQnJzMgw8+aNCm98/pdGK32/2+HSNIVv+QrP4hWf2jvaxTp04t2OcdWutu3YBTgKfavp4CLGjve3JycnRPWLJkSY9sxwiS1T8kq39IVv9oLyuwSu+jphpxjHwycJpSaicwHzhGKfWKAeMKIYTogG4Xcq31nVrrNK31QOB8YLHWema3kwkhhOgQmUcuhBBBzqiTnQBorZcCS40cUwghxIHJHrkQQgQ5KeRCCBHkpJALIQKK9mpq/7cDd0WT2VGChhRyIUTA0FpT+/F26r8spim/0uw4QUMKuRAiYNQvLcL5TSn2ySnYj0g1O07QkEIuhAgIzu92UbeogMjx/YidPhillNmRgoYUciGE6Ro3lFPz/lbCRyYQf/YwlEWKeGdIIRdCmKr5x2qq5m/GlhFDwoUjUSFSljpLfmJCCNO4iuqpfHkj1qRIHJeNxmILMTtSUJJCLoQwRWtZIxXP52Kx23DMysYSYeiF5n2KFHIhRI9z1zRTMW8DWBRJs7MJibGZHSmoSSEXQvQoj9NFxbxcvC0eHLOyCU2MMDtS0JNCLoToMd4WNxUv5OGubsFx2WhsKcGxck+gk0IuhOgR2u2l8qWNtJY6SZw5irCBsWZH6jWkkAsh/E57NVWvb6JlWy3x54wgYmSC2ZF6FSnkQgi/0lpT895WmvIqiT1lMFHj+5kdqdeRQi6E8Ku6RTtp+H430cekEy39U/xCCrkQwm/qvyqmfmkxUYf2J+b4TLPj9FpSyIUQftGwag+1C3cQMdZB3Iyh0gTLj6SQCyEM17Sxkup3txA2LI6Ec0dIEyw/k0IuhDBUy/YaKl/Lx5oaTeLMLFSolBl/k5+wEMIwrhInFS9uJDQh3NcEK0yaYPUEKeRCCEO0VjT5mmBFhOKYPYaQKKvZkfoMKeRCiG7z1LX4mmBpjWN2NqGxYWZH6lOkkAshusXb2Er5vFy8DW4cl2djTYo0O1KfI4VcCNFlXpfH1wSroonES7OwpUWbHalPkkIuhOgS7fZS+Uo+rqJ6Ei8YSfiQOLMj9VlSyIUQnaeh6q0ttGypJv7MYURkO8xO1KfJ2kpCiE7RWuPIVzQVlhN70kCiDulvdqQ+T/bIhRCdUv9FIXGFFuxHpRJ9dLrZcQRSyIUQneD8tpS6zwupS/USe9Igs+OINlLIhRAd0riujJoPtxGelUjZaC1NsAKIFHIhRLuaN1dR9eYWbINiSbxgpFSOACO/DiHEAbUU1FH5Sj7W/pE4Ls1CWaVsBJpu/0aUUuFKqZVKqR+UUnlKqb8YEUwIYb7W3Q1UPJ9HSGwYjlnZWMJlolsgMuK30gIco7V2KqWswDKl1P+01isMGFsIYRJ3VTPl83JRNguOWdmE2G1mRxL70e09cu3jbPurte2muzuuEMI8nnoXFfM2oN1ekmZnE5oQbnYkcQBK6+7XXKVUCLAaGAo8qbW+fR+PmQPMAUhOTs6ZP39+t7fbHqfTid1u9/t2jCBZ/UOydp6lFVJXWrA2QskhXlr2ceV9oGTtiN6UderUqau11gf/5g6ttWE3IA5YAmQf6HE5OTm6JyxZsqRHtmMEyeofkrVzvC633vPvdbroD1/rps1V+31cIGTtqN6UFVil91FTDT39rLWuAZYC04wcVwjhf9qjqXxtE66ddSScO5zw4fFmRxIdZMSslSSlVFzb1xHAccCm7o4rhOg52qupfmcLzflVxM0YQuRB/cyOJDrBiFkrA4AX246TW4A3tdYLDBhXCNEDtNbULtxB45oyYo7LwD4pxexIopO6Xci11uuB8QZkEUKYoP7LYpzLSrAfnkL0sRlmxxFdIJdoCdGHOVfuou6TnUSMSyL2lMHSPyVISSEXoo9q3FBBzXtbCR8RT8I5w1EWKeLBSgq5EH1Q89ZqquZvwpYRQ8JFo1AhUgqCmfz2hOhjXMX1VL6UT6gjAselWVhsIWZHEt0khVyIPqS1rJGK53Ox2K0kzc7GEmk1O5IwgBRyIfoId00LFfNyQSmSZmUTEhNmdiRhECnkQvQBnoZWKuZtwNvsxjErm1BHhNmRhIGkkAvRy3lb3FQ8n4u7ugXHpaOxpQRHAynRcVLIhejFtNtL5cv5tJY6SbxwJGGDY82OJPxACrkQvZT2aqre2EzL1hrizxpORFai2ZGEn0ghF6IX0lpT8/5WmjZUEDt9MFE5yWZHEn4khVyIXqju0wIaVu4memo60Uemmh1H+JkUciF6mfqvS6hfUkTUxP7EnJBpdhzRA6SQC9GLNKzeQ+3H24kY4yDu9KHSBKuPkEIuRC/RtLGS6ne2EDY0joTzRkgTrD5ECrkQvUDL9loqX9uENcVO4sWjUKHy0u5L5LctRJBzlTqpeDGP0PgwHJdnYwkzYuEvEUykkAsRxNwVTVQ8l4slPBTH7DGEREkTrL5ICrkQQcpT56L8uVzwahyzswmNkyZYfZUUciGCkLexlYrnNuB1tuK4PBtrv0izIwkTSSEXIsh4XR4qXtxIa3kTiZeMwpYebXYkYTIp5EIEEe3xUvVqPq7COhLOH0n40HizI4kAIIVciCChvZqqt7bQvLmauDOGEjnGYXYkESCkkAsRBLTW1C7YTtO6cmKmDcQ+cYDZkUQAkUIuRBCoX1yE89tS7EekEn10mtlxRICRQi5EgHMuL6XuswIiJ/Qj9uRB0j9F/IYUciECWOMPZdR8uI3wUQnEnzVc+qeIfZJCLkSAat5STdWbW7ANjCHxwpGoECniYt+kkAsRgFoK66h8eSPWfpE4Lh2NsoaYHUkEMCnkQgSY1j0NVDyfR0iMDcesbCzh0gRLHJgUciECiLuqmfJ5uahQi68JVrTN7EgiCEghFyJAhLRAxXO5aJeXpNnZhCaEmx1JBAn5zCZEAPA2uxmw2oKnqQXHFWOw9o8yO5IIIrJHLoTJdKuXihc3ElYPCTNHEZYZY3YkEWS6XciVUulKqSVKqXylVJ5S6nojggnRF2iPpvL1Tbh21rJnjCZiRILZkUQQMmKP3A3crLUeBUwCrlFKZRkwrhC9mtaa6nd/pHljJXGnDsGZos2OJIJUtwu51nqX1npN29f1QD6Q2t1xhejtav+3k8bVe4g+NgP74SlmxxFBTGlt3F6AUmog8BWQrbWu+9V9c4A5AMnJyTnz5883bLv743Q6sdvtft+OESSrfwRq1rjtCscWCzUZXipGaVCBm3VfJKt/tJd16tSpq7XWB//mDq21ITfADqwGzmzvsTk5ObonLFmypEe2YwTJ6h+BmNW5cpcuuv0rXfFavvZ6vP/374GYdX8kq3+0lxVYpfdRUw2ZtaKUsgLvAK9qrd81YkwheqOm3Aqq3/2RsOHxJJwjTbCEMYyYtaKAeUC+1vrh7kcSondq3lZD5eubsKVHkzhzFCpUZv8KYxjxTJoMXAwco5Ra13Y72YBxheg1XMX1VL64kVBHBI7LRmOxSRMsYZxuX9mptV4GyOdDIfajtbyRiudzsUSFkjQrG0uk1exIopeRz3ZC+JG7toWKebmglK8JVmyY2ZFELySFXAg/8TS0UjFvA94mN47Ls7E6IsyOJHopKeRC+IG3xUPFC3m4q5pxXJqFLTU45jGL4CSFXAiDabeXylc20lpST+KFowgbHGd2JNHLSSEXwkDaq6l6czMtP9YQf+ZwIrISzY4k+gAp5EIYRGtNzQdbaVpfQezJg4g6ONnsSKKPkEIuhEHqPiug4bvdRE9JI/qoNLPjiD5ECrkQBqhfVkL94iKiDulPzIkDzY4j+hgp5EJ0U8PaMmoXbCdidCJxZwzF17VCiJ4jhVyIbmjaVEX1W5sJGxJLwvkjpQmWMIUUciG6qGVnLZWv5GNNsZN4SRbKKi8nYQ555gnRBa5dDVS8kEdofJivCVZYt9sWCdFlUsiF6CR3ZRMVz23AEhaCY3Y2IXab2ZFEHyeFXIhO8NS5KJ+XCx6NY/YYQuPCzY4khBRyITrK2+Sm4rlcvE6XrwlWv0izIwkBSCEXokO8Lg8VL+bRWt5I4sVZ2NKjzY4kxP+RQi5EO7THS9Vrm3AV1JFw/gjCh8WbHUmIX5BCLsQBaK+m+u0fad5URdzpQ4kck2R2JCF+Qwq5EPuhtab24+00ri0j5sRM7IcOMDuSEPskhVyI/ahfUoTzm1Lsk1OInpJudhwh9ksKuRD74Fyxi7pPC4gc34/Y6YOlf4oIaFLIhfiVxvXl1HywlfCRCcSfPUz6p4iAJ4VciL00/1hN1RubsWXGkHjRSFSIvERE4JNnqRBtWgrrqHx5I9Z+kTguHY2yhpgdSYgOkUIuBNC6p4HKF/KwRNtwzMrGEiFNsETwkEIu+jx3dTMV83IhRJE0K5uQaGmCJYKLFHLRp3mcLirm5eJ1eXHMGkNoYoTZkYToNCnkos/yNrupeD4PT20LjsuysA2IMjuSEF0ihVz0SbrVS+VLG2nd1UDCRaMIGxhrdiQhukwKuehztEdTOX8TLdtrSThnOBEjE8yOJES3SCEXfYrWmur3fqQ5r5LYUwcTOb6f2ZGE6DYp5KJPqftkJ42r9hB9TDrRk1PNjiOEIaSQiz6j/qti6r8sJmrSAGKOzzQ7jhCGkUIu+oSGVbupXbiDiLEO4k4bIk2wRK9iSCFXSj2nlCpTSuUaMZ4QRmrKq6D6nR8JGxZHwrkjpAmW6HWM2iN/AZhm0FhCGCaiEipf34QtLZrEi7NQofIhVPQ+hjyrtdZfAVVGjCWEUVwlTgassRCaEEHiZaOx2KQJluidlNbamIGUGggs0Fpn7+f+OcAcgOTk5Jz58+cbst0DcTqd2O12v2/HCJLVWNYGSP3OgldpSg7TeMLNTtS+YPi5/kSy+kd7WadOnbpaa33wb+7QWhtyAwYCuR15bE5Oju4JS5Ys6ZHtGEGyGsdd06xL7/9Ol/x1uV62YInZcTos0H+ue5Os/tFeVmCV3kdNlQOGolfxNrZSPi8Xb5Mbx+WjaZX2KaIPkEIueg2vy0PFC3m4q5pIvCQLW1q02ZGE6BFGTT98HVgOjFBKFSulZhsxrhAdpd1eKl/Jx1VUT+IFIwkfEmd2JCF6jCHLoGitLzBiHCG6Qns1VW9upmVLNfFnDSNitMPsSEL0KDm0IoKa1pqaD7fRtL6C2JMGEXVIf7MjCdHjpJCLoFb3eSENK3ZhPzqN6KPTzI4jhClkhdkeVNPoYmuZk8KqRoqqmqhqaKHB5aHJ5aGsvJkF5T8QZQshISqM9IQIMhIiGdYvmthIq9nRA5LzmxLqvygk8uBkYqcNNDuO6CytobYIKrdC9U6oKYTmWnA1QmsjhFjBFgU2O0QPgPhMiB8ISSMhNMzs9AFFCrkf1TS6WLypjG+2VrK2qJrt5Q2/uD82wkqkLYQIWwiNjV6Kt1bQ6PJQ29T6i8cNTopiQkY8k4cmcsyIZCnsQOPaMmo+2k54ViLxZwyTJljBonIbbF4IBcuh+HtoKPv5PksohMf6irc1EjwuX1F3OX23n4TYYMBBkDYRhh0HA4/0Ff0+TAq5wZpbPSxYv4v31hazYnsVHq8mIcrGhIw4zpqQRtaAGDISI0mNiyDc+vMl40uXLmXKlCn/N0ZxdRNFVY1s3FXH2sJqFm8q4+3VxYRaFJMGJ3LG+FSmjx3wizH6iqZNVVS9tYWwwbEkXjASFSJFPKDV74Y1L0Huu1Ce7/u3hCEw9FhIzYF+o3x72tEDwLKf53NTDdQU+N4IStf63gRWzYMVT/qK//BpMOESyJwMffBNXQq5QfbUNTNv2Q7eXFVETWMrAxMjmXPUYE4c3Z+xqbFYOtFxL9wawtB+dob2szN1pG8FG69X80NxDYvy9vBJ7i5ufusH7vl4I+cenM7sIwbRLyYIrkE3QMvOWqpezcc6IIrES7JQVjnNE7CKV8G3T8CmBeB1+4rstL/DyOkQl9G5sSLifLcBB0H2mb5/a22CbYth08e+bax/w3fYZeKVMP7iPnX4RQp5N1U3uPj3l9t44duduL2aE0cnM3NSJocNTjT0477FohifEc/4jHhunzaC5dsqeXlFAc8u28GLy3dy6eEDueqoIcRH2QzbZqBp3d1AxQsbCYkNw3H5aCzh8vQNSLtzYfE9sOV/EB4Hh14FB8+CxCHGbsca4XtTGDkdXA9C3ruw8hn4+GZY9igcfTsc1DdmRssroYu8Xs0bq4q4f2E+9S1uzhifyg3HDicjMdLv21ZKcfhQB4cPdVBQ2cCjn//If7/azuvfFfLH6aM49+D0XnfM2F3ZRPm8XCw2C47Z2YTYe+8bVtBqroMv/grfPwthMXDM3b4iHtYDDatskTB+Joy7CLYvgS/+Bh9eCyueJjrtUmCK/zOYSAp5FxRUNnDr2+tZuaOKQwcl8LfTsxmebM7l4JmJUTxy3jiuOnoId3+Qy+3vbOC9tSU8cNZBPfKm0hM89S7Kn8sFjxfH78YSGt83DiMFlS2L4KMboH4XHPo7mHIHRMT3fA6lYMgxMHgq5H8I/7uDCWtuB+sOOO7Pvr34XkgOMHbSJxt2ccrjy8jfVcc/zhrD/DmTTCviexvRP5r5V07i/jPHkFdax/QnvuazjXvMjtVt3iY3Fc/l4q13kXjZaKzJ0gUrkGh3K3z2J3jtXN8x7Cs+h5P+YU4R35tSkDUDrvmO0pST4Lun4dnjfSdLeyEp5B3k8Xj56KrbWXf7nxiUFMXC647kvEMyAuoQhsWiuGBiBguvO5LMxEiufGkVD3yyCa/XmJ7zPU23eqh4MY/WskYSZ2YRlhFjdiSxl8bdO3nz2tPZvugl3zHwK5dA2m9bZZsqPIYfh/8OLnzLN2f9v1N8nx56GSnkHdDq8XLL2+tZX1TDjO3LmJdQRHpC4B62SE+I5O2rDueCiek8tXQbt7z9A60er9mxOkV7vFS+tglXQR0J540gfLjJe3jiF8Kay7G+cQ6upkY+3jOeynE3gjWAD3kNPwGu+hoSBsHrF8C618xOZCgp5O1obvUw56VVvLe2hLibbiLqyCOouO8+Gr5baXa0Awq3hnDfGWO46fjhvLumhKteXk2L22N2rA7RXk31Oz/SnF9F3IwhRI5NMjuS2FvlNsavvQNr425m3HgroRFRvPfAX2mqrzM72YHFZcBlH8PAI+D9q2H5U2YnMowU8gPweDU3zF/Hks3l3HtGNtccO4LUhx/GlplJyXXX4SosNDviASmluO7YYfzt9Gy+2FTGTW/8gCfAD7NoralduIPGNWXEHJ+JfVKK2ZHE3upK4aXTsXhdcNnHxIw7mRm33IWzqpKPHr4fj7u1/THMFBYNF70Fo06DRXfCmpfNTmQIKeT7obXmrvdz+SRvN3efksVFh2YCEBIdTfpTTwJQdPVcPE7ngYYJCBdPyuSu6aP4eMMu/vRB7k9L8wWk+qXFOJeVYJ+cQvQx6WbHEXtrqoaXz4SmataP/TMMGAtAyvCRnPC76yjauIHFz/0noJ9fgO9CobPmwZBj4aPrIH+B2Ym6TQr5frz47U5eX1nI3ClDmH3EoF/cZ8vMJPWxx3AVFFBy881oT+AfsrjiyMFcdfQQXv2ukFdWFJgdZ5+c3+2ibtFOIsf3I3b64IA6kdzneb3wzhW+Blfnv4ozeugv7s46cioTZ5zN+i8+Ye0nQVAYQ21w3suQMgHevRLK8s1O1C1SyPdhTWE19y7M57hR/bjlhBH7fEzUpEPpf9ddNHz5FWUPPtTDCbvmthNHMHVEEn9dsJF1RTVmx/mFxg3l1Ly/lfCRCcSfPQzViZYGogd8/RBs/RxO+jsMPnqfDzni/EsYcvAklr74DDt/WNPDAbvAFgXnv+rrrvjGxdBSb3aiLpNC/ivOFje/f20tyTHhPHTOuAP2SIk//zziL7qIquefp+add3swZddYLIpHzhtHv+hwrnl1DQ0tbrMjAdD8YzVV8zdjy4gh4cKRqBB5WgaUguWw5F4Ycy4cvP9VHJXFwsnX3kRiegYLHv0HVaXFPRiyi6L7w9nPQdU2WHib2Wm6TF4xv/LPTzZRWtvEY+eP71C72OQ77yDq8MPY9f/+H42rV/dAwu6Ji7Tx2PnjKKlp4sFPN5sdB1dRPZUvb8SaFIHjstFYbH2vm2NAc7f4jiPHpsMpj7TbWdAWEckZt/0JS2go7z/wV5qD4BwSg46EI2+GH16DrV+YnaZLpJDvZU1hNS+tKOCSSZnkZHZs3rIKDSX1kUewpaZSfO3vcRWX+Dll9x08MIGZkzJ44dudph5iaS1rpOL5XCx2G45ZY7BESMeIgPP1w1CxxVfEO9gzJSapH6fd/Adqy8r46NG/43EHxie/AzryFkgcBgtu9PVADzJSyNtorfnbgo30iw7j1mkjO/W9IbGxpD39FNrjoXjuXDzOhva/yWS3TxtJkj2MexZsNGWWgbummYp5G8CiSJqdTUiMNMEKOHWlsOwRyD7Lt4BDJ6SNHM3xc66lcMM6lr70rJ8CGsgaDqc+6ut5/t3TZqfpNCnkbRZvKmNtYQ03HDcce1jn9wzDBg0i9ZGHadm2jdLbbkN7A/tKyuhwK9cdO4xVBdUs3VLeo9v2OF1UzMvF2+LBMSub0MTe2cgo6H31T9BeOPZPXfr27CnHkXPKGaxbtIAfPltocDg/GHiEb4GKbx7zLWQRRKSQ49sbf+jTLQxMjOTsnK4v4GufPJnkO+7AuXgx5Y88amBC/zj34HQyEiJ56NPNPbZX7m1xU/FCHu7qFhyXjsaW0gMtTkXnVRf4VvXJudS3ek8XHXXRZQwafzCLn/8PhbnrjcvnL8fc5Vs3dPmTZifpFCnkwHc7qti4q465U4Zi7eaMifiZFxF33nlUPvMMtR9+aFBC/7CFWpg7ZQi5JXWsKqj2+/a020vlSxtpLXWSeNFIwgbF+n2boou+f9a3OPIRN3VrGIslhOnX3Ub8gFQ+evg+qneXGhTQT/qPgZGn+JaRc7eYnabDpJADLy8vIDbCymnjun85uFKK/nf9kciJE9l11900rVtnQEL/OW1cCtHhoby03L8XCWmvpur1TbRsqyX+7OFEjEr06/ZEN7Q2wdqXfSvvxKZ2e7iwyEhOv/VuUIr3H/gbLY0Bfg7pkCugsRLy3jc7SYf1+UJe4WxhUd5uzs5JM2whY2W1kvrYo4QmJ1N07e9p3bXLkHH9IdIWytk5aXySu4uqBpdftqG1pub9rTTlVRJ7ymCiJiT7ZTvCIPkLfJfjH7L/OeOdFdd/AKfddCc1u0tZ8NgDeL0BfDX0oKMhcSisft7sJB3W5wv5F/l7cHs1Z4zv/p7H3kLj40l/+il0czNFc6/B2xi4U5rOGJ9Kq0fzRb5/FqKoW1RAw8rdRE9NJ/oIY3/Owg82vu9b0X7gUYYOmz56LMfOupqd61bz1SvPGTq2oSwW38VPhSugPjgWZ+nzhXxR3h5S4yIYnWL8ogVhQ4eS+vBDtGzeTOkddwbsTJYxqbGkxIazKM/4J239V8XULy0i6tD+xJyQafj4wmCuRt9FMSOn+wqawcYeN43xJ53K6o8/YMPiTw0f3zCjTgE0bA6C2Tb08ULe3Oph2dYKThid7LcGTfajjqLfrbdS/+mnVPzrX37ZRncppThhdH++/rHc0J7lDav3ULtwBxFjHMTNGCpNsILBzq/B3eQr5H4y5eIryBw7ns+ffYrijbl+20639MuC+EGw5ROzk3RIny7kuSW1uNxeDhvs3xNvCZddSuxZZ1Lx1NPULQzMd/hJgxNocXvJKzVmcYCmjZVUv7OFsGFxJJw3QppgBYvCFWAJhfRJftuEJSSEU264ndh+yXz48H3Ulu3227a6TCkYOBmKVvpm7wS4Pl3I1xT6ptyNz/DvMmJKKfr/+c9E5ORQeucfaNoQeHshP/0M1hgwDbFlew2Vr+VjTY0mcWYWKrRPP82CS/H3kJwNNv8uZRgeZef02/6E1+vh/Qf+hqspAM8hpU2EpqqgWLC5T7/C1hfXkhoXQVJ0mN+3ZbHZSHvicUITEym+5hpa9wTWSZTkmHBSYsPZUFLbrXFcJU4qXtxIaEK4rwlWmDTBChpaQ+naHltAOSEllVNvuJPKkiI+fuLBwJvJknaI78/Stebm6IA+XcgLqxoZnBTVY9sLTUgg7emn8TqdFF9zLd7m5h7bdkcMTrJTUNn1PaPWiiZfE6yIUByzxxAS1X73SBFAGirA5fQ1j+ohmWPHMfWyOWxfvZJl8wNs2bWEwb4/q3eYm6MD+nwhz0jw70fIXwsfMZyUB/9Jc14eu/7wx4BaFis9IZKiqq4Vck9di68JltY4ZmcTGuv/TznCYNU7fX9245L8rhh/4ikcdPzJfP/B2+R9GUBtZK3hEJ3y888lgBlSyJVS05RSm5VSW5VSdxgxpr81utzUNLaSGt/zDZuijzmGpJtupG7hQir//e8e3/7+pMVHUNngorm1cx9xvY2tlM/LxdvgxnF5Ntaknn1zFAapLfL9Gdfza6VOvWwO6aPH8tl/n6B0SwAtuxaXDjWBvcg6GFDIlVIhwJPASUAWcIFSKqu74/qbs211nOhwcz7+J15xBbEzTqP8scepWxQY82l/6vrYmZWDvC6PrwlWRROJl2RhS4v2Vzzhb662RSDCjL+moj0hoaGceuMdRCcm8cGD91JXUdbjGU47rioAABUgSURBVPYpLPrnn0sAM6KT/0Rgq9Z6O4BSaj4wA9howNh+0+TyEJb8EW8UvcriGnMOA4Qc4eX8dXZct97Im5dn8GLzi6bk+ElpTRNhyXFUOI8i0d7+z0S7vVS+ko+rqJ7Ei0YRPjSuB1IKf/n66ygqKv8Gz5SC5cCtjWtqvFSvNn5dTnu/synZ+F9evPUPpIyajSWk+6/NbmUtOQkHeRzZ7RT+pbp7jFYpdTYwTWt9RdvfLwYO1Vpf+6vHzQHmACQnJ+fMnz+/W9vtCKfTid2+7zapxfVe/rb1TdITdxNtM2+Os73ezdxndgLw1JUDcUabt0pOidNLbf0AZvU/m0MH7D+H0+nEHmUneb0iepeFstFe6tID51j/3g70HAg0Zmd1frmF1opmauOy0O18WPd4PISE+GdGkqthB86S97Hah2AfcGq3LyTrTta4mlwc1h14TpvRrQwd1d5zYOrUqau11r+dVqS17tYNOAd4dq+/Xww8caDvycnJ0T1hyZIl+71vZ4VTZ96+QL+9qqhHshxIU16ezhszVu849zztaW42LcfTS7fqzNsX6M276w74uCWLl+iq93/URbd/pWuXFPZQuq450HMg0JiedfnTWv85RuuGynYf6u+sqxa8rx88d7peNv+lbo/Vraz/OtT3M+kh7WUFVul91FQjTnYWA3ufHUkDArzpMES0LfLb2MkTe/4QnpVF7eWX0fTDD+z+059Mm8kS0db9MTHqwMuuxW9TNCzfhf2oVKKP7vpCHCLA/HQRkMv8NrMTTj6N7KknsOLdN8j/5kvzgsSmQsp487bfQUYU8u+BYUqpQUopG3A+ENgrKgDRYb6TnHVNrSYn8WkZP56k66+j9oMPqXzWnDUOfzoBHHWApe6cy0tJ3GohMieZ2JMGSf+U3iSs7UR1c/cuCjOCUorjrria1JGj+fTpx9i1dbM5QVrqwRb4h+a6Xci11m7gWmARkA+8qbXO6+64/hZhCyE+0kpJTZPZUf5P4lVXEXPyyZQ//Aj1ixf3+PaLqxtx2G377cveuK6Mmg+34eyniT9zmBTx3ia27YP1T9MQTRYSauW0m/9AZFw8Hzx4L/VVFT0foqYQ4jJ6frudZMg8cq31Qq31cK31EK31vUaM2RO6cwGMPyilGHDfvYSPHk3pLbfSvHlLj26/sKqRtPh9zwFv3lxF1ZtbsA2MZc9BXlSIFPFeJ66tzXC1f1eL6ozImFjOuO1uXE1NfPDPe2ht6cGroVuboX5Xj18g1RV9+srOzMQotpebfzxwb5bwcNKefBKL3U7x1VfjrqrqsW1vL28gM/G3hbyloI7KV/Kx9o/EcWkWWtqn9E5RDrBFQ0XP7kC0x5ExkOnX3cqeHdv45OnHeu4cUlVbs6z4QT2zvW7o04V8bGosJTVNlNUHVs8Ta3I/0p58EndlJcW/vw6vyz9LsO1td20zu2qbGZv2y7ngrbsbqHghj5DYMByzsrGEmzc9UviZUpAyDkpWmZ3kN4bkTOSoCy9jy/KvWfGO/6cuA75OkACpE3pme93Qpwv5hExf0VpbWGNykt+KGJNNyv330bR6Nbv/31/8vheytq2l74SMnwu5u6qZ8nm5KKsFx6xsQuwHns0ieoH0ibA7NyBmrvzawaeeSdZRx/DtW6+yZcUy/2+w6HuITPy5eVYA69OFfHRKLLZQC8u3VZodZZ9iTj4Zx9yrqX33Xape8O9Vnyu2VxIWaiGrbck7T72Linkb0G4vSbOyCU0I9+v2RYBInwTa41tgIsAopTj+ymsZMHwk/3vyEfZs3+q/jWntWy0pbaLvk0qA69OFPNwawpFDHXy2cU9AdSHcm+Paa4k+4QTK/vlPnF/6Zz6t1ppPN+7hqOFJhIWG4G12U/FcLp46F47LRmPt33OtfoXJBh0J1kjY9LHZSfYp1GZjxs1/JCI6hvcfvAdntZ/OIe3Jg5oCGDHNP+MbrE8XcoATR/enpKbJsCXOjKYsFlL+fj9hI0dQcvMttGw1fi9kfXEtu2qbOXF0f3Srh4oX82gtayTx4izCMnu+gZIwkTUChh7rK+QBulh4VFw8p992N83Oej588F7c/jiHtGkBoGDEycaP7Qd9vpAfO6ofoRbFu2tKzI6yX5bISNKffBIVHk7R3GtwV3d/Oba9vbe2BGuI4pjhSVS+tgnXzjoSzh1O+HD/LoEnAlTW6eDcDTuWmp1kv/oNHMzJ197Mrq2b+fQ/jxv7idrrhfVvQObhYO9n3Lh+1OcLeaI9jGnZ/Xl7dRFNLvMv198f64ABpP/rCdy7d1Ny/Q3oVmOuSG10uXlndTEnZ/eHTwpozq8i7rQhRB4UHE9g4QcjT/Gd5Pt+ntlJDmjYxMOZfN7F5C9bysoP3jZu4O1LoGo75Fxu3Jh+1ucLOcDFkzKpa3bzwbrA3SsHiBg3jgH3/I3GlSvZfc+9huyFvL+2lPoWN3MJp3H1HmKOy8B+WIoBaUXQsobD+JmweSHUBMZVnvtz6BnnMnLy0Syb/xJbvzfoBO33z0KkA7JOM2a8HiCFHJg4KIHRKTE8tXQbrZ7APC74k9jTTiPxyiupeeMNql99rVtjudxenlq6lZtjYoj6oRL74SlEHxv4lyOLHnDIFaBCYNnDZic5IKUUJ1x1Hf0HD2XhEw9SXtDN9TV3/eB7AztkNoQGz3KFUsjxPRluOWEEhVWNvLkqsPdAAJJuvAH7scey5/77cX7zTZfHeeP7QsZXuzmjDiLGJRF7ymDpnyJ84jIg51JY8xJUBfbiw1ZbGDNuvZuwqCjee+CvNNZ247qQxfdAeBwcdo1xAXuAFPI2U0YkkZMZz2Of/0h9c2B0RNwfZbGQ8o9/EDZkCCU33EjL9s6/0GqbWlm5aDu3EkHY8HgSzh6OskgRF3s58hawhMIXfzE7Sbvs8QmcfuvdNNXV8cFD9+HuyjmkHV/Bj5/CETdAeKzxIf1ICnkbpRR3n5JFubOFBz4xqWVmJ4TYo0h/+imU1Urx3Ll4ajvXevSVNzZwY3Monv6RJM4chQqVp4L4lZgBcMRNkPcebAmMdWUPJHnwUKbNvZHSzRv5/JknO3cOqbUJPrrB11dl4u/8F9JP5NW7l3HpcVx2+EBe+a6AVTt7rllVV1lTU0l74nFcJSWU3Hgj2t2xRZPXrizh+E1OGiJDyZwzFotNumCJ/TjiRkgaCQtu9PXmDnAjDjuCw86+gLwvP2f1gvc6/o1fPuBrknXqoz8vsBFEpJD/yi0njCA1LoLr56+jusH/zaq6KzInhwF/+QsN3y5nz9//0e7jKwpqCXtvG04LDJo7DkuktQdSiqAVaoNTH4f6Uvjoet+l6wHusLMuYPikI/jy1efZvvb79r9h+1L45lEYNxMGT/FzOv+QQv4rUWGhPHnhBMrqm7npzXV4vYH/xI078wwSLr+c6ldeoXr+G/t9XGtVM7ufWY9ba8IuHIndEXx7HsIEGYfCMXdB7juw8hmz07RLWSxMm3sD/QYO5uPHHqDpQAtS1JXC27PBMRxOan9HKFBJId+Hg9Lj+NMpWSzZXM4/Fm0yO06H9LvlZqKOPord99xDw4rvfnO/p6GVzU+sIdSt2T41hdHZcsGP6ITJN8LwabDoD7Ct51ev6ixrWDin33o31rBwtv7vPRrr9nEOqaUe5l/kOz5+7ksQFvhLuu2PFPL9mDkpk5mTMvjPl9t55qvtZsdplwoJIfWhh7ANzKTk+utxFfy8you3xU3+46uJaHKzeHQ0M04YamJSEZQsFjjjP5A0AubPJLruR7MTtSs60cGMW+6itcHJR4/cj8e910wWdwu8MdM3b/zseb7/VxCTQr4fSin+clo208cO4N6F+by8fKfZkdoVYreT/tRTABTNvQZPfT3a7WXDE2uw17pYMDCCqy4aJ3PFRddExMHMdyDKwdj1f4HStWYnateAYSMYOGUaxRtz+eK5f/tmsrQ2w9uzfMfGZ/wLRpxkdsxuk0J+ACEWxcPnHsRxo/px9wd5PPr5loBtd/sTW0YGqY8/jquggJKbbmHFY9+TWNHCRyk2rr4iB4vMFRfdEd0fLnkfT0gEvHCqb+51gEsYPopDzziXDV8sYu1Hb8GrZ/u6G570Txh3odnxDCGFvB1hoSH8e2YOZ01I49HPf+SOdzbQ3Bq4zbUAog6dSMKdf8Rdl056uYuv0sK4+pqJ2GSuuDBCwmDWTPg7xKbBK2fB2lfNTtSuyefOZOhBY1n66ovs3LgRznwGDp1jdizDyCu7A0JDLDx4zliunTqUN1YVcdbT31JY2Wh2rP0qqGxgVnkq2+xJ7HFt5/y5BxMaIr9qYRxXWCJcvhAyJsEHc+GDa30nDQOU2vIJJ+nXSIlqxHvUnTD2XLMjGUpe3R2klOKWE0fw7CUHU1TVyPTHv+aVFQUBNT3R69W8vHwn0x9fRnFtEyE3n8aEhy7BYpFfs/CDyAS4+H3fpfxrX4b/HA2Fv50xZaqmGoZvfgrmX4DNkcF5Dz7H4FOuNDuV4eQV3knHZSXz8XVHMiYtlrvez+W8/y4nf5f5qwttLK3jnP8s5+4P8hiXHseC3x/BMVkD5MSm8C9LCBx7t+8kaGsjPHcifHwzNJp8ZbTX65v3/uShDNj1GRx2Lcz6FJUwyNxcfiKFvAvSEyJ59YpD+efZY/mxzMnJj3/Nda+vZUdFz688vq3cybWvreHkx79mW7mTh845iJdnTyQ9QS72ET1o6HEwdwVMuhpWPQePjoWlf4fmHt7J0drXF+a/R/tmptj7sTrnn3Divb4+671UqNkBgpVSinMOTueErP789+ttPLdsJwvWl3LsqGQuOSyTyUMcfpsh4vVqlm2t4OUVBXyRv4dwawjXTh3KlUcNJjZCLrkXJgmzw7T7YcKlsOReWHo/LH8Kxl3g62/uGOa/bbsaYMNbvkUhdm+AuEzfvPcx5+D86mv/bTdASCHvpthIK7eeOJLLDh/E89/s4I3vi/hs4x7SEyKYNro/J47uz4SM+G4XdY9Xs6awmk/zdvNJ3m6KqppIjLJx9ZQhXD55EA578DTBF71cv5Fw3su+eebLn/QtGffdvyH9UN8ycqNOgYTB3d9OixO2fQH5C2DLJ9BSB8nZvt4wB13g6xPTR0ghN0hSdBi3TRvJ9ccN438bdvPe2hJe+HYnz3y9g9gIK+PS45iQEU9WSgwZCZGkJ0QQadv3j7/R5aaoqonCqkbySmtZW1jD2sJq6prd2EIsHD40kVtOGMG07P6EhUrnQhGgUsbDWc/Ciff5TobmvQef3e27xWVA2kRIO8R3VWX8QN90xpB9fKLU2nfMvXqnby3N0jVQtBJ2rwePCyISYNRpMOESSJ8IffC8kBRyg4WFhnD6+FROH59KXXMrSzaVsXxbJWsLa3j0iy2/aB5nDwslwhZCpC2EpqYm9Def0+Ty4Gz5uR2tUjC8XzTTxw7gsCEOpo5IIjpcDp+IIGLvB0fe7LtVF8Dm/0Hht1DwLeTutWiyCvEdnrFGgTUCPK3Q2uA7bOJu/vlxoRG+N4lJV8PQ4yHjMAjp26Wsb//v/Swm3MqMcanMGJcKQF1zK9vLGyisaqSwsoGqhlaaWt00tHgoK2thUHo/IqyhJNptpCdEkpEQyZCkKCncoveIz4RJV/lu4Os+WLnNt7ddUwDNteBq9M2ACbH5eoNbIyEm1bfXHj/Qd6x9X3vufZgU8h4UE+47xDIuPe439y1dupQpU8aakEoIE8Wk+G6DjjQ7SVCT6YdCCBHkpJALIUSQk0IuhBBBrluFXCl1jlIqTynlVUodbFQoIYQQHdfdPfJc4Ewg8JsSCyFEL9WtWSta63xAGjMJIYSJ5Bi5EEIEOdXe0mVKqc+B/vu4649a6w/aHrMUuEVrveoA48wB5gAkJyfnzJ8/v6uZO8zpdGK3B8fK2JLVPySrf0hW/2gv69SpU1drrX9zPrLdQt4RHSnkv3p8OVDQ7gO7zwFU9MB2jCBZ/UOy+odk9Y/2smZqrZN+/Y+mXNm5ryD+oJRata93r0AkWf1DsvqHZPWPrmbt7vTDM5RSxcBhwMdKqUXdGU8IIUTndXfWynvAewZlEUII0QW9fdbKf80O0AmS1T8kq39IVv/oUlZDTnYKIYQwT2/fIxdCiF5PCrkQQgS5Xl/IA72xl1JqmlJqs1Jqq1LqDrPzHIhS6jmlVJlSKtfsLO1RSqUrpZYopfLbfv/Xm51pf5RS4UqplUqpH9qy/sXsTAeilApRSq1VSi0wO0t7lFI7lVIblFLrlFIdus7FLEqpOKXU20qpTW3P28M6+r29vpATwI29lFIhwJPASUAWcIFSKsvcVAf0AjDN7BAd5AZu1lqPAiYB1wTwz7YFOEZrfRAwDpimlJpkcqYDuR7INztEJ0zVWo8LgrnkjwGfaK1HAgfRiZ9xry/kWut8rfVms3Psx0Rgq9Z6u9baBcwHZpicab+01l8BVWbn6Ait9S6t9Zq2r+vxvShSzU21b9rH2fZXa9stIGchKKXSgOnAs2Zn6U2UUjHAUcA8AK21S2td09Hv7/WFPMClAkV7/b2YAC02wUwpNRAYD3xnbpL9aztcsQ4oAz7TWgdq1keB2wCv2UE6SAOfKqVWt/V7ClSDgXLg+bbDVs8qpaI6+s29opArpT5XSuXu4xawe7dt9tX/NyD3xIKVUsoOvAPcoLWuMzvP/mitPVrrcUAaMFEplW12pl9TSp0ClGmtV5udpRMma60n4Dt8eY1S6iizA+1HKDABeFprPR5oADp8zsyUXitG01ofZ3aGLioG0vf6expQalKWXkcpZcVXxF/VWr9rdp6O0FrXtDWhm4bv/E4gmQycppQ6GQgHYpRSr2itZ5qca7+01qVtf5Yppd7Ddzgz4M6X4asFxXt9EnubThTyXrFHHsS+B4YppQYppWzA+cCHJmfqFZRvtZN5QL7W+mGz8xyIUipJKRXX9nUEcBywydxUv6W1vlNrnaa1Hojvubo4kIu4UipKKRX909fACQTemyMAWuvdQJFSakTbPx0LbOzo9/f6Qh7Ijb201m7gWmARvpNxb2qt88xNtX9KqdeB5cAIpVSxUmq22ZkOYDJwMXBM29SzdW17koFoALBEKbUe35v7Z1rrgJ/aFwSSgWVKqR+AlcDHWutPTM50IL8HXm17HowD7uvoN8ol+kIIEeR6/R65EEL0dlLIhRAiyEkhF0KIICeFXAghgpwUciGECHJSyIUQIshJIRdCiCD3/wE7SewBmbh9nwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "L = 1.6\n",
    "GCR = 0.35\n",
    "P = L/GCR\n",
    "print(f'pitch = {P}')\n",
    "\n",
    "# tracker 1\n",
    "pts1 = np.radians(np.arange(360))\n",
    "pts1 = np.stack((L/2*np.cos(pts1), L/2*np.sin(pts1)), axis=1)\n",
    "circle1 = LinearRing(pts1)\n",
    "plt.plot(*circle1.xy)\n",
    "\n",
    "# tracker 2\n",
    "pts2 = np.radians(np.arange(360))\n",
    "pts2 = np.stack((P + L/2*np.cos(pts2), L/2*np.sin(pts2)), axis=1)\n",
    "circle2 = LinearRing(pts2)\n",
    "plt.plot(*circle2.xy)\n",
    "\n",
    "# tracker 1 surface\n",
    "tracker1 = LineString([(-L/2, 0), (L/2, 0)])\n",
    "plt.plot(*tracker1.xy)\n",
    "tracker1rot = affinity.rotate(\n",
    "    tracker1, tracker_theta, use_radians=True)\n",
    "plt.plot(*tracker1rot.xy)\n",
    "\n",
    "# tracker 2 surface\n",
    "tracker2 = LineString([(P-L/2, 0), (P+L/2, 0)])\n",
    "plt.plot(*tracker2.xy)\n",
    "center2 = shapely.geometry.Point((P, 0))\n",
    "tracker2rot = affinity.rotate(\n",
    "    tracker2, angle=tracker_theta, use_radians=True, origin=center2)\n",
    "plt.plot(*tracker2rot.xy)\n",
    "\n",
    "# sunray\n",
    "a, b = tracker1rot.coords\n",
    "c = P * np.tan(tracker_theta-np.pi/2)\n",
    "sunray = LineString([a, (P, c)])\n",
    "plt.plot(*sunray.xy)\n",
    "\n",
    "plt.gca().axis('equal')\n",
    "plt.grid()\n",
    "list(sunray.coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'tracker_theta': array([129.68948615]),\n",
       " 'aoi': array([61.35300178]),\n",
       " 'surface_azimuth': array([292.53615425]),\n",
       " 'surface_tilt': array([56.42233095])}"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result_noback180"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However if the trackers are bifacial then it may be advantageous to backtrack\n",
    "so that the back of the panel is facing the sun.\n",
    "\n",
    "In this particular situation the trackers are tilted 30-degrees. At a (ze, az)\n",
    "of (80, 338) the solar vector is `[[-0.37], [0.9131], [0.17365]]` which means\n",
    "that the sun in the y-z plane is only 10.77-degrees above the horizon, but the\n",
    "tracker is tilted 30-degrees, so **the sun is coming from below the trackers**.\n",
    "\n",
    "This means that there's no chance of shading the bottom of the next row, but\n",
    "it might shade the top. The backtrack condition is different because if the\n",
    "tracker rotation R > 90, then cos(R) < 0, and one tracker will shade the top\n",
    "of the next if\n",
    "\n",
    "    LR = -L/cos(R) > x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.63862662])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.cos(np.pi-tracker_theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'tracker_theta': array([129.68948615]),\n",
       " 'aoi': array([61.35300178]),\n",
       " 'surface_azimuth': array([292.53615425]),\n",
       " 'surface_tilt': array([56.42233095])}"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result_noback180"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(result_noback180)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

gh656.py

import pvlib
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# in Brazil so facing north
axis_azimuth = 0.0 
axis_tilt = 20
max_angle = 75.0
gcr = 0.35

# Brazil, timezone is UTC-3[hrs]
starttime = '2017-01-01T00:30:00-0300'
stoptime = '2017-12-31T23:59:59-0300'
lat, lon = -27.597300, -48.549610
times = pd.DatetimeIndex(pd.date_range(
    starttime, stoptime, freq='5T'))
solpos = pvlib.solarposition.get_solarposition(
    times, lat, lon)

# get the early times
ts0 = '2017-01-01 05:30:00-03:00'
ts1 = '2017-01-01 12:30:00-03:00'
apparent_zenith = solpos['apparent_zenith'][ts0:ts1]
azimuth = solpos['azimuth'][ts0:ts1]

# current implementation
sat = pvlib.tracking.singleaxis(
    apparent_zenith, azimuth, axis_tilt, axis_azimuth, max_angle, True, gcr)

# turn off backtracking and set max angle to 180[deg]
sat180no = pvlib.tracking.singleaxis(
    apparent_zenith, azimuth, axis_tilt, axis_azimuth, max_angle=180, gcr=gcr, backtrack=False)

# calculate cos(R)
# cos(R) = L / Lx, R is rotation, L is surface length,
# Lx is shadow on ground, tracker shades when Lx > x
# x is row spacing related to GCR, x = L/GCR
lrot = np.cos(np.radians(sat180no.tracker_theta))

# proposed backtracking algorithm for sun below trackers
# Note: if tr_rot > 90[deg] then lrot < 0 
# which *can* happen at low angles if axis tilt > 0
# tracker should never backtrack more than 90[deg], when lrot = 0
# if sun below trackers then use abs() to reverse direction of trackers
cos_rot = np.minimum(np.abs(lrot) / gcr, 1.0)
backtrack_rot = np.degrees(np.arccos(cos_rot))

# combine backtracking correction with the true-tracked rotation
# Note: arccosine always positive between [-90, 90] so change
# sign of backtrack correction depending on which way tracker is rotating
tracker_wbacktrack = sat180no.tracker_theta - np.sign(sat180no.tracker_theta) * backtrack_rot

# plot figure
df = pd.DataFrame({
'sat': sat.tracker_theta,
'sat180no': sat180no.tracker_theta,
'lrot': lrot,
'cos_rot': cos_rot,
'backtrack_rot': backtrack_rot,
'tracker_wbacktrack': tracker_wbacktrack})
plt.ion()
df[['sat', 'sat180no', 'tracker_wbacktrack']].iloc[:25].plot()
plt.title('proposed backtracking for sun below tracker')
plt.ylabel('tracker rotation [degrees]')
plt.yticks(np.arange(-30,200,15))
plt.grid()


from shapely.geometry.polygon import LinearRing
from shapely import affinity
from shapely.geometry import LineString

L = 1.6  # length of trackers
P = L/gcr  # distance between rows

f = plt.figure('trackers')  # new figure

# true track position at 5:30AM
tracker_theta = -np.radians(df.sat180no.values[0])

# tracker 1 circle
pts1 = np.radians(np.arange(360))
pts1 = np.stack((L/2*np.cos(pts1), L/2*np.sin(pts1)), axis=1)
circle1 = LinearRing(pts1)
plt.plot(*circle1.xy, ':')

# tracker 2 circle
pts2 = np.radians(np.arange(360))
pts2 = np.stack((P + L/2*np.cos(pts2), L/2*np.sin(pts2)), axis=1)
circle2 = LinearRing(pts2)
plt.plot(*circle2.xy, ':')

# tracker 1 surface
tracker1 = LineString([(-L/2, 0), (L/2, 0)])
plt.plot(*tracker1.xy, '-.')
tracker1rot = affinity.rotate(
    tracker1, tracker_theta, use_radians=True)
plt.plot(*tracker1rot.xy)

# tracker 2 surface
tracker2 = LineString([(P-L/2, 0), (P+L/2, 0)])
plt.plot(*tracker2.xy, '-.')
center2 = shapely.geometry.Point((P, 0))
tracker2rot = affinity.rotate(
    tracker2, angle=tracker_theta, use_radians=True, origin=center2)
plt.plot(*tracker2rot.xy)

# sunray
a, b = tracker2rot.coords
d0 = b[0] - P
d1 = b[1] - P * np.tan(tracker_theta-np.pi/2)
sunray2 = LineString([b, (d0, d1)])
plt.plot(*sunray2.xy, '--')

# backtracking
tracker_theta = -np.radians(df.tracker_wbacktrack.values[0])

# backtrack tracker 1 surface
tracker1 = LineString([(-L/2, 0), (L/2, 0)])
tracker1rot = affinity.rotate(
    tracker1, tracker_theta, use_radians=True)
plt.plot(*tracker1rot.xy)

# tracker 2 surface
tracker2 = LineString([(P-L/2, 0), (P+L/2, 0)])
center2 = shapely.geometry.Point((P, 0))
tracker2rot = affinity.rotate(
    tracker2, angle=tracker_theta, use_radians=True, origin=center2)
plt.plot(*tracker2rot.xy)

# parallel sunrays
sun_angle1 = np.arctan2(*reversed(np.diff(sunray1.xy)))
# sun_angle2 = np.arctan2(*reversed(np.diff(sunray2.xy)))
a, b = tracker1rot.coords
c0 = a[0] + P + L
c1 = a[1] + (P+L) * np.tan(sun_angle1)
sunray1 = LineString([a, (c0, c1)])
plt.plot(*sunray1.xy, '--')

# alternate backtracking
tracker_theta = -np.radians(df.sat.values[0])

# backtrack tracker 1 surface
tracker1 = LineString([(-L/2, 0), (L/2, 0)])
tracker1rot = affinity.rotate(
    tracker1, tracker_theta, use_radians=True)
plt.plot(*tracker1rot.xy)

# tracker 2 surface
tracker2 = LineString([(P-L/2, 0), (P+L/2, 0)])
center2 = shapely.geometry.Point((P, 0))
tracker2rot = affinity.rotate(
    tracker2, angle=tracker_theta, use_radians=True, origin=center2)
plt.plot(*tracker2rot.xy)

plt.gca().axis('equal')
plt.ylim([-2,6])
plt.xlim([-2,6])
plt.grid()
plt.title('Backtracking with sun below trackers')
plt.xlabel('distance between rows')
plt.ylabel('height above "system" plane')
plt.legend([
    'tracker 1',
    'tracker 2',
    'tracker 1: system plane',
    'tracker 1: true track 98.3[deg]',
    'tracker 2: system plane',
    'tracker 2: true track 98.3[deg]',
    'sunray',
    'tracker 1: backtrack 32.5[deg]',
    'tracker 2: backtrack 32.5[deg]',
    'parallel sunray',
    'tracker 1: alt backtrack -16[deg] or 164[deg]',
    'tracker 2: alt backtrack -16[deg] or 164[deg]'])