bugfix

Extraterrestrial radiation should depend on latitude

tree-water-demand/src/main/java/de/hft/stuttgart/water/TreeWaterDemand.java

@@ -49,12 +49,15 @@ public class TreeWaterDemand {
 			}
 		}
 
-		double longitude = 9.933333;
+		// Würzburg
+		double longitude = 9.9;
+		double latitude = 49.8;
 		// germany = 345
+		// Longitude of middle of time zone.
 		double lz = 345;
-		// Würzburg =
+		// West_positive_longitude
 		double lm = 360 - longitude;
-		double height = 177;
+		double altitude = 177;
 
 		double svf = 1;
 		double advection = 1.6;
@@ -101,10 +104,10 @@ public class TreeWaterDemand {
 				int dayOfYear = counter.getDayOfYear();
 				double hourOfDay = counter.getHour() + 0.5;
 
-				double extraterrestrialRadiation = ra(dayOfYear, hourOfDay, lz, lm);
-				double netLongwaveRadiation = rnl(temp, actualVaporPressure, incomingRadiation, height, extraterrestrialRadiation);
+				double extraterrestrialRadiation = ra(dayOfYear, hourOfDay, lz, lm, latitude);
+				double netLongwaveRadiation = rnl(temp, actualVaporPressure, incomingRadiation, altitude, extraterrestrialRadiation);
 				double netRadiation = rn(incomingRadiation, netLongwaveRadiation);
-				double et0 = et0(temp, actualVaporPressure, incomingRadiation, height, windSpeed, humidity, netRadiation);
+				double et0 = et0(temp, actualVaporPressure, incomingRadiation, altitude, windSpeed, humidity, netRadiation);
 				double et0u = et0u(svf, advection, et0);
 
 
@@ -112,7 +115,7 @@ public class TreeWaterDemand {
 				bw.write(String.format(";%.1f", windSpeed));
 				bw.write(String.format(";%.5f", saturationVaporPressure));
 				bw.write(String.format(";%.5f", actualVaporPressure));
-				bw.write(String.format(";%.5f", gamma(height)));
+				bw.write(String.format(";%.5f", gamma(altitude)));
 				bw.write(String.format(";%.5f", extraterrestrialRadiation));
 				bw.write(String.format(";%.5f", netRadiation));
 				bw.write(String.format(";%.3f", et0));
@@ -158,9 +161,9 @@ public class TreeWaterDemand {
 		return svf * advection * et0;
 	}
 
-	private static double et0(double temperature, double actualVaporPressure, double shortWaveRadiation, double height,
+	private static double et0(double temperature, double actualVaporPressure, double shortWaveRadiation, double altitude,
 			double windSpeed, double humidity, double netRadiation) {
-		double gamma = gamma(height);
+		double gamma = gamma(altitude);
 		double saturationVaporPressure = saturationVaporPressure(temperature);
 		double delta = delta(temperature, saturationVaporPressure);
 		double g = g(netRadiation, shortWaveRadiation > 0);
@@ -177,8 +180,8 @@ public class TreeWaterDemand {
 		return 0.6108 * Math.exp((17.27 * temperature) / (temperature + 237.3));
 	}
 
-	private static double gamma(double height) {
-		double pressure = 101.3 * Math.pow((293 - 0.0065 * height) / 293, 5.26);
+	private static double gamma(double altitude) {
+		double pressure = 101.3 * Math.pow((293 - 0.0065 * altitude) / 293, 5.26);
 		return 0.665 * 0.001 * pressure;
 	}
 
@@ -199,12 +202,12 @@ public class TreeWaterDemand {
 
 	// NetLongwaveRadiation
 	// unit is apparently [MJ / (m².h)]
-	private static double rnl(double temperature, double actualVaporPressure, double shortWaveRadiation, double height,
+	private static double rnl(double temperature, double actualVaporPressure, double shortWaveRadiation, double altitude,
 			double extraterrestrialRadiation){
 		// divide by 24 for hourly
 		double sigma = 4.903E-9 / 24;
 		double rs = shortWaveRadiation;
-		double rso = (0.75 + 2 * 10E-5 * height) * extraterrestrialRadiation;
+		double rso = (0.75 + 2 * 10E-5 * altitude) * extraterrestrialRadiation;
 		double rsRsoQuotient = 0.5;
 		if (rso != 0) {
 			rsRsoQuotient = rs / rso;
@@ -217,12 +220,12 @@ public class TreeWaterDemand {
 	}
 
 	// ExtraterrestrialRadiation
-	private static double ra(int dayOfYear, double hourOfDay, double lz, double lm) {
+	private static double ra(int dayOfYear, double hourOfDay, double lz, double lm, double latitude) {
 		double dr = 1 + 0.033 * Math.cos(2 * Math.PI * dayOfYear / 365d);
 		double b = (2 * Math.PI * (dayOfYear - 81)) / 364;
 		double sc = 0.1645 * Math.sin(2 * b) - 0.1255 * Math.cos(b) - 0.025 * Math.sin(b);
 		double omega = Math.PI / 12 * ((hourOfDay + 0.06667 * (lz - lm) + sc) - 12);
-		double phi = lm * Math.PI / 180;
+		double phi = latitude * Math.PI / 180;
 		double delta = 0.409 * Math.sin(2 * Math.PI / 365 * dayOfYear - 1.39);
 		double omegas = Math.acos(-Math.tan(phi) * Math.tan(delta));
 		double ra = 0;