/*
 * To change this license header, choose License Headers in Project Properties.
 * To change this template file, choose Tools | Templates
 * and open the template in the editor.
 */
package d.dataNodes;

import d.util.RangeMap;
import d.util.SOLbuilder;
import java.io.File;
import java.text.DecimalFormat;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import soils.Component;
import soils.Horizon;
import soils.db.tables.TableHorizon;
import soils.db.tables.TableHorizonCalculations;
import soils.utils.SoilUtils;

/**
 *
 * @author sidereus
 */
public class DSSATData {

  DecimalFormat df1 = new DecimalFormat("#.###");
  private static String SLASH = "/";

  // Romero, C. C., Hoogenboom, G., Baigorria, G. A., Koo, J., Gijsman, A. J., & Wood, S. (2012).
  // Reanalysis of a global soil database for crop and environmental modeling.
  // Environmental Modelling & Software, 35, 163-170.
  // Tab 4, columns Hydrological Group and Runoff curve number (at different slope angles)
  private final Map<String, RangeMap<Double, Integer>> table_runoffCurveNumber = new HashMap() {
    {
      put("A", new RangeMap() {
        {
          // slope angle start, slope angle end, runoff curve number value
          put(0d, 2d, 61);
          put(2d, 5d, 64);
          put(5d, 10d, 68);
          put(10d, 100d, 71);
        }
      });
      put("B", new RangeMap() {
        {
          put(0d, 2d, 73);
          put(2d, 5d, 76);
          put(5d, 10d, 80);
          put(10d, 100d, 83);
        }
      });
      put("C", new RangeMap() {
        {
          put(0d, 2d, 81);
          put(2d, 5d, 84);
          put(5d, 10d, 88);
          put(10d, 100d, 91);
        }
      });
      put("D", new RangeMap() {
        {
          put(0d, 2d, 84);
          put(2d, 5d, 87);
          put(5d, 10d, 91);
          put(10d, 100d, 94);
        }
      });
    }
  };

  public void executeComputation(Component comp, String fileName, String id_soil, double longitude, double latitude) {

    List<String> firstSet = new ArrayList();
    List<String> secondSet = new ArrayList();

    double salb = albedo(comp);
    double slu1 = evaporationLimit(comp);
    double sldr = drainageRate(comp);
    double slro = runoffCurveNumber(comp);
    double slnf = 1.0;
    double slpf = 1.0;
    String smhb = "IB001";
    String smpx = "IB001";
    String smke = "IB001";
    String scom = "-99"; // Color, moist, Munsell hue

    horizons(comp, firstSet, secondSet);

    SOLbuilder fileBuilder = new SOLbuilder();
    fileBuilder.getSolFile(fileName,
        "USA", id_soil, latitude, longitude, "-99",
        scom, salb, slu1, sldr, slro, slnf, slpf, smhb, smpx, smke,
        firstSet.iterator(), secondSet.iterator());

  }

  private double albedo(Component comp) {
    double salb;
    if (Double.isNaN(comp.albedodry_r())) {
      comp.albedodry_r(SoilUtils.averageMissingValue(comp.albedodry_h(),
          comp.albedodry_l(),
          comp.albedodry_r()));
    }
    if (Double.isNaN(comp.albedodry_r())) {
      comp.albedodry_r(-99d);
    }
    salb = comp.albedodry_r();

    if (salb < 0.09 && salb > 0.17) {
      String msg = "Albedo value ranges from 0.09 for a black soil"
          + " to 0.17 for a yellow soil. Current value is: " + salb;
      throw new IllegalArgumentException(msg);
    }
    return salb;
  }

  // modified (upper limit 12) from https://github.com/ldemaz/rcropmod to adhere to
  // Romero, C. C., Hoogenboom, G., Baigorria, G. A., Koo, J., Gijsman, A. J., & Wood, S. (2012).
  // Reanalysis of a global soil database for crop and environmental modeling.
  // Environmental Modelling & Software, 35, 163-170.
  private double evaporationLimit(Component comp) {
    double slu1 = Double.NaN;
    for (Map.Entry<String, Horizon> horizon : comp.horizons.entrySet()) {
      TableHorizon h = horizon.getValue().getTableHorizon();
      computeAverageMissingValue(h);
      double clay = h.claytotal_r();
      slu1 = 5 + 0.15 * (100 - (98 - clay));
      if (slu1 > 12) {
        slu1 = 12;
      }
      break;
    }
    return slu1;
  }

  private double drainageRate(Component comp) {
    return parseDrainClass(comp.getTableComponent().drainagecl());
  }

  private double parseDrainClass(String drainClass) {
    String actualClass = drainClass.trim().toLowerCase().replace(" ", "");
    switch (actualClass) {
      case "verypoorlydrained":
        return 0.01;
      case "poorlydrained":
        return 0.05;
      case "somewhatpoorlydrained":
        return 0.25;
      case "moderatelywelldrained":
        return 0.40;
      case "welldrained":
        return 0.60;
      case "somewhatexcessivelydrained":
        return 0.75;
      case "excessivelydrained":
        return 0.85;
      default:
        throw new IllegalArgumentException();
    }
  }

  // if still NAN exit program
  private double runoffCurveNumber(Component comp) {
    if (Double.isNaN(comp.slope_r())) {
      comp.slope_r(SoilUtils.averageMissingValue(comp.slope_h(), comp.slope_l(), comp.slope_r()));
    }
    if (Double.isNaN(comp.slope_r())) {
      throw new IllegalArgumentException();
    }
    double slope = comp.slope_r();
    String hd = comp.hydgrp();
    // if soils in dual groups (A/D, B/D, C/D) get first soil type (drained condition)
    if (hd.contains(SLASH)) {
      hd = hd.split(SLASH)[0];
    }
    return table_runoffCurveNumber.get(hd).get(slope);
  }

  private void horizons(Component comp, List<String> firstSet, List<String> secondSet) {
    double baseLayer_depth = 0;

    for (Map.Entry<String, Horizon> horizon : comp.horizons.entrySet()) {
      TableHorizon th = horizon.getValue().getTableHorizon();
      TableHorizonCalculations cth = horizon.getValue().getTableHorizonCalculations();
      computeAverageMissingValue(th);
      double clayPerc = th.claytotal_r();
      double siltPerc = th.silttotal_r();
      double organic_matterPerc = th.om_r();
      double bulkDensity = th.dbthirdbar_r();

      if (!isMandatorysSoilPropertiesValid(clayPerc, siltPerc, organic_matterPerc, bulkDensity)) {
        break;
      }

      // depth, base of layer [cm]
      double slb = baseLayer_depth + th.hzdepb_r();
      firstSet.add(firstDataSet(th, cth, clayPerc, siltPerc, organic_matterPerc, slb, bulkDensity, baseLayer_depth));
      secondSet.add(secondDataSet(th, slb));
      baseLayer_depth = slb;

    }

  }

  private String firstDataSet(TableHorizon h, TableHorizonCalculations cth, double clayPerc, double siltPerc,
      double organic_matterPerc, double slb, double bulkDensity, double baseLayer_depth) {

    double clay = clayPerc / 100;
    double silt = siltPerc / 100;
    double organic_matter = organic_matterPerc / 100;
    StringBuilder tmpFirstString = new StringBuilder();
    tmpFirstString.append(spaceFormat(slb));

    // Master horizon
    double slmh = -99d;
    tmpFirstString.append(spaceFormat(slmh));

    // lower limit, or wilitng point [cm3/cm3]
    double slll = computeWiltingPoint(clay, silt, organic_matter);
    tmpFirstString.append(spaceFormat(slll));

    // drained upper limit, or field capacity [cm3/cm3]
    double sdul = computeFieldCapacity(clay, silt, organic_matter);
    tmpFirstString.append(spaceFormat(sdul));

    // upper limit, saturated [cm3/cm3]
    double ssat = computeSaturatedUpperLimit(sdul, clay, silt, organic_matter);
    tmpFirstString.append(spaceFormat(ssat));

    // root growth factor [0-1]
    double srgf = computeRootGrowth(baseLayer_depth, slb);
    tmpFirstString.append(spaceFormat(srgf));

    // Sat. hydraulic conductivity, macropore, [cm/h]
    double ssks = computeSatHydraulicConductivity(h);
    tmpFirstString.append(spaceFormat(ssks));

    // bulk density [g/cm3]
    double sbdm = bulkDensity;
    tmpFirstString.append(spaceFormat(sbdm));

    // soil organic carbon [%]
    // Saxton K., Rawls W. (2006) Soil water characteristic estimates by texture and organic matter for
    // hydrologic solutions. Soil Science Society of America Journal 70:1569-1578
    double sloc = organic_matterPerc / 2;
    tmpFirstString.append(spaceFormat(sloc));

    // Clay (<0.002 mm), [%]
    double slcl = clayPerc;
    tmpFirstString.append(spaceFormat(slcl));

    // Silt (0.05 to 0.002 mm), [%]
    double slsi = siltPerc;
    tmpFirstString.append(spaceFormat(slsi));

    // Coarse fraction (>2 mm), [%]
    double slcf = cth.fragvol_r();
    tmpFirstString.append(spaceFormat(slcf));

    // Total nitrogen, [%]
    double slni = -99d;
    tmpFirstString.append(spaceFormat(slni));

    // pH in water
    double slhw = h.ph1to1h2o_r();
    tmpFirstString.append(spaceFormat(slhw));

    // pH in buffer
    double slhb = h.ph01mcacl2_r();
    tmpFirstString.append(spaceFormat(slhb));

    // Cation exchange capacity, [cmol/kg]
    double scec = h.cec7_r();
    tmpFirstString.append(spaceFormat(scec));

    // ???
    double sadc = -99d;
    tmpFirstString.append(spaceFormat(sadc));

    return tmpFirstString.toString();

  }

  private String secondDataSet(TableHorizon tb, double slb) {
    StringBuilder tmpSecondString = new StringBuilder();
    tmpSecondString.append(spaceFormat(slb));

    // Phosphorus, extractable, [mg/kg]
    double slpx = tb.pbray1_r();
    tmpSecondString.append(spaceFormat(slpx));

    // Phosphorus, total, [mg/kg]
    double slpt = tb.ptotal_r();
    tmpSecondString.append(spaceFormat(slpt));

    // Phosphorus, organic, [mg/kg]
    double slpo = -99d;
    tmpSecondString.append(spaceFormat(slpo));

    // CaCO3 content, [g/kg]
    double caco3 = tb.caco3_r();
    tmpSecondString.append(spaceFormat(caco3));

    // Aluminum
    double slal = tb.extral_r();
    tmpSecondString.append(spaceFormat(slal));

    // Iron
    double slfe = tb.freeiron_r();
    tmpSecondString.append(spaceFormat(slfe));

    // Manganese
    double slmn = -99d;
    tmpSecondString.append(spaceFormat(slmn));

    // Base saturation, [cmol/kg]
    double slbs = -99d;
    tmpSecondString.append(spaceFormat(slbs));

    // Phosphorus isotherm A, [mmol/kg]
    double slpa = -99d;
    tmpSecondString.append(spaceFormat(slpa));

    // Phosphorus iostherm B, [mmol/l]
    double slpb = -99d;
    tmpSecondString.append(spaceFormat(slpb));

    // Potassium, exchangeable, [cmol/kg]
    double slke = -99d;
    tmpSecondString.append(spaceFormat(slke));

    // Magnesium, [cmol/kg]
    double slmg = -99d;
    tmpSecondString.append(spaceFormat(slmg));

    // Sodium, [cmol/kg]
    double slna = -99d;
    tmpSecondString.append(spaceFormat(slna));

    // Sulphur
    double slsu = -99d;
    tmpSecondString.append(spaceFormat(slsu));

    // Electric conductivity, seimen
    double slec = tb.ec_r();
    tmpSecondString.append(spaceFormat(slec));

    // Calcium, exchangeable, [cmol/kg]
    double slca = -99d;
    tmpSecondString.append(spaceFormat(slca));

    return tmpSecondString.toString();
  }

  private String spaceFormat(double value) {
    return String.format("%8s", df1.format(value));
  }

  private boolean isMandatorysSoilPropertiesValid(double clay, double silt,
      double organic_matter, double bulkDensity) {
    if (Double.compare(clay, -99d) == 0) {
      return false;
    } else if (Double.compare(silt, -99d) == 0) {
      return false;
    } else if (Double.compare(organic_matter, -99d) == 0) {
      return false;
    } else if (Double.compare(bulkDensity, -99d) == 0) {
      return false;
    } else {
      return true;
    }
  }

  // Saxton K., Rawls W. (2006) Soil water characteristic estimates by texture and organic matter for
  // hydrologic solutions. Soil Science Society of America Journal 70:1569-1578
  // table 1, eq. 1
  private double computeWiltingPoint(double clay, double silt, double organic_matter) {
    double tmp_slll = -0.024 * silt + 0.487 * clay + 0.006 * organic_matter
        + +0.005 * (silt * organic_matter) - 0.013 * (clay * organic_matter)
        + +0.068 * (silt * clay) + 0.031;

    return tmp_slll + (0.14 * tmp_slll - 0.02);
  }

  // Saxton K., Rawls W. (2006) Soil water characteristic estimates by texture and organic matter for
  // hydrologic solutions. Soil Science Society of America Journal 70:1569-1578
  // table 1, eq. 2
  private double computeFieldCapacity(double clay, double silt, double organic_matter) {
    double tmp_sdul = -0.251 * silt + 0.195 * clay + 0.011 * organic_matter
        + +0.006 * (silt * organic_matter) - 0.027 * (clay * organic_matter)
        + +0.452 * (silt * clay) + 0.299;

    return tmp_sdul + (1.283 * Math.pow(tmp_sdul, 2) - 0.374 * tmp_sdul - 0.015);
  }

  // Saxton K., Rawls W. (2006) Soil water characteristic estimates by texture and organic matter for
  // hydrologic solutions. Soil Science Society of America Journal 70:1569-1578
  // table 1, eq. 3 & 5
  private double computeSaturatedUpperLimit(double fieldCapacity, double clay, double silt, double organic_matter) {
    double ttmp_ssat = 0.278 * silt + 0.034 * clay + 0.022 * organic_matter
        + -0.018 * (silt * organic_matter) - 0.027 * (clay * organic_matter)
        + -0.584 * (silt * clay) + 0.078;
    double tmp_ssat = ttmp_ssat + (0.636 * ttmp_ssat - 0.107);

    return fieldCapacity + tmp_ssat - 0.097 * silt + 0.043;
  }

  private double computeRootGrowth(double top, double bottom) {
    if (bottom < 15) {
      return 1;
    } else {
      double midlayer = 0.5 * (bottom - top);
      return Math.exp(-0.02 * midlayer);
    }
  }

  private double computeSatHydraulicConductivity(TableHorizon h) {
    return h.ksat_r() * 0.36; // why 0.36???
  }

  private void computeAverageMissingValue(TableHorizon h) {
    if (Double.isNaN(h.claytotal_r())) {
      h.claytotal_r(SoilUtils.averageMissingValue(h.claytotal_h(), h.claytotal_l(), h.claytotal_r()));
    }
    if (Double.isNaN(h.claytotal_r())) {
      h.claytotal_r(-99d);
    }

    if (Double.isNaN(h.sandtotal_r())) {
      h.sandtotal_r(SoilUtils.averageMissingValue(h.sandtotal_h(), h.sandtotal_l(), h.sandtotal_r()));
    }
    if (Double.isNaN(h.sandtotal_r())) {
      h.sandtotal_r(-99d);
    }

    if (Double.isNaN(h.silttotal_r())) {
      h.silttotal_r(SoilUtils.averageMissingValue(h.silttotal_h(), h.silttotal_l(), h.silttotal_r()));
    }
    if (Double.isNaN(h.silttotal_r())) {
      h.silttotal_r(-99d);
    }

    if (Double.isNaN(h.om_r())) {
      h.om_r(SoilUtils.averageMissingValue(h.om_h(), h.om_l(), h.om_r()));
    }
    if (Double.isNaN(h.om_r())) {
      h.om_r(-99d);
    }

    if (Double.isNaN(h.ksat_r())) {
      h.ksat_r(SoilUtils.averageMissingValue(h.ksat_h(), h.ksat_l(), h.ksat_r()));
    }
    if (Double.isNaN(h.ksat_r())) {
      h.ksat_r(-99d);
    }

    if (Double.isNaN(h.hzdepb_r())) {
      h.hzdepb_r(SoilUtils.averageMissingValue(h.hzdepb_h(), h.hzdepb_l(), h.hzdepb_r()));
    }
    if (Double.isNaN(h.hzdepb_r())) {
      h.hzdepb_r(-99d);
    }

    if (Double.isNaN(h.dbthirdbar_r())) {
      h.dbthirdbar_r(SoilUtils.averageMissingValue(h.dbthirdbar_h(), h.dbthirdbar_l(), h.dbthirdbar_r()));
    }
    if (Double.isNaN(h.dbthirdbar_r())) {
      h.dbthirdbar_r(-99d);
    }

    if (Double.isNaN(h.caco3_r())) {
      h.caco3_r(SoilUtils.averageMissingValue(h.caco3_h(), h.caco3_l(), h.caco3_r()));
    }
    if (Double.isNaN(h.caco3_r())) {
      h.caco3_r(-99d);
    }

    if (Double.isNaN(h.ph1to1h2o_r())) {
      h.ph1to1h2o_r(SoilUtils.averageMissingValue(h.ph1to1h2o_h(), h.ph1to1h2o_l(), h.ph1to1h2o_r()));
    }
    if (Double.isNaN(h.ph1to1h2o_r())) {
      h.ph1to1h2o_r(-99d);
    }

    if (Double.isNaN(h.ph01mcacl2_r())) {
      h.ph01mcacl2_r(SoilUtils.averageMissingValue(h.ph01mcacl2_h(), h.ph01mcacl2_l(), h.ph01mcacl2_r()));
    }
    if (Double.isNaN(h.ph01mcacl2_r())) {
      h.ph01mcacl2_r(-99d);
    }

    if (Double.isNaN(h.cec7_r())) {
      h.cec7_r(SoilUtils.averageMissingValue(h.cec7_h(), h.cec7_l(), h.cec7_r()));
    }
    if (Double.isNaN(h.cec7_r())) {
      h.cec7_r(-99d);
    }

    if (Double.isNaN(h.pbray1_r())) {
      h.pbray1_r(SoilUtils.averageMissingValue(h.pbray1_h(), h.pbray1_l(), h.pbray1_r()));
    }
    if (Double.isNaN(h.pbray1_r())) {
      h.pbray1_r(-99d);
    }

    if (Double.isNaN(h.ptotal_r())) {
      h.ptotal_r(SoilUtils.averageMissingValue(h.ptotal_h(), h.ptotal_l(), h.ptotal_r()));
    }
    if (Double.isNaN(h.ptotal_r())) {
      h.ptotal_r(-99d);
    }

    if (Double.isNaN(h.extral_r())) {
      h.extral_r(SoilUtils.averageMissingValue(h.extral_h(), h.extral_l(), h.extral_r()));
    }
    if (Double.isNaN(h.extral_r())) {
      h.extral_r(-99d);
    }

    if (Double.isNaN(h.freeiron_r())) {
      h.freeiron_r(SoilUtils.averageMissingValue(h.freeiron_h(), h.freeiron_l(), h.freeiron_r()));
    }
    if (Double.isNaN(h.freeiron_r())) {
      h.freeiron_r(-99d);
    }
  }

}