diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..4a2c337 --- /dev/null +++ b/.gitignore @@ -0,0 +1,89 @@ +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] + +# C extensions +*.so + +# Distribution / packaging +.Python +env/ +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +*.egg-info/ +.installed.cfg +*.egg + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover + +# Translations +*.mo +*.pot + +# Django stuff: +*.log + +# Sphinx documentation +docs/_build/ + +# PyBuilder +target/ + +# DotEnv configuration +.env + +# Database +*.db +*.rdb + +# Pycharm +.idea + +# VS Code +.vscode/ + +# Spyder +.spyproject/ + +# Jupyter NB Checkpoints +.ipynb_checkpoints/ + +# exclude data from source control by default +/data/ + +# Mac OS-specific storage files +.DS_Store + +# vim +*.swp +*.swo + +# Mypy cache +.mypy_cache/ \ No newline at end of file diff --git a/Examples.ipynb b/Examples.ipynb index 2170a81..af77686 100644 --- a/Examples.ipynb +++ b/Examples.ipynb @@ -1,98 +1,224 @@ { "cells": [ { - "cell_type": "code", - "execution_count": 3, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "from solarlib import *\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np" + "# Using solarlib" ] }, { "cell_type": "code", "execution_count": 4, + "metadata": { + "ExecuteTime": { + "end_time": "2021-02-21T06:57:42.864211Z", + "start_time": "2021-02-21T06:57:42.861210Z" + } + }, + "outputs": [], + "source": [ + "from solarlib.location import Location\n", + "import datetime\n", + "import pytz" + ] + }, + { + "cell_type": "markdown", "metadata": {}, + "source": [ + "## Define a location" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "ExecuteTime": { + "end_time": "2021-02-21T06:57:15.310798Z", + "start_time": "2021-02-21T06:57:15.308798Z" + } + }, "outputs": [], "source": [ - "day=0 # for January 1st\n", - "h=13.5 # 1:30 pm\n", - "lat=51.5074 # Latitude\n" + "latitude = -31.9505\n", + "longitude = 115.8605\n", + "timezone = 'Australia/Perth'\n", + "perth = Location(latitude,longitude,timezone)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "ExecuteTime": { + "end_time": "2021-02-21T06:58:39.006792Z", + "start_time": "2021-02-21T06:58:39.003791Z" + } + }, + "source": [ + "## Sunrise\n", + "The output is a python datetime object." ] }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 20, "metadata": { - "scrolled": true + "ExecuteTime": { + "end_time": "2021-02-21T07:05:29.130067Z", + "start_time": "2021-02-21T07:05:29.127067Z" + } }, "outputs": [ { - "data": { - "text/plain": [ - "(15.799012631177902, '15:47:56')" - ] - }, - "execution_count": 5, - "metadata": {}, - "output_type": "execute_result" + "name": "stdout", + "output_type": "stream", + "text": [ + "2020-06-27 07:17:33.495057+08:00\n" + ] } ], "source": [ - "Sunset(lat,day)" + "print(perth.sunrise('2020-06-27'))" ] }, { - "cell_type": "code", - "execution_count": 6, + "cell_type": "markdown", "metadata": {}, + "source": [ + "## Sunset\n", + "The output is a python datetime object." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "ExecuteTime": { + "end_time": "2021-02-21T06:59:27.398504Z", + "start_time": "2021-02-21T06:59:27.395502Z" + } + }, "outputs": [ { - "data": { - "text/plain": [ - "(8.200987368822098, '08:12:03')" - ] - }, - "execution_count": 6, - "metadata": {}, - "output_type": "execute_result" + "name": "stdout", + "output_type": "stream", + "text": [ + "2020-06-27 17:21:50.791503+08:00\n" + ] } ], "source": [ - "Sunrise(lat,day)" + "print(perth.sunset('2020-06-27'))" ] }, { - "cell_type": "code", - "execution_count": 7, + "cell_type": "markdown", "metadata": {}, + "source": [ + "## Day length\n", + "The output is a python timedelta object." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "ExecuteTime": { + "end_time": "2021-02-21T07:00:03.861630Z", + "start_time": "2021-02-21T07:00:03.857624Z" + } + }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "252.5009632059124\n" + "10:04:17.296446\n" ] } ], "source": [ - "h=np.linspace(0,24,49)\n", - "irr=SolarIrradiance(lat,day,13)\n", - "print(irr)" + "print(perth.day_length('2020-06-27'))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Solar irradiance\n", + "The output is irradiance in $kW/m^2$." ] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 11, + "metadata": { + "ExecuteTime": { + "end_time": "2021-02-21T07:00:51.740736Z", + "start_time": "2021-02-21T07:00:51.735734Z" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "0.7360521598862617" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "perth.solar_irradiance('2020-06-27 14:27:00')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Irradiance throughout a day\n" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": { + "ExecuteTime": { + "end_time": "2021-02-21T07:08:04.096221Z", + "start_time": "2021-02-21T07:08:04.094221Z" + } + }, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 22, "metadata": { - "scrolled": false + "ExecuteTime": { + "end_time": "2021-02-21T07:08:21.739676Z", + "start_time": "2021-02-21T07:08:21.637643Z" + } }, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD8CAYAAACW/ATfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHyBJREFUeJzt3XlwnPd93/H3F/d9L0CQxEFCoHgpImmYki0fkmXZktpG9thupKa2mnGrTCt3bNdtYrdN7cxEqSetHY87tjJy7EbO+Kgb2ZZcK7ZlRZajODpIieJ9gQRAECBO4gZx/voHnqVACudez+4+n9cMZhcPdvf5Plzisw9+z/f5PeacQ0RE0luG3wWIiEj8KexFRAJAYS8iEgAKexGRAFDYi4gEgMJeRCQAFPYiIgGgsBcRCQCFvYhIAGT5XQBAVVWVa2xs9LsMEZGUcvDgwX7nXGgtj02KsG9sbOTAgQN+lyEiklLMrH2tj9UwjohIACjsRUQCQGEvIhIACnsRkQBQ2IuIBIDCXkQkABT2IiIBkBR99iLJ5sLgBL9p7adnZIqygmz21Zeza2MJZuZ3aSIRUdiLLHK+f5xHfnqcX57ofdPPfmtzKf/l3h3csrXSh8pEoqOwF/E89XoXf/g3h8nKND793m3805trqSsvYGB8il+e6OUvftXK/d94kX//nmY+/d5m7eVLSlHYiwB//WI7f/Tjo+xvrOB//Yu91JTkXf1ZbWk+H721gQ/t28R/e/IYX332DH2jU/zpB3cr8CVlKOwl8H56uJs/+vFR7txezdf/5T5yszKXfFxBThb/48O/Rag4l0d/1UpFYTb/6f3bE1ytSGQU9hJopy6N8h//7+u8paGcr/3u8kEfZmb8wftv5PL4NF97rpXdG0u556baBFUrEjm1XkpgzczN8+n/c4jC3Cwe/d195GWvHPRhZsYf37eLvfVl/METh+kenoxzpSLRU9hLYP3l35/nePcIf/KBXVQvGqNfi9ysTL7yO3uYnXN89okjOOfiVKVIbCjsJZDO94/zlV+e5u5dG7h7d2TDMA2VhXzmfdt4/nQff3fyza2aIslEYS+B9Mc/OUZOVgZ/fN+uqF7nwbc3sjVUyCM/PcH07HyMqhOJPYW9BM6BtkF+daqPh++44ZoWy0hkZ2bwX//JDs71j/PXL675okEiCaewl8D50i9OU1WUy8fe1hCT17vjxmretS3EV589w9jUbExeUyTWFPYSKL85288/nhvg4TuaKMiJTeexmfEf7trG8OQM33upIyavKRJrCnsJlK/88gy1pXk8sL8+pq+7p66MtzdV8pcvnGNqdi6mry0SC6uGvZnVmdlzZnbCzI6Z2Se95V8ws4tmdsj7unfRcz5nZmfN7JSZvT+eGyCyVkcvDvNy2yAff8eWNffUr8e/vb2JnpEpnnytK+avLRKttezZzwKfcc7tAG4FHjaznd7P/tw5t8f7ehrA+9n9wC7gbuDrZhb73yyRdfrf/9BGQU4mH2mpi8vrv+OGKnZvKuEvnm9lfl5995JcVg1751y3c+5V7/4ocALYtMJT7gO+75ybcs6dB84C+2NRrEik+sem+MnrXXz4LZspzc+OyzrMjH/zzq2c6x/nH1r747IOkUita8zezBqBvcBL3qJPmNlhM/uWmZV7yzYBFxY9rZOVPxxE4u67L3UwPTfPx97WGNf13L17A+UF2XzvZR2oleSy5rA3syLgCeBTzrkR4FGgCdgDdANfCj90iae/6W9aM3vIzA6Y2YG+vr51Fy6yVnPzju+81M47m6u4obooruvKzcrkw2/ZzC+O9dA3OhXXdYmsx5rC3syyWQj67zjnfgjgnOtxzs055+aBb/DGUE0nsHhQdDPwpiNWzrnHnHMtzrmWUCgUzTaIrOiFswuXF4x1B85y7t9fz+y8428OdiZkfSJrsZZuHAO+CZxwzn150fLFE4p8EDjq3X8KuN/Mcs1sC9AMvBy7kkXW54mDnZTmZ3PnjuqErK8pVMQtWyr4/isdOlArSWMte/a3AR8F3nNdm+WfmdkRMzsM3AF8GsA5dwz4AXAc+BnwsHNOjcfii5ErM/z82CX+2c21q85VH0sP7K+nfWCCV9oGE7ZOkZWsegqhc+4Flh6Hf3qF5zwCPBJFXSIx8fThbqZm5/nQvs0JXe/7dtWQn53JU6936QLlkhR0Bq2ktSde7WRrqJA9dWUJXW9BThZ37qjm6SPdzMxpNkzxn8Je0tbFoUleabvMh/Zt9uXC4L9980YuT8zwwln13Iv/FPaStn529BIA9/p0jdh33xiiJC+LnxzS9AniP4W9pK2/PdLN9g3FbKkq9GX9uVmZ3L17A7843sOVGfUoiL8U9pKWekaucLDjMvdEeMnBWPntmzcxNjXL86d14qD4S2Evaennxy7hHNx70wZf67hlawUleVk8c7zH1zpEFPaSlv72yCWaQoU01xT7Wkd2ZgZ3bK/m7072MqcTrMRHCntJO4Pj07x0fsD3IZywu3bWMDg+zcH2y36XIgGmsJe086tTvcy7hRObksG7t4XIzjSeOX7J71IkwBT2knaeO9VHVVEuuzeW+l0KAMV52dy6tZJnjvfgnIZyxB8Ke0krs3PzPH+ql9tvDJGRkfgTqZbzvp01tA1M0No35ncpElAKe0krr3YMMXJllvdsT8wMl2v13p0LQ0rPnuj1uRIJKoW9pJXnTvWSlWG8o7nK71KuUVuaz7aaIn59Rv324g+FvaSV50720tJYTklefK4zG413NYd45fxlJqd1Nq0knsJe0kbX0CQnL40m3RBO2Lu2hZiem+fF8wN+lyIBpLCXtPHCmYXZJd+9LTnDfv+WCnKzMvi1pk4QHyjsJW28cLafUHEu22rie1HxSOVlZ3LL1kqFvfhCYS9pwTnHb1r7ua2p0pe569fqXc1VtPaNc3Fo0u9SJGAU9pIWTvWM0j82zW03JFcXzvXevS0EoL17STiFvaSF8Hh9sof9DdVF1JTk8ptWHaSVxFLYS1r4TesAW6sK2ViW73cpKzIzbt1ayYvnBjR1giSUwl5S3szcPC+dG+DtN1T6XcqavG1rJX2jU7T2jftdigSIwl5S3usXhhifnuMdST6EE3br1oUPpRfPaShHEkdhLykvHJq3bEmNPfuGygJqS/P4R4W9JJDCXlLeS+cH2b6hmPLCHL9LWZPwuP1LGreXBFLYS0qbnZvn1fbL7N9S4Xcp6/K2rZX0j01ztldTHktiKOwlpR3vHmF8ei7lwl7j9pJoq4a9mdWZ2XNmdsLMjpnZJ73lFWb2jJmd8W7LveVmZl81s7NmdtjM9sV7IyS4Xj4/CMD+xtQK+7qKfDaW5vHiuUG/S5GAWMue/SzwGefcDuBW4GEz2wl8FnjWOdcMPOt9D3AP0Ox9PQQ8GvOqRTwvnR+ksbKA6pI8v0tZFzOjpbGCA+2DGreXhFg17J1z3c65V737o8AJYBNwH/C497DHgQ949+8Dvu0WvAiUmVltzCuXwJufd7zSNphyQzhhLY3l9IxM0XlZ8+RI/K1rzN7MGoG9wEtAjXOuGxY+EIDwvLKbgAuLntbpLROJqTO9YwxNzLA/RVour9fSsPAhdbD9ss+VSBCsOezNrAh4AviUc25kpYcusexNf6ea2UNmdsDMDvT1aVIoWb+X21JzvD7sxg3FFOdm8Uqbxu0l/tYU9maWzULQf8c590NvcU94eMa7DV9JuROoW/T0zUDX9a/pnHvMOdfinGsJhUKR1i8B9lr7ZULFudRVJPd8OMvJzDD21Jdpz14SYi3dOAZ8EzjhnPvyoh89BTzo3X8QeHLR8o95XTm3AsPh4R6RWHrtwhB768qSev761by1sYJTPaMMT874XYqkubXs2d8GfBR4j5kd8r7uBb4I3GVmZ4C7vO8BngbOAWeBbwD/LvZlS9BdHp/mfP84e+vL/S4lKi0N5TgHr3Zo717iK2u1BzjnXmDpcXiAO5d4vAMejrIukRW9dmEhHPfVl/lcSXT21JeRmWEcbLvMHTcm57VzJT3oDFpJSa91DJGZYdy0udTvUqJSkJPFztoSDrTrIK3El8JeUtJrHUNs31BMQc6qf5wmvb31ZRzpHGZuXidXSfwo7CXlzM07Dl0YYm+KD+GE7akrY3x6jjO9o36XImlMYS8pp7VvjLGpWfbWpfbB2bA9dQsfWoc6hnyuRNKZwl5Szmte50q67NlvqSqkND+bQxcU9hI/CntJOa91DFGan82WqkK/S4kJM+PmujKFvcSVwl5SzmsdC+P1qXwy1fX21pVxumeU8alZv0uRNKWwl5QyemWG072jaTNeH7anvox5B4c7h/0uRdKUwl5SyuHOYZxLn/H6sD2bvYO0GsqROFHYS0p51Zs07Oa69Ar78sIcGisLOHRB0yZIfCjsJaW83jlEU2iheyXd3FxXxusXNIwj8aGwl5Ry9OIIN21K7SkSlnPTplIujVyhb3TK71IkDSnsJWX0jU5xaeQKu9M07MPbdbRLe/cSewp7SRnhEEzXsN+5sQSAYxcV9hJ7CntJGUe9tsRdXiimm5K8hRPFjijsJQ4U9pIyjnYNs6WqkOK89Ds4G7Z7UylHL650iWeRyCjsJWUcvTiStkM4Ybs3lnBxaJLL49N+lyJpRmEvKWFwfJqLQ5PctCk9h3DCbtJBWokThb2khKPeOPbujem9Z7/L2z6N20usKewlJYTDb1eaD+OUFmRTX1Fw9cNNJFYU9pISjnUNU19RkJZnzl7vJh2klThQ2EtKOHJxOG3PnL3erk0ldAxOMDwx43cpkkYU9pL0hidmuDA4mfadOGHhD7VjOkgrMaSwl6T3xpmz6d2JE7ZbB2klDhT2kvSC0okTVl6Yw6ayfI52adxeYkdhL0nvyMVhNpXlU16Y43cpCbN7U4k6ciSmFPaS9I53jQRmCCds98ZSzvePM6Zr0kqMKOwlqU1Oz3F+YJwdtcEK+/D2nro06nMlki5WDXsz+5aZ9ZrZ0UXLvmBmF83skPd176Kffc7MzprZKTN7f7wKl2A43TOKc7B9Q7DCfnttMQAnL2ncXmJjLXv2fwXcvcTyP3fO7fG+ngYws53A/cAu7zlfN7PMWBUrwXOieyHsdgZsz35TWT7FeVlXt18kWquGvXPu18DgGl/vPuD7zrkp59x54CywP4r6JOBOXhqlMCeTzeX5fpeSUGbGjg0lnOzWMI7ERjRj9p8ws8PeME+5t2wTcGHRYzq9ZW9iZg+Z2QEzO9DX1xdFGZLOTnSPcOOGYjIyzO9SEm57bTEnL43inPO7FEkDkYb9o0ATsAfoBr7kLV/qN3LJ/6nOuceccy3OuZZQKBRhGZLOnHOc6B5he8CGcMK2byhhbGqWzsuTfpciaSCisHfO9Tjn5pxz88A3eGOophOoW/TQzUBXdCVKUHUPX2HkymzgOnHCdngHaTVuL7EQUdibWe2ibz8IhDt1ngLuN7NcM9sCNAMvR1eiBFW4E2XHhmKfK/HHtppizBaOW4hEK2u1B5jZ94DbgSoz6wQ+D9xuZntYGKJpA34fwDl3zMx+ABwHZoGHnXNz8Sld0t0J7+DktoCGfWFuFg0VBWq/lJhYNeydcw8ssfibKzz+EeCRaIoSgYXhi83l+ZSk8QXGV7N9Q8nVDz2RaOgMWklaJy+NBna8PmxHbQltA+NMTGvaBImOwl6S0pWZOc71jQV2vD5se20xzsHpnjG/S5EUp7CXpHS2d4x5R2DbLsN2eNNEnFRHjkRJYS9J6bgXbtsDvme/uTyfwpxMtV9K1BT2kpROdo+Sn51JQ2Wh36X4KiPD2F5bwgm1X0qUFPaSlE5eGmHbhmIyAzhNwvW2byjmZPeIpk2QqCjsJek45zh5aZTtNcEewgnbXlvCyJVZuoav+F2KpDCFvSSd/rFpBsenuTHg4/VhN3ofeqd7NJQjkVPYS9I544XaNu3ZA7Ctpgh4499FJBIKe0k6p6+GfZHPlSSHsoIcQsW56rWXqCjsJemc6R2jND+bUHGu36UkjW01RZzpVdhL5BT2knTO9IzRXF2EmTpxwpqriznbowuZSOQU9pJUnHOc7h2lWeP112iuKWJ8eo6LQ7qQiURGYS9JpW9siqGJGY3XXyd8sPqMxu0lQgp7SSrhMFMnzrW2Vav9UqKjsJekEm4vbNae/TVKC7KpLs7VQVqJmMJeksrpcCdOkTpxrtdcU6Ree4mYwl6SypmeUbbVqBNnKc3VxZzpHWN+Xh05sn4Ke0kazjlO94ypE2cZ22qKmVBHjkRIYS9Jo290iuHJGbZVa7x+KVenTejVUI6sn8Jekkb44KP27JfWXK32S4mcwl6Sxml14qwo3JGjOXIkEgp7SRqne8YoK1Anzkq21RRrGEciorCXpHGmZ5Rt1cXqxFnBQvulOnJk/RT2khScc5zpHeMGDeGsqLm6mMkZdeTI+insJSmoE2dt1JEjkVLYS1I426dOnLW4wfswbO0d97kSSTWrhr2ZfcvMes3s6KJlFWb2jJmd8W7LveVmZl81s7NmdtjM9sWzeEkfrX0L4dUU0p79SsoKcqgqyqG1Tx05sj5r2bP/K+Du65Z9FnjWOdcMPOt9D3AP0Ox9PQQ8GpsyJd219o5RmJNJTYk6cVazNVSksJd1WzXsnXO/BgavW3wf8Lh3/3HgA4uWf9steBEoM7PaWBUr6au1b4wmXZ1qTZpCRVf/EhJZq0jH7Gucc90A3m21t3wTcGHR4zq9ZW9iZg+Z2QEzO9DX1xdhGZIuzvWNawhnjZpChQyOTzM4Pu13KZJCYn2AdqndsiUbgp1zjznnWpxzLaFQKMZlSCqZmJ7l4tAkW6sK/S4lJYQ/FM9pKEfWIdKw7wkPz3i3vd7yTqBu0eM2A12RlydBcC58cFZtl2sSDnuN28t6RBr2TwEPevcfBJ5ctPxjXlfOrcBweLhHZDnh0NIwztpsKs8nJytD4/ayLlmrPcDMvgfcDlSZWSfweeCLwA/M7ONAB/AR7+FPA/cCZ4EJ4PfiULOkmda+cTIMGioL/C4lJWRmGFurCmnVJQplHVYNe+fcA8v86M4lHuuAh6MtSoKltW+MuooC8rIz/S4lZTSFijjWNex3GZJCdAat+K61d0xDOOvUFCqkY3CCqdk5v0uRFKGwF1/NzTvO94/TFFInzno0VRcx76B9YMLvUiRFKOzFV11Dk0zNzmvPfp2uduRo3F7WSGEvvgpPgKa2y/XZ4p2ToPZLWSuFvfgqvGeqPfv1KczNYmNpntovZc0U9uKr1r5xyguyqSjM8buUlNNUrQnRZO0U9uKr1j514kSqKVREa+8YCx3PIitT2IuvzinsI9YUKmR8eo6ekSm/S5EUoLAX3wxNTNM/Nk1TtdouI6E5cmQ9FPbiG12dKjrhDiaFvayFwl58E56id4umNo5IdXEuhTmZV2cNFVmJwl580zYwTmaGUVehCdAiYWY0VhVyvl9hL6tT2ItvzvePU19RQHam/htGaovCXtZIv2Xim3N94xrCidLWqkI6L08wPTvvdymS5BT24ov5eUf7wASNlQr7aDRWFTLvoGNQE6LJyhT24oue0StMzsyxRbNdRiX8l5GGcmQ1CnvxRTicdJHx6ITDvk1hL6tQ2IsvwmHfqLCPSllBDuUF2ZxT2MsqFPbii/N94+RmZVBbkud3KSlvoSNHJ1bJyhT24ou2gYVOnIwM87uUlLelqoi2fh2glZUp7MUX5/rH1YkTI1uqCrg0coXxqVm/S5EkprCXhJudm6djYEKdODGypWphjpy2AY3by/IU9pJwF4cmmZ13OqEqRt7oyNFQjixPYS8JF+4cUdjHRmPVwtxCOkgrK1HYS8K1KexjqiAniw0leZzXnr2sQGEvCXe+f5zivCwqdd3ZmFH7paxGYS8Jd75/oe3STG2XsaKpjmU1UYW9mbWZ2REzO2RmB7xlFWb2jJmd8W7LY1OqpItw2EvsbK0q5PLEDEMT036XIkkqFnv2dzjn9jjnWrzvPws865xrBp71vhcB4MrMHBeHJhX2MaYJ0WQ18RjGuQ943Lv/OPCBOKxDUlTH4ATO6eBsrDUq7GUV0Ya9A35hZgfN7CFvWY1zrhvAu62Och2SRs6rEycu6isKyDDNfinLy4ry+bc557rMrBp4xsxOrvWJ3ofDQwD19fVRliGpQrNdxkdOVgabyws0+6UsK6o9e+dcl3fbC/wI2A/0mFktgHfbu8xzH3POtTjnWkKhUDRlSAppHxinsjCHkrxsv0tJO41VhbQPqNdelhZx2JtZoZkVh+8D7wOOAk8BD3oPexB4MtoiJX20D0xQX1ngdxlpqaGigHbNjyPLiGYYpwb4kdcrnQV81zn3MzN7BfiBmX0c6AA+En2Zki7aByZ4a6O6ceOhobKAkSuzDE1MU1agE9bkWhGHvXPuHHDzEssHgDujKUrS09TsHF3DkzRUbva7lLTU4E0Z3TYwwR6FvVxHZ9BKwnRensS5hT1Qib3wv6uGcmQpCntJmHAINeiiJXFRXxEOex2klTdT2EvChENIe/bxkZedyYaSPIW9LElhLwnTPjBBYU6mZruMo/rKAjoGNYwjb6awl4RpHxinoVKzXcZTY2UBbdqzlyUo7CVh2gcnNIQTZw2VhfSNTjExrYuPy7UU9pIQc/OOC4MTOjgbZ2905GjvXq6lsJeE6B6eZGbOac8+zhoqFj5MFfZyPYW9JMTVTpwKhX081avXXpahsJeEuBr2mu0yrkrzsykvyKZ9UHv2ci2FvSRE++A4OZkZbCjJ87uUtFdfWUiHhnHkOgp7SYj2/gnqKvLJzFDbZbwttF9qGEeupbCXhGhXJ07CNFQU0DU0yfTsvN+lSBJR2EvcOedoHxi/OneLxFd9ZSHzDjovayhH3qCwl7jrH5tmYnqORrVdJkT431kHaWUxhb3EXXiuFg3jJEa4/VIHaWUxhb3EXVu/ZrtMpFBRLgU5mTpIK9dQ2EvctQ9OkGGwuVxhnwhmRn1Fgfbs5RoKe4m79oFxNpblk5Ol/26J0lhZqD17uYZ++yTu2gc022WiNVQWcOHyJPPzzu9SJEko7CXuOgYnqK/QwdlEqq8sYHp2nksjV/wuRZKEwl7iauTKDIPj02q7TLBGr/NJQzkSprCXuOrQdWd9ET6BTQdpJUxhL3EV3rNUj31ibSzLJzvTdIlCuUphL3EVntpYUyUkVmaGUVeui4/LGxT2ElftA+NUFeVSmJvldymBU19ZcPWENhGFvcTV6Z4xtoY0hOOHplAR5/rHmJnT7JcSx7A3s7vN7JSZnTWzz8ZrPZK8rszMcaxrmH315X6XEkh768u4MjPPie4Rv0uRJBCXsDezTOBrwD3ATuABM9sZj3VJ8jrcOczMnKOlQWHvh7d4/+4H2i77XIkkg3jt2e8HzjrnzjnnpoHvA/fFaV2SpA62L4TMPoW9L2pL89lUls/BDoW9QLyOmm0CLiz6vhO4JdYref50H3/y/47H+mUlRnpGrrC1qpCKwhy/SwmsfQ3l/PzYJe768vN+lyLL+J231vGv37k17uuJV9gvdaHRaybpMLOHgIcA6uvrI1pJUW4WzTVFET1X4q+5poi7d9f6XUag/d5tjcw7h3OaIydZVRXlJmQ9Fo//BGb2NuALzrn3e99/DsA599+XenxLS4s7cOBAzOsQEUlnZnbQOdeylsfGa8z+FaDZzLaYWQ5wP/BUnNYlIiKriMswjnNu1sw+AfwcyAS+5Zw7Fo91iYjI6uJ2WqNz7mng6Xi9voiIrJ3OoBURCQCFvYhIACjsRUQCQGEvIhIACnsRkQCIy0lV6y7CrA9oj/DpVUB/DMtJNUHe/iBvOwR7+7XtCxqcc6G1PCkpwj4aZnZgrWeQpaMgb3+Qtx2Cvf3a9vVvu4ZxREQCQGEvIhIA6RD2j/ldgM+CvP1B3nYI9vZr29cp5cfsRURkdemwZy8iIqtI6bAP8kXNzazNzI6Y2SEzS/uLAZjZt8ys18yOLlpWYWbPmNkZ7zYtr3+4zLZ/wcwueu//ITO7188a48XM6szsOTM7YWbHzOyT3vKgvPfLbf+63/+UHcbxLmp+GriLhcsevgI84JwLxHUKzawNaHHOBaLX2MzeBYwB33bO7faW/Rkw6Jz7ovdhX+6c+0M/64yHZbb9C8CYc+5/+llbvJlZLVDrnHvVzIqBg8AHgH9FMN775bb/n7PO9z+V9+x1UfMAcc79Ghi8bvF9wOPe/cdZ+CVIO8tseyA457qdc69690eBEyxc4zoo7/1y279uqRz2S13UPKJ/hBTlgF+Y2UHver5BVOOc64aFXwqg2ud6Eu0TZnbYG+ZJy2GMxcysEdgLvEQA3/vrth/W+f6nctivelHzNHebc24fcA/wsPenvgTHo0ATsAfoBr7kbznxZWZFwBPAp5xzI37Xk2hLbP+63/9UDvtOoG7R95uBLp9qSTjnXJd32wv8iIVhraDp8cY0w2ObvT7XkzDOuR7n3Jxzbh74Bmn8/ptZNgtB9x3n3A+9xYF575fa/kje/1QO+8Be1NzMCr2DNZhZIfA+4OjKz0pLTwEPevcfBJ70sZaECged54Ok6ftvZgZ8EzjhnPvyoh8F4r1fbvsjef9TthsHwGs3+gpvXNT8EZ9LSggz28rC3jwsXEf4u+m+7Wb2PeB2Fmb86wE+D/wY+AFQD3QAH3HOpd2BzGW2/XYW/oR3QBvw++Ex7HRiZu8A/h44Asx7i/8zC+PWQXjvl9v+B1jn+5/SYS8iImuTysM4IiKyRgp7EZEAUNiLiASAwl5EJAAU9iIiAaCwFxEJAIW9iEgAKOxFRALg/wP/2SgqoqJLygAAAABJRU5ErkJggg==\n", + "text/plain": [ + "[]" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD4CAYAAAANbUbJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXRc5X3/8fdXo32zbEmWvK8SxiRgQBhCgARSUmd1aNOGpc3WHOImNO2vK/2lTZqmS5b+Uk4aUtehtElLQkgawBQHEpITIAlgG2wDtvEmb7JsWYstWbIlWdL398cdhUFI1kgezZ3l8zpHB82dOzMfRne+fua5z30ec3dERCSz5YQdQEREpp6KvYhIFlCxFxHJAir2IiJZQMVeRCQL5Ib1wlVVVb5w4cKwXl5EJC09//zzbe5ePdHHhVbsFy5cyObNm8N6eRGRtGRmByfzOHXjiIhkARV7EZEsoGIvIpIFVOxFRLKAir2ISBZQsRcRyQJxFXszW2Vmu8xsr5ndOcr908zsETPbZmbbzewjiY8qIiKTNe44ezOLAHcDNwJNwCYzW+/uO2J2+ySww93fY2bVwC4zu8/d+6cktcgUOTs4xEtHOuk8fZYzZwc53T/ImbOD9Eb/W1maz6XzplNfU0puRF+MJX3Ec1HVSmCvuzcCmNn9wGogttg7UGZmBpQCHcBAgrOKTInO02f52e7jPLHzOD/bdZxTveMfusX5ES6eO43L5k/n0vnTuXzBdGaU5CchrcjkxFPs5wCHY243AVeO2OdrwHqgGSgDPuDuQyOfyMxuB24HmD9//mTyiiTEiZ5+frDlCE/saGHjgQ4Gh5zKknxWXVTLDctmUjutkOL8XIryIhTlBz+FuTk0n+xly+ETbDl0ki2HTrDuqUYGhpy8iPH+y+fyibcuZd6M4rD/90ReJ55ib6NsG7m81a8DW4EbgCXAj83saXfves2D3NcB6wAaGhq0RJYk3eCQ8+3nDvJPP9pN55mz1NeU8vHrFvO2C2tYMa+CSM5oh/ur5lcWM7+ymNUr5gDQe3aQl4908vDWZr676TDf29zEb142l09ev5T5lSr6kjriKfZNwLyY23MJWvCxPgJ8wYM1Dvea2X5gGbAxISlFEmDj/g4+u347O492cfWSSv763cu5cFb5eT1nYV6EhoUzaFg4g09cv4R/e7KRb288xPdfaOKmS+dwx/VLWVhVkqD/A5HJs/HWoDWzXGA38DbgCLAJuNXdt8fs869Ai7v/jZnVAC8Al7h721jP29DQ4JoITZLhWGcv//jDnTy8tZk5FUV8+l0X8o431BKcYkq8lq5e1j65j28/d4iBIeePb6znE29dMmWvJ9nFzJ5394YJPy6eBcfN7J3AXUAEuNfd/97M1gC4+1ozmw38JzCLoNvnC+7+3+d6ThV7mWruzjd/eYAvPb6LgSFnzXWL+f23LqUoP5KU1z/e1cvnH93JI9uaedcbZ/Hl37qY4vzQJpqVDDGlxX4qqNjLVHJ3/mHDTr7x9H6uv6Caz733DaH0obs7655q5AuPvcKFteWs++DlzJ2uvnyZvMkWew0UlowzMDjEn3//Rb7x9H4+fPVC/v1DV4R2stTM+PhblnDvh6/g8InTvPdrv+C5xvZQskh2U7GXjNJ7dpBP3PcC33u+if/za/V89j3LyRlnhE0yXH/BTB765JupKM7jtnue47+endT6EyKTpmIvGaO7b4CP/McmfrSjhc+99yL+8NfqUuqk6JLqUh765Ju5tq6Kv37oZf5hw86wI0kWUbGXjNDe3cet33iWjQc6uOsDK/jQ1QvDjjSq8sI87vnQFfzuVQtY91Qj9z2nFr4kh4YGSNo70dPPb//bMzSdOMM3Png5NyyrCTvSOUVyjL9570U0nTjNZx/ezqKqEq5eUhV2LMlwatlLWnN3/vR72zjccYZvfXRlyhf6YZEc46u3XMqiqhI+cd8LHGzvCTuSZDgVe0lr//7z/fzkleP833cu48rFlWHHmZCywjzu+VAwgu73vrmZrt6zISeSTKZiL2lr2+GTfPGxV3j78pqU7aMfz4LKEr5+22UcaOvhU9/ZwuCQpoySqaFiL2mpq/csd3znBWaWFfKl91+cUqNuJurqJVV8bvVF/GxXK/+oEToyRXSCVtKOu/OXP3iJ5pO9PPDxq6goTv955G+7cgF7Wrq55+f7qa8p47evmDf+g0QmQC17STvf2XiYR188yp+8vZ7LF8wIO07C/NW7LuSapVV8+qGXaGztDjuOZBgVe0krrxzr4nOPbOfauirWXLck7DgJlRvJ4Z8/sIKC3Ah//6i6cySxVOwlbZzuH+CT971AeVEeX/ntFSkxDUKiVZcV8Ac3LOUnrwRLJIokioq9pI2v/Gg3jW093PWBFVSXFYQdZ8p85M2LWFRVwuf/dwdnB1+3uqfIpKjYS1o41tnLt549yG9eNpc3L83sq03zc3P4q3ddyL7WHr71jKZTkMSIq9ib2Soz22Vme83szlHu/zMz2xr9ednMBs0sc86cSej+5ad7cHf+8G11YUdJihuWzeS6+mruemI37d19YceRDDBusTezCHA38A5gOXCLmS2P3cfdv+zuK9x9BfCXwJPu3jEVgSX7HO44zXc3HebmK+Yzb0Z2LPxhZnzm3Rdyun+Q//fj3WHHkQwQT8t+JbDX3RvdvR+4H1h9jv1vAb6TiHAiAHc9sYdIjnHHDUvDjpJUS2eW8cE3LeA7Gw+xvbkz7DiS5uIp9nOAwzG3m6LbXsfMioFVwP+cfzQR2Hu8mwe3NPHBNy2gprww7DhJ90dvq6eiKI+/fWQHYS0hKpkhnmI/2vi2sY669wC/GKsLx8xuN7PNZra5tbU13oySxf75id0U5UVY85bMGlMfr2nFefzpr1/Ac/s72PDSsbDjSBqLp9g3AbHXbs8FmsfY92bO0YXj7uvcvcHdG6qrq+NPKVlpe3Mnj754lI9es4jK0swdajmem6+Yz7LaMv5hw056zw6GHUfSVDzFfhNQZ2aLzCyfoKCvH7mTmU0D3gI8nNiIkq3++ce7KS/M5WPXLg47SqgiOcZn33MRR06e4dvPHQo7jqSpcYu9uw8AdwCPAzuBB9x9u5mtMbM1MbveBPzI3bUKg5y3Fw6d4Imdx/n4W5YwrSgv7Dihe9OSSi6dX8F9zx1U371MSlzj7N19g7vXu/sSd//76La17r42Zp//dPebpyqoZJev/Gg3lSX5fDhN56mfCreunM++1h6e269RzTJxuoJWUs4z+9r5+d42fv+tSygp0Czcw9598WzKC3PVlSOTomIvKeeuJ3ZTU17A71y1IOwoKaUoP8JvXDaXH758VFfVyoSp2EtKOdgedFN8+OpFFOZFwo6Tcm67cj5nB53vP98UdhRJMyr2klIe2tKMGaxeMTvsKCmprqaMlQtn8O2NhxjSerUyASr2kjLcnYe2HuGqRZXMrigKO07KuvXK+RxsP80v97WHHUXSiIq9pIxtTZ3sb+vhpktHnY1Dola9oZbpxXnc95ymP5b4qdhLynhoyxHyc3NY9cbasKOktMK8CO+/fC4/3tHC8a7esONImlCxl5RwdnCIR7Y1c+OFNZQX6iKq8dyycj4DQ84Dmw+Pv7MIKvaSIn6+p432nn7epy6cuCyuLuXqJZV8Z+NhBnWiVuKgYi8p4cEtR6gozuMt9ZogL163XjmfIyfP8NQezSAr41Oxl9B19w3wox3HePfFs8jP1SEZr7cvr6WqNJ/7ntUVtTI+fbIkdI+/fIzes0MahTNB+bk5/FbDPH76SgtHO8+EHUdSnIq9hO6hrUeYN6OIy+ZPDztK2rnlivk4cP9GnaiVc1Oxl1C1dPXyi71t3LRiDmajLYom5zK/spirl1Tyvy+OtZ6QSEDFXkL1yLZmhhxWqwtn0q6/YCb7WntoPqmuHBmbir2E6sEtR7hk7jSWVJeGHSVtXVsXjGB6WqNy5BziKvZmtsrMdpnZXjO7c4x93mpmW81su5k9mdiYkol2t5xie3OXxtafp/qaUmaWFfDUnrawo0gKG3dlCDOLAHcDNxIsPr7JzNa7+46YfSqArwOr3P2Qmc2cqsCSOR7acoRIjvHuizXD5fkwM66tq+Ynr7QwOOREcnTuQ14vnpb9SmCvuze6ez9wP7B6xD63Aj9w90MA7n48sTEl0wwNOQ9vbebauiqqywrCjpP2rquv4uTps7x8pDPsKJKi4in2c4DYcV1N0W2x6oHpZvYzM3vezD442hOZ2e1mttnMNre2qn8xm21rOsmRk2c0b32CXLO0ClC/vYwtnmI/2nfCkZNx5AKXA+8Cfh34azOrf92D3Ne5e4O7N1RX67L4bDa8aPY1S3UcJEJlaQFvmFPOU7vVby+ji6fYNwHzYm7PBUYO6m0CHnP3HndvA54CLklMRMlEm/Z3sLi6RF04CXRtXTUvHDrBqd6zYUeRFBRPsd8E1JnZIjPLB24G1o/Y52HgWjPLNbNi4EpgZ2KjSqYYHHI2HujgykUzwo6SUa6tq2JgyHm2sSPsKJKCxi327j4A3AE8TlDAH3D37Wa2xszWRPfZCTwGvAhsBO5x95enLraks13HTnGqd4ArFqrYJ9LlC6ZTlBdRv72MatyhlwDuvgHYMGLb2hG3vwx8OXHRJFNtOhC0PFeqZZ9QBbkR3rSkkqd2q9jL6+kKWkm6jfs7mFNRxNzpxWFHyTjX1lVxoP00h9pPhx1FUoyKvSSVu/Pc/g6uWKgZLqfCr6ZO2KvWvbyWir0k1YH207R197FyUWXYUTLSkuoSZk8r5GkNwZQRVOwlqTbubwdg5SK17KeCmXFdfTW/2NfGwOBQ2HEkhajYS1Jt3H+CGSX5muVyCl1bV82p3gG2NZ0MO4qkEBV7SaqNB9pZuXCGFiqZQm9eWokZuppWXkPFXpLmaOcZDnec4QoNuZxSFcX5XDy3QuPt5TVU7CVpNkbnw9GVs1Pvuroqth4+SecZTZ0gARV7SZqN+zsoLcjlwlnlYUfJeNfVVzPk8Mw+deVIQMVekmbTgQ4uXzBdi2skwYp5FZQW5PKk+u0lSsVekqKjp5/dLd2aIiFJ8iI5v5o6wX3kjOSSjVTsJSk0H07yvXlJJUdOnqG5szfsKJICVOwlKTbt7yA/N4eL504LO0rWuGhO8F7vOtYVchJJBSr2khQbD3SwYl4FBbmRsKNkjfqaMgBeOXYq5CSSClTsZcp19w2wvblLQy6TbFpRHrOnFfLKURV7UbGXJHjh4AkGh1z99SFYNqucXWrZC3EWezNbZWa7zGyvmd05yv1vNbNOM9sa/flM4qNKutq4v4NIjnHZfE1+lmwX1Jaxr7Wb/gFNipbtxl2pyswiwN3AjQQLi28ys/XuvmPErk+7+7unIKOkuY0HOnjD7HJKCuJaGE0SaFltGQNDTmNbN8tqdTFbNounZb8S2Ovuje7eD9wPrJ7aWJIp+gYG2Xr4pLpwQjJc4NVvL/EU+znA4ZjbTdFtI73JzLaZ2Q/N7KLRnsjMbjezzWa2ubVVkzRlgxebOukfGNLi4iFZXF1CXsQ0IkfiKvajXds+8pK8F4AF7n4J8C/AQ6M9kbuvc/cGd2+orq6eWFJJS8OTn6nYhyMvksOS6lKNtZe4in0TMC/m9lygOXYHd+9y9+7o7xuAPDOrSlhKSVu7jp1i7vQippfkhx0lay2rLVPLXuIq9puAOjNbZGb5wM3A+tgdzKzWoqtRmNnK6PO2JzqspJ/Gtm4Wa1WqUF1QW87Rzl46T2u642w2brF39wHgDuBxYCfwgLtvN7M1ZrYmutv7gZfNbBvwVeBm1+xLWc/d2d/aw+KqkrCjZLVls4IraXe1qHWfzeIaCxftmtkwYtvamN+/BnwtsdEk3bV09dHTP8iSahX7MC2rHZ42oUujorKYrqCVKdPY2g2gbpyQ1ZYXUl6Yq377LKdiL1NmX1sPEAz/k/CYmaZNEBV7mTqNrd0U50eoLS8MO0rWW1Zbxq5jp7SQSRZTsZcp09jaw6KqEqIDtSREF9SW0d03QNOJM2FHkZCo2MuU0bDL1DE8bYK6crKXir1Mid6zgzSdOKNhlynigpgROZKdVOxlShxsP427Ts6mitKCXOZOL9KInCymYi9TYl902OUSdeOkjGW1GpGTzVTsZUoMj7FfpG6clLGstozGth76BgbDjiIhULGXKdHY2kNteaEWLEkhF9SWMTjk7D3eHXYUCYGKvUyJfW096q9PMRcOz5GjrpyspGIvCefuNLZ2q9inmIWVJeTn5qjYZykVe0m4tu5+TvUO6ORsismN5LC0upSdKvZZScVeEk4ToKWuZbPKtGpVllKxl4RrHJ4ATSNxUs6y2jJauvo40dMfdhRJsriKvZmtMrNdZrbXzO48x35XmNmgmb0/cREl3TS2dlOQm8OciqKwo8gIF0SnTdDFVdln3GJvZhHgbuAdwHLgFjNbPsZ+XyRY0Uqy2PAEaDk5mgAt1VxYOzwiR1052Saelv1KYK+7N7p7P3A/sHqU/f4A+B/geALzSRpq1LDLlFVdVsD04jwtUZiF4in2c4DDMbebott+xczmADcBazkHM7vdzDab2ebW1taJZpU00D8wxKGO0yyu0snZVGRmXFBbxs6jKvbZJp5iP9p38ZErINwF/IW7n/M6bHdf5+4N7t5QXV0db0ZJI4c6TjM45GrZp7BlteXsbjnF0JAWMskm8VzL3gTMi7k9F2gesU8DcH90kYoq4J1mNuDuDyUkpaQNDbtMfctqyzjdH0xBPb+yOOw4kiTxtOw3AXVmtsjM8oGbgfWxO7j7Indf6O4Lge8Dn1Chz06NWnc25Q3Pbb9TJ2mzyrgte3cfMLM7CEbZRIB73X27ma2J3n/OfnrJLo2t3VSVFlBemBd2FBnD0pnBt6790X+YJTvENSWhu28ANozYNmqRd/cPn38sSVeNrRqJk+rKCvMoK8jlWGdv2FEkiXQFrSRUY1sPS1TsU96sikKOdmrx8WyiYi8Jc6Knn46efg27TAO104o4qpZ9VlGxl4RpbBseiaOWfaqbPa2Q5pMq9tlExV4SZl/r8EgctexTXe20Qtq6++gfGAo7iiSJir0kTGNrD3kRY950TYCW6mZPC/5GLV1q3WcLFXtJmMbWbubPKCY3osMq1c2qKARQv30W0adSEiYYiaMunHQwa9pwsdeInGyhYi8JMTA4xMH2HvXXp4naaDeOWvbZQ8VeEqLpxBnODmoCtHRRWpBLWWEuR0+qZZ8tVOwlIYaHXeqCqvQxW2Pts4qKvSRE4/CwS11QlTZqpxWq2GcRFXtJiH2tPUwvzmN6SX7YUSROsytU7LOJir0kRGNrt07Oppna8iLauvvoGzjnmkOSIVTsJSEOtp9mYaX669PJ8Fj7ls6+kJNIMqjYy3kbHHJau/uYHS0ekh401j67qNjLeWvv7mNwyJlZrmKfTmZprH1WiavYm9kqM9tlZnvN7M5R7l9tZi+a2VYz22xm1yQ+qqSqlq6gG6CmrCDkJDIRr7bsVeyzwbgrVZlZBLgbuJFg8fFNZrbe3XfE7PYTYL27u5ldDDwALJuKwJJ6hifTqlHLPq2UFORSXpirbpwsEU/LfiWw190b3b0fuB9YHbuDu3e7u0dvlgCOZI2WUyr26WqWLqzKGvEU+znA4ZjbTdFtr2FmN5nZK8CjwEdHeyIzuz3azbO5tbV1MnklBbV09pJjUFWqMfbpRssTZo94ir2Nsu11LXd3f9DdlwHvAz4/2hO5+zp3b3D3hurq6okllZTV0tVHVWmBpjZOQ7OmFWnh8SwRz6ezCZgXc3su0DzWzu7+FLDEzKrOM5ukiZZTverCSVOzphXS1t2vC6uyQDzFfhNQZ2aLzCwfuBlYH7uDmS01M4v+fhmQD7QnOqykppauPmrKNRInHQ2PyFHrPvONOxrH3QfM7A7gcSAC3Ovu281sTfT+tcBvAh80s7PAGeADMSdsJcMd7+rl0vkVYceQSYgda79AV0BntHGLPYC7bwA2jNi2Nub3LwJfTGw0SQd9A4O09/RTq26ctPTq8oQ6SZvpdEZNzkvrqegFVerGSUu6sCp7qNjLeRm+elZTJaSn4vxcphXlcfSkin2mU7GX83J8+OrZMhX7dDVLi5hkBRV7OS/HosW+dpqKfboKir367DOdir2cl5auPvIixvTivLCjyCTV6sKqrKBiL+fleFcvM8sKiV5mIWlo9rRC2nv66T2rC6symYq9nJfg6lmNxElnw11ww7OXSmZSsZfzElw9q/76dDa7IriwqlkjcjKair2cl5ZOzYuT7mq1PGFWULGXSevpG+BU34CKfZqbreUJs4KKvUzacV09mxGK8iNUFOepZZ/hVOxl0rQcYeaoLS/U8MsMp2Ivk/ZqsVfLPt3NrijSCdoMp2Ivk6aWfeaonVb4q6uhJTOp2MuktXT1UZwfobQgrpmyJYXNnlZIhy6symgq9jJpLV3BsEtdPZv+aqMjctRvn7niKvZmtsrMdpnZXjO7c5T7bzOzF6M/vzSzSxIfVVLN8a4+Zpapvz4TzI6OtW/WiJyMNW6xN7MIcDfwDmA5cIuZLR+x237gLe5+MfB5YF2ig0rqaTnVq9kuM0St1qLNePG07FcCe9290d37gfuB1bE7uPsv3f1E9OazwNzExpRU4+4c09WzGWOWLqzKePEU+znA4ZjbTdFtY/k94Iej3WFmt5vZZjPb3NraGn9KSTldZwboGxhSN06GGL6wqvmkunEyVTzFfrSzbz7qjmbXExT7vxjtfndf5+4N7t5QXV0df0pJOS2nNOwy08zSvPYZLZ4xc03AvJjbc4HmkTuZ2cXAPcA73L09MfEkVWmMfeaZNa2QZhX7jBVPy34TUGdmi8wsH7gZWB+7g5nNB34A/K677058TEk1wwuN16rYZ4xZ0wo5ptE4GWvclr27D5jZHcDjQAS41923m9ma6P1rgc8AlcDXo2OuB9y9YepiS9iGW/YzNVVCxpg1rZATp89ypn+QovxI2HEkweK69NHdNwAbRmxbG/P7x4CPJTaapLKWrl6mFeVRmKeikCmGR+Qc6+plUVVJyGkk0XQFrUxKcPWsWvWZZNbwIiYakZORVOxlUrQcYeaZVaGx9plMxV4mZXheHMkcs7Q8YUZTsZcJGxpyjp/qUzdOhinMizC9OE/DLzOUir1MWHtPP4NDrpZ9Bpo1rUh99hlKxV4m7FfDLstU7DPNvBlFHGw/HXYMmQIq9jJhx6NTJWjGy8xTX1PGgfYeLWKSgVTsZcKOdQZXz6rPPvPU1ZQx5LC/rSfsKJJgKvYyYS1dvZhBVamKfaaprykFYHfLqZCTSKKp2MuEHT/VS2VJAXkRHT6ZZlFVCZEcY09Ld9hRJMH0aZUJCy6oUqs+ExXkRlhQWayWfQZSsZcJa+nq1WyXGax+Zhl7jqtln2lU7GXCWrp6malin7Hqa0o5qBE5GUfFXibk7OAQbd396sbJYMMjchpbNSInk6jYy4S0nhoedqmWfaaqrykDYM9x9dtnkriKvZmtMrNdZrbXzO4c5f5lZvaMmfWZ2Z8mPqakileXI1TLPlMNj8jRSdrMMu7iJWYWAe4GbiRYj3aTma139x0xu3UAnwLeNyUpJWUML0eoln3mys/NYWFlMbs1/DKjxNOyXwnsdfdGd+8H7gdWx+7g7sfdfRNwdgoySgrRQuPZob6mjD1q2WeUeIr9HOBwzO2m6LYJM7PbzWyzmW1ubW2dzFNIyFq6esnNMWYU54cdRaZQXU0ZBztOa0ROBomn2Nso23wyL+bu69y9wd0bqqurJ/MUErKWrj5mlhWQkzPaYSGZor6mFHfYq/H2GSOeYt8EzIu5PRdonpo4kuqOn+qlRrNdZrzhETkq9pkjnmK/Cagzs0Vmlg/cDKyf2liSqlq6eqnRPPYZb2FlCbkakZNRxh2N4+4DZnYH8DgQAe519+1mtiZ6/1ozqwU2A+XAkJn9EbDc3bumMLuE4FhnL29aXBl2DJli+bk5LKoq0YicDDJusQdw9w3AhhHb1sb8foyge0cy2Jn+Qbp6BzRVQpaorynj5ebOsGNIgugKWonb8ApVGnaZHZbOLOVQx2nO9GtETiZQsZe4bW8OeuXmzygOOYkkQ31NGe6wr1VdOZlAxV7i9uCWI1SXFXD5gulhR5Ek0KpVmUXFXuJyoqefn+06zupLZhPRGPussLCqhLyIaW77DKFiL3F59KWjnB10brpsUhdPSxrKiwQjcjRtQmZQsZe4PLjlCPU1pSyfVR52FEmiupoyDb/MECr2Mq5D7ad5/uAJbrp0Lmbqwskm9TPLOHxCI3IygYq9jOuhrUcAWL1idshJJNnqNEdOxlCxl3Nydx7ccoSrFs9gdkVR2HEkyTQiJ3Oo2Ms5bWvqZH9bD79xqS6QzkYLKoMRObu1RGHaU7GXc3rwhSYKcnNY9cbasKNICPIiOSyuKmWvTtKmPRV7GdPZwSEeefEov7a8hvLCvLDjSEjqakrVss8AKvYypqf3tNLR089NKzS2PpvV15RxuOMMp/sHwo4i50HFXsb0gxeOML04j+vqtapYNhs+SasROelNxV5Gdar3LD/e0cJ7LplNfq4Ok2xWF121ShdXpTd9imVUj718jL6BId53qbpwst2CGcXkR3I0bUKai6vYm9kqM9tlZnvN7M5R7jcz+2r0/hfN7LLER5VkenDLERZWFnPpvIqwo0jIciM5LK4u0YRoaW7cYm9mEeBu4B3AcuAWM1s+Yrd3AHXRn9uBf01wTkmio51neKaxnfddOkfTIwgwPEeOWvbpLJ5lCVcCe929EcDM7gdWAzti9lkNfMvdHXjWzCrMbJa7H0104Cd3t/J3/7tj/B1l0rr7BnCH92kUjkTVzyzlkW3N3PiVJ8OOkhE+cMU8Pnbt4qS+ZjzFfg5wOOZ2E3BlHPvMAV5T7M3sdoKWP/Pnz59oVgBKC3Kpi44OkKlzS205C6tKwo4hKeK9K2az53g3A0NDYUfJCFWlBUl/zXiK/Wjf430S++Du64B1AA0NDa+7Px6XL5jO5Qsun8xDRWSSFlSW8NVbLg07hpyHeE7QNgHzYm7PBZonsY+IiIQknmK/Cagzs0Vmlg/cDGWG+7QAAAbdSURBVKwfsc964IPRUTlXAZ1T0V8vIiKTM243jrsPmNkdwONABLjX3beb2Zro/WuBDcA7gb3AaeAjUxdZREQmKp4+e9x9A0FBj922NuZ3Bz6Z2GgiIpIouoJWRCQLqNiLiGQBFXsRkSygYi8ikgUsOLcawgubtQIHJ/nwKqAtgXESJVVzQepmU66JUa6JycRcC9x9wotMhFbsz4eZbXb3hrBzjJSquSB1synXxCjXxCjXq9SNIyKSBVTsRUSyQLoW+3VhBxhDquaC1M2mXBOjXBOjXFFp2WcvIiITk64texERmQAVexGRbODuU/oDrAJ2EcyIeeeI+/4get924EtjPP7LwCvAi8CDQEV0+23A1pifIWDFKI9fBDwH7AG+C+TH5DoB9BOsqHVZiuT6G6AX6CNYJ+CaFMm1CjgUzdUCPJkiud4PnIrmOgy8Icm57iA4th2oitn+xZi/40HgkhTJ9efAYDRXM/CZFMn1G0B3zPH1kSTnui/6Gi8D9wJ50e0fA85E8z6awPqVB3wTeAnYCfzlGI8f67g34KvR9/JFYurXWD9TXegjwD5gMZAPbAOWR++7HngCKIjenjnGc7wdyI35AH1xlH3eCDSO8fgHgJujv68Ffj+a6yjws2iuPcCLKZKrMeb92jX8+BTItT/6Pi2J/h2vSZFcJ4G7ou/XK8AzSc51KbAQOEC0eEVzNQErorkaQzi+xsp1BPgp4X0ex8rVHv275hMU3JPR35OV650EBdSA78QcXwcI1tj+R4J/HBPyfgG3AvdHfy+Ovs7CeI77mLw/jOa9CnhutNeP/ZnqbpxfLVbu7v3A8GLlELyZX3D3PgB3Pz7aE7j7j9x9IHrzWYJVsEa6heAP9BpmZsANwPejm74JvC+aawD4t2iue4FZZjYrBXLtiXm/fgoML7gbdq4zwHfdfR/B3/HaFMl1Fvh+9P36JlBvZjXJyBV9/BZ3PzBi80pgu7tvjeb6L2BB9L6wcx0BTofxeRwnVztBq/sswXTqgwSf0WTl2uBRwMbo41cCu9z9YYJvHFtJ3PvlQImZ5QJFBD0MXbGPPcdxTzTHt6KRnwUqovVrTFNd7MdaiBygHrjWzJ4zsyfN7Io4nu+jBP+ajfQBRv8jVgInY97s4defQ/BmH47Z3hvdHnauw2Z2k5m9Avwu8FR0n7BzDQHTzexnBIvG35AiuY4QdAFA0BKcTvCBSkausYw87uuimUmBXK3Am8xsG8Gqc29MkVzPABcStJ4/BfzS3YeSncvM8gg+d4/x+r/jSRJXv74P9BD0MBwC/sndO0bsP9ZxzyjZYu8b1VQX+3MtRJ5L8MG8Cvgz4IHov2SjP5HZpwn+pb9vxPYrCVoqL0/g9cd6HU+FXO7+oLsvI+ieuDJ6f9i5DLgceBfwJeAKM6tPgVzbCP4R2grcSFDMBpKUa8yni3n89cBbgF9GN4Wdq41gbpVLgB8TdAekQq45BC3n2cCngevMrDyEXF8HnnL3p5na+rWS4NvLbIJ++T8xs8UjH3aO1z/XfaOa6mJ/roXIm4AfRL+GbCRoOVaZ2X+Y2VYz+9XKWGb2IeDdwG3Rr1mxbmbsf63bCL7eDK/INfz6TQRv1ryY7YUx94WZK/b96iH4qleVArnygMfcvQeoIOi/vyQFctW6+0fcfQXwKEHrfn+Sco2lCZhnZhcD9xB0lzTG3Bdmrlp3747ebgc8icfXuXJdMvz6BH/DNmBZMnOZ2WeBauCPY3LFfh4rSFz9upXg83Q22gX0C2DkXDljHfejZYu9b3Q+Tqf++fwQ/OvXSPAv1/AJoYui960B/jb6ez3BVxIb5TlWATuA6lHuy4n+Ty8+R4bv8doTHJ+I5hp5gvalFMl1KOb92k0wMsFSINdhggOyiODs/x7gDSmQ60D0dfOjr/1wMv+OMfse4NUTjrkEI3AOANcRwnE/Tq7Y4+toMo+vc+Tq4tUT7dsJjvuqJB5fHyP49lUUsy22fn2eoJgm5O8I/AXwH9H3viS6z8XxHPfR39/Fa0/Qbhz3PR9vh/P9IfiauJtgVM6nY7bnA/9NcOb9BeCGMR6/N/pGDg+dWhtz31uBZ8d5/cUEJ1z2Rt+4gphcJwlOCB0DGlIk170EJ4P6CAr/NSmS650EXST9BB/EP0qRXH8czdRPUCSmJznXpwgKyQBBMbgnuv0xXh3ieBTYnCK5vj7i+Lo6RXL9DsE32eGhl7+T5FwDBDVq+PGfiW6/jaBGDBEMUmgCys83F8HAi+8RHLM7gD+b4HFvwN3RzC8RrV/n+tF0CSIiWUBX0IqIZAEVexGRLKBiLyKSBVTsRUSygIq9iEgWULEXEckCKvYiIlng/wM4iquA+vh4GAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] @@ -104,7 +230,9 @@ } ], "source": [ - "irr=DailyIrradiance(lat,day,True)" + "output = perth.daily_irradiance('2020-06-27')\n", + "time, irradiance = list(zip(*output))\n", + "plt.plot(time,irradiance)" ] }, { @@ -116,6 +244,9 @@ } ], "metadata": { + "jupytext": { + "formats": "ipynb,md" + }, "kernelspec": { "display_name": "Python 3", "language": "python", @@ -131,7 +262,49 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.3" + "version": "3.8.3" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": false + }, + "varInspector": { + "cols": { + "lenName": 16, + "lenType": 16, + "lenVar": 40 + }, + "kernels_config": { + "python": { + "delete_cmd_postfix": "", + "delete_cmd_prefix": "del ", + "library": "var_list.py", + "varRefreshCmd": "print(var_dic_list())" + }, + "r": { + "delete_cmd_postfix": ") ", + "delete_cmd_prefix": "rm(", + "library": "var_list.r", + "varRefreshCmd": "cat(var_dic_list()) " + } + }, + "types_to_exclude": [ + "module", + "function", + "builtin_function_or_method", + "instance", + "_Feature" + ], + "window_display": false } }, "nbformat": 4, diff --git a/Examples_files/Examples_14_1.png b/Examples_files/Examples_14_1.png new file mode 100644 index 0000000..f002297 Binary files /dev/null and b/Examples_files/Examples_14_1.png differ diff --git a/LICENSE.txt b/LICENSE.txt new file mode 100644 index 0000000..3754756 --- /dev/null +++ b/LICENSE.txt @@ -0,0 +1,17 @@ +MIT License +Copyright (c) 2018 YOUR NAME +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. \ No newline at end of file diff --git a/README.md b/README.md index dd3a92e..b9de1b5 100644 --- a/README.md +++ b/README.md @@ -1 +1,85 @@ -# SolarLib +# Solarlib + + +```python +from solarlib.location import Location +import datetime +import pytz +``` + +## Define a location + + +```python +latitude = -31.9505 +longitude = 115.8605 +timezone = 'Australia/Perth' +perth = Location(latitude,longitude,timezone) +``` + +## Sunrise +The output is a python datetime object. + + +```python +print(perth.sunrise('2020-06-27')) +``` + + 2020-06-27 07:17:33.495057+08:00 + + +## Sunset +The output is a python datetime object. + + +```python +print(perth.sunset('2020-06-27')) +``` + + 2020-06-27 17:21:50.791503+08:00 + + +## Day length +The output is a python timedelta object. + + +```python +print(perth.day_length('2020-06-27')) +``` + + 10:04:17.296446 + + +## Solar irradiance +The output is irradiance in $kW/m^2$. + + +```python +perth.solar_irradiance('2020-06-27 14:27:00') +``` + + + + + 0.7360521598862617 + + + +## Irradiance throughout a day + + + +```python +import matplotlib.pyplot as plt +``` + + +```python +output = perth.daily_irradiance('2020-06-27') +time, irradiance = list(zip(*output)) +plt.plot(time,irradiance) +``` + +![png](Examples_files/Examples_14_1.png) + + diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..af44f19 --- /dev/null +++ b/requirements.txt @@ -0,0 +1 @@ +pytz diff --git a/setup.cfg b/setup.cfg new file mode 100644 index 0000000..224a779 --- /dev/null +++ b/setup.cfg @@ -0,0 +1,2 @@ +[metadata] +description-file = README.md \ No newline at end of file diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..0f96352 --- /dev/null +++ b/setup.py @@ -0,0 +1,12 @@ +from setuptools import find_packages, setup + +setup( + name='solarlib', + packages=find_packages(), + version='0.1.0', + description='A library for calculating solar irradiance, sunrise, sunset, and more. ', + author='Pooya Darvehei', + url='https://github.com/pooyad359/solarlib', + license='MIT', + install_requires=['pytz'] +) diff --git a/solarlib/__init__.py b/solarlib/__init__.py index 9684c9c..8b13789 100644 --- a/solarlib/__init__.py +++ b/solarlib/__init__.py @@ -1,2 +1 @@ -from .solarlib import * diff --git a/solarlib/calendar.py b/solarlib/calendar.py new file mode 100644 index 0000000..7b0fad3 --- /dev/null +++ b/solarlib/calendar.py @@ -0,0 +1,23 @@ +import datetime +import pytz + +JULIAN_DAYS_INIT = 1721424.5 +CENTURY21 = 2451545 +CENTURY_DAYS = 36525 + +def time2seconds(time): + hour = time.hour + minute = time.minute + seconds = time.second + return hour * 3600 + minute * 60 + seconds + + +def julian_day(time): + utc = pytz.timezone('UTC') + utc_time = time.astimezone(utc) + time_of_day = time2seconds(utc_time)/3600/24 + return JULIAN_DAYS_INIT + utc_time.toordinal() + time_of_day + +def julian_century(time): + return (julian_day(time)-CENTURY21)/CENTURY_DAYS + \ No newline at end of file diff --git a/solarlib/error_handling.py b/solarlib/error_handling.py new file mode 100644 index 0000000..6fe3ef5 --- /dev/null +++ b/solarlib/error_handling.py @@ -0,0 +1,14 @@ +class ValueOutOfRange(Exception): + def __init__(self, name, rng, value): + self.message = f'{name} = {value} ' + self.message += 'but is expected to be in {rng} range.' + super().__init__(self.message) + + +def type_check(variable, types_list): + checks = [isinstance(variable, o) for o in types_list] + return any(checks) + + +def raise_type_error(variable_name, input_type): + raise TypeError(f'Unrecognised type {input_type} for {variable_name}') diff --git a/solarlib/location.py b/solarlib/location.py new file mode 100644 index 0000000..7cf3bab --- /dev/null +++ b/solarlib/location.py @@ -0,0 +1,140 @@ +from solarlib.error_handling import ( + ValueOutOfRange, + raise_type_error, + type_check) +from solarlib.solar import ( + sunrise, + sunset, + estimated_irradiance, + solar_noon + ) +import pytz +import datetime + + +class Location: + def __init__(self, latitude, longitude, timezone=None): + self.set_timezone(timezone) + self.set_latitude(latitude) + self.set_longitude(longitude) + + def set_latitude(self, latitude): + if not type_check(latitude, [int, float]): + raise_type_error('latitude', type(latitude)) + elif -90 <= latitude <= 90: + self.__latitude = latitude + else: + raise ValueError( + 'latitude must be between -90 and 90 degrees ' + + f'but recieved {latitude}' + ) + + def set_longitude(self, longitude): + if not type_check(longitude, [int, float]): + raise_type_error('longitude', type(longitude)) + elif -180 <= longitude <= 180: + self.__longitude = longitude + else: + raise ValueOutOfRange( + 'longitude must be between -180 and 180 ' + + f'degrees but recieved {longitude}' + ) + + def set_timezone(self, timezone): + if timezone is None: + self.__timezone = pytz.timezone('UTC') + elif isinstance(timezone, str): + self.__timezone = pytz.timezone(timezone) + elif isinstance(timezone, float) or isinstance(timezone, int): + timezone_hours = round(timezone*2)/2 + timedelta = datetime.timedelta(hours=timezone_hours) + self.__timezone = datetime.timezone(timedelta) + else: + raise_type_error('timezone', type(timezone)) + + @property + def latitude(self): + return self.__latitude + + @property + def longitude(self): + return self.__longitude + + @property + def timezone(self): + return self.__timezone + + def __repr__(self): + lat = self.__latitude + lon = self.__longitude + tz = self.__timezone + output = f'''Location(latitude={lat}, + longitude={lon}, + timezone={str(tz)})''' + return output.replace(' '*4, '').replace('\n', '') + + def __str__(self): + return repr(self) + + def sunrise(self, date, fmt='%Y-%m-%d'): + date = self.__parse_date(date, fmt='%Y-%m-%d') + lat = self.latitude + lon = self.longitude + base = sunrise(date, lat, lon) + return sunrise(base, lat, lon) + + def sunset(self, date, fmt='%Y-%m-%d'): + date = self.__parse_date(date, fmt=fmt) + lat = self.latitude + lon = self.longitude + base = sunset(date, lat, lon) + return sunset(base, lat, lon) + + def solar_irradiance(self, time, fmt='%Y-%m-%d %H:%M:%S'): + time = self.__parse_time(time, fmt=fmt) + lat = self.latitude + lon = self.longitude + return estimated_irradiance(time, lat, lon) + + def daily_irradiance(self, date, fmt='%Y-%m-%d', freq_min=30): + date = self.__parse_date(date, fmt=fmt) + delta = datetime.timedelta(minutes=freq_min) + n = 1440//freq_min + 1 + output = [ + ( + date+i*delta, + self.solar_irradiance(date+i*delta) + ) for i in range(n) + ] + return output + + def day_length(self, date, fmt='%Y-%m-%d'): + date = self.__parse_date(date, fmt=fmt) + rise_time = self.sunrise(date, fmt) + set_time = self.sunset(date, fmt) + return set_time-rise_time + + def solar_noon(self, date, fmt='%Y-%m-%d'): + date = self.__parse_date(date, fmt=fmt) + lon = self.longitude + base = solar_noon(date, lon) + return solar_noon(base, lon) + + def __parse_date(self, date, fmt='%Y-%m-%d'): + if isinstance(date, str): + return datetime.datetime.strptime(date, fmt) + elif isinstance(date, datetime.datetime): + return date + elif isinstance(date, datetime.date): + time = datetime.time() + return datetime.datetime.combine(date, time) + else: + raise_type_error('date', datetime.datetime) + + def __parse_time(self, date, fmt='%Y-%m-%d %H:%M:%S'): + if isinstance(date, str): + return datetime.datetime.strptime(date, fmt) + elif isinstance(date, datetime.datetime): + return date + else: + raise_type_error('date', datetime.datetime) diff --git a/solarlib/solar.py b/solarlib/solar.py new file mode 100644 index 0000000..6b76324 --- /dev/null +++ b/solarlib/solar.py @@ -0,0 +1,270 @@ +import datetime +import pytz +from math import degrees, radians, atan2, asin, acos, sin, cos, tan +from solarlib.calendar import julian_century +from solarlib.utils import get_minutes + + +def geom_mean_long(time): + ''' + Sun mean longitude (degrees) + ''' + julian_cent = julian_century(time) + angle = 280.46646 + julian_cent * (julian_cent * 3.032e-4 + 36000.76983) + return angle % 360 + + +def geom_mean_anom(time): + ''' + Sun Mean Anomaly (degrees) + ''' + julian_cent = julian_century(time) + angle = 357.52911 + julian_cent * (35999.05029 - julian_cent * 1.537e-4) + return angle + + +def eccent_earth_orbit(time): + ''' + Eccentricity of Earth orbit + ''' + julian_cent = julian_century(time) + return 0.016708634 - julian_cent * (4.2037e-5 + 1.267e-7 * julian_cent) + + +def equation_of_center(time): + ''' + Sun Equation of Center + ''' + julian_cent = julian_century(time) + anom = geom_mean_anom(time) + anom = radians(anom) + t1 = sin(anom) * (1.914602 - julian_cent * (0.004817 + 1.4e-5 * julian_cent)) + t2 = sin(2 * anom) * (0.019993 - 0.000101 * julian_cent) + t3 = sin(3 * anom) * 0.000289 + return t1 + t2 + t3 + + +def true_longitude(time): + return geom_mean_long(time) + equation_of_center(time) + + +def true_anomaly(time): + return geom_mean_anom(time) + equation_of_center(time) + + +def rad_vector(time): + ''' + in Astronomical Unit (AUs) + ''' + ecc = eccent_earth_orbit(time) + anom = true_anomaly(time) + return (1.000001018 * (1 - ecc * ecc)) / (1 + ecc * cos(radians(anom))) + + +def app_long(time): + julian_cent = julian_century(time) + tlon = true_longitude(time) + return tlon - 0.00569 - 0.00478 * sin( + radians(125.04 - 1934.136 * julian_cent)) + + +def mean_obliq_ecliptic(time): + julian_cent = julian_century(time) + res = 46.815 + julian_cent * (0.00059 - julian_cent * 0.001813) + res = 21.448 - julian_cent * (res) + res = 23 + (26 + res / 60) / 60 + return res + + +def obeliq_corr(time): + julian_cent = julian_century(time) + moe = mean_obliq_ecliptic(time) + return moe + 0.00256 * cos(radians(125.04 - 1934.136 * julian_cent)) + + +def solar_ascention(time): + lon = app_long(time) + oblq = obeliq_corr(time) + a = cos(radians(oblq)) * sin(radians(lon)) + b = cos(radians(lon)) + ascn = atan2(a, b) + return degrees(ascn) + + +def solar_declination(time): + lon = app_long(time) + oblq = obeliq_corr(time) + decn = asin(sin(radians(oblq)) * sin(radians(lon))) + return degrees(decn) + + +def equation_of_time(time): + oblq = obeliq_corr(time) + lon = geom_mean_long(time) + anom = geom_mean_anom(time) + ecc = eccent_earth_orbit(time) + var_y = tan(radians(oblq / 2)) * tan(radians(oblq / 2)) + result = [ + var_y * sin(2 * radians(lon)), -2 * ecc * sin(radians(anom)), + 4 * ecc * var_y * sin(radians(anom)) * cos(2 * radians(lon)), + -0.5 * var_y * var_y * sin(4 * radians(lon)), + -1.25 * ecc * ecc * sin(2 * radians(anom)) + ] + return 4 * degrees(sum(result)) + + +def sunrise_hour_angle(time, latitude): + decn = radians(solar_declination(time)) + a1 = cos(radians(90.833)) / (cos(radians(latitude)) * cos(decn)) + a2 = tan(radians(latitude)) * tan(decn) + return degrees(acos(a1 - a2)) + + +def solar_noon(time, longitude): + ''' + Solar noon time + ''' + tz = time.tzinfo + utc = pytz.timezone('UTC') + julian_cent = julian_century(time) + eot = equation_of_time(time) + noon_offset = (720 - 4 * longitude - eot) / 1440 + base_time = datetime.time() + midnight = datetime.datetime.combine(time.date(), base_time) + midnight = utc.localize(midnight) + noon_utc = midnight + datetime.timedelta(days=noon_offset) + return noon_utc.astimezone(tz) + + +def sunrise(time, latitude, longitude): + ''' + Sunrise time + ''' + noon = solar_noon(time, longitude) + julian_cent = julian_century(time) + hour_angle = sunrise_hour_angle(time, latitude) + delta = datetime.timedelta(minutes=4 * hour_angle) + return noon - delta + + +def sunset(time, latitude, longitude): + ''' + Sunset time + ''' + noon = solar_noon(time, longitude) + julian_cent = julian_century(time) + hour_angle = sunrise_hour_angle(time, latitude) + delta = datetime.timedelta(minutes=4 * hour_angle) + return noon + delta + + +def daytime_length(time, latitude): + ''' + Length of day time (Sunlight duration) in hours. + ''' + julian_cent = julian_century(time) + hour_angle = sunrise_hour_angle(time, latitude) + return 8 * hour_angle / 60 + + +def true_solar_time(time, longitude): + ''' + True solar time in minutes from midnight + ''' + utc = pytz.timezone('UTC') + julian_cent = julian_century(time) + eot = equation_of_time(time) + delta = datetime.timedelta(minutes=eot + longitude * 4) + true_time = (time + delta).astimezone(utc) + hour = true_time.hour + minute = true_time.minute + second = true_time.second + return get_minutes(true_time) + + +def hour_angle(time, longitude): + solar_time = true_solar_time(time, longitude) + return solar_time / 4 - 180 + + + +def solar_zenith(time, latitude, longitude): + ''' + Solar Zenith Angle in degrees + ''' + julian_cent = julian_century(time) + decn = solar_declination(time) + h_angle = hour_angle(time, longitude) + a1 = sin(radians(latitude)) * sin(radians(decn)) + a2 = cos(radians(latitude)) * cos(radians(decn)) * cos(radians(h_angle)) + return degrees(acos(a1 + a2)) + + +def solar_elevation_angle(time, latitude, longitude): + ''' + Solar Elevation Angle in degrees + ''' + return 90 - solar_zenith(time, latitude, longitude) + + +def atmospheric_refraction(time, latitude, longitude): + ''' + Approximate atmospheric refraction in degrees + ''' + elev = solar_elevation_angle(time, latitude, longitude) + if elev > 85: + refraction = 0 + elif elev > 5: + refraction = 58.1 / tan(radians(elev)) - 0.07 / pow( + tan(radians(elev)), 3) + 0.000086 / pow(tan(radians(elev)), 5) + elif elev > -.575: + refraction = 1735 + elev * (-518.2 + elev * (103.4 + elev * + (-12.79 + elev * 0.711))) + else: + refraction = -20.772 / tan(radians(elev)) + + return refraction / 3600 + + +def corrected_solar_elevation(time, latitude, longitude): + ''' + Solar Elevation corrected for atmospheric refraction + ''' + elev = solar_elevation_angle(time, latitude, longitude) + refr = atmospheric_refraction(time, latitude, longitude) + return elev + refr + + +def solar_azimuth(time, latitude, longitude): + ''' + Azimuth angle (measured Clockwise from North) + ''' + julian_cent = julian_century(time) + zenith = radians(solar_zenith(time, latitude, longitude)) + decln = radians(solar_declination(time)) + h_angle = hour_angle(time, longitude) + lat = radians(latitude) + if h_angle > 0: + a1 = sin(lat) * cos(zenith) - sin(decln) + a2 = cos(lat) * sin(zenith) + azimuth = degrees(acos( a1 / a2)) + 180 + return azimuth%360 + + else: + a1 = sin(lat) * cos(zenith) - sin(decln) + a2 = cos(lat) * sin(zenith) + azimuth = 540 - degrees(acos(a1 /a2)) + return azimuth%360 + +def estimated_irradiance(time, latitude, longitude): + ''' + Estimated solar irradiance in kW/m^2 + ''' + zenith = solar_zenith(time, latitude, longitude) + zenith = min(zenith,92) + if zenith<92: + amf = 1/cos(radians(zenith)) # Air Mass Factor + amf = 1/(cos(radians(zenith))+0.50572*(96.07995-zenith)**(-1.6364)) + return 1.353*0.7**(amf**0.678) + else: + return 0 diff --git a/solarlib/solarlib.py b/solarlib/solarlib.py deleted file mode 100644 index fedaaa7..0000000 --- a/solarlib/solarlib.py +++ /dev/null @@ -1,110 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -import datetime -deg2rad=np.pi/180 -rad2deg=180/np.pi -def SolarDeclination(DayOfYear): - ''' - calculate solar declination angle (deg) - DayOfYear=0 ---> January 1st - ''' - sin_delta=np.sin(-23.44*deg2rad)*np.cos(2*np.pi/365.24*(DayOfYear+2)+2*0.0167*(2*np.pi/365.24*(DayOfYear-2))) - delta=np.arcsin(sin_delta) - return delta -def ZenithAngle(Latitude,DayOfYear,HourOfDay): - ''' - calculate solar zenith angle (deg) - Latitude: (+) North, (-) South - DayOfYear: 0 for January 1st - HourOfDay: hours since midnight - ''' - delta=SolarDeclination(DayOfYear) - HourAngle=15*(HourOfDay-12)*deg2rad - phi=Latitude*deg2rad - cosz=np.sin(phi)*np.sin(delta)+np.cos(phi)*np.cos(delta)*np.cos(HourAngle) - zenith=np.arccos(cosz)*rad2deg - return zenith -def SolarIrradiance(Latitude,DayOfYear,HourOfDay): - ''' - calculate solar irradiance (kW) - Latitude: (+) North, (-) South - DayOfYear: 0 for January 1st - HourOfDay: hours since midnight - ''' - z=ZenithAngle(Latitude,DayOfYear,HourOfDay) - irr=1050.*np.cos(z*deg2rad) - if irr<0: - irr=0 - return irr -def DailyIrradiance(Latitude,DayOfYear,display=False): - ''' - return solar irradiance for every 5 min - Latitude: (+) North, (-) South - DayOfYear: 0 for January 1st - ''' - h=np.linspace(0,24,24*12+1) - irr=np.array([SolarIrradiance(Latitude,DayOfYear,hr) for hr in h]) - if display: - plt.plot(h,irr) - return h,irr -def Sunrise(Latitude,DayOfYear): - ''' - calculate sunrise time - Latitude: (+) North, (-) South - DayOfYear: 0 for January 1st - ''' - f=lambda h:np.cos(ZenithAngle(Latitude,DayOfYear,h)*deg2rad) - hl,hr=0,12 - h_sunrise=Bisection(f,hl,hr) - sr_time=hour2time(h_sunrise) - return h_sunrise,sr_time - -def Sunset(Latitude,DayOfYear): - ''' - calculate sunset time - Latitude: (+) North, (-) South - DayOfYear: 0 for January 1st - ''' - f=lambda h:np.cos(ZenithAngle(Latitude,DayOfYear,h)*deg2rad) - hl,hr=12,24 - h_sunset=Bisection(f,hl,hr) - ss_time=hour2time(h_sunset) - return h_sunset,ss_time - -def Bisection(func,a,b,max_iter=10000,epsilon=1e-9): - ''' - bisection method for finding root of a function - func: a single variable function - a,b: upper and lower bounds - ''' - if func(a)*func(b)>0: - return None - elif abs(func(a))