EFH2HydrologyModel.java [src/java/m/efh2] Revision:   Date:
package m.efh2;

import m.efh2.EFH2RainfallCoefficients.RfC;
import java.util.List;

/**
 * Storm runoff model based on conventions in Engineering Field Handbook. This
 * uses a lookup method to compute total runoff, based on a number of measures
 * of the drained field.
 */
public class EFH2HydrologyModel {

  /**
   * storm type parameter, from values in EFH2RainfallCoefficients
   */
  private String stormType = "II";
  /**
   * total precipitation parameter, inches
   */
  private double precip = 1.0;
  /**
   * watershed length parameter, feet
   */
  private double watershedLength = 200.0;
  /**
   * watershed slope parameter, percent
   */
  private double watershedSlope = 0.1;
  /**
   * runoff curve number parameter, for lookups
   */
  private int runoffCurveNumber = 40;
  // outputs
  /**
   * runoff depth output, in
   */
  private double runoffQ = Double.NaN;
  /**
   * time of concentration, hours
   */
  private double Tc;
  /**
   * unit peak discharge output, cfs/ac-in
   */
  private double qu;


  public void simulate() {
    // lookup Q based on P and CN, using Table 2.2
    runoffQ = EFH2RunoffTable.lookupQ(precip, runoffCurveNumber);
    computeTc();
    computeQu();
  }


  private void computeTc() {
    // compute Tc using Eq. 2-5
    Tc = Math.pow(watershedLength, 0.8) * Math.pow((1000.0 / runoffCurveNumber) - 9.0, 0.7)
        / (1140.0 * Math.sqrt(watershedSlope));
    if (Tc < 0.1) {
      Tc = 0.1;
    } else if (Tc > 10) {
      throw new IllegalArgumentException(
          "Computed Tc > 10 hours; EFH2 model is not valid for this range.\n" + "Please estimate peak discharge using NEH-4 procedure, TR-20 program");
    }
  }


  private void computeQu() {
    // lookup rainfall coefficients for qu computation
    double Ia = EFH2IaTable.lookupIa(runoffCurveNumber);
    double SAPG = Ia / precip;    // Ia/P
    List<RfC> rfcs = EFH2RainfallCoefficients.getRainfallCoeff(stormType, SAPG);
    if (rfcs.size() < 1) {
      //Logger logger = Logger.getLogger(this.getClass());
      System.err.println("Failed to lookup rainfall coefficients for storm " + stormType + " and Ia/P = " + SAPG);
      qu = Double.NaN;
      return;
    }

    // compute qu, per algorithm in Basic source
    double logTc = Math.log10(Tc);
    RfC rfc = rfcs.get(0);
    double logLower = rfc.getRfC2() + (rfc.getRfC3() + rfc.getRfC4() * logTc) * logTc;
    double lowerQug = Math.pow(10.0, logLower);
    if (rfcs.size() == 1) {

      // Ia/P matched one of the standard values, use just one set of coeffs
      qu = lowerQug / 640.0;

    } else {
      // Ia/P fell between two standard values, read second RfC coeffs
      double lowerSAPG = rfc.getSAPG();
      rfc = rfcs.get(1);
      double upperSAPG = rfc.getSAPG();

      // compute second Qug result and interpolate
      double logUpper = rfc.getRfC2() + (rfc.getRfC3() + rfc.getRfC4() * logTc) * logTc;
      double upperQug = Math.pow(10.0, logUpper);
      double Qug = lowerQug + (upperQug - lowerQug) * (SAPG - lowerSAPG) / (upperSAPG - lowerSAPG);
      qu = Qug / 640.0;
    }
  }


  /**
   * get storm type parameter, see enumerated values
   */
  public String getStormType() {
    return stormType;
  }


  /**
   * set storm type parameter
   *
   * @param stormType
   */
  public void setStormType(String stormType) {
    this.stormType = stormType;
  }


  public List<String> getStormTypes() {
    return EFH2RainfallCoefficients.getStormTypes();
  }


  /**
   * get total precipitation parameter
   *
   * @return precip in inches
   */
  public double getPrecip() {
    return precip;
  }


  /**
   * set total precipitation parameter
   *
   * @param precip total storm precip in inches
   */
  public void setPrecip(double precip) {
    if (precip < 1.0 || precip > 15) {
      throw new IllegalArgumentException("Precip out of range: (1.0 .. 15.0) : " + precip);
    }
    this.precip = Math.abs(precip);
  }


  /**
   * get watershed length parameter, feet
   */
  public double getWatershedLength() {
    return watershedLength;
  }


  /**
   * set watershed length parameter
   *
   * @param watershedLength feet
   */
  public void setWatershedLength(double watershedLength) {
    if (this.watershedLength < 200 || this.watershedLength > 26000) {
      throw new IllegalArgumentException(
          "Watershed flow length " + watershedLength
          + " out of bounds; please use TR-55 program to compute Tc and qu");
    }
    this.watershedLength = Math.abs(watershedLength);
  }


  /**
   * get watershed slope parameter
   *
   * @return slope, ft/ft
   */
  public double getWatershedSlope() {
    return watershedSlope;
  }


  /**
   * set watershed slope parameter
   *
   * @param watershedSlope ft/ft
   */
  public void setWatershedSlope(double watershedSlope) {
    if (Math.abs(watershedSlope) > 64.0 || Math.abs(watershedSlope) < 0.1) {
      throw new IllegalArgumentException(
          "Watershed slope must be in percent: 0.1% .. 64%");
    }
    this.watershedSlope = Math.abs(watershedSlope);
  }


  /**
   * get runoff curve number parameter
   *
   * @return curve number
   */
  public int getRunoffCurveNumber() {
    return runoffCurveNumber;
  }


  /**
   * set runoff curve number parameter
   *
   * @param runoffCurveNumber curve number
   */
  public void setRunoffCurveNumber(int runoffCurveNumber) {
    if (runoffCurveNumber < 40 || runoffCurveNumber > 95) {
      throw new IllegalArgumentException("Runoff Curve Number: 40 .. 95");
    }
    this.runoffCurveNumber = runoffCurveNumber;
  }


  /**
   * get computed runoff, in inches
   */
  public double getRunoffQ() {
    return runoffQ;
  }


  /**
   * get computed Tc time of concentration, hours
   */
  public double getTimeOfConcentration() {
    return Tc;
  }


  /**
   * get computed peak discharge rate, cfs/ac-in??
   */
  public double getUnitPeakDischarge() {
    return qu;
  }
}