diff --git a/.github/workflows/quality_checks.yml b/.github/workflows/quality_checks.yml new file mode 100644 index 0000000..4fdb121 --- /dev/null +++ b/.github/workflows/quality_checks.yml @@ -0,0 +1,78 @@ +name: Quality Checks + +on: + - push + - pull_request + +jobs: + format: + name: Check formatting + runs-on: ubuntu-22.04 + steps: + - uses: actions/checkout@v4 + + - uses: actions/setup-python@v4 + with: + python-version: "3.10" + + - name: Install tox + run: python -m pip install tox + + - name: Format with black + run: tox -e format + + + lint: + name: Lint + runs-on: ubuntu-22.04 + steps: + - uses: actions/checkout@v4 + + - uses: actions/setup-python@v4 + with: + python-version: "3.10" + + - name: Install tox + run: python -m pip install tox + + - name: Run linter with flake8 + run: tox -e lint + + typecheck: + name: Type check + runs-on: ubuntu-22.04 + steps: + - uses: actions/checkout@v4 + + - uses: actions/setup-python@v4 + with: + python-version: "3.10" + + - name: Install tox + run: python -m pip install tox + + - name: Type check with mypy + run: python -m tox -e typecheck + + test: + name: Test + runs-on: ubuntu-22.04 + strategy: + fail-fast: false # Whether to cancel all jobs if any matrix job fails + matrix: + include: + - {python-version: "3.11", toxenv: "py311"} + - {python-version: "3.10", toxenv: "py310"} + - {python-version: "3.9", toxenv: "py39"} + steps: + - uses: actions/checkout@v4 + + - uses: actions/setup-python@v4 + with: + python-version: ${{ matrix.python-version }} + + - name: Install tox + run: python -m pip install tox + + - name: Execute tests with pytest + run: tox -e ${{ matrix.toxenv }} diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 71c9e92..657ffda 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -7,7 +7,7 @@ repos: - id: trailing-whitespace - id: requirements-txt-fixer - repo: https://github.com/psf/black - rev: 23.7.0 + rev: 23.11.0 hooks: - id: black diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..261eeb9 --- /dev/null +++ b/LICENSE @@ -0,0 +1,201 @@ + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright [yyyy] [name of copyright owner] + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000..c5c8468 --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1,2 @@ +graft src +recursive-exclude __pycache__ *.py[cod] diff --git a/dev-requirements.txt b/dev-requirements.txt new file mode 100644 index 0000000..bc6bd27 --- /dev/null +++ b/dev-requirements.txt @@ -0,0 +1,4 @@ +build>=0.10.0 +pre-commit>=3.3.3 +scipy>=1.11.4 +tox>=4.11.1 diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..c066da7 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,7 @@ +[build-system] +requires = ["setuptools", "wheel"] +build-backend = "setuptools.build_meta" + +[tool.black] +line-length = 88 +target-version = ['py310'] diff --git a/requirements.txt b/requirements.txt deleted file mode 100644 index c9d8775..0000000 --- a/requirements.txt +++ /dev/null @@ -1 +0,0 @@ -pre-commit>=3.3.3 diff --git a/setup.cfg b/setup.cfg new file mode 100644 index 0000000..287b8ec --- /dev/null +++ b/setup.cfg @@ -0,0 +1,90 @@ +[metadata] +name = pysolorie +version = 1.0.0 +description = Orientation Analysis of Solar Panel +long_description = file: README.md +long_description_content_type = text/markdown +url = https://github.com/aaghamohammadi/pysolorie +author = Alireza Aghamohammadi +author_email = "Alireza Aghamohammadi" +license = {file = "LICENSE" } +classifiers = + License :: OSI Approved :: Apache Software License + +[options] +package_dir = + =src +packages = find: +include_package_data = True +install_requires = + scipy>=1.11.4 + +[options.packages.find] +where = src +exclude = + tests* + +[mypy] +python_version = 3.10 +warn_unused_configs = True +show_error_context = True +pretty = True +namespace_packages = True +check_untyped_defs = True + +[flake8] +max-line-length = 88 + +[tool:pytest] +testpaths = tests +addopts = --cov --strict-markers +xfail_strict = True + +[coverage:run] +source = pysolorie +branch = True + +[coverage:report] +show_missing = True +skip_covered = True + +[coverage:paths] +source = + src/pysolorie + */site-packages/pysolorie + +[tox:tox] +envlist = py39, py310, py311 +isolated_build = True + + +[testenv] +deps = + pytest + pytest-cov +commands = + pytest {posargs} + +[testenv:typecheck] +deps = + mypy + pytest +commands = + mypy {posargs:src tests} + +[testenv:format] +skip_install = True +deps = + black + isort +commands = + black {posargs:--check --diff src tests} + isort {posargs:--check --diff src tests} --profile=black + +[testenv:lint] +skip_install = True +deps = + flake8 + flake8-bugbear +commands = + flake8 {posargs:src tests} diff --git a/src/pysolorie/__init__.py b/src/pysolorie/__init__.py new file mode 100644 index 0000000..aad9c11 --- /dev/null +++ b/src/pysolorie/__init__.py @@ -0,0 +1,29 @@ +# Copyright 2023 Alireza Aghamohammadi + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from .atmospheric_transmission import AtmosphericTransmission +from .irradiance import SolarIrradiance +from .model import HottelModel +from .numerical_integration import IrradiationCalculator +from .observer import Observer +from .sun_position import SunPosition + +__all__ = [ + "HottelModel", + "SunPosition", + "SolarIrradiance", + "Observer", + "AtmosphericTransmission", + "IrradiationCalculator", +] diff --git a/src/pysolorie/atmospheric_transmission.py b/src/pysolorie/atmospheric_transmission.py new file mode 100644 index 0000000..cbec258 --- /dev/null +++ b/src/pysolorie/atmospheric_transmission.py @@ -0,0 +1,81 @@ +# Copyright 2023 Alireza Aghamohammadi + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +import math + +from .model import HottelModel +from .observer import Observer + + +class AtmosphericTransmission: + """ + This class calculates the effective atmospheric transmission coefficient + of the direct beam. + """ + + def __init__( + self, + climate_type: str, + observer_altitude: int, + observer_latitude: float, + ): + r""" + Initialize the AtmosphericTransmission class. + + :param climate_type: The type of climate. + :type climate_type: str + :param observer_altitude: The altitude of the observer in meters. + :type observer_altitude: int + :param observer_latitude: The latitude of the observer in degrees. + :type observer_latitude: float + """ + self.hottel_model = HottelModel() + ( + self.a0, + self.a1, + self.k, + ) = self.hottel_model.calculate_transmittance_components( + climate_type, observer_altitude + ) + self.observer = Observer(observer_latitude) + + def calculate_transmittance(self, day_of_year: int, solar_time: float) -> float: + r""" + Calculate the effective atmospheric transmission coefficient of the direct beam. + + The effective atmospheric transmission coefficient of the direct beam + is calculated using the formula: + + .. math:: + \tau_b = a_0 + a_1 \cdot \exp\left(-\frac{k}{\cos(\theta_z)}\right) + + where: + \(\tau_b\) is the effective atmospheric transmission coefficient + of the direct beam, \(a_0\), \(a_1\), and \(k\) + are the components of clear-sky beam radiation transmittance, + and \(\theta_z\) is the zenith angle. + + :param day_of_year: The day of the year. + :type day_of_year: int + :param solar_time: The solar time in seconds. + :type solar_time: float + :return: The effective atmospheric transmission coefficient of the direct beam. + :rtype: float + """ + zenith_angle = self.observer.calculate_zenith_angle(day_of_year, solar_time) + cos_zenith_angle = math.cos(zenith_angle) + EPSILON = 1e-8 # Small constant to prevent division by zero + if abs(cos_zenith_angle) < EPSILON: + return self.a0 + else: + return self.a0 + self.a1 * math.exp(-self.k / cos_zenith_angle) diff --git a/src/pysolorie/irradiance.py b/src/pysolorie/irradiance.py new file mode 100644 index 0000000..1b04a10 --- /dev/null +++ b/src/pysolorie/irradiance.py @@ -0,0 +1,71 @@ +# Copyright 2023 Alireza Aghamohammadi + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import math + +from .sun_position import SunPosition + + +class SolarIrradiance: + def __init__(self, sun_position: SunPosition): + r""" + Initialize the SolarIrradiance class. + + :param sun_position: An instance of the SunPosition class. + :type sun_position: SunPosition + """ + self.sun_position = sun_position + + def calculate_extraterrestrial_irradiance(self, day_of_year: int) -> float: + r""" + Calculate the extraterrestrial solar irradiance for a given day of the year. + + The extraterrestrial solar irradiance, \(E\), + is the amount of solar energy received per unit area on + a surface perpendicular to the Sun's rays outside Earth's atmosphere. + + The formula used is: + + .. math:: + + E = \text{{SOLAR\_CONSTANT}} + \times (1 + 0.33 \times \cos (2 + \times \pi \times \frac{{\text{{day\_of\_year}}}}{365})) + + where: + \(\text{{SOLAR\_CONSTANT}}\) is the average solar radiation arriving outside + of the Earth's atmosphere, which is approximately 1367 Watts per square meter. + This is also known as the solar constant.The factor 0.033 accounts + for the variation in the Earth-Sun distance + due to the Earth's elliptical orbit. + + :param day\_of\_year: The day of the year, ranging from 1 to 365. + :type day\_of\_year: int + + :return: The extraterrestrial solar irradiance in Watts per square meter. + :rtype: float + """ + + # Solar constant (W/m^2) + SOLAR_CONSTANT = 1367 + + # Factor to account for the Earth's orbital eccentricity + earth_orbital_eccentricity = 0.033 + + # Calculate the extraterrestrial solar irradiance + extraterrestrial_irradiance = SOLAR_CONSTANT * ( + 1 + earth_orbital_eccentricity * math.cos(2 * math.pi * day_of_year / 365) + ) + + return extraterrestrial_irradiance diff --git a/src/pysolorie/model.py b/src/pysolorie/model.py new file mode 100644 index 0000000..dec37a8 --- /dev/null +++ b/src/pysolorie/model.py @@ -0,0 +1,133 @@ +# Copyright 2023 Alireza Aghamohammadi + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from typing import Dict, Tuple + + +class HottelModel: + r""" + Hottel Model for estimating clear-sky beam radiation transmittance + based on climate type, and observer altitude. + + :ivar CLIMATE_CONSTANTS: Correction factors for different climate types + (\(r_0\), \(r_1\), and \(r_k\)). + """ + + CLIMATE_CONSTANTS: Dict[str, Tuple[float, float, float]] = { + "TROPICAL": (0.95, 0.98, 1.02), + "MIDLATITUDE SUMMER": (0.97, 0.99, 1.02), + "SUBARCTIC SUMMER": (0.99, 0.99, 1.01), + "MIDLATITUDE WINTER": (1.03, 1.01, 1.00), + } + + def _convert_to_km(self, observer_altitude: int) -> float: + r""" + Convert altitude from meters to kilometers. + + :param observer_altitude: Altitude of the observer in meters. + :type observer_altitude: float + :return: Altitude in kilometers. + :rtype: float + """ + return observer_altitude / 1000.0 + + def _calculate_a0_star(self, observer_altitude: int) -> float: + r""" + Calculate \(a_0^*\) based on observer altitude. + + The formula used to calculate \(a_0^*\) is: + + .. math:: + a_0^* = 0.4237 - 0.00821 \cdot (6 - A)^2 + + :param observer_altitude: Altitude of the observer in meters. + :type observer_altitude: float + :return: \(a_0^*\) value. + :rtype: float + """ + observer_altitude_km = self._convert_to_km(observer_altitude) + observer_altitude_diff = 6.0 - observer_altitude_km + return 0.4237 - 0.00821 * observer_altitude_diff**2 + + def _calculate_a1_star(self, observer_altitude: int) -> float: + r""" + Calculate \(a_1^*\) based on observer altitude. + + The formula used to calculate \(a_1^*\) is: + + .. math:: + a_1^* = 0.5055 + 0.00595 \cdot (6.5 - A)^2 + + :param observer_altitude: Altitude of the observer in meters. + :type observer_altitude: float + :return: \(a_1^*\) value. + :rtype: float + """ + observer_altitude_km = self._convert_to_km(observer_altitude) + observer_altitude_diff = 6.5 - observer_altitude_km + return 0.5055 + 0.00595 * observer_altitude_diff**2 + + def _calculate_k_star(self, observer_altitude: int) -> float: + r""" + Calculate \(k^*\) based on observer altitude. + + The formula used to calculate \(k^*\) is: + + .. math:: + k^* = 0.2711 + 0.01858 \cdot (2.5 - A)^2 + + :param observer_altitude: Altitude of the observer in meters. + :type observer_altitude: float + :return: \(k^*\) value. + :rtype: float + """ + observer_altitude_km = self._convert_to_km(observer_altitude) + observer_altitude_diff = 2.5 - observer_altitude_km + return 0.2711 + 0.01858 * observer_altitude_diff**2 + + def calculate_transmittance_components( + self, climate_type: str, observer_altitude: int + ) -> Tuple[float, float, float]: + r""" + Calculate the components of clear-sky beam radiation transmittance + (\(a_0\), \(a_1\), and \(k\)) based on climate type and observer altitude. + + Correction factors adjust the clear-sky beam radiation transmittance components. + + The formulas used to calculate the components are: + + .. math:: + a_0 = r_0 \cdot a_0^* + a_1 = r_1 \cdot a_1^* + k = r_k \cdot k^* + + :param climate_type: Climate type (e.g., "TROPICAL"). + :type climate_type: str + :param observer_altitude: Altitude of the observer in meters. + :type observer_altitude: float + :return: Components of clear-sky beam radiation transmittance + (\(a_0\), \(a_1\), \(k\)). + :rtype: tuple of floats + :raises ValueError: If an invalid climate type is provided. + """ + if climate_type.upper() not in self.CLIMATE_CONSTANTS: + raise ValueError("Invalid climate type") + + r0, r1, rk = self.CLIMATE_CONSTANTS[climate_type.upper()] + + a0_star = self._calculate_a0_star(observer_altitude) + a1_star = self._calculate_a1_star(observer_altitude) + k_star = self._calculate_k_star(observer_altitude) + + return r0 * a0_star, r1 * a1_star, rk * k_star diff --git a/src/pysolorie/numerical_integration.py b/src/pysolorie/numerical_integration.py new file mode 100644 index 0000000..0ab1b55 --- /dev/null +++ b/src/pysolorie/numerical_integration.py @@ -0,0 +1,156 @@ +# Copyright 2023 Alireza Aghamohammadi + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import math + +import numpy as np +from scipy import integrate # type: ignore + +from .atmospheric_transmission import AtmosphericTransmission +from .irradiance import SolarIrradiance +from .observer import Observer +from .sun_position import SunPosition + + +class IrradiationCalculator: + OMEGA = 7.15 * 1e-5 + + def __init__( + self, climate_type: str, observer_altitude: int, observer_latitude: float + ): + """ + Initialize the IrradiationCalculator class. + + :param climate_type: The type of climate. + :type climate_type: str + :param observer_altitude: The altitude of the observer in meters. + :type observer_altitude: int + :param observer_latitude: The latitude of the observer in degrees. + :type observer_latitude: float + """ + self._observer = Observer(observer_latitude=observer_latitude) + self._sun_position = SunPosition() + self._solar_irradiance = SolarIrradiance(self._sun_position) + self._atmospheric_transmission = AtmosphericTransmission( + climate_type, observer_altitude, observer_latitude + ) + + @staticmethod + def _heaviside(x: float) -> int: + """ + Heaviside step function. + + :param x: Input to the function. + :type x: float + :return: 1 if x >= 0, else 0. + :rtype: int + """ + return 1 if x >= 0 else 0 + + def _calculate_irradiance_component( + self, hour_angle: float, panel_orientation: float, day_of_year: int + ) -> float: + r""" + Calculate a component of the irradiance integral. + + :param hour_angle: The hour angle. + :type hour_angle: float + :param panel_orientation: The orientation of the solar panel in radians. + :type panel_orientation: float + :param day_of_year: The day of the year. + :type day_of_year: int + :return: A component of the irradiance integral. + :rtype: float + """ + observer_latitude = self._observer._validate_latitude() + solar_declination = self._sun_position.solar_declination(day_of_year) + cos_theta = math.sin(solar_declination) * math.sin( + observer_latitude - panel_orientation + ) + math.cos(solar_declination) * math.cos(hour_angle) * math.cos( + observer_latitude - panel_orientation + ) + transmittance = self._atmospheric_transmission.calculate_transmittance( + day_of_year, self._sun_position.solar_time(hour_angle) + ) + irradiance = self._solar_irradiance.calculate_extraterrestrial_irradiance( + day_of_year + ) + return ( + irradiance + * transmittance + * cos_theta + * self._heaviside(cos_theta) + / self.OMEGA + ) + + def calculate_direct_irradiation( + self, panel_orientation: float, day_of_year: int + ) -> float: + r""" + Calculate the total direct irradiation + for a given solar panel orientation (beta). + + The total direct irradiation is calculated using the formula: + + .. math:: + E_b(n,\phi) = \frac{E}{\Omega} \int_{\omega_s}^{\omega_t} + \cos(θ) H(\cos(θ)) \times \tau_b d\omega + + + where: + n is the day of year + \(\phi\) is the latitude of the observer, and + \(E\) is the amount of solar energy received + \Omega is a constant equal to 7.15 * 1e-5 + theta is incidence angle + \omega_s is the sunrise hour angle + \omega_t is the sunset hour angle + H is the heaviside step function + + + :param panel_orientation: The orientation of the solar panel in radians. + :type panel_orientation: float + :param day_of_year: The day of the year. + :type day_of_year: int + :return: The total direct irradiation. + :rtype: float + """ + sunrise_hour_angle, sunset_hour_angle = self._observer.calculate_sunrise_sunset( + day_of_year + ) + irradiance_components = [ + self._calculate_irradiance_component( + hour_angle, panel_orientation, day_of_year + ) + for hour_angle in np.arange(sunrise_hour_angle, sunset_hour_angle, 0.01) + ] + return integrate.simpson(irradiance_components, dx=0.01) + + def find_optimal_orientation(self, day_of_year: int) -> float: + """ + Find the optimal orientation (beta) that maximizes the total direct irradiation. + + :param day_of_year: The day of the year. + :type day_of_year: int + :return: The optimal orientation (beta) in degrees. + :rtype: float + """ + betas = np.arange(-math.pi / 2, math.pi / 2, 0.005) # Discretize beta + irradiations = [ + self.calculate_direct_irradiation(beta, day_of_year) for beta in betas + ] + optimal_beta = betas[ + np.argmax(irradiations) + ] # Find beta that gives max irradiation + return math.degrees(optimal_beta) diff --git a/src/pysolorie/observer.py b/src/pysolorie/observer.py new file mode 100644 index 0000000..051eb70 --- /dev/null +++ b/src/pysolorie/observer.py @@ -0,0 +1,111 @@ +# Copyright 2023 Alireza Aghamohammadi + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import math +from typing import Optional + +from .sun_position import SunPosition + + +class Observer: + def __init__( + self, + observer_latitude: Optional[float] = None, + observer_longitude: Optional[float] = None, + ): + """ + Initialize the Observer class. + + :param observer_latitude: The latitude of the observer in degrees (optional). + :type observer_latitude: Optional[float] + :param observer_longitude: The longitude of the observer in degrees (optional). + :type observer_longitude: Optional[float] + """ + self.observer_latitude = ( + math.radians(observer_latitude) if observer_latitude is not None else None + ) + self.observer_longitude = ( + math.radians(observer_longitude) if observer_longitude is not None else None + ) + self.sun_position = SunPosition() + + def calculate_zenith_angle(self, day_of_year: int, solar_time: float) -> float: + r""" + Calculate the zenith angle. + + The zenith angle is calculated using the formula: + + .. math:: + \cos(\theta_z) = \sin(\phi) \cdot \sin(\delta) + + \cos(\phi) \cdot \cos(\delta) \cdot \cos(\omega) + + where: + \(\theta_z\) is the zenith angle, + \(\phi\) is the latitude of the observer, + \(\delta\) is the solar declination, and + \(\omega\) is the hour angle. + + :param day_of_year: The day of the year. + :type day_of_year: int + :param solar_time: The solar time in seconds. + :type solar_time: float + :return: The zenith angle in radians. + :rtype: float + """ + observer_latitude = self._validate_latitude() + + solar_declination = self.sun_position.solar_declination(day_of_year) + hour_angle = self.sun_position.hour_angle(solar_time) + return math.acos( + math.sin(observer_latitude) * math.sin(solar_declination) + + math.cos(observer_latitude) + * math.cos(solar_declination) + * math.cos(hour_angle) + ) + + def calculate_sunrise_sunset(self, day_of_year: int) -> tuple: + r""" + Calculate the hour angle at sunrise and sunset. + + The hour angle at sunrise and sunset is calculated using the formula: + + .. math:: + \cos(\omega) = -\tan(\phi) \cdot \tan(\delta) + + where: + \(\omega\) is the hour angle, + \(\phi\) is the latitude of the observer, and + \(\delta\) is the solar declination. + + :param day_of_year: The day of the year. + :type day_of_year: int + :return: The hour angle at sunrise and sunset in radians. + :rtype: tuple + """ + observer_latitude = self._validate_latitude() + + solar_declination = self.sun_position.solar_declination(day_of_year) + hour_angle = math.acos( + -math.tan(observer_latitude) * math.tan(solar_declination) + ) + + sunrise = -hour_angle + sunset = hour_angle + + return sunrise, sunset + + def _validate_latitude(self) -> float: + if self.observer_latitude is None: + raise ValueError("Observer latitude must be provided.") + return self.observer_latitude diff --git a/src/pysolorie/sun_position.py b/src/pysolorie/sun_position.py new file mode 100644 index 0000000..62af9ac --- /dev/null +++ b/src/pysolorie/sun_position.py @@ -0,0 +1,109 @@ +# Copyright 2023 Alireza Aghamohammadi + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +import math + + +class SunPosition: + def solar_declination(self, day_of_year: int) -> float: + r""" + Calculate the solar declination angle in radians. + + The solar declination angle is the angle between the rays of the sun + and the plane of the Earth's equator. + + The formula used to calculate the solar declination angle is: + + .. math:: + \delta = \sin \left( \frac{2 \pi}{365} + \times (284 + \text{{day\_of\_year}}) \right) + \times \left(\frac{23.45 \pi}{180}\right) + + :param day_of_year: The day of the year. + :type day_of_year: int + :return: The solar declination angle in radians. + :rtype: float + """ + # tilt of the Earth's axis (in degrees) + earth_tilt_degrees = 23.45 + + # Convert the tilt to radians + earth_tilt_radians = math.radians(earth_tilt_degrees) + + # Offset to ensure declination angle is zero at the March equinox + equinox_offset_days = 284 + + # Calculate the solar declination angle + solar_declination = ( + math.sin((2 * math.pi) * (equinox_offset_days + day_of_year) / 365) + * earth_tilt_radians + ) + + return solar_declination + + def hour_angle(self, solar_time: float) -> float: + r""" + Calculate the hour angle based on the solar time. + + The hour angle is a measure of time, expressed in angular terms, + from solar noon. + + The formula used to calculate the hour angle is: + + .. math:: + \omega = (t - \text{{seconds\_in\_half\_day}}) + \times \frac{\pi}{\text{{seconds\_in\_half\_day}}} + + :param solar_time: The solar time in seconds. + :type solar_time: float + :return: The hour angle in radians. + :rtype: float + """ + # The number of seconds in half a day (12 hours) + seconds_in_half_day = 12 * 60 * 60 + + # The Earth rotates by pi/seconds_in_half_day radians per second + earth_rotation_rate = math.pi / seconds_in_half_day + + # Calculate the hour angle + hour_angle = (solar_time - seconds_in_half_day) * earth_rotation_rate + + return hour_angle + + def solar_time(self, hour_angle: float) -> float: + r""" + Calculate the solar time based on the hour angle. + + The solar time is a measure of time, expressed in seconds, from solar noon. + + The formula used to calculate the solar time is: + + .. math:: + t = \omega \times \frac{\text{{seconds\_in\_half\_day}}}{\pi} + + \text{{seconds\_in\_half\_day}} + + :param hour_angle: The hour angle in radians. + :type hour_angle: float + :return: The solar time in seconds. + :rtype: float + """ + # The number of seconds in half a day (12 hours) + seconds_in_half_day = 12 * 60 * 60 + + # The Earth rotates by pi/seconds_in_half_day radians per second + earth_rotation_rate = math.pi / seconds_in_half_day + + # Calculate the solar time + solar_time = hour_angle / earth_rotation_rate + seconds_in_half_day + + return solar_time diff --git a/tests/__init__.py b/tests/__init__.py new file mode 100644 index 0000000..7150904 --- /dev/null +++ b/tests/__init__.py @@ -0,0 +1,13 @@ +# Copyright 2023 Alireza Aghamohammadi + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. diff --git a/tests/test_pysolorie.py b/tests/test_pysolorie.py new file mode 100644 index 0000000..b8a2142 --- /dev/null +++ b/tests/test_pysolorie.py @@ -0,0 +1,272 @@ +# Copyright 2023 Alireza Aghamohammadi + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import math +from typing import Tuple + +import pytest + +from pysolorie import ( + AtmosphericTransmission, + HottelModel, + IrradiationCalculator, + Observer, + SolarIrradiance, + SunPosition, +) + + +@pytest.mark.parametrize( + "climate_type, observer_altitude, expected_result", + [ + ("MIDLATITUDE SUMMER", 1200, (0.228, 0.666, 0.309)), # Tehran Summer + ("MIDLATITUDE WINTER", 1200, (0.242, 0.679, 0.303)), # Tehran Winter + ("TROPICAL", 26, (0.124, 0.739, 0.392)), # Medan + ("SUBARCTIC SUMMER", 136, (0.140, 0.739, 0.379)), # Fairbanks + ], +) +def test_calculate_transmittance_components( + climate_type: str, + observer_altitude: int, + expected_result: Tuple[float, float, float], +) -> None: + hottel_model: HottelModel = HottelModel() + result: Tuple[ + float, float, float + ] = hottel_model.calculate_transmittance_components(climate_type, observer_altitude) + assert pytest.approx(result, abs=1e-3) == expected_result + + +def test_invalid_climate_type() -> None: + with pytest.raises(ValueError, match="Invalid climate type"): + hottel_model: HottelModel = HottelModel() + hottel_model.calculate_transmittance_components("INVALID", 1000) + + +@pytest.mark.parametrize( + "day_of_year, solar_time, expected_declination, expected_hour_angle", + [ + (1, 12 * 60 * 60, -0.4014257279586958, 0), # January 1st at noon + (81, 10 * 60 * 60, 0, -math.pi / 6), # March 22nd at 10am (equinox) + (81, 12 * 60 * 60, 0, 0), # March 22nd at noon (equinox) + (1, 13 * 60 * 60, -0.4014257279586958, math.pi / 12), # January 1st at 1pm + ], +) +def test_sun_position( + day_of_year: int, + solar_time: int, + expected_declination: float, + expected_hour_angle: float, +) -> None: + sun_position: SunPosition = SunPosition() + declination: float = sun_position.solar_declination(day_of_year) + hour_angle: float = sun_position.hour_angle(solar_time) + assert declination == pytest.approx(expected_declination, abs=1e-3) + assert hour_angle == pytest.approx(expected_hour_angle, abs=1e-3) + + +@pytest.mark.parametrize( + "hour_angle, expected_solar_time", + [ + (0, 12 * 60 * 60), # solar noon + (-math.pi / 6, 10 * 60 * 60), # 10am + (math.pi / 12, 13 * 60 * 60), # 1pm + ], +) +def test_solar_time( + hour_angle: float, + expected_solar_time: int, +) -> None: + sun_position: SunPosition = SunPosition() + solar_time: float = sun_position.solar_time(hour_angle) + assert solar_time == pytest.approx(expected_solar_time, abs=1e-3) + + +@pytest.mark.parametrize( + "day_of_year, expected_irradiance", + [ + (1, 1412.104), # January 1st + (81, 1374.918), # March 22nd (equinox) + (172, 1322.623), # June 21st (summer solstice) + (264, 1359.464), # September 23rd (equinox) + (355, 1411.444), # December 21st (winter solstice) + ], +) +def test_calculate_extraterrestrial_irradiance( + day_of_year: int, expected_irradiance: float +) -> None: + sun_position = SunPosition() + solar_irradiance = SolarIrradiance(sun_position) + irradiance: float = solar_irradiance.calculate_extraterrestrial_irradiance( + day_of_year + ) + assert irradiance == pytest.approx(expected_irradiance, abs=1e-3) + + +@pytest.mark.parametrize( + "observer_latitude," + + "observer_longitude," + + "day_of_year," + + "solar_time," + + "expected_zenith_angle", + [ + ( + 35.69, + 51.39, + 81, + 12 * 60 * 60, + 0.623, + ), # Test at Tehran on the 81st day of the year at noon + ( + 35.69, + 51.39, + 355, + 12 * 60 * 60, + 1.032, + ), # Test at Tehran on the 355th day of the year at noon + ( + 3.59, + 98.67, + 81, + 10 * 60 * 60, + 0.527, + ), # Test at Medan on the 81st day of the year at 10am + ( + 64.84, + -147.72, + 1, + 13 * 60 * 60, + 1.547, + ), # Test at Fairbanks on the 1st day of the year at 1pm + (90, 0, 1, 13 * 60 * 60, 1.972), # Test at North Pole on January 1st at 1pm + ], +) +def test_calculate_zenith_angle( + observer_latitude: float, + observer_longitude: float, + day_of_year: int, + solar_time: int, + expected_zenith_angle: float, +) -> None: + observer: Observer = Observer(observer_latitude, observer_longitude) + zenith_angle: float = observer.calculate_zenith_angle(day_of_year, solar_time) + assert zenith_angle == pytest.approx(expected_zenith_angle, abs=1e-3) + + +def test_calculate_zenith_angle_without_latitude(): + observer = Observer(None, 0) + with pytest.raises( + ValueError, + match="Observer latitude must be provided.", + ): + observer.calculate_zenith_angle(1, 12 * 60 * 60) + + +@pytest.mark.parametrize( + "climate_type," + + "observer_altitude," + + "observer_latitude," + + "day_of_year," + + "solar_time," + + "expected_transmittance", + [ + ("MIDLATITUDE SUMMER", 1200, 35.69, 81, 12 * 60 * 60, 0.683), # Tehran Summer + ("MIDLATITUDE WINTER", 1200, 35.69, 355, 12 * 60 * 60, 0.618), # Tehran Winter + ("TROPICAL", 63, 3.59, 81, 10 * 60 * 60, 0.597), # Medan + ("SUBARCTIC SUMMER", 136, 64.84, 1, 13 * 60 * 60, 0.140), # Fairbanks + ], +) +def test_calculate_transmittance( + climate_type: str, + observer_altitude: int, + observer_latitude: float, + day_of_year: int, + solar_time: float, + expected_transmittance: float, +) -> None: + atmospheric_transmission: AtmosphericTransmission = AtmosphericTransmission( + climate_type, observer_altitude, observer_latitude + ) + result: float = atmospheric_transmission.calculate_transmittance( + day_of_year, solar_time + ) + assert pytest.approx(result, abs=1e-3) == expected_transmittance + + +@pytest.mark.parametrize( + "day_of_year, expected_sunrise_hour_angle, expected_sunset_hour_angle", + [ + (1, -1.261, 1.261), # January 1st + (81, -math.pi / 2, math.pi / 2), # March 22nd (equinox) + (172, -1.888, 1.888), # June 21st (solstice) + ], +) +def test_calculate_sunrise_sunset( + day_of_year: int, + expected_sunrise_hour_angle: float, + expected_sunset_hour_angle: float, +) -> None: + observer: Observer = Observer(observer_latitude=35.69) # Tehran + sunrise_hour_angle, sunset_hour_angle = observer.calculate_sunrise_sunset( + day_of_year + ) + assert sunrise_hour_angle == pytest.approx(expected_sunrise_hour_angle, abs=1e-3) + assert sunset_hour_angle == pytest.approx(expected_sunset_hour_angle, abs=1e-3) + + +@pytest.mark.parametrize( + "climate_type, observer_altitude, observer_latitude, day_of_year, expected_result", + [ + ( + "MIDLATITUDE SUMMER", + 1200, + 35.6892, + 172, + 0.241, + ), # Tehran Summer, day_of_year=172 (June 21) + ( + "MIDLATITUDE WINTER", + 1200, + 35.6892, + 355, + 63.839, + ), # Tehran Winter, day_of_year=355 (Dec 21) + ( + "TROPICAL", + 26, + 3.5952, + 100, + -6.635, + ), # Medan, day_of_year=100 (April 10) + ( + "SUBARCTIC SUMMER", + 132, + 64.84361, + 200, + 32.613, + ), # Fairbanks Summer, day_of_year=200 (July 19) + ], +) +def test_find_optimal_orientation( + climate_type: str, + observer_altitude: int, + observer_latitude: float, + day_of_year: int, + expected_result: float, +) -> None: + irradiation_calculator = IrradiationCalculator( + climate_type, observer_altitude, observer_latitude + ) + result = irradiation_calculator.find_optimal_orientation(day_of_year) + assert pytest.approx(result, abs=1e-3) == expected_result