mikofski/vfsky.ipynb

## vfsky.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e93c6f41-c341-4d90-b0c8-e9ba9469358f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "attachments": {
    "8a27b5e2-3f7d-4258-b2e6-b4c253987a29.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAABbCAYAAABd59MPAAAHc0lEQVR4nO3daY/kKAyAYWq0//8v13wYRU2zHAZsYuB9pNX2dFeOymFzhYQAAAAAAAAAAABwns/bOwDgVd/k38QEAEDTN/n5G/6fUAAAG7IM6rl1WicQEtRh/ry9AwBesbqp6kkcs9slCTny39s7ADhVC1KfxmfoR9BVOs6cIwCuSUq8p5SKrZuvtNaTruumc+QKTViAjt1LtLskjxDGj/Xu58gdEghQ12oK+YZzAtOu3+Omc+QKCQTAigBLAD8QnejAuNHAG5eU49Jzbl21Zpvc33Lrbu1L+jmthJI+ZyLdp9q6epen9mGIBAKMmU0ekmaX1jokv6vt4zf5v5b4O6Y/rwzoJA9jNGEB/bSSR/pz+rnaOkrrzdUmpOup7ZNEPNKp9h2lI6ZmEhvJYwFqIEAfi8AUl9Sff4eQD6C5v6X7k66vtE1Nkv1v7ZdGk9yzDMljAWogQJtlsD3NTO0lV4OSro+E8QISCCCnWbJNm2h2DoCzzWGaCZXax0IkEKBNax6n3v6O3ZSaniS/05gjS2M96EACAdYqJZFdE8lsR3cIspFiJAaHSCDAerm2/hD2TSIlpaDf0xkOx0ggQJ12+3yslEh21Nt8pVmzING8hAQCyFgG+lOSiMRI7YORWE6RQIC1bikt52okPMdxGBIIUGYV7C2bc95Smhal9d1y72Xv6R+6JSG7tPMFC1jTfk4jDo65/pDcdnu89Ra+kafeJfsoOf6nPEuzJc0DfkIpClg9Bchpepqfakk0/nvvXFqco0W0Ht7JrXOmaskFAOyHQuRlZidTrFW7JZO+xWjLBICNaHWiS0oerVIJpRZgf9zHF/E4CosLENgPLQgX8phAAACXqE0EJ5kkjpILAGyIGggAYEjvKCzLB6sAABuRJhCtdxUDAA4hfd1kaxqBmXcb3Jx8Vj1sST8TgBHVONOqgWgEnhtenDPq5uQJYHO1TnTLaQkInABwKOnwW4bvAsClRofxEvwBwD/TWJ3rA1k1oybNWADQpzUF/tL11Tq4R1/gop2A3pwi2nrbTHkPQKr0QrJarG7NgN47a/ovIw8SfsLaJqza9qwDfOu7zpYGSAIAJNJkMBuHVdYnGYX1/Nz7rmJLkk76FfvwCTqvIj0VxwOoe+t5OZX1tTrR40BdSh7WATSXuNLAbaWWNNPsjR+1AgeAH6M1CUlXg/n9l0sguY2WdiRuP0sPgodaghUCYxmvNQX6PElEGi8t42rXfVvqA+l9e6B2sDi54/xkHDtgTJxEJPH3ude0ahpHFPbTLDzzMKP1tnPL3+zkGiewSm9NJPfZmSaxLr2jsCxZtNm1DkrcBDez7dsDJzWPH9rj9HGXnpqIluFteUogFlbevLcGituTZyw3rp7BBPAsd32Kr1kvbyQsfYkdtk2A+Of2Y8DIPGjJPdxXamIfDv7hdyEn/k/MSwIpefMmHH1Xyk0l8pu+awsJAxZqo1utmt3NnkRfJc20K29O6ba/yf8BwIo0Bmp/zn4lStK5WFoPMOaMfp+ebde2P7MPb5g5jj1z59zaye69eVP7PkLe6PU/ND/VSi53CktIR6hJlpcmea/Xm6RA0BtsvSdNzfNfWpYENXcdkEDgUq0DLiatgUn7gbxfb5KbXRoQPNc+Zs+/ZF29nzmVVgIZXYcp753o0PcENsmUNaP9O59Q3sYJdk8e1uc/x+vxWOHY704CQermKelbTS+SxLD7yLzR8z9zzE41e97dHzevo7Bgx/1FubEdRuZx/qGGBIJeXoLjG8O7Je38ObcGbe0pgnpHB0qXy3VWa11fmufeXW2OBIIaVxfri6QzpJ5mNvjPDqv/JL+TPJdVWi79fWlbGgUkyRPjR6APBKhzV+rbgHbyyP29JF0uTSatbc0mzdKUNl5q7qpIIEgdXWLqdGPymD3/GsdMMkIs3WbrM7PbrynVcnYeTCFCAgF+uy1haNA6ZpLaxci2arUQyXZrSv0npe0fhQSCGLWPH9Q+xpZffcy0SvXpTLQ3D2cXI4HgQfL459bjoPHE9G7HrGeevZzrR94xCgsxbxd962Y+8QaeKVFb9D2cLp489dGbEC2byFwjgSAEv801rY5Tj/s8643v5O1YSp650dzfUt+FdNjwtWjCgrfg8aYbg4L2w37WNKcHyXV2Wz48eNx9RgK5m+RmvGZESeSG7xjC2Pkv0Q680uZJ6TKj+9HzWevmK3fXJU1Y9+qZhhvn2eH8a46wCqHcab4iMLsL/hpIIHfKTbVQsvPU5T3eDpYraZ1/i2OW69RO/15bJr0epbWs3DJaz5zMPpXvFgnkPisuytrU3o8Tks6O3AelUE4itWsm7gjvaXadffjvST6zI7m2dPSXg6m3k4FGzWe29nUjjtkP62Px9j3W5HKnsAX3FzewkEVTrvt7jFFYGGUx+gXYkfX17zJ5hEACAQCPtiiUkUAww23JCFjM6l5wfY+RQDDr6BfmAA0W1/029xIJBAB8cl37ADTlxt8D6MN9hGtx8QPjuH9wPW4CoB/3DRDhZgBkSB5ABjcFUMc9AgAAAAAAAAAA4MxfVZP8Zx7BsawAAAAASUVORK5CYII="
    },
    "a2e1b105-bbbd-4836-8851-25cd9305ea37.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUcAAADlCAYAAADA6jU6AAAKGUlEQVR4nO3d63LbNhAGUCbT939l9UerhqV5J0DsAufMZJrUtkQT4KfFRdQ0AQAAAAAAAAAAAAAAAAAAAAAAAEBKn9YHwLh+tz4AgIiEI9F8pp8V49G/obi/Wh8AzHymafo1+zs0o3IkinkwTou/AwxrrVL8rPxdRckrVI4AENyyUrQQA/CveSj2Eo5rQU9wJr3JYrlgk8EyELMd/9Bs5YHy1qpEwZiMBiODLFXj0dA5w+/AvzQWPHdUKX42/j+BGVbDfWeGzxZikhKOcN3dOUVVYyLCEc6z0DIQ4QjHSoSiEE1Gg5HFN6De7LMlQtFiTFIqRzJ4e1Gj1PDZYkxiwpFMaldeNUNR1ZiMcCS6N6ovocgPGo7Ias/XCUU2qRzJoHTQCEUOCUeiWL5/usZwWihymgYliq0wLNFHhSKXqRzpmVDkNg1MFCVv9yUUeUxDE8XZOca9PisUKUaDE8XVBZijxRuhyCPmHMnoG1hCkWqEI9n8moQiL9ARiODJnkahSBUqR7ISilQlHMlmK8z2PqFQKHKZcCSLoy08a18XitwmHInuKMzWglEo8phwJCpzijSl49Dak0D7Vo1CkeJ0oP7sLUxEUuIuPEKRanSkvmQIxr09jVcqxrs/C6foUP2IHoylP+b07mPAKRZkqE0okpJw7EPEqrFmKPo8aKoTjpRWu1KM+EJAh4RjflHCwvCZrgjH3CIEo1CkS8KRu4QiXdMRuSpCKJasmLfem+3aGJzKkbMihCK8RjhypMdQ/B7P05tbCHYY0GflT8nHeLpX8ckx7T3GZ+X71n5263vohMqRpTcqxRKhsvYJhFfv/Xh1Q7lKcSDCka+3hs+lq6210L0bYvPw/DX9DFMGIhzzqHWhvjmnWHsYWuL8nD1Gwdk54ZhDjQvx7YWWCPNzy2pwOf+4Fei2+gxIOI6nxepzhGD8ms8z3rmJRcZgzHjMzTlh8ZXq2K225JS4uW0U2ULm6Rzs0FSO/Wu5TzFSxVjK8neKGDxCsQAnL7YnlUrrzdtHwajvlScUC1I59qd1KG79PPUIxQqczLiuVo0RQnHrMdboe4SmcswvWyhCCl6984oSikePt0XfIzQdNJ+IoXjm8Ws+X++ybSHqghOeR4ZQPHquN563JxZaGjLnGF+2UFx7HvOR1wjFAIRjXNlDces5BeU2oRiIRohlKziyhiLnCMWANEYcPVWKkJ4LJ4ann0siFKEwF1B7T4JRKMZnG05SFmTaurs4IRTjM4+YnHCM4+4twFx8sQjFTmjAdq5+at6StotFKHZG5fiu7/zTlQ9xWnLxxSIUO6VB37H8zJLl55hsff+ctoIXqRzr2grArcpRKEIQwrGOs0OttY8IXX4NaEA4lnVl/mlr7lEoxmKf4qA0ehlnQ7GnjyntnYWWwWn4Z65cQKVuKkFdQpFpmgyr7ypxAbn4YhGK/I9wvEYo9kcowgOfyU1aAf7zNBQFKtzXtCgxlFj39P6K38dwfuGctWtuGYxn3k1W7Jr7XeqBOvF9pRJq/VC95xNiGksI/GM5Kf90kl7AtmehJY/SQVikzUdfrd4KRfISirmtvXPszJ2s5t9TpA+MGo5rJ69UMKoa2xCKY9gqZIoXNqOFY81QpA2hmNfWYspW9bj8d9W3444SjldvHfbkeVyk73K+c3oafmuK9oXew/HqDWXnXHRQ3pWtN1tfe+VuVr0GwJNQ/Or13EALpQKt2r7GpVED4ExAjnpuojFVkV+JQHstFF97gsAMq2Oz0JJfjVC8+ziXjd7x3Hw2HqGYX4lAa36X/N4XZM44s22A+oRifunmFfeM3BGXc1klbjbBdUKxDynnFfeoHP+4uq/KQkEZzmFuqecV9zQ/gEaeBptgZHRdDaHXqByvE4yMrrsh9JoRw1G4vcv57ke3Q+g1IQ+qIsPpd1jc6ksXW3OuGrFypA6B2J/u5xX3pDjIQlSN5QnEfg0xr7hH5XiOYFznnPRnqHnFPaOEo3CDfUMPodekPfALogZj1ONiPMMPodeMUjlyjwDvmyH0ji5+iR1RL+7oCxne79y3IbfmXNXVL7MQORh/rfw3AqHYN/OKFxhWv2stCCMEpFDsn3nFi3oNx7fCJvvNcoVi/8wr3tRrOJZy96NbtzrOWpXYsmrsvoMPzBD6oR5/0RJhM3+MOwF59PzRF2T4Y60/tZ4GOWIIXUBvlWPpYJymex8ufiTCPCNt1XiBNIQuqLdwfGorsGoE5BsE8D1P52K3zvtWIJZ4viVbcx7qLRxrNuaVgGwdShZa7juaUjnagnVnkW7+fFfazLxiRb2F4xNnOuaZgFxeWEcXRElC8ZmnUyrz4LzjynSLecXKhOM/rrxin12o+Xb0N6pIofiOeVuuBVmL7WPmFSv53foAHlg28N1X67vhtTU3ufbvWvOV8/DVud/xmf2pYavyXHvOUtWivrNC5fjM2WF2rQpSpy5rWQ3O23ZvjrFGu+7NfRpCvyBz5VjK005y5udrV5CU8+Z0yNnjmP+/O6vQgvGG7Ccp4s0b9mQ5Tp552s625gSgcnzXndVPxmFeMZAeTli2auzMyrLV5/EcheKZfm74XJAFmfftVY9CcTxnA+3Mi+nZ7+eE7CcwW9U4t7YamfV34TrvbglO5diWUByTrTkJZD6hd6rGaJVmtOOhLu9uSSRz5Zg9GBmHrTkJZQvHuwEnGGnBvGJimcIxY8AdHXOmDexcY14xuUzheNfb4VPjDs/kYV6xE1nCMXp1JRAxhO5MlnC8q2aoCkS+DKE7lKEBRliEyXSs/GEI3bHeK0eowdacAURvjB4rqit32untd8/OvOJAIjdK62B8a75yS+S2GZF5xcEYVv/k/c7MmVccVNRwbFE1CkXmDKEHFzUc39QiFI/uCO4CassQmpDh+FbVqFJkyRCa/0QMx7uuhqoOy5etOfwQrfFG2PA954Jqy7wim3qoHLMGI22ZV2RXpHAsHXIZQnO5MBP9eHtgXpFTIoXjHWsBaKGFNYbQXBIlHEtUeUKRLYbQXBalgZ98WFYPoZhhCiCjK+9jv0JbDSBCI98NxjPOPG6p9zk/eZy1z7BueTzZH6eWtePzPnmqudPZW14gxPeZ/dn7OoR1poPqxJQiFEnhqJPqyJSiL5HKXmfVkcdQu52FIre1mijeWoTpYeWZ82psj7GpniIi7XOcJp15mn6+cIy0zedJP9CHKKpFR1rbtqJD/998D2fP56bENhh9iCpaVY46NGdtDZP1Iapq9VEEIzt7zp9Wjb2fa6FIVW92MK/05xlWb+v5fBDIG8NqoUgJ+g+vqtnhhOI9I61We18yYdXoeEKRs3wCI2GVHlb3XOXwDv2HEEqFo2qRJ/QbuuO9q9yh39AtnZs79BuAGaEIAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA9Oxv6fKxmTJrm1gAAAAASUVORK5CYII="
    },
    "baf72685-fd9a-47a4-8b3c-297ce6a4c2c2.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAABDIAAAIXCAIAAADheunIAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAADsMAAA7DAcdvqGQAADw3SURBVHhe7d0xjtvM1ibgWchEFwP8y5h0ggsvpTPHdwUOHYy30NEswGg4c+DYgBMnBhyMMYk38E11q1SiSlUkJVFkkXwenIAWSxRJldTntdit//YPAADAosQSAABgYWIJAACwMLEEAABYmFgCAAAsTCwBAAAWJpYAAAALE0sAAICFiSUAAMDCxBIAAGBhYgkAALAwsQQAAFiYWAIAACxMLAEAABYmlgAAAAsTSwAAgIWJJQAAwMLEEgAAYGFiCQAAsDCxBAAAWJhYAgAALEwsAQAAFiaWAAAACxNLAACAhYklAADAwsQSAABgYWIJAACwMLEEAABYmFgCAAAsTCwBAAAWJpYAAAALE0sAAICFiSUAAMDCxBIAAGBhYgkAW/P527//9X/yev87ru1RuuOHz3ElAA8klgCwOb8/XKSLf//r29e4turr++wuoYbvBcAExBIAtqcUMAY/9yiFmTGfsQBwP7EEgA36+eMpCxih+jOGK7gAFiSWALBFf5/f5Rmj/4qswgcs7378iisBeDCxBIBN+vXpJY8ZfZ9+FK7gevr0N64E4NHEEgC26arruApXcL08/4wrAXg4sQSAjbriOi5XcAEsTCwBYLPG/ha7K7gAliaWALBd4/7mryu4ABYnlgCwYWO+IbEwxteVAMxMLAFgy4av4yp8ouLrSgDmJpYAsGmlX3zvfhhSyC19X28CwEOIJQBsW+kLTE7BwxVcAE0QSwDYuNIXmBwv03IFF0AbxBIAtq5+HZcruAAaIZYAsH2V+HF5BZevKwFYhlgCwA6ULtb69OPiRl9XArAQsQSAPSh9gclFvfvxKw4HYF5iCQC7UPoCk6xcwQWwGLEEgH0o/eL7WbmCC2A5YgkAO1H6ApNOuYILYEFiCQB7UfoCk1S+rgRgSWIJALvRcx2XrysBWJRYAsCO1H7x/fDtigAsRSwBYE8KX2ASyhVcAAsTSwDYk2IscQUXwNLEEgB2pHgRlyu4ABYnlgCwG+VfeXcFF8DyxBIAdqLyvSWu4AJogFgCwEb97Xw94t+vle9SfPr0Nw4BYEFiCQCbNPCd7rF8VALQBrEEgG2qfUVJp3xUAtAKsQSAbRqMJf4AF0A7xBIAtunnj6csh5zq5YPPSQCaIpYAsFG/n9+9nCeTl6d3354//e78KjwAbRBLAACAhYklAADAwsQSAABgYWIJAACwMLEEAABY2AZjyf/+8X+VUkoppdRKK7Z07MzWYkk2rZVSSiml1Bor9nbshliilFJKKaVarNjesQ+biiXZVFZKKaWUUuut2OGxDz4tUUoppZRSLVZs79gHsUQppcr1r//z7+wWpZRSc1Zs79iHrcWSIJvQSil1W4klSim1bMXe7ii8LccltmiDseQgm9ZKKfX//vv/mLOyR1dKKXVVxZbuSCbZvM3GEoAkCwwzV9wJAO4glmyeWAJsXBYSFqm4KwDcRCbZA7EE2LgsIYyvf/7zX9ktN1fcFQCuJ5PshFgCbFkWD64qsQRgcTLJfoglwMZlCWF8iSUAixNL9kMsATYuSwjjSywBWJZMsitiCbB9WUgYWVPFkrgTAFxDJtkbsQSg4j//FRf8dASYnTfevRFLACo6sSTwAxJgNt5yd0gsAagQSwCW4P12n8QSgIrzWBL4SQlc5fO3f4f3jfP68DmuvN2DNtuIcDhxiZ0RSwAqxBLgPmLJDcLhxCV2RiwBqBBLgPuIJdcKxxKX2B+xBKDiIpYEfmR2bfu/bOF+YslVwoHEJXZJLAGoEEuGiCXQTyy5SjiQuMQuiSUAFaVYEvjBmYgl0E8sGS8cRVxir8QSgAqxZIhYAv28RkYKZyYusWNiCUBFJZYEfoIeaLmgn9fIGOG0xCX2TSwBqBBLhmi52Iu/Xz99e3p3mudP7749f/4bV/YYeo38+vzjw7uXp9Pal7Fb3pJw4HGJfRNLACrEEtiLX59eUmyI9f53XPX5Wyc2nNdghKjHkrdAkq86Vdjyz7ctFG3pfwTCzscldk8sAaiox5LAj1LYkFosKdx+UU+f6smkkh/GbPbf/3qpJo3NxJKw53EJxBKAqmZjya3Xk3T9fN3I69Uj+f/XhltePnz68bXnf2phc4qxZFx4eK2r8kN4iWW31Ovb17ihc2IJmySWAFT0xpLgYT9Qa/9x+7rqnutJDn7+7rt0pFt3XKAC6zI+gVTqpXzNVek1clWVP4rZxksv7HZcgjdiCUBFc7FkTOfUdz1J0JNqKnXDBSqwOr0vrpcPn//+igP/fq2NPP7fwZm+WPLy4dPv42br/+Pw7kcac7KBl17Y57gER2IJQMVQLAke85O1GEvG/29utTv5+ePaTHKo6gbFErai/vr6VrigsRw2Sh+YVGNJabPlV+jYza7ppRd2OC5Bh1gCUNFSLLmuiteT/H0eee1WobZ+gTu7V3vR1eZzcXzho8VKLKlt9uv7fKRYwo6IJQAVbcaSW68nKf1H7NP7H19/dhupv7U/WrrlC9yh9qIrXkAV/f6QDQ51Ob4cSyo5v7IbhdfUql96YW/jEpwTSwAqRsSS4AE/Yuux5PbrSQrbLF4K/6r4uUqpkRJL2Irii67/N7VKn2xcvEyKL8/qS2/0a2q9L72wq3EJLoglABXNxZJa2zGmo7psoYZ+l/3126Y/vP/x/Pn3r5/p85lzYglbUXwR9U/m0l0uLrgqvUaGXnr5eLGEvRBLACrGxZJg6h+05Vhy3/UkxWvWu38I6BZiCVsxKmNkxsx/saQr7GdcghKxBKCirVhy5/Uk5ahzqHffnkM+6e/AisQStqL0Aqn+Bkh0ayzpe41sOJaEnYxLUCGWAFSMjiXBpD9xixGiv+cY/r/ecX8d+PXCrfERRSxhK8SSxwp7GJegTiwBqGgqlkxwPUnpE5Weenl6P5RPVtcbQcWDXnRXv0ZGjl/dSy/sYVyCOrEEoKKpWDLFf9wWfwVlqF7DSbz7hdX1RlAhljxQ2L24BL3EEoCKa2JJMN2P3ofFkn/++fm7+LUkQ1Xpz9bVG0Fd6UU3MJlHJRmxJOxbXIIhYglARVOxZIr/uD2qfWdiXxX/DphYwlYUY0nfF4zc8b0lYgmUiSUAFVfGkmCiH8CPjiVHf79++vY0Op8U/hSYWMJWlGNJ36eUt3/L+45iSdixuAQjiCUAFU3FkoGe45Yk0/Hr8+/n90MR5fJ/jsUStqISS6p/mLs4foLovqVYEvYqLsE4YglAxfWxJJjiJ3G5Q7r/epJxftV+/+T+/wmGVtViSTnel//cdmnkbmNJ2KW4BKOJJQAVjcWSO68n+fk35I2vn398eP1U5CU0Vb1NzN/nQjK5+7p5aFU9loR6ef0ynzjw7dLHfMBbFf/jQCyB8cQSgIrWYskd15OUmpjej19KscSnJWxXbywZU5VrJvcZS8L+xCW4hlgCUHFTLAnu/pFc75Buvp6k/I0lHz6Xc055m363hO0qvehexv9BiOq032EsCTsTl+BKYglARYOxJNSN15PUtvn0/sfXn51w8rO6zZX1RnCNYix5Di+Hwm9t5VWN94FYAuOJJQDTu+8Hc28sGVPF60lu+Yr3U/neEjatEksOqyrhP9S7b+Vrt5K9xZKwJ3EJrieWAEzvvp/NxQ5pgutJypd7janKb9uLJWxFTyx5dfgUsfMafAqBpOdDkmRXsSTsRlyCm4glAA9xx0/oR11PEtT++G9Pvfv2tfb/wWIJWzEQSxgUzlhcgluJJQAPcccP6UddTxLVfyMlr+7vsZSIJWyFWHKvcMbiEtxKLAF4iDt+SD/qepJzx292f/0Ok7Spt629PL3/8ZALVKBVYsldwumKS3AHsQTgUW79Ua1Dgpl50d0unKu4BPcRSwAe5daf1jokmJkX3e3CuYpLcB+xBOCBbvqBrUOCmXnR3SicqLgEdxNLAB7opp/ZOiSYmRfdLcJZikswBbEE4LGu/8mtQ4KZedFdLZyiuAQTEUsAHuv6H946JJiZF93VwimKSzARsQTgsa7/4a1Dgpl50V0nnJ+4BNMRSwAe7sof4TokoF3hHSkuwaTEEoCHu/KnuFgCtCu8I8UlmJRYAjCHa36QiyVAo8LbUVyCqYklAHPwsxxYO+9jPJRYAjATP9GB9fIOxqOJJQAz8UMdWC/vYDyaWAIwEz/UgZXy9sUMxBKA+fjRDqyONy7mIZYAzMdPd2B1vHExD7EEYFZ+wAMr4i2L2YglALPyMx5YC+9XzEksAZibn/RA+7xTMTOxBGBuftgD7fNOxczEEoC5+WEPNM7bFPMTSwAW4Ec+0CxvUCxCLAFYgJ/6QLO8QbEIsQRgGX7wAw3y1sRSxBKAZfjZD7TG+xILEksAFqMDANrhHYlliSUAi9EEAO3wjsSyxBKAxWgCgEZ4O2JxYgnAkrQCwOK8EdECsQRgSboBYHHeiGiBWAKwMA0BsCBvQTRCLAFYmJ4AWIr3H9ohlgAsT2cAzM87D00RSwCWpzkA5uedh6aIJQDL0xwAM/O2Q2vEEoAmaBGA2XjDoUFiCUATdAnAbLzh0CCxBKAVGgVgBt5qaJNYAtAKvQLwaN5naJZYAtAQHQPwON5haJlYAtAQTQPwON5haJlYAtAQTQPwIN5eaJxYAtAWrQMwOW8stE8sAWiL7gGYnDcW2ieWADRHAwFMyFsKqyCWADRHDwFMxfsJayGWALRIJwHczzsJKyKWALRIMwHczzsJKyKWALRIMwHcydsI6yKWADRKSwHczBsIqyOWADRKVwHczBsIqyOWALRLYwHcwFsHaySWALRLbwFcy/sGKyWWADRNhwGM5x2D9RJLAJqmyQDG847BeoklAE3TZAAjebtg1cQSgNZpNYBB3ihYO7EEoHW6DWCQNwrWTiwBWAENB9DDWwQbIJYArICeA6jx/sA2iCUA66DzAC55Z2AzxBKAddB8AJe8M7AZYgnAOmg+gIy3BbZELFnQny8f//Pm+Xu8ZdtmPN4/3798/xOXF/D9+XCgH788eCcamkJpV8buzNV3aEjc9yX2WwsCJN4Q2BixZEEN9ZSzmOd4/3x/fn2YhyeCPnuMJaejDoZ3Z82h5Hioi+y4LgRIvCGwMWLJgsSSB5gtEfTZZyy5IpisO5QsmUpeaUSAwFsB2yOWMBuxZPtGBZOVh5J4kMs9uXoRwPsAmySWMBuxZA8Gg8kplKz09CydSl7pSGDPvAOwVWIJsxFL9uEUTArHv/pQEo9g4b3XlMCeeQdgq5qNJal7qfz471v/5/W3nj8eu5/g48fnL9//9LUR+V2G71E1/tGLbXp245/vX8LWDrcEh229DSyb5kBiY/l2bsMepO2Frb09+m3PTvF4k7v3/NQNn7l4sPBAhXN6xQN1vW0rbiidnoFYcsWRjn4i4lH2n+JX1adumplzUA8maU351FRVntyKERuvnoej2nN4uGO69bRjl6c8PUZ35a1n4JymBHbLy58N21ws6XQCl8rb6rnLsc0c67pHT6O7DU3nxvrWyk3ndAcSW6ewx3k7eDiI9EjlM1pZ3zm0eMvRJHue72p09mCvrXe8PXftc13b69dDPu7K5em58khHPxHHo6yf44Py+iv3aoTT7nbPQbr18sQMyA+/35jNp2MuDz4+YL423u10+jrnLjvnp33urrj9HGS0JrBDXvhs28ZiyXDzkndrnaaiIr9H3bWPnh67e/PwDr25OC1THsjxSM4+9HkTt1E++yfl9cXjnW7PK6e/c+fBJ+iKU9S318/P8YGy03P9kY5+ItIdKyc5Kq29fq9GOZ3s41lIt1RmTZ/Bp+7MqAdIx10afXy8fF2809kZOZ3A7ujTHp+fvnvOwhndCeyQFz7btqlY0mkPns8uP3m9/CWuOd9ep9l5u2Ql3ny4oCWuGNuVXf/o6R7dBzht5k13W2c7dX5epjyQoLO510e6+P/y0tnvKq8vHu/Ee562d7lnlUc6uwxr9AOdbayzrc6mgrOduOVIO/d53Vr9iTjdrXiWo8K6W/ZqpNOmX+9enhNL6tmj465frDmsyE/H6UiPd7i85RE0KLArXvJs3qZiSWoFSk1UoSM73aH8IAOrM1c/+sCNQWFLp9XdlZMeSJDuULlH2ovKBsvrS8c79Z6ne+R3OJ240qb6117oG35ad772piNNt1buVTqlpdsO0tZOa27aq9E6u3/8vOe2DT3Icf/yM3U8hRc7e7jDxYntHunrys5xP/Rw9SiwH17v7MFWPy25/H/lgmPzUG8dhkecXPvonXt025zuZkpbKR34tAcSpLaqMr60E13l9YXjnXzPa+NPp7XQU74aOuSugcM/PVZn/W1HOrRXhVN6ujG/T9rYI89/Jj3mwY1beZjj7p1PitoJjMPLB5GOtHPF3QyHq1OBPfBKZye29bslWQ/0+seEnr9UM8JxdE/vUOz6aq579KC49RsOfOoD6RxKZfgNOxlc7sT0e17b5NARXfVAN5yeG4909COdrT/eev5ghV24ca+ukg6h92GWcty77tGVT19wGFw7inSmjmY5XM0K7IFXOjuxsV957/ZAmbcL5882VR97aVyHcc2jB8WO74YDf9yB1Ebf9uxcHu/0e17Z9aEdDkYMORo6PYUBNx7p0CNdntI3x5u7dytt6sa9uk5lH69zzZ5eta/HDZ92rnTy3hxW9Gw7Hemrm8/XdTQrsHle5uzH1mJJWHH2G8wXTsPPeoghY3uM0Y/+Ku3BfbHkAQcyuhsevZOv0q3H433Anld2/eKhC4YOOhk6+sKmbj3S0U9EdlzH20/3K23p1r26Tm0fr3Lc/XGu2dfjlvNJebGNw4reTXdP6D1Hex0tC2yYFzi7sr1YEoV8UAsI8R5DW7jH8KO/SnvQ7WCuP/AHHMjobnj0Tr5Kt+Yd4IR7Xtn1EY90xc4MnZ7LAbce6egnIm+CjyuOdyxu6Na9uk51H69x3P9xrjqe46YPu5efuCSu6DmIdKAHDz2pXboW2DAvcHZls7Hk5M+f71lIiHdJnc497dKQ6qO/rjseQ3cHbjjw6Q+k2MR23PbsXB7vA56Cyq4PP9IV+zI0dLrnqHI0yeUpPTquOayobOfGvbpOfR8bcTwLr/t33NmLEx5X9BxDOpnpV95rz9r0NC6wSV7a7E2zsWSgIUutTmf9YPuTNWGnbUzRL1376EHxLunGSk9TWD/tgbzqP/nB9c9OcHm80+95bc8GHyk9PdVjPhnYWOnwbzzSoScibfZyq8dVr2tqm7lxr67Ts4+NiHsYzs5xXy/P9+EU1g+hO39mOa1n9C6wPV7X7FD7saT0o/3PaWW3geh2BvGmM2l93OKpewi3lO6RBlQ2eObaR6/0a0OPWVo/7YEEx32tDk8Hc7bvB5Vnp7MXp/tMvufVXT89UmlL/Wsv9A4vH/5tRzr0RKQ7XT4Nad3z9+NWCs/VTXt1nb59bETcxY9fvhzO1OXBHk5h9QjScx7veZoDsx20Dga2xCuafWo3lnQaprPv0e77Iu1TM/B2l05rcXa3U6fQuUN2j7M7jGvJrn70Yr821AeW1096IKfNVcff8uyUj3fiPT9t7+IeZ8/P6Q83n/+VgmvPUNDZ67Ct+uHfcqTH+9R2q3hKj7pP0qvSmJv26jq9+9iIuI/xK0cuD/Zwlgafhc4hns7rTIetiYEt8Ypmn9qNJecNU+70P8DnrULffaLzLmHEHcZ3ZFc+erFfSzdWHra2/iEH0nOPvserPDvF4514zzstYtR5sMFHOt+xfvWNpf90z/f7+iMtnseOyimNzs9F7eCmPf+X+vexEd1TdXm0h1NUOQunu54NqNx8Ot13ndRLmhjYDC9ndqvlWPL6E7zbViVvP89rDdt5K5b5+FxoBSqPcnBl73DVo6fB3X4t3Vh55J710x3IUDf85upnp3i8bx75FJzvQv2ROh/6jFTe2Ovj1c/flUc69ETUT+mb7rnoCwVTnv8LA/vYiNOpujzew7ryeei53/HJO1+Vbr33vF7SysAGeCGzZ23HkuDtSpLUM308XX3T17CFO33p3Cs43LHeBhzu0LnH0B16jH70Yr+Wbqx0LQPrpzmQoW746Lpnp3i8yYRPwZ+wpbiRsJn80fLnp7Pb1zvbVtjSYUP95++KIx16IvpPaWf9cCaY8PyfG9rHRtRP9eEAyk/B8V7loyutTbfVntPb6WZgA7yQ2bPmYwnAww0FwJXQ0MCqeQmzc2IJsHfpE52VpxI9DayY1y+IJcCu/Tn9/bS2rzMbR2cDa+SVC4FYAuzS6Tc/orV/VHKguYE18sqFQCwBdun09wDebCOUaG5ghbxs4UAsAfbp9HeRr//r0E3T4sCKeMFCIpYAbIouB1bECxYSsQRgazQ6sApeqtAllgBsjV4H2ud1ChmxBGCDdDzQMq9QuCSWAGyQpgda5hUKl8QSgA3S9ECzvDyhSCwB2CatDzTICxNqxBKAbdL9QIO8MKFGLAHYLA0QNMVLEnqIJQCbpQeCdng9Qj+xBGDLdELQAq9EGCSWAGyZZgha4JUIg8QSgC3TDMHivAxhDLEEYOO0RLAgL0AYSSwB2DhdESzICxBGEksAtk9jBIvw0oPxxBKA7dMbwfy87uAqYgnALuiQYE5ecXAtsQRgFzRJMCevOLiWWAKwC5okmI2XG9xALAHYC60SzMALDW4jlgDshW4JZuCFBrcRSwB2RMMED+UlBjcTSwB2RM8Ej+P1BfcQSwD2RecEj+CVBXcSSwD2RfMEj+CVBXcSSwD2RfMEk/OygvuJJQC7o4WCCXlBwSTEEoDd0UXBhLygYBJiCcAeaaRgEl5KMBWxBGCP9FJwP68jmJBYArBTOiq4x0SvoF+fXv4dNnVRHz7HAbmfP54uBod6+vQ3DoCVEksAdiq0MnEJuN5kr6C/z+/yjPFW377GAV2Vwe9+/IoDYLXEEoCdCq1MXAKuNPHL5/eHFDC69f53XH9U+Wjl5flnHAArJpYA7FdoaOISMNojXjifv12Ejdc6v5SrnF5cvsVGiCUA+xV6mrgEjPaYF87X93neeKvTpVzlAS7fYjPEEoBdC21NXAJGeOBLpvdSrvLHKS7fYkPEEoBdC51NXAKGPPr1Us0eLt9iB8QSgL0L/U1cAupmeaVULuUqlcu32BixBGDvQn8Tl4C6mV4plUu58nL5FpsjlgDsXWhx4hJQMefLpPJXubrl8i02SCwBQDKBPrO/QAYu5XL5FpsklgAglkCfBV4gPZdyuXyLjRJLAHgV2p24BHQs9dKoXcp18dXvsBGbiSV/vnz8z5vn7/GWycRNT7/hHsXDWWJHprSH52hOaQc+fvkTbxqhuNtLnD/aE9qduAQcLfa6+Pv87iKQHOv8q99hK8SSIcs0bJXDWXf3uI/naD5pB+6PJUudQ5oTOp64BCz5ivj16SWLIud1+up32A6xpN9SzVr1cOKKqxrRVuzlOZpL2oEpYklascqpxWRCuxOXgOVeET9/PJ0SSKVcysX2+N2SPk02at+fX/dJ9xjtuJm+LZb0MLUQS+BksZdD3+Vb3XIpF1sjltQdu76l/jO8ptX9WsIKzsVxF6dv9iePJaYWr0K7E5dgx5Z7IQxdvtUtl3KxLWJJVcP/dex/taM1nIg1xRJTiyD0OnEJdmypF0Ll8q0Pnyt/L9gXmLAljceSP3++P3/8GJuvVx8/Pn/5/qfQMqUWrfsfvdmNf75/CVs73BIctvU28FJPM3mx1c5GwzaLu3dw5+GcHFeP/l/t2G2OdFVT6jnqfxZ69rLjlh1O94mbPrtT9dQN7PZxdf9BsXGh14lLsEuLvQQql28dfo2k8veCfd0729FwLEntU8lFj1dstjo31rdW7L+Ow0u9ZPehypvtvVfJuMPp6Nu/kgfFkgkOqnNjfWurfI6ivt08uWWH07iwKkSSw3ImJJo4+mhot8ftL9sWep24BPuz3PyvXL51ulKr8tXvvl2RrWg2lgz30eddVbHZKvd3Fy4bsOOjF1uz00N9qW8/u+ckh9PVu4eXhh+/a7KN7vs5OjgO6z+pt+xwus/ZpzsXst0b3O0rpxYbFXqduAR7suDMr1y+df5hiEu52LRGY0lqnV7/t7d7/cqf76em7axtKjZbp8286W7r9YKXePNlA3a8X7kxK2z1uKa70c6OTHQ4ZxZvHic6qMLZPN5l9c/Rm/4dPbp+h4v3Od2ps+78kYd3Wy4hCI1OXII9WWzmVy7fuswbLuViwxqNJek/rkutU7Gv6r8xKGzptPp85UDj1t3qZed2WntaN9HhnBke8WCeo8KNaSujdPbs+h0+v88pxiSn1WeHVNrtc8Mj2IHQ6MQl2I3lpn3l8q3y1Vku5WKz2v+0pP4bz2eKrVR3M6WtpPXnq4+9ZflOg1tNrWlaPdHhnBvYy4eb6KCGzmZaf7660ecobXeUzr5dv8PD90nrz1aXdjuz9NSiEaHRiUuwAwtO+NqXJ1a/MNGlXGxUq79bkrd3Hz8+P3/paxeLzVa5Lzu5rW8b2upp59MGpjmccyOGPJbn6HIf8o306+zcDTs8eJ/y+tJuZ0YMYQ9ClxOXYAcWm/C1L0/s/U4Sl3KxSc3+yntPg1f8g63FTuoxfdtg01YaMMXhnBsx5OTmdrmX5+hN/xNwHNZ/Um/Z4XTTpKcuGDGEnQiNTlyCTTPVoQXtxpLQHL1+hcWhPSo578WKndRj+rbBpq2y2bsP51z5USoeE0umOKihw6iczIHzM7S+ttmJn6NXx2GV44tu2eHyMZyU1w8+0vCG2Q+9GntgnkMjWo4lUWgVa71ip2sqNluP6duG1vcPuONwzo0YcvKoWBJ5juqOw/pP6i07nG6a9NQFI4awHzo2ts0Mh3asIJac/PnzPesXT/1WsZN6TN82tNVTBujv6q4+nHMjhizAc5Q7DqvtycEtOzx0n/L6Ebs9Ygj7oWlj28xwaEebsWSwLbrs0Yp3KfdlJ5X1x60P3auyexc7N9XhnBvYy0eb6qDSjVt8jl4dh/U/T1fvcJDuc9WpG7HbAyeXfdG0sWGmNzSl0U9LUhNWaYxG9oi39W1DjVtaHVyOOK09rZvocM4Mj3gwz1H/Lhwdh1W2G12/w51brzp16dbqbg+PYGe0bmySiQ2tafUirtQEhp4q+9Lt7rden/qmYit1W992ur18t3SvN29/pCmuqH4h9zSHc+a4ycqhzcBzFG+sPkdXuGGHq6fmqLx+eLeXn1o0RvfGJpnY0Jp2f7ek0yXWdPuqYrN1W98W9DZmp3t9jAsFec83xeF0NdE6eo4ubrzVLTt826kb3O0mphat0cCxMaY0NKjlX3lP/VPJx+fzpqnYbN3Wt706rindsftQxU72Y7HhS3crGXU4Ha20jvcfVLpxc8/RVW7Z4XSfq07d0G5LJZTo4dgS8xna1HIsefXn7ZskYh/15uPbBS4XHVOx2bqtb3tzXFW4Z/ZQfzp/qCnsXM/Xgoexdx3OybF1rHSW8/IcTfAs3LLDPafmTXn9wG43NbVohzaOzTCZoVmtx5IFxf6s3g0v17lVd21vGn6OrtTKDptaVGnm2AYzGZolltQdO8W8UVy6g6zt1x61+hxdrY0dNrXo+p//61u3NHNsgGkMLRNL+sQuLfuv44U7SP+ffabJ5+h6TeywqcVJlkkOpaVjpdIEPizEW4HGtBhL3lojrhPP3Vzio7ID8SlnTw6t22WJJaxRmr1pJoeK64CWiCUbEc/dXOKjsgPxKWdPut1bt8IqyYR1SbNXLIH2uYhrWLy+Jl1ds9j1NvmOkDTzHN1q4R02tTjptm5ZhbViCetymLpZJgkVVwMtEUsAiLLWLaswQCxhRQ7z9jKThIojgJaIJQC8yvq2rOIgyYSVSFNXLIG1EEsAeJX1bVnFQWIJa5DmrUwCKyKWADA2kxxIJjTuMG+LmSRUHAQ0RiwB2LusacsqDuoQS2jZYd7KJLA6YgnA3mV9W7fiiAuSCW1KU7cYS+IgoEliCcCuZX1bVnHQBbGENh3mrUwCaySWAOxX1rdlFQeViCU06DBvi5kkVBwEtEosAdiprGnLKg6qk0xoSpq6PiqBlRJLAHYq69uyioPqxBLakeatTALrJZYA7FHWt2UVBw2RTGjEYd4WM0moOAhom1gCsDtZ05ZVHDSCWEILDvNWJoG1E0sAdifr27oVR4wmmbCsNHWLsSQOAtZALKn7/vyf6Pl7vCn58+VjXFdYyd6k6fDxy59407njXKqthzllfVtWcdBoYgnLOsxbmQQ2QCzpUc8efYmFHRrIJVIJ7cj6tqzioGuIJSzoMG+LmSRUHASshFjS6xRMug2lUEKuN5dIJbQia9qyioOuJ5mwiDR1fVQC2yCWDDhFkGNPeXkL9OUSqYRmZH1bVnHQ9cQS5pfmrUwCmyGWDDr/bEQooeI4NfKP0I6BxYRhYVnfllUcdCvJhJkd5m0xk4SKg4BVEUtGSFHk48fjNV16TC6Uc4lUQhOypi2rOOgOYglzOsxbmQQ2RiwZ4/QrJpEWk4JSLpFKaEDWtGUVB91NMmEeaeoWY0kcBKyQWDLOWTDRYVJ2mUukElqQ9W1ZxUF3E0uYQZq3Mglsj1gyUjeX5L88AFGeS6QSlpf1bVnFQVMQS5jBYd4WM0moOAhYJ7FklLMPSzSZVJ3nEqmEpWVNW1Zx0HQkEx4qTV0flcAmiSVjHHvN/6RfeddnUtbNJVIJi8v6tm7FEZMSS3icNHVlEtgqsWTYKZR8+XP62MSVXBTFKdKZLFIJC8n6tqzioKlJJjzIYd7KJLBhYsmQbig5+7dgQlHKJV8OU0UqYRFZ35ZVHPQAYgmPcJi3xUwSKg4CVk4s6Vf6dEQwodcxlxwu+ZNKWEDWtGUVBz2MZMK00tT1UQlsm1jS5xRKznrLys1wcJogpgjLyPq2rOKghxFLmFCatzIJbJ5YUteTPk4fmOg6uSS4sqSsb8sqDnoksYQJHeZtMZOEioOATRBLqvqv1XIlFz2O00MqYW5Z05ZVHPR4kgmTSFPXRyWwB2IJTE4qYTFZ39atOGIWYgn3S1NXJoGdEEtgYukSLqmEmWV9W1Zx0FwkE+50mLcyCeyHWAJT+vM9/V6Jy/uYVda3ZRUHzUgs4R6HeVvMJKHiIGBbxBKYwumXjSIflTCnrGnLKg6anWTCbdLU9VEJ7IpYAlPo/lHgQChhXlnfllUcNDuxhBukeSuTwN6IJTCJ78/pjwI/f5dJmFPWt2UVBy1BLOEGh3lbzCSh4iBgi8QSgBXLmras4qDlSCZcJU1dH5XADoklACuW9W3diiMWJZYwXpq6Mgnsk1gCsFZZ35ZVHLQ0yYSRDvNWJoHdEksAVinr27KKgxogljDGYd4WM0moOAjYNLEEYH2ypi2rOKgZkgn90tT1UQnsmVgCsD5Z35ZVHNQMsYQead7KJLBzYgnAymR9W1ZxUEvEEnoc5m0xk4SKg4AdEEsA1iRr2rKKg9ojmVCUpq6PSgCxBGBNsr6tW3FEk8QSLqWpK5MAgVgCsBpZ35ZVHNQqyYTMYd7KJMCBWAKwDlnfllUc1DCxhK7DvC1mklBxELAnYgnACmRNW1ZxUPMkEw7S1PVRCZCIJQArkPVtWcVBzRNLCNK8lUmALrEEoHVZ35ZVHLQGYgnBYd4WM0moOAjYH7EEoGlZ05ZVHLQeksnOpanroxIgI5YANC3r27oVR6yKWLJnaerKJMAlsQSgXVnfllUctDaSyW4d5q1MAhSJJQCNyvq2rOKgFRJL9ukwb4uZJFQcBOyYWALQoqxpyyoOWi3JZG/S1PVRCVAjlgC0KOvbuhVHrJlYsitp6sokQA+xBKA5Wd+WVRy0ZmLJrhzmrUwC9BNLANqS9W1ZxUHrJ5nsRJq6YgnQTywBaEjWtGUVB22CWLIHaerKJMAgsQSgIVnfllUctBWSybaleSuTAGOIJQCtyPq2rOKgDRFLtu0wb4uZJFQcBHAklgA0IWvasoqDNkcy2ao0dX1UAowklgA0IevbuhVHbJFYsklp6sokwHhiCcDysr4tqzhoi8SSTTrMW5kEuIpYArCwrG/LKg7aLslkY9LUFUuAq4glAEvKmras4qBNE0u2JE1dmQS4llgCsKSsb8sqDto6yWQb0ryVSYAbiCUAi8n6tqzioB0QS7bhMG+LmSRUHARQIZYALCNr2rKKg3ZDMlm7NHV9VALcRiwBWEbWt3UrjtgTsWTV0tSVSYCbiSUAC8j6tqzioD0RS1btMG9lEuAeYgnA3LK+Las4aH8kk5VKU1csAe4hlgDMKmvasoqDdkksWaM0dWUS4E5iCcCssr4tqzhorySTdUnzViYB7ieWAMwn69uyioN2TCxZl8O8LWaSUHEQwDhiCcBMsqYtqzho9ySTtUhT10clwCTEEoCZZH1bt+IIxJKVSFNXJgGmIpYAzCHr27KKgxBLVuIwb2USYEJiCcDDZX1bVnEQR5JJ49LUFUuACYklAI+VNW1ZxUF0iCUtS1NXJgGmJZYAPFbWt2UVB3FOMmlTmrcyCTA5sQTggbK+Las4iAtiSZsO87aYSULFQQA3EUsAHiVr2rKKg6iQTFqTpq6PSoBHEEsAHiXr27oVR1AnljQlTV2ZBHgQsQTgIbK+Las4iDqxpCmHeSuTAI8jlgBML+vbsoqDGCKZNCJNXbEEeByxBGBiWdOWVRzECGJJC9LUlUmAhxJLACaW9W1ZxUGMI5ksK81bmQR4NLEEYEpZ35ZVHMRoYsmyDvO2mElCxUEAUxBLACaTNW1ZxUFcSTJZSpq6PioBZiCWAEwja9qyioO4nliyiDR1ZRJgHmIJwDSyvi2rOIjriSXzS/NWJgFmI5YATCDr27KKg7iVZDKzNHXFEmA2YgnAvbKmLas4iDuIJXNKU1cmAeYklgDcK+vbuhVHcDfJZB5p6sokwMzEEoC7ZH1bVnEQdxNL5nGYtzIJMD+xBOB2Wd+WVRzERCSTR0tTVywB5ieWANwoa9qyioOYjljyUGnqyiTAIsQSgBtlfVtWcRDTEUseJ81bmQRYilgCcIusb8sqDmJqksmDpKkrlgBLEUsArpY1bVnFQTyAWPIIaerKJMCCxBKAq2V9W7fiCB5GMplWmroyCbAssQTgOlnfllUcxMOIJdM6zFuZBFicWAJwhaxvyyoO4sEkk6mkqSuWAIsTSwDGypq2rOIgHk8smUSaujIJ0AKxBGCsrG/LKg7i8cSS+6V5K5MAjRBLAEbJ+ras4iDmIpncKU1dsQRohFgCMCxr2rKKg5iRWHKPNHVlEqAdYgnAsKxv61Ycwewkk9ukqSuTAE0RSwAGZH1bVnEQsxNLbnOYtzIJ0BqxBKBP1rdlFQexEMnkWmnqiiVAa8QSgKqsacsqDmI5YslV0tSVSYAGiSUAVVnfllUcxHLEkvHSvJVJgDaJJQBlWd+WVRzE0iSTkdLUFUuANoklAAVZ05ZVHEQDxJIx0tSVSYBmiSUABVnf1q04gmZIJv3S1JVJgJaJJQC5rG/LKg6iGWJJv8O8lUmAxoklAGeyvi2rOIjGSCY1aeq2E0s+f/t32Jnz+vA5rrzdgzYLzEYsATjJmras4iDaE3rQuERHmrrh/KTlVHHQ7MQSoEgsATjJ+ras4iDaE3rQuMRRmrfh5KTlVHHQEsQSoEgsAYiyvi2rOIhWhTY0LvEmTd1wZtJyqjhoCWIJUCSWALzKmras4iAaFtrQuETDmSQQS4AisQRAJtmI0InGpX1LUzeckLScKg5ajlgCFIklAGLJRoRONC7tWJq34Wyk5VRx0KLkB6BILAH2LuvbsoqDWInQ48alvUpTN5yKtJwqDlqUWAIUiSXArmVNW1ZxEOsRety4tEtp6obzkJZTxUFT+fv107end6do8fTu2/Pnv3Flj6FY8uvzjw/vXp5Oa1/GbhlYNbEE2LWsb+tWHMGqhC42Lu1PmrrhJKTlVHHQFX59ekmxIdb733HV52+d2HBegxGiHkveAkm+6lRhyz/ftlDkQxhYO7EE2K+sb8sqDmJtQj8al3YmTd1wBtLyoeKI69RiSeH2i3r6VE8mlfwwZrP//tdLNWmIJbB2YgmwU1nfllUcxAqFfjQu7UmauuHw03KqOOg6xVgyLjy81lX54end2M3++1/fvsYNnRNLYO3EEmCPsqYtqziI1QotaVzahzR1w4Gn5VRx0NXGJ5BKvZSvuSrlh6uq/FGMWAJrJ5YAe5T1bVnFQaxWaEnj0g6keRuOOi2nioNu0RtLXj58/vsrDvz7tTby+LsoZ/piycuHT7+Pm63/Bsu7H2nMiVgCayeWALuT9W1ZxUGsXOhK49LWpakbDjktp4qDblGPJd++Xn4MUg4bpQ9MqrGktNmfP0rJZOxmxRJYE7EE2JesacsqDmL9QlcalzYtTd1wvGk5VRx0o1osqfX6xfGFC64qsaS22a/v85FiCWyTWALsS9a3dSuOYBNCVxqXtitN3XCwaTlVHHS7ciwpXkAV/f6QDQ51Ob4cSyq/yF7ZjULeEEtg7cQSYEeyvi2rOIitCI1pXNqoNHXDkablQ8URdxn76UdH6ZONi7xRjCXF30I5GJk3xBJYO7EE2Iusb8sqDmJDQmMal7YoTd1wmGk5VRx0l7EfU3SU7nJxwVUpP9zwPSc5sQTWTiwBdiFr2rKKg9ic0JvGpW1JUzccYFpOFQfda1TGyIzJBmIJUCSWALuQ9W1ZxUFsTuhN49KGpHkbji4tp4qDJlCKJdXfAIlujSV9+UEsgZ0QS4Dty/q2rOIgNiq0p3FpK9LUDYeWllPFQRMQS4BZiSXAxmVNW1ZxENsV2tO4tAlp6objSsup4qBpzHkRl1gCiCXA1mV9W7fiCB5s2X4xPFxcWr80dcNBpeVUcdBkxBJgVmIJsGVZ35ZVHMSDLd4vhkeMSyuXpm44orR8qDhiSqVYMvDEjUoyYglQJJYAm5X1bVnFQTyeWDKJNHXD4aTlVHHQlIqxpO8LRu743hKxBBBLgI3Kmras4iBm0UK/GB40Lq1TmrrhQNJyqjhoYuVY0vdb77d/y7tYAoglwEZlfVtWcRCzEEvulOZtOIq0nCoOml4lllS/Y6Q4vjBYLAGKxBJgg7K+Las4iLk00i+Gx41La5OmbjiEtJwqDppeLZaUf/H954+nfFio0kixBCgSS4CtyZq2rOIg9if0qXFpVdLUDfufllPFQQ9RjyWhXj58+n28Ouvv10/fSpmk8osoYglQJJYAW5P1bd2KI9il0KfGpfVIUzfsfFpOFQc9Sm8sGVOVvyYslgBFYgmwKVnfllUctBNv/4f97tSiPb379vy5/FsBVT9fN/Lh3Ut3O28Vbnn58OnH1/5vsWhN2PO4tBJp6oY9T8uHiiMeqBRLLmdCtaqRQCwBisQSYDuyvi2rOGg7Cl3j8ZqZX58rF9WEGhlOfv7+MLIBHdxgO/1ieOi4tAZp6obdTsup4qAHKsaS55BUC38FOK8PPVNCLAGKxBJgI7KmLas4aFNqsWTMtTe1P6YU9aSaSvVtsKl+MTx6XGpbmrphh9NyqjjosSqx5LCqN/f2f4YmlgBFYgmwEVnfllUctCnFWDImkxyq2rGV/6TScFU3KJZcK83bsLdpOVUc9HA9seTV2wV+t1wlKJYARWIJsAVZ35ZVHLQ14xNIpYq/kfz3efQvD1xU5Yv2WusXww7EpValqRt2NS2nioMebiCWAExLLAFWL2vasoqDNqg3lrx8+Py38/dbKyMv/35r6aOSp/c/vv7s/i/431+ffxR/86R8KZdYcpU0dcN+puVUcdAcxBJgVmIJsG5Z05ZVHLRN9VjyrfAHskrZ4LLLLGyz+NUTr4qfq5Q+MBFLxktTN+xkWk4VB81ELAFmJZYA65b1bVnFQdtUiyW1jr84Pvt84/KPLA39LvvL07tvH97/eP78+9fP9PnMudZiSRD2IS61pDt1wx52/xkqDpqPWALMSiwBVizr27KKgzarHEve/Shng1e/P2SDQ52PL/3t1+73ed9ELBkpTd2we2k5VRw0H7EEmJVYAqxV1rRlFQdt2ZhPPzKl1HF22VU56hzq3bfnkE9u6EobjCVB2I241IY0dcOOpeVUcdCsxBJgVmIJsFZZ39atOGLjihGiv90fbjTH/XXg1wu3xkcUsWRQmrphr9JyqjhobmIJMCuxBFilrG/LKg7auFu6xhEJYcx3eHfq5en9UD5pM5YEYU/i0tLS1A27lJYPFUcsQCwBZiWWAOuT9W1ZxUHbV+oaK98ckoxKCKVfQRmq13AS735BLOmXpm7Yn7ScKg5agFgCzEosAVYma9qyioN24WGx5J9/fv4ufi3JUFV6VrGkR5q6YWfScqo4aBliCTArsQRYmaxvyyoO2oVbusZrEkLtOxP7qvh3wJqNJUHYmbi0hO7UDXvS/WeoOAhgH8QSYE2yvi2rOGgvHh1Ljv5+/fTtaXQ+KfwpsIfEkv/8V1y4T9iZuLSENHXDbqTlVHEQwD6IJcBqZE1bVnHQjpRiyUC7f+dlOb8+/35+PxRRLr8V/lGflqw8maSpG3YgLaeKgwB2QywBViPr27oVR+xLMZYUUkHH4PeWjPer9vsnl9dxPfAirpBM7g4nYX/i0ozS1A2PnpZTxUEAeyKWAOuQ9W1ZxUH7Uo4lfTFj+Fve//n5N+SNr59/fHj9VOTlaSA//H0uJJOLHXj475asMJmkqRseOi0fKo4A2BmxBFiBrG/LKg7anUosqX7Re3H82eBSfuj9+KUUS2b9tCS5L5mEXYpLs0hTNzxuWk4VBwHsjFgCtC5r2rKKg/aoFkvKvy5S/vr2bGT5G0s+fC7nnPI25/vdkswdySTsUlx6vDR1w4Om5VRxEMD+iCVA67K+Las4aI/qsSTUy4dP6fsN3/6UVj7grS4iRPUTmPc/vv7shJOf1W0W8sZMsSRoPpl0p254xO4/Q8VBALsklgBNy/q2rOKgneqNJWOq+De4bvmK91Mt/70lIZncFE7CXsWlR0pTNzxcWk4VBwHsklgCtCtr2rKKg/arFEtexn/BSDUYlC/3GlOV37afNZYcNJlM0tQND5SWU8VBAHsllgDtyvq2bsURu1aMJc8//5b+CnBe1V8XOaj98d+eevfta+37TxaIJcH1ySTsWFx6gDR1w6Ok5VRxEMCOiSVAo7K+Las4aNcqseSwqvLLJKHefRv3/Yn130jJq/t7LCXLxJKgpWSSpm54iLR8qDgCYN/EEqBFWd+WVRy0dz2x5NXht9I7H3o8hUDS/yFJyfGb3V+/wyRt6m1rL0/vf4za4GKxJLjyV03CvsWlSaWpG7afllPFQQD7JpYAzcmatqziIIZiCSejk0k4h3FpOmnqho2n5VRxEMDuiSVAc7K+Las4CLHkKgslk+7UDVvu/jNUHASAWAK0JuvbsoqDeCWWXGlcMgmnMS5NIU3dsNm0nCoOAkAsAZqSNW1ZxUFEYsn1xv2qSTiTcek+aeqGDablVHEQAG/EEqAhWd/WrTiCE7HkVkPJJJzJuHSHNHXD1tJyqjgIgCOxBGhF1rdlFQdxIpbc4fHJJE3dsKm0fKg4AoAOsQRoQta3ZRUHcUYsuU9vMgknMy7V/b///j+6FW99k6Zu2E5aThUHAdAhlgDLy5q2rOIgcmLJ3erJJJzMuFSRZZJDHValqRs2kpZTHcYAkBFLgOVlfVtWcRA5sWQKIZlUwkk4n3HpQpZGUoVV3akbttD9Z6jD3QG4JJYAC8v6tqziIArEkumUkkk4n3HpQpZGUoVVaeqGu6flVIe7A3BJLAGWlDVtWcVBMIPRySSLIt1KU1cmAbiWWAIsJmvasoqDYDYXyeQylmQ5JKvD1JVJAG4glgCLyfq2rOIgmNPFr5p0k0kWQrJKU/cylsT7A1AnlgDLyPq2rOIgWEQnmVwbS3xUAnAbsQRYQNa0ZRUHwYKOySTFkiyEZHWYujIJwM3EEmABWd/WrTgCFtdJJlkIySrN3stYctgCAIPEEmBuWd+WVRwELXj7VZMshGSVpq5MAnAPsQSYVda3ZRUHQVN6k8lh6rp8C+BOYgkwn6xpyyoOgpYcskctmRymrkwCcD+xBJhP1rdlFQdBM7oJ5DKZpKnr8i2A+4klwEyyvi2rOAhakuWQbjJJU9dHJQCTEEuAmWR9W7fiCGhJSiDdCsnkEE4OU1cmAZiKWALMJGvdUsXV0JJuFLmskEwOs9flWwBTEUuAmWTdW6q4GpqRhZCswoDXefuf/5JJACYklgAzyRq4Q8V10JIsh3QrRJEwIExdl28BTEssAeajh6N9WQ7J6hBLZBKAyYklABBlISSrw5iQSS5jyWEVADcTSwAgynJIt+KIt1gSlwCYjlgCAFEWRVLF1QA8jFgCAFGWRlLF1QA8jFgCAFGWRg4V1wHwSGIJAJzIJACLEEsAAICFiSUAAMDCxBIAAGBhYgkAALAwsQQAAFiYWAIAACxMLAEAABYmlgAAAAsTSwAAgIWJJQAAwMLEEgAAYGFiCQAAsDCxBAAAWJhYAgAALEwsAQAAFiaWAAAACxNLAACAhYklAADAwsQSAABgYWIJAACwMLEEAABYmFgCAAAsTCwBAAAWJpYAAAALE0sAAICFiSUAAMDCxBIAAGBhYgkAALAwsQQAAFiYWAIAACxMLAEAABYmlgAAAIv655//D8aPKoYKT45aAAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "id": "431d06ac-7e68-45ae-b18b-e40a7698b4f7",
   "metadata": {},
   "source": [
    "# Calculate the sky view factor\n",
    "Use [A-10 from Siegel's & Howell's \"Thermal Radiation\"](http://www.thermalradiation.net/sectiona/A-10.html):\n",
    "\n",
    "[![image.png](attachment:a2e1b105-bbbd-4836-8851-25cd9305ea37.png)](http://www.thermalradiation.net/sectiona/A-10.html)\n",
    "[![image.png](attachment:8a27b5e2-3f7d-4258-b2e6-b4c253987a29.png)](http://www.thermalradiation.net/sectiona/A-10.html)\n",
    "\n",
    "with $Y = y/x$ in the image above.\n",
    "\n",
    "$$\n",
    "dF_{d1-d2} = \\frac{Y\\sin^2\\phi \\mathit{dY}}{2\\left(1+Y^2-2Y\\cos\\phi\\right)^{3/2}}\n",
    "$$\n",
    "\n",
    "Just want to integrate from zero to $Y_{\\psi} = \\frac{sin\\left(\\psi+\\phi\\right)}{\\sin\\psi}$ where $y$ and $x$ are related using the sine rule. ",
    "Use [trig identities](https://en.wikipedia.org/wiki/List_of_trigonometric_identities) to replace $\\sin\\left(\\pi-\\psi-\\phi\\right)$, the interior angle facing `y`, ",
    "with $sin\\left(\\psi+\\phi\\right)$.\n",
    "![vfsky.png](attachment:baf72685-fd9a-47a4-8b3c-297ce6a4c2c2.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6ad5c2b6-ab59-4885-8437-4b30b95b67e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_vf(y, phi):\n",
    "    \"\"\"calculate view factor from d1 -> d2\"\"\"\n",
    "    phirad = np.radians(phi)\n",
    "    sinphi = np.sin(phirad)\n",
    "    cosphi = np.cos(phirad)\n",
    "    numerator = y * sinphi**2\n",
    "    denominator = 1 + y**2 - 2*y*cosphi\n",
    "    denominator *= np.sqrt(denominator)\n",
    "    return numerator / (2*denominator)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a6be2fd-638d-41a4-936f-66c68205fc9d",
   "metadata": {},
   "source": [
    "## tests\n",
    "Try several examples and compare to these formulations:\n",
    "1. $\\mathit{VF}_{sky} = \\frac{\\cos\\psi + \\cos\\phi}{2}$\n",
    "2. $\\mathit{VF}_{sky} = \\frac{1 + \\cos\\left(\\psi + \\phi\\right)}{2}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a3e99cf2-e237-478e-a65f-2dab6ee71042",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Front side tilted 20-deg\n",
    "phi = 20.0\n",
    "psi = 10.0\n",
    "ypsi = np.sin(np.radians(phi+psi))/np.sin(np.radians(psi))\n",
    "y = np.linspace(0, ypsi, 1000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "06c10a6a-b61f-4c05-8029-d12c909dcd9f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9330127018921618"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vf = np.trapz(calc_vf(y, phi), y)\n",
    "vf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "84f94ee7-6b49-4cd8-8e89-c745d9aafdb8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9622501868990583"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test 1\n",
    "(np.cos(np.radians(phi)) + np.cos(np.radians(psi)))/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "da8a6d0b-7445-4a50-b2ae-b057cb64c5af",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9330127018922194"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test 2\n",
    "(1+np.cos(np.radians(phi+psi)))/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9cbb2592-feec-49c7-b6f0-28dce1814056",
   "metadata": {},
   "outputs": [],
   "source": [
    "psi = 5.0\n",
    "ypsi = np.sin(np.radians(phi+psi))/np.sin(np.radians(psi))\n",
    "y = np.linspace(0, ypsi, 1000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7f5b3b51-8d0c-4987-89f4-e3fc8f8d0786",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9531538935182058"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vf = np.trapz(calc_vf(y, phi), y)\n",
    "vf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b89618c5-2777-40db-b92e-e5e488260931",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.967943659438827"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test 1\n",
    "(np.cos(np.radians(phi)) + np.cos(np.radians(psi)))/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a2beac94-8388-4234-b95f-67d2ed449bfe",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.953153893518325"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test 2\n",
    "(1+np.cos(np.radians(phi+psi)))/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "794e71de-e34e-43f6-944d-d2372ffcdd48",
   "metadata": {},
   "outputs": [],
   "source": [
    "# test close to zero, should compare to (1+cos(phi))/2\n",
    "psi = 0.001\n",
    "ypsi = np.sin(np.radians(phi+psi))/np.sin(np.radians(psi))\n",
    "y = np.linspace(0, ypsi, 1000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "73d922d0-13ed-4a9b-8dfc-f23ab448a4aa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9698463103929542"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test psi=0\n",
    "(1+np.cos(np.radians(phi)))/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "d69c716d-6c2c-4bf7-8f73-afa9fc690ad0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9698414540979982"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vf = np.trapz(calc_vf(y, phi), y)\n",
    "vf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "7b832383-5001-4309-ab18-b6701a873da7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9698463103167998"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test 1\n",
    "(np.cos(np.radians(phi)) + np.cos(np.radians(psi)))/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "7bb0cc53-8748-465b-94f4-a7af09438012",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.969843325632588"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test 2\n",
    "(1+np.cos(np.radians(phi+psi)))/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "ed508c2f-fb7b-4613-8970-14c7c1e1b190",
   "metadata": {},
   "outputs": [],
   "source": [
    "# backside of 20-deg tilted module\n",
    "phi = 160.0\n",
    "psi = 10.0\n",
    "ypsi = np.sin(np.radians(phi+psi))/np.sin(np.radians(psi))\n",
    "y = np.linspace(0, ypsi, 1000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "87acea9c-e4cb-46af-81e9-ee312693ece2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0075961234938907805"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vf = np.trapz(calc_vf(y, phi), y)\n",
    "vf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "863c7b2c-bdaf-4cd1-bd00-d45a78faf294",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.02255756611314985"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test 1\n",
    "(np.cos(np.radians(phi)) + np.cos(np.radians(psi)))/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "22025b71-6522-4dc7-b7ad-dae88253e4e4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.00759612349389599"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test 2\n",
    "(1+np.cos(np.radians(phi+psi)))/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "dd030411-1e33-4709-9294-6149464215a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "psi = 5.0\n",
    "ypsi = np.sin(np.radians(phi+psi))/np.sin(np.radians(psi))\n",
    "y = np.linspace(0, ypsi, 1000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "68d2e6a6-ce36-4b0e-b174-6169b2eeb182",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.017037086855421994"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vf = np.trapz(calc_vf(y, phi), y)\n",
    "vf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "2e2fc9c9-250a-4ac7-b72d-64669fe19c10",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.028251038652918614"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test 1\n",
    "(np.cos(np.radians(phi)) + np.cos(np.radians(psi)))/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "9c57a4f1-abe9-4f46-88ae-19d1415ec990",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0170370868554659"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test 2\n",
    "(1+np.cos(np.radians(phi+psi)))/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "74c0ecbe-346b-447f-b636-068fe0c1171b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# test close to zero, should compare to (1+cos(phi))/2\n",
    "psi = 0.001\n",
    "ypsi = np.sin(np.radians(phi+psi))/np.sin(np.radians(psi))\n",
    "y = np.linspace(0, ypsi, 1000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "1ae9ef1b-de7b-4d65-9c62-e1fe01a69703",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.03015368960704584"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test psi=0\n",
    "(1+np.cos(np.radians(phi)))/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "120d0826-51c6-4335-8150-bfd038f88d44",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.030148833814086378"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vf = np.trapz(calc_vf(y, phi), y)\n",
    "vf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "2a45c267-8534-4323-a985-7d2e6e33a538",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.03015368953089148"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test 1\n",
    "(np.cos(np.radians(phi)) + np.cos(np.radians(psi)))/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "863cfc35-1d31-4a8f-9821-5ab3bcdd9d47",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.030150704989803045"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test 2\n",
    "(1+np.cos(np.radians(phi+psi)))/2"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}