Evapotranspiration.java

public class Evapotranspiration {
   static public double calculate ( //
      int dayOfYear, // 1-366
      double tempMinC, //
      double tempMaxC, //
      double relativeHumidityMin, //
      double relativeHumidityMax, //
      double windSpeedAvgMs, // At 2 meters above the ground.
      double solarRadiationAvgWM2, //
      double albedo, // 0.23 for grass.
      double elevationM, //
      double latitudeDegrees //
   ) {
      // Calculations from:
      // Step by Step Calculation of the Penman-Monteith Evapotranspiration (FAO-56 Method)
      // Lincoln Zotarelli, Michael D. Dukes, Consuelo C. Romero, Kati W. Migliaccio, and Kelly T. Morgan
      // http://n4te.com/x/7163-ae45900.pdf (was https://edis.ifas.ufl.edu/pdffiles/AE/AE45900.pdf)
      // https://github.com/hwestbrook/evapotranspiration_calculator

      // Step 1: average daily air temperature C
      var dailyAirTemperatureAvgC = (tempMaxC + tempMinC) / 2;

      // Step 2: average solar radiation MJ
      var solarRadiationAvgMJM2 = WtoMJ(solarRadiationAvgWM2);

      // Step 3: avg wind speed m/s at 2m height

      // Step 4: slope of saturation vapor pressure curve
      var saturationVaporPressureCurveSlope = saturationVaporPressureCurveSlope(dailyAirTemperatureAvgC);

      // Step 5: atmospheric pressure
      var atmosphericPressue = atmosphericPressue(elevationM);

      // Step 6: psycometric constant
      var psychrometricConstant = psychrometricConstant(atmosphericPressue);

      // Step 7: delta term (DT)
      var deltaTerm = deltaTerm(saturationVaporPressureCurveSlope, psychrometricConstant, windSpeedAvgMs);

      // Step 8: psi term (PT)
      var psiTerm = psiTerm(saturationVaporPressureCurveSlope, psychrometricConstant, windSpeedAvgMs);

      // Step 9: temperature term (TT)
      var temperatureTerm = temperatureTerm(dailyAirTemperatureAvgC, windSpeedAvgMs);

      // Step 10: average saturation vapor pressure curve
      var saturationVaporPressureAvg = (saturationVapor(tempMaxC) + saturationVapor(tempMinC)) / 2;

      // Step 11: actual vapor pressure
      var saturationVaporPressureActual = saturationVaporPressureActual(tempMinC, tempMaxC, relativeHumidityMin,
         relativeHumidityMax);

      // Step 11.1: vapor pressure deficit
      var saturationVaporPressureDeficit = saturationVaporPressureAvg - saturationVaporPressureActual;

      // Step 12.1: relative sun earth difference
      var relativeEarthSunDifference = relativeEarthSunDifference(dayOfYear);

      // Step 12.2: relative sun earth difference
      var solarDeclination = solarDeclination(dayOfYear);

      // Step 13: latitude degrees to radians
      var latitudeRadians = latitudeDegrees * Math.PI / 180;

      // Step 14: sunset hour angle
      var sunsetHourAngle = sunsetHourAngle(latitudeRadians, solarDeclination);

      // Step 15: extraterrestrial radiation
      var extraterrestrialRadiation = extraterrestrialRadiation(relativeEarthSunDifference, sunsetHourAngle, latitudeRadians,
         solarDeclination);

      // Step 16: clear sky solar radiation
      var clearSkySolarRadiation = clearSkySolarRadiation(elevationM, extraterrestrialRadiation);

      // Step 17: clear sky solar radiation
      var netSolarRadiation = (1 - albedo) * solarRadiationAvgMJM2;

      // Step 18: net outgoing long wave solar radiation
      var netOutgoingLongWaveSolarRadiation = netOutgoingLongWaveSolarRadiation(tempMinC, tempMaxC, saturationVaporPressureActual,
         solarRadiationAvgMJM2, clearSkySolarRadiation);

      // Step 19: net radiation in mm
      var netRadiationMM = (netSolarRadiation - netOutgoingLongWaveSolarRadiation) * 0.408;

      // Final step: evapotranspiration in mm
      var radiationTerm = deltaTerm * netRadiationMM;
      var windTerm = psiTerm * temperatureTerm * (saturationVaporPressureAvg - saturationVaporPressureActual);
      return radiationTerm + windTerm;
   }

   static private final double CtoK (double c) {
      return c + 273.15;
   }

   static private final double WtoMJ (double mj) {
      return mj * 0.0864;
   }

   static private final double saturationVapor (double temperature) {
      return 0.6108 * Math.exp(17.27 * temperature / (temperature + 237.3));
   }

   static private final double saturationVaporPressureCurveSlope (double temperature) {
      var top = 4098 * saturationVapor(temperature);
      var bottom = Math.pow(temperature + 237.3, 2);
      return top / bottom;
   }

   static private final double atmosphericPressue (double elevationM) {
      var inner = (293 - 0.0065 * elevationM) / 293;
      return 101.3 * Math.pow(inner, 5.26);
   }

   static private final double psychrometricConstant (double atmosphericPressue) {
      return 0.000665 * atmosphericPressue;
   }

   static private final double deltaTerm (double saturationVaporPressureCurveSlope, double psychrometricConstant,
      double avgWindSpeedMs) {
      var top = saturationVaporPressureCurveSlope;
      var bottom = saturationVaporPressureCurveSlope + psychrometricConstant * (1 + 0.34 * avgWindSpeedMs);
      return top / bottom;
   }

   static private final double psiTerm (double saturationVaporPressureCurveSlope, double psychrometricConstant,
      double avgWindSpeedMs) {
      var top = psychrometricConstant;
      var bottom = saturationVaporPressureCurveSlope + psychrometricConstant * (1 + 0.34 * avgWindSpeedMs);
      return top / bottom;
   }

   static private final double temperatureTerm (double dailyAirTemperatureAvgC, double avgWindSpeedMs) {
      return 900 / CtoK(dailyAirTemperatureAvgC) * avgWindSpeedMs;
   }

   static private final double saturationVaporPressureActual (double minTempC, double maxTempC, double minRelativeHumidity,
      double maxRelativeHumidity) {
      var rel1 = saturationVapor(minTempC) * (maxRelativeHumidity / 100);
      var rel2 = saturationVapor(maxTempC) * (minRelativeHumidity / 100);
      return (rel1 + rel2) / 2;
   }

   static private final double relativeEarthSunDifference (int dayOfYear) {
      return 1 + 0.033 * Math.cos(2 * Math.PI / 365 * dayOfYear);
   }

   static private final double solarDeclination (int dayOfYear) {
      return 0.409 * Math.sin(2 * Math.PI / 365 * dayOfYear - 1.39);
   }

   static private final double sunsetHourAngle (double latitudeRadians, double solarDeclination) {
      return Math.acos(-1 * Math.tan(latitudeRadians) * Math.tan(solarDeclination));
   }

   static private final double extraterrestrialRadiation (double relativeEarthSunDifference, double sunsetHourAngle,
      double latitudeRadians, double solarDeclination) {
      var rel1 = 24 * 60 / Math.PI;
      var rel2 = 0.0820 * relativeEarthSunDifference;
      var rel3 = sunsetHourAngle * Math.sin(latitudeRadians) * Math.sin(solarDeclination)
         + Math.cos(latitudeRadians) * Math.cos(solarDeclination) * Math.sin(sunsetHourAngle);
      return rel1 * rel2 * rel3;
   }

   static private final double clearSkySolarRadiation (double elevationM, double extraterrestrialRadiation) {
      return (0.75 + 2 * Math.pow(10, -5) * elevationM) * extraterrestrialRadiation;
   }

   static private final double netOutgoingLongWaveSolarRadiation (double minTempC, double maxTempC,
      double saturationVaporPressureActual, double solarRadiationAvgMJM2, double clearSkySolarRadiation) {
      var rel1 = 4.903 * Math.pow(10, -9);
      var rel2 = (Math.pow(CtoK(maxTempC), 4) + Math.pow(CtoK(minTempC), 4)) / 2;
      var rel3 = 0.34 - 0.14 * Math.sqrt(saturationVaporPressureActual);
      var rel4 = 1.35 * solarRadiationAvgMJM2 / clearSkySolarRadiation - 0.35;
      return rel1 * rel2 * rel3 * rel4;
   }

   static public void main (String[] args) throws Throwable {
      System.out.println(Evapotranspiration.calculate(141, 75, 84, 0.63, 0.74, 3.61, 296, 0.23, 4.26, 18.450574));
   }
}