TreeWaterDemand.java

package de.hft.stuttgart.water;

import java.io.IOException;
import java.nio.file.Files;
import java.nio.file.Path;
import java.nio.file.Paths;
import java.time.Duration;
import java.time.LocalDateTime;
import java.time.Month;
import java.util.ArrayList;
import java.util.List;
import java.util.stream.Collectors;

import com.opencsv.CSVReader;
import com.opencsv.exceptions.CsvException;

public class TreeWaterDemand {

	/**
	 * Rs/Rso für Nacht = 0.5
	 * @param args
	 * @throws IOException
	 * @throws CsvException
	 */
	public static void main(String[] args) throws IOException, CsvException {
		List<String[]> csvValues;
		Path csvPath = Paths.get("Wuerzburg-hour-Meteonorm8.csv");
		try (CSVReader csvReader = new CSVReader(Files.newBufferedReader(csvPath))) {
			csvValues = csvReader.readAll();
		}

		int tempIndex = -1;
		int irrIndex = -1;
		int humidityIndex = -1;
		int windSpeedIndex = -1;
		String[] header = csvValues.get(1);
		for (int i = 0; i < header.length; i++) {
			if (header[i].startsWith("GHI (W/m^2)")) {
				irrIndex = i;
			} else if (header[i].startsWith("Dry-bulb (C)")) {
				tempIndex = i;
			} else if (header[i].startsWith("RHum (%)")) {
				humidityIndex = i;
			} else if (header[i].startsWith("Wspd (m/s)")) {
				windSpeedIndex = i;
			}
		}

		double longitude = 9.933333;
		// germany = 345
		double lz = 345;
		// Würzburg =
		double lm = 360 - longitude;
		double height = 177;

		double svf = 1;
		double advection = 1.6;
		double tr = 0.9;

		LocalDateTime startOfYear = LocalDateTime.of(2005, 1, 1, 0, 0);
		
		LocalDateTime start = LocalDateTime.of(2005, 1, 1, 0, 0);
		LocalDateTime end = LocalDateTime.of(2006, 1, 1, 0, 0);
		
//		LocalDateTime start = LocalDateTime.of(2005, 4, 1, 0, 0);
//		LocalDateTime end = LocalDateTime.of(2005, 10, 1, 0, 0);
		
		List<Double> et0us = new ArrayList<>();
		List<Double> et0s = new ArrayList<>();
		List<Double> demands = new ArrayList<>();
		List<Double> rain = new ArrayList<>();
		LocalDateTime counter = start;
		System.out.println("Date,ETIa (L/m²/d),ET0 (L/m²/d)");
		int startHours = (int) Duration.between(startOfYear, start).toHours();
		int hours = (int) Duration.between(start, end).toHours();
		for (int i = startHours; i < startHours + hours; i++) {
			String[] row = csvValues.get(i + 2);
			double temp = Double.parseDouble(row[tempIndex]);
			double humidity = Double.parseDouble(row[humidityIndex]);
			double incomingRadiation = Double.parseDouble(row[irrIndex]);
			// convert to MJ/h - 1 W = 0.0036 MJ/h
			incomingRadiation *= 0.0036;
			double windSpeed = Double.parseDouble(row[windSpeedIndex]);
			double saturationVaporPressure = e0(temp);
			double actualVaporPressure = ea(saturationVaporPressure, humidity);
			int dayOfYear = counter.getDayOfYear();
			double hourOfDay = counter.getHour() + 0.5;
			double et0 = et0(temp, actualVaporPressure, incomingRadiation, height, windSpeed, humidity, dayOfYear,
					hourOfDay, lz, lm);
			double et0u = et0u(svf, advection, et0);
			
			et0s.add(et0);
//			System.out.println("et0: " + et0);
//			System.out.println("etia: " + etia);
			et0us.add(et0u);
			counter = counter.plusHours(1);
		}
		
//		System.out.println(rain.stream().collect(Collectors.summarizingDouble(Double::doubleValue)));
		
		Double et0Sum = et0s.stream().collect(Collectors.summingDouble(Double::doubleValue));
		System.out.println("Summe aller ET0: " + et0Sum);
		
		
		Double et0uSum = et0us.stream().collect(Collectors.summingDouble(Double::doubleValue));
		System.out.println("Summe aller ET0u: " + et0uSum);
		double et0uFactor = (0.865 * Math.log10(1 / et0uSum) + 3.36);
		
		for (double et0u : et0us) {
			double etia = 0;
			if (et0u > 0) {
				etia = etia(tr, et0u, et0uFactor);
			}
			demands.add(etia);
		}

		System.out.println("Summe aller ETIa: " + demands.stream().collect(Collectors.summingDouble(Double::doubleValue)));

//		printHourly(start, et0s, demands);
		
//		printDaily(start, et0s, demands);

//		test();
	}

	private static void printHourly(LocalDateTime start, List<Double> et0s, List<Double> demands) {
		LocalDateTime date = start;
		for (int i = 0; i < demands.size(); i++) {
			System.out.println(date + "," + demands.get(i) + "," + et0s.get(i));
			date = date.plusHours(1);
		}
	}

	private static void printDaily(LocalDateTime start, List<Double> et0s, List<Double> demands) {
		LocalDateTime date = start;
		for (int i = 0; i < demands.size(); i = i + 24) {
			double demandSum = 0;
			double et0Sum = 0;
			for (int j = i; j < i + 24; j++) {
				demandSum += demands.get(j);
				et0Sum += et0s.get(j);
			}
			System.out.println(date + "," + demandSum + "," + et0Sum);
			date = date.plusHours(24);
		}
	}

	private static void test() {
		double lz = 15;
		double lm = 16.25;
		double height = 8;

		double temp = 38;
		double humidity = 52;
		double windSpeed = 3.3;
		double incomingRadiation = 2.45;

		double saturationVaporPressure = e0(temp);
		double actualVaporPressure = ea(saturationVaporPressure, humidity);
		int dayOfYear = 274;
		double hourOfDay = 14 + 0.5;
		double et0 = et0(temp, actualVaporPressure, incomingRadiation, height, windSpeed, humidity, dayOfYear,
				hourOfDay, lz, lm);
		System.out.println(et0);
	}

	private static double etia(double tr, double et0u, double et0uFactor) {
//		return 1 * 1.6 * et0u; 
		return tr * 1.25 * et0u * et0uFactor;
//		return tr * (1.61 * Math.log10(800) - 3.39) * et0u * (0.865 * Math.log10(1 / et0u) + 3.36);
	}

	private static double et0u(double svf, double advection, double et0) {
		return svf * advection * et0;
	}

	private static double et0(double temperature, double actualVaporPressure, double shortWaveRadiation, double height,
			double windSpeed, double humidity, int dayOfYear, double hourOfDay, double lz, double lm) {
		double nominalRadiation = rn(temperature, actualVaporPressure, shortWaveRadiation, height, dayOfYear, hourOfDay,
				lz, lm);
		double gamma = gamma(height);
		double saturationVaporPressure = e0(temperature);
		double delta = delta(temperature, saturationVaporPressure);
		double g = g(nominalRadiation, shortWaveRadiation > 0);
		double n = 0.408 * delta * (nominalRadiation - g)
				+ gamma * (37 / (temperature + 273)) * windSpeed * (saturationVaporPressure - actualVaporPressure);
		return n / (delta + gamma * (1 + 0.34  * windSpeed));
	}

	private static double ea(double saturationVaporPressure, double humidity) {
		return saturationVaporPressure * humidity / 100;
	}

	private static double e0(double temperature) {
		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);
		return 0.665 * 0.001 * pressure;
	}

	private static double g(double nominalRadiation, boolean daylight) {
		if (daylight) {
			return 0.1 * nominalRadiation;
		} else {
			return 0.5 * nominalRadiation;
		}
	}

	private static double rn(double temperature, double actualVaporPressure, double shortWaveRadiation, double height,
			int dayOfYear, double hourOfDay, double lz, double lm) {
		// 1 - 0.23: albedo effect of grass
		return (1 - 0.23) * shortWaveRadiation
				- rnl(temperature, actualVaporPressure, shortWaveRadiation, height, dayOfYear, hourOfDay, lz, lm);
	}

	private static double rnl(double temperature, double actualVaporPressure, double shortWaveRadiation, double height,
			int dayOfYear, double hourOfDay, double lz, double lm) {
		// divide by 24 for hourly
		double sigma = 4.903E-9 / 24;
		double rs = shortWaveRadiation;
		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 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;
		if (!(omega < -omegas || omega > omegas)) {
			double omega1 = omega - Math.PI * 1 / 24;
			double omega2 = omega + Math.PI * 1 / 24;
			ra = 12 * 60 / Math.PI * 0.0820 * dr * ((omega2 - omega1) * Math.sin(phi) * Math.sin(delta)
					+ Math.cos(phi) * Math.cos(delta) * (Math.sin(omega2) - Math.sin(omega1)));
		}
		double rso = (0.75 + 2 * 10E-5 * height) * ra;
		double rsRsoQuotient = 0.5;
		if (rso != 0) {
			rsRsoQuotient = rs / rso;
		}
		if (rsRsoQuotient > 1) {
			rsRsoQuotient = 1;
		}
		return sigma * Math.pow(temperature + 273.16, 4) * (0.34 - 0.14 * Math.sqrt(actualVaporPressure))
				* (1.35 * rsRsoQuotient - 0.35);
	}

	private static double delta(double temperatur, double e0) {
		return (4098 * e0) / Math.pow(temperatur + 237.3, 2);
	}

}