Some fixes in PvModel and test

docs/readthedocs/_static/bibliography/bibtexAll.bib

@@ -12,12 +12,6 @@ @Article{Maleki.2017
 DOI = {10.3390/en10010134}
 }
 
-@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/}}
-}
-
 @article{Spencer.1971,
   added-at = {2018-06-18T21:23:34.000+0200},
   author = {Spencer, J. W.},
@@ -73,7 +67,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 +86,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 +117,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

@@ -127,7 +127,7 @@ $$
 **References:**
 
 * {cite:cts}`Maleki.2017`
-* {cite:cts}`Itaca_Sun`
+* {cite:cts}`Duffie.2013` p. 17 (formula 1.6.10)
 
 
 ### Solar Altitude Angle
@@ -146,7 +146,7 @@ $$
 **References:**
 
 * {cite:cts}`Maleki.2017` p. 5
-* {cite:cts}`Itaca_Sun`
+* {cite:cts}`Duffie.2013` p. 15 (formula 1.6.5) with $\sin (\alpha_s) = \cos (\theta_z)$
 
 
 ### Zenith Angle
@@ -198,11 +198,11 @@ $$
 Calculating the air mass ratio by dividing the radius of the earth with approx. effective height of the atmosphere (each in kilometer)
 
 $$
-airmassratio = (\frac{6371 km}{9 km}) = 707.8\overline{8}
+\mathrm{airmassratio} = (\frac{6371 km}{9 km}) = 707.8\overline{8}
 $$
 
 $$
-airmass = \sqrt{(707.8\overline{8} \cdot \cos({\theta_z}))^2 +2 \cdot 707.8\overline{8} +1)} - 707.8\overline{8} \cdot \cos{(\theta_z)})
+\mathrm{airmass} = \sqrt{(707.8\overline{8} \cdot \cos({\theta_z}))^2 +2 \cdot 707.8\overline{8} +1)} - 707.8\overline{8} \cdot \cos{(\theta_z)})
 $$
 
 **References:**
@@ -234,6 +234,8 @@ $$
 
 * {cite:cts}`Zheng.2017` p. 53, formula 2.3b
 * {cite:cts}`Iqbal.1983`
+* {cite:cts}`Spencer.1971`
+* {cite:cts}`Duffie.2013` (the fifth ed. seems to have a typo in formula (1.4.1b): factor $0.000719$ is missing a zero)
 
 
 ### Beam Radiation on Sloped Surface
@@ -274,7 +276,7 @@ b = (\cos(\phi) \cdot \cos(\delta)) \cdot (\sin(\omega_{2}) - \sin(\omega_{1}))
 $$
 
 $$
-E_{beam,S} = E_{beam,H} \cdot \frac{a}{b}
+E_{\mathrm{beam},S} = E_{\mathrm{beam},H} \cdot \frac{a}{b}
 $$
 
 **Please note:** $\frac{1}{180}\pi$ is omitted from these formulas, as we are already working with data in *radians*.
@@ -286,7 +288,7 @@ $$
 **$\omega_1$** = hour angle $\omega$\
 **$\omega_2$** = hour angle $\omega$ + 1 hour\
 **$\alpha_e$** = surface azimuth angle\
-**$E_{beam,H}$** = beam radiation (horizontal surface)
+**$E_{\mathrm{beam},H}$** = beam radiation (horizontal surface)
 
 **Reference:**
 
@@ -300,13 +302,15 @@ The diffuse radiation is computed using the Perez model, which divides the radia
 A cloud index is defined by
 
 $$
-\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}
+\epsilon = \frac{\frac{E_{\mathrm{dif},H} + E_{\mathrm{beam},N}}{E_{\mathrm{dif},H}} + 5.535 \cdot 10^{-6} \cdot \theta_{z}^3}{1 + 5.535 \cdot 10^{-6} \cdot \theta_{z}^3}
 $$
 
+with angle $\theta_z$ values in **degrees** ({cite:cts}`Duffie.2013` p. 94) and $E_{\mathrm{beam},N} = \frac{E_{\mathrm{beam},H}}{\cos (\theta_z)}$ ({cite:cts}`Duffie.2013` p. 95).
+
 Calculating a brightness index
 
 $$
-\Delta = m \cdot \frac{E_{dif,H}}{I_{0}}
+\Delta = m \cdot \frac{E_{\mathrm{dif},H}}{I_{0}}
 $$
 
 **Perez Fij coefficients (Myers 2017):**
@@ -384,7 +388,7 @@ $$
 the diffuse radiation can be calculated:
 
 $$
-E_{dif,S} = E_{dif,H} \cdot (\frac{1}{2} \cdot (1 +
+E_{\mathrm{dif},S} = E_{\mathrm{dif},H} \cdot (\frac{1}{2} \cdot (1 +
 cos(\gamma_{e})) \cdot (1- F_{1}) + \frac{a}{b} \cdot F_{1} +
 F_{2} \cdot \sin(\gamma_{e}))
 $$
@@ -396,25 +400,27 @@ $$
 **$\gamma_{e}$** = slope angle of the surface\
 **$I_{0}$** = Extraterrestrial Radiation\
 **$m$** = air mass\
-**$E_{beam,H}$** = beam radiation (horizontal surface)\
-**$E_{dif,H}$** = diffuse radiation (horizontal surface)
+**$E_{\mathrm{beam},H}$** = beam radiation (horizontal surface)\
+**$E_{\mathrm{beam},N}$** = beam radiation (normal incidence, thus radiation on a plane normal to the direction of the beam)\
+**$E_{\mathrm{dif},H}$** = diffuse radiation (horizontal surface)
 
 **References:**
 
 * {cite:cts}`Perez.1987`
 * {cite:cts}`Perez.1990`
 * {cite:cts}`Myers.2017` p. 96f
+* {cite:cts}`Duffie.2013` p. 95f
 
 
 ### Reflected Radiation on Sloped Surface
 
 $$
-E_{ref,S} = E_{Ges,H} \cdot \frac{\rho}{2} \cdot (1-
+E_{\mathrm{ref},S} = E_{\mathrm{Ges},H} \cdot \frac{\rho}{2} \cdot (1-
 \cos(\gamma_{e}))
 $$
 
 *with*\
-**$E_{Ges,H}$** = total horizontal radiation ($E_{beam,H} + E_{dif,H})$\
+**$E_{\mathrm{Ges},H}$** = total horizontal radiation ($E_{\mathrm{beam},H} + E_{\mathrm{dif},H})$\
 **$\gamma_e$** = slope angle of the surface\
 **$\rho$** = albedo
 
@@ -428,13 +434,13 @@ $$
 Received energy is calculated as the sum of all three types of irradiation.
 
 $$
-E_{total} = E_{beam,S} + E_{dif,S} + E_{ref,S}
+E_{\mathrm{total}} = E_{\mathrm{beam},S} + E_{\mathrm{dif},S} + E_{\mathrm{ref},S}
 $$
 
 *with*\
-**$E_{beam,S}$** = Beam radiation\
-**$E_{dif,S}$** = Diffuse radiation\
-**$E_{ref,S}$** = Reflected radiation
+**$E_{\mathrm{beam},S}$** = Beam radiation\
+**$E_{\mathrm{dif},S}$** = Diffuse radiation\
+**$E_{\mathrm{ref},S}$** = Reflected radiation
 
 A generator correction factor (depending on month surface slope $\gamma_{e}$) and a temperature correction factor (depending on month) multiplied on top.
 

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

@@ -24,7 +24,6 @@ import tech.units.indriya.unit.Units._
 
 import java.time.ZonedDateTime
 import java.util.UUID
-import java.util.stream.IntStream
 import scala.math._
 
 final case class PvModel private (
@@ -545,20 +544,21 @@ final case class PvModel private (
     val gammaEInRad = gammaE.toRadians
 
     // == brightness index beta  ==//
-    val beta = eDifH * airMass / extraterrestrialRadiationI0
+    val delta = eDifH * airMass / extraterrestrialRadiationI0
 
     // == cloud index epsilon  ==//
-    // if we have no clouds,  the epsilon bin is 8, as epsilon bin for an epsilon in [6.2, inf.[ = 8
-    var x = 8
+    val x = if (eDifH.value.doubleValue > 0) {
+      // if we have diffuse radiation on horizontal surface we have to consider
+      // the clearness parameter epsilon, which then gives us an epsilon bin x
 
-    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 +
+      // Beam radiation is required on a plane normal to the beam direction (normal incidence),
+      // thus dividing by cos theta_z
+      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,
       ))
 
@@ -579,18 +579,23 @@ final case class PvModel private (
         // get the corresponding bin
         val finalEpsilon = epsilon
 
-        x = IntStream
-          .range(0, discreteSkyClearnessCategories.length)
-          .filter((i: Int) =>
-            (finalEpsilon - discreteSkyClearnessCategories(i)(
-              0
-            ) >= 0) && (finalEpsilon - discreteSkyClearnessCategories(
-              i
-            )(1) < 0)
-          )
-          .findFirst
-          .getAsInt + 1
+        discreteSkyClearnessCategories.indices
+          .find { i =>
+            (finalEpsilon -
+              discreteSkyClearnessCategories(i)(0) >= 0) &&
+            (finalEpsilon -
+              discreteSkyClearnessCategories(i)(1) < 0)
+          }
+          .map(_ + 1)
+          .getOrElse(8)
+      } else {
+        // epsilon in [6.2, inf.[
+        8
       }
+    } else {
+      // if we have no clouds, the epsilon bin is 8,
+      // as the epsilon bin for an epsilon in [6.2, inf.[ is 8
+      8
     }
 
     // calculate the f_ij components based on the epsilon bin
@@ -604,19 +609,17 @@ final case class PvModel private (
     val f22 = 0.0012 * pow(x, 3) - 0.0067 * pow(x, 2) + 0.0091 * x - 0.0269
     val f23 = 0.0052 * pow(x, 3) - 0.0971 * pow(x, 2) + 0.2856 * x - 0.1389
 
-    // calculate circuumsolar brightness coefficient f1 and horizon brightness coefficient f2
-    val f1 = max(0, f11 + f12 * beta + f13 * thetaZInRad)
-    val f2 = f21 + f22 * beta + f23 * thetaZInRad
+    // calculate circumsolar brightness coefficient f1 and horizon brightness coefficient f2
+    val f1 = max(0, f11 + f12 * delta + f13 * thetaZInRad)
+    val f2 = f21 + f22 * delta + f23 * thetaZInRad
     val aPerez = max(0, cos(thetaGInRad))
     val bPerez = max(cos(1.4835298641951802), cos(thetaZInRad))
 
     // finally calculate the diffuse radiation on an inclined surface
     eDifH * (
-      ((1 + cos(
-        gammaEInRad
-      )) / 2) * (1 - f1) + (f1 * (aPerez / bPerez)) + (f2 * sin(
-        gammaEInRad
-      ))
+      ((1 + cos(gammaEInRad)) / 2) * (1 - f1) +
+        (f1 * (aPerez / bPerez)) +
+        (f2 * sin(gammaEInRad))
     )
   }
 

src/test/scala/edu/ie3/simona/agent/em/EmAgentIT.scala

@@ -259,11 +259,11 @@ class EmAgentIT
         scheduler.expectMessage(Completion(storageAgent))
 
         /* TICK 0
-         LOAD: 0.000269 MW
-         PV:  -0.005685 MW
+         LOAD: 0.269 kW
+         PV:  -5.617 kW
          STORAGE: SOC 0 %
          -> charge with 5 kW
-         -> remaining -0.0004161 MW
+         -> remaining -0.348 kW
          */
 
         emAgentActivation ! Activation(0)
@@ -272,7 +272,7 @@ class EmAgentIT
           0,
           weatherService.ref.toClassic,
           WeatherData(
-            WattsPerSquareMeter(400d),
+            WattsPerSquareMeter(540d),
             WattsPerSquareMeter(200d),
             Celsius(0d),
             MetersPerSecond(0d),
@@ -285,19 +285,19 @@ class EmAgentIT
             emResult.getInputModel shouldBe emInput.getUuid
             emResult.getTime shouldBe 0L.toDateTime
             emResult.getP should equalWithTolerance(
-              -0.000416087825.asMegaWatt
+              -0.00034885012.asMegaWatt
             )
-            emResult.getQ should equalWithTolerance(0.0000882855367.asMegaVar)
+            emResult.getQ should equalWithTolerance(0.00008828554.asMegaVar)
         }
 
         scheduler.expectMessage(Completion(emAgentActivation, Some(7200)))
 
         /* TICK 7200
-         LOAD: 0.000269 MW (unchanged)
-         PV:  -0.003797 MW
+         LOAD: 0.269 kW (unchanged)
+         PV:  -3.651 kW
          STORAGE: SOC 63.3 %
-         -> charge with 3.5282 kW
-         -> remaining 0 MW
+         -> charge with 3.382 kW
+         -> remaining 0 kW
          */
 
         emAgentActivation ! Activation(7200)
@@ -322,36 +322,36 @@ class EmAgentIT
             emResult.getQ should equalWithTolerance(0.0000882855367.asMegaVar)
         }
 
-        scheduler.expectMessage(Completion(emAgentActivation, Some(13107)))
+        scheduler.expectMessage(Completion(emAgentActivation, Some(13362)))
 
-        /* TICK 13107
-         LOAD: 0.000269 MW (unchanged)
-         PV:  -0.003797 MW (unchanged)
+        /* TICK 13362
+         LOAD: 0.269 kW (unchanged)
+         PV:  -3.651 kW (unchanged)
          STORAGE: SOC 100 %
          -> charge with 0 kW
-         -> remaining -0.003528 MW
+         -> remaining -3.382 kW
          */
 
-        emAgentActivation ! Activation(13107)
+        emAgentActivation ! Activation(13362)
 
         resultListener.expectMessageType[ParticipantResultEvent] match {
           case ParticipantResultEvent(emResult: EmResult) =>
             emResult.getInputModel shouldBe emInput.getUuid
-            emResult.getTime shouldBe 13107L.toDateTime
+            emResult.getTime shouldBe 13362L.toDateTime
             emResult.getP should equalWithTolerance(
-              -0.0035281545552.asMegaWatt
+              -0.003382375474.asMegaWatt
             )
             emResult.getQ should equalWithTolerance(0.0000882855367.asMegaVar)
         }
 
         scheduler.expectMessage(Completion(emAgentActivation, Some(14400)))
 
         /* TICK 14400
-         LOAD: 0.000269 MW (unchanged)
-         PV:  -0.000066 MW
+         LOAD: 0.269 kW (unchanged)
+         PV:  -0.066 kW
          STORAGE: SOC 100 %
          -> charge with -0.202956 kW
-         -> remaining 0 MW
+         -> remaining 0 kW
          */
 
         // send weather data before activation, which can happen