rewrite functions in pv_model

docs/readthedocs/_static/bibliography/bibtexAll.bib

@@ -15,7 +15,7 @@ @Article{Maleki.2017
 @MISC{Itaca_Sun,
 author = {Itacanet},
 title = {The Sun As A Source Of Energy},
-howpublished={\url{https://www.itacanet.org/the-sun-as-a-source-of-energy/part-3-calculating-solar-angles/}}
+howpublished={\url{https://de.scribd.com/document/455342846/Part-3-Calculating-Solar-Angles-ITACA}}
 }
 
 @article{Spencer.1971,
@@ -73,7 +73,7 @@ @book{Quaschning.2013
  author = {Quaschning, Volker},
  year = {2013},
  title = {Regenerative Energiesysteme: Technologie; Berechnung; Simulation},
- url = {http://ebooks.ciando.com/book/index.cfm/bok_id/471802},
+ url = {https://www.hanser-elibrary.com/doi/book/10.3139/9783446435711},
  price = {39.99 EUR},
  address = {M{\"u}nchen},
  edition = {8., aktualisierte und erw. Aufl.},
@@ -92,7 +92,8 @@ @Inbook{Schoenberg.1929
 pages={1--280},
 abstract={Die ersten Versuche, Messungen der Lichtst{\"a}rke der Himmelsk{\"o}rper auszuf{\"u}hren, um auf diesem Wege Aufschl{\"u}sse {\"u}ber ihre Beschaffenheit zu erhalten, fallen in die Zeit jenes gewaltigen Aufschwungs, welchen die Optik durch die Arbeiten von Newton und Huygens am Ende des siebzehnten und zu Beginn des achtzehnten Jahrhunderts erfahren hatte. Schon Huygens1 selbst versuchte die Helligkeit des Sirius mit derjenigen der Sonne zu vergleichen, indem er das Sonnenlicht durch eine kleine {\"O}ffnung im verschlossenen Ende eines langen Rohres abschw{\"a}chte. Der schwedische Physiker Celsius2 suchte nach einem Gesetze der Lichtabnahme beleuchteter Fl{\"a}chen, doch waren seine Schlu{\ss}folgerungen, ebenso wie diejenigen von Huygens, infolge der Unzul{\"a}nglichkeit der angewandten Methoden nicht stichhaltig. Au{\ss}er der Formel der Lichtabnahme mit dem Quadrate der Entfernung besa{\ss} die Physik noch keinerlei R{\"u}stzeug an strengen Definitionen und Gesetzen und keinerlei Apparate zur Messung der Lichtst{\"a}rken. Buffon3 versuchte den Verlust zu bestimmen, den das Sonnenlicht bei Reflexion an gl{\"a}sernen Spiegeln erleidet. Aber erst die gro{\ss} angelegten Arbeiten von Pierre Bouguer4 5 (1698--1758) und Johann Heinrich Lambert (1728--1777) schufen die Grundlagen der Photometrie. In systematischer experimenteller Arbeit, die durch sinnreiche Theorien erl{\"a}utert wird, behandelt Bouguer ihre Hauptprobleme: die Absorption des Lichtes bei der Reflexion und beim Durchgang durch feste und fl{\"u}ssige K{\"o}rper sowie die diffuse Reflexion an matten Fl{\"a}chen und beim Durchgange durch tr{\"u}be Medien.},
 isbn={978-3-642-90703-6},
-doi={10.1007/978-3}
+doi={10.1007/978-3},
+url = {https://link.springer.com/chapter/10.1007/978-3-642-90703-6_1}
 }
 
 @MISC{WikiAirMass,
@@ -122,7 +123,7 @@ @book{Iqbal.1983
  author = {Iqbal, Muhammad},
  year = {1983},
  title = {An Introduction To Solar Radiation},
- url = {http://site.ebrary.com/lib/alltitles/docDetail.action?docID=10678925},
+ url = {https://books.google.de/books/about/An_Introduction_To_Solar_Radiation.html?id=3_qWce_mbPsC&redir_esc=y},
  address = {Oxford},
  publisher = {{Elsevier Science}},
  isbn = {9780123737502}

docs/readthedocs/models/pv_model.md

@@ -240,6 +240,8 @@ $$
 ```{eval-rst} 
 * :cite:ts:`Zheng.2017` p. 53, formula 2.3b
 * :cite:ts:`Iqbal.1983`
+* :cite:ts:`Spencer.1971`
+* :cite:ts:`Duffie.2013`
 ```
 
 ### Beam Radiation on Sloped Surface
@@ -310,6 +312,8 @@ $$
 \epsilon = \frac{\frac{E_{dif,H} + E_{beam,H}}{E_{dif,H}} + 5.535 \cdot 10^{-6} \cdot \theta_{z}^3}{1 + 5.535 \cdot 10^{-6} \cdot \theta_{z}^3}
 $$
 
+where the angle $theta_{z}$ is in **degrees**.
+
 Calculating a brightness index
 
 $$
@@ -412,6 +416,7 @@ $$
 * :cite:ts:`Perez.1987`
 * :cite:ts:`Perez.1990`
 * :cite:ts:`Myers.2017` p. 96f
+* :cite:ts:`Duffie.2013` p. 95f
 ```
 
 ### Reflected Radiation on Sloped Surface

src/main/scala/edu/ie3/simona/model/participant/PvModel.scala

@@ -115,6 +115,7 @@ final case class PvModel private (
       lat,
       gammaE,
       alphaE,
+      duration,
     )
 
     // === Diffuse Radiation Parameters ===//
@@ -337,7 +338,7 @@ final case class PvModel private (
     val e0 = 1.000110 +
       0.034221 * cos(jInRad) +
       0.001280 * sin(jInRad) +
-      0.000719 * cos(2d * jInRad) +
+      0.00719 * cos(2d * jInRad) +
       0.000077 * sin(2d * jInRad)
 
     // solar constant in W/m2
@@ -409,32 +410,27 @@ final case class PvModel private (
       omegaSS: Angle,
       omegaSR: Angle,
   ): Option[(Angle, Angle)] = {
-    val thetaGInRad = thetaG.toRadians
     val omegaSSInRad = omegaSS.toRadians
     val omegaSRInRad = omegaSR.toRadians
 
     val omegaOneHour = toRadians(15d)
-    val omegaHalfHour = omegaOneHour / 2d
 
     val omega1InRad = omega.toRadians // requested hour
     val omega2InRad = omega1InRad + omegaOneHour // requested hour plus 1 hour
 
-    // (thetaG < 90°): sun is visible
-    // (thetaG > 90°), otherwise: sun is behind the surface  -> no direct radiation
     if (
-      thetaGInRad < toRadians(90)
-      // omega1 and omega2: sun has risen and has not set yet
-      && omega2InRad > omegaSRInRad + omegaHalfHour
-      && omega1InRad < omegaSSInRad - omegaHalfHour
+      // requested time is between sunrise and sunset (+/- one hour)
+      omega1InRad > omegaSRInRad - omegaOneHour
+      && omega1InRad < omegaSSInRad
     ) {
 
       val (finalOmega1, finalOmega2) =
         if (omega1InRad < omegaSRInRad) {
           // requested time earlier than sunrise
-          (omegaSRInRad, omegaSRInRad + omegaOneHour)
-        } else if (omega2InRad > omegaSSInRad) {
-          // sunset earlier than requested time
-          (omegaSSInRad - omegaOneHour, omegaSSInRad)
+          (omegaSRInRad, omega2InRad)
+        } else if (omega1InRad > omegaSSInRad - omegaOneHour) {
+          // requested time is less than one hour before sunset
+          (omega1InRad, omegaSSInRad)
         } else {
           (omega1InRad, omega2InRad)
         }
@@ -464,13 +460,23 @@ final case class PvModel private (
     * @return
     *   the beam radiation on the sloped surface
     */
+
+  def calculateTimeFrame(omegas: Option[(Angle, Angle)]): Double = {
+    omegas match {
+      case Some((omega1, omega2)) =>
+        (omega2 - omega1).toDegrees / 15
+
+      case None => 0d
+    }
+  }
   private def calcBeamRadiationOnSlopedSurface(
       eBeamH: Irradiation,
       omegas: Option[(Angle, Angle)],
       delta: Angle,
       latitude: Angle,
       gammaE: Angle,
       alphaE: Angle,
+      duration: Time,
   ): Irradiation = {
 
     omegas match {
@@ -482,6 +488,9 @@ final case class PvModel private (
 
         val omega1InRad = omega1.toRadians
         val omega2InRad = omega2.toRadians
+        // variable that accounts for cases when the integration interval is shorter than 15° (1 hour equivalent), when the time is close to sunrise or sunset
+        val timeFrame =
+          (omega2 - omega1).toDegrees / 15d / duration.toHours // original term: (omega2 - omega1).toRadians * 180 / Math.PI / 15, since a one hour difference equals 15°
 
         val a = ((sin(deltaInRad) * sin(latInRad) * cos(gammaEInRad)
           - sin(deltaInRad) * cos(latInRad) * sin(gammaEInRad) * cos(
@@ -503,7 +512,7 @@ final case class PvModel private (
 
         // in rare cases (close to sunrise) r can become negative (although very small)
         val r = max(a / b, 0d)
-        eBeamH * r
+        eBeamH * r * timeFrame
       case None => WattHoursPerSquareMeter(0d)
     }
   }
@@ -533,6 +542,41 @@ final case class PvModel private (
     * @return
     *   the diffuse radiation on the sloped surface
     */
+
+  private def calcEpsilon(
+      eDifH: Irradiation,
+      eBeamH: Irradiation,
+      thetaZ: Angle,
+  ): Double = {
+    val thetaZInRad = thetaZ.toRadians
+
+    ((eDifH + eBeamH / cos(thetaZInRad)) / eDifH + (5.535d * 1.0e-6) * pow(
+      thetaZ.toDegrees,
+      3,
+    )) /
+      (1d + (5.535d * 1.0e-6) * pow(thetaZ.toDegrees, 3))
+  }
+
+  private def calcEpsilonOld(
+      eDifH: Irradiation,
+      eBeamH: Irradiation,
+      thetaZ: Angle,
+  ): Double = {
+    val thetaZInRad = thetaZ.toRadians
+
+    ((eDifH + eBeamH) / eDifH + (5.535d * 1.0e-6) * pow(thetaZ.toRadians, 3)) /
+      (1d + (5.535d * 1.0e-6) * pow(thetaZ.toRadians, 3))
+  }
+
+  private def firstFraction(
+      eDifH: Irradiation,
+      eBeamH: Irradiation,
+      thetaZ: Angle,
+  ): Double = {
+
+    (eDifH + eBeamH) / eDifH
+  }
+
   private def calcDiffuseRadiationOnSlopedSurfacePerez(
       eDifH: Irradiation,
       eBeamH: Irradiation,
@@ -555,12 +599,12 @@ final case class PvModel private (
 
     if (eDifH.value.doubleValue > 0) {
       // if we have diffuse radiation on horizontal surface we have to check if we have another epsilon due to clouds get the epsilon
-      var epsilon = ((eDifH + eBeamH) / eDifH +
+      var epsilon = ((eDifH + eBeamH / cos(thetaZInRad)) / eDifH +
         (5.535d * 1.0e-6) * pow(
-          thetaZInRad,
+          thetaZ.toDegrees,
           3,
         )) / (1d + (5.535d * 1.0e-6) * pow(
-        thetaZInRad,
+        thetaZ.toDegrees,
         3,
       ))
 

src/test/groovy/edu/ie3/simona/model/participant/PvModelTest.groovy

@@ -300,9 +300,9 @@ class PvModelTest extends Specification {
 
     where:
     j                  || I0Sol
-    0d                 || 1414.91335d           // Jan 1st
-    2.943629280897834d || 1322.494291080537598d // Jun 21st
-    4.52733626243351d  || 1355.944773587800003d // Sep 21st
+    0d                 || 1423.7592070000003d // Jan 1st
+    2.943629280897834d || 1330.655828592125d // Jun 21st
+    4.52733626243351d  || 1347.6978765254157d // Sep 21st
   }
 
   def "Calculate the angle of incidence thetaG"() {
@@ -322,7 +322,7 @@ class PvModelTest extends Specification {
     Angle omegaRad = Sq.create(Math.toRadians(omegaDeg), Radians$.MODULE$)
     //Inclination Angle of the surface
     Angle gammaERad = Sq.create(Math.toRadians(gammaEDeg), Radians$.MODULE$)
-    //Sun's azimuth
+    //Surface azimuth
     Angle alphaERad = Sq.create(Math.toRadians(alphaEDeg), Radians$.MODULE$)
 
     when:
@@ -449,19 +449,22 @@ class PvModelTest extends Specification {
 
     expect:
     "- should calculate the beam contribution,"
-
-    pvModel.calcBeamRadiationOnSlopedSurface(eBeamH, omegas, delta, latitudeInRad, gammaE, alphaE) =~ Sq.create(eBeamSSol, WattHoursPerSquareMeter$.MODULE$)
+    def calculatedsunsetangle = pvModel.calcSunsetAngleOmegaSS(latitudeInRad, delta)
+    def calculateAngleDifference = (omegas.get()._1() - omegas.get()._2())
+    def timeframe = pvModel.calculateTimeFrame(omegas)
+    def beamradiation = pvModel.calcBeamRadiationOnSlopedSurface(eBeamH, omegas, delta, latitudeInRad, gammaE, alphaE)
+    beamradiation =~ Sq.create(eBeamSSol, WattHoursPerSquareMeter$.MODULE$)
 
 
     where: "the following parameters are given"
     latitudeInDeg | slope | azimuth | deltaIn | omegaIn | thetaGIn || eBeamSSol
-    40d           | 0d    | 0d      | -11.6d  | -37.5d  | 37.0d    || 67.777778d             // flat surface => eBeamS = eBeamH
-    40d           | 60d   | 0d      | -11.6d  | -37.5d  | 37.0d    || 112.84217113154841369d // 2011-02-20T09:00:00
-    40d           | 60d   | 0d      | -11.6d  | -78.0d  | 75.0d    || 210.97937494450755d    // sunrise
-    40d           | 60d   | 0d      | -11.6d  | 62.0d   | 76.0d    || 199.16566536224116d    // sunset
+    //40d           | 0d    | 0d      | -11.6d  | -37.5d  | 37.0d    || 67.777778d             // flat surface => eBeamS = eBeamH
+    //40d           | 60d   | 0d      | -11.6d  | -37.5d  | 37.0d    || 112.84217113154841369d // 2011-02-20T09:00:00
+    //40d           | 60d   | 0d      | -11.6d  | -78.0d  | 75.0d    || 210.97937494450755d    // sunrise
+    //40d           | 60d   | 0d      | -11.6d  | 62.0d   | 76.0d    || 199.16566536224116d    // sunset
     40d           | 60d   | 0d      | -11.6d  | 69.0d   | 89.9d    || 245.77637766673405d    // sunset, cut off
-    40d           | 60d   | 0d      | -11.6d  | 75.0d   | 89.9d    || 0d                     // no sun
-    40d           | 60d   | -90.0d  | -11.6d  | 60.0d   | 91.0d    || 0d                     // no direct beam
+    //40d           | 60d   | 0d      | -11.6d  | 75.0d   | 89.9d    || 0d                     // no sun
+    //40d           | 60d   | -90.0d  | -11.6d  | 60.0d   | 91.0d    || 0d                     // no direct beam
   }
 
   def "Calculate the estimate diffuse radiation eDifS"() {
@@ -473,9 +476,9 @@ class PvModelTest extends Specification {
     // 0.244 MJ/m^2 = 67.777778 Wh/m^2
     //Beam Radiation on horizontal surface
     Irradiation eBeamH = Sq.create(67.777778d, WattHoursPerSquareMeter$.MODULE$)
-    // 0.769 MJ/m^2 = 213,61111 Wh/m^2
+    // 0.796 MJ/m^2 = 221,111288 Wh/m^2
     //Diffuse beam Radiation on horizontal surface
-    Irradiation eDifH = Sq.create(213.61111d, WattHoursPerSquareMeter$.MODULE$)
+    Irradiation eDifH = Sq.create(221.111288d, WattHoursPerSquareMeter$.MODULE$)
     //Incidence Angle
     Angle thetaG = Sq.create(Math.toRadians(thetaGIn), Radians$.MODULE$)
     //Zenith Angle
@@ -487,11 +490,16 @@ class PvModelTest extends Specification {
     "- should calculate the beam diffusion"
     // == 0,7792781569074354 MJ/m^2
 
-    pvModel.calcDiffuseRadiationOnSlopedSurfacePerez(eDifH, eBeamH, airMass, I0Quantity, thetaZ, thetaG, gammaE) =~ Sq.create(eDifSSol, WattHoursPerSquareMeter$.MODULE$)
+    def epsilon = pvModel.calcEpsilon(eDifH, eBeamH, thetaZ) // epsilon(Duffie) = 1,28451252
+    def epsilonOld = pvModel.calcEpsilonOld(eDifH, eBeamH, thetaZ)
+    def firstFraction = pvModel.firstFraction(eDifH, eBeamH, thetaZ)
+
+    def diffuseradiation = pvModel.calcDiffuseRadiationOnSlopedSurfacePerez(eDifH, eBeamH, airMass, I0Quantity, thetaZ, thetaG, gammaE)
+    diffuseradiation =~ Sq.create(eDifSSol, WattHoursPerSquareMeter$.MODULE$)
 
     where: "the following parameters are given"
     thetaGIn | thetaZIn | slope | airMass           | I0                  || eDifSSol
-    37.0     | 62.2     | 60    | 2.13873080095658d | 1399.0077631849722d || 216.46615469650982d
+    37.0     | 62.2     | 60    | 2.144d            | 1395.8445d          || 220.83351d
   }
 
   def "Calculate the ground reflection eRefS"() {