/*
 * $Id$
 *
 * This file is part of the Cloud Services Integration Platform (CSIP),
 * a Model-as-a-Service framework, API, and application suite.
 *
 * 2012-2017, OMSLab, Colorado State University.
 *
 * OMSLab licenses this file to you under the MIT license.
 * See the LICENSE file in the project root for more information.
 */
package soils.utils;

/**
 *
 * @author joelevin
 */
public class WEPSUtils {

  //  table for settled bulk density as a function of sand (across the top) and clay (down the side)
  static final double[][] mbdtv = {
    {1.48d, 1.25d, 1.00d, 1.06d, 1.16d, 1.22d, 1.30d, 1.39d, 1.45d, 1.51d, 1.52d},
    {1.52d, 1.40d, 1.19d, 1.25d, 1.32d, 1.40d, 1.52d, 1.58d, 1.63d, 1.65d},
    {1.52d, 1.40d, 1.25d, 1.35d, 1.45d, 1.53d, 1.60d, 1.64d, 1.72d},
    {1.52d, 1.40d, 1.29d, 1.41d, 1.50d, 1.57d, 1.63d, 1.68d},
    {1.50d, 1.40d, 1.35d, 1.43d, 1.53d, 1.61d, 1.64d},
    {1.46d, 1.40d, 1.40d, 1.43d, 1.53d, 1.62d},
    {1.45d, 1.40d, 1.38d, 1.42d, 1.50d},
    {1.42d, 1.37d, 1.33d, 1.33d},
    {1.33d, 1.32d, 1.20d},
    {1.23d, 1.18d},
    {1.15d}
  };

  static final double averageOrganicMatterBulkDensity = 0.224d;


  /**
   *
   * @param clay
   * @param sand
   * @return
   */
  public static double estimateMineralBulkDensity(double clay, double sand) {

    //The fractions can't sum to more than 1.0d.
    double sum = sand + clay;
    if (sum > 1.0d) {
      sand /= sum;
      clay /= sum;
      //System.out.println("essbd: sand plus clay fractions greater than 1.0."+
      //"  Values adjusted by averaging the difference.");
    }

    /* li           = index into table less than clay content
     * lj           = index into table less than sand content
     * hi           = index into table greater than clay content
     * hj           = index into table greater than sand content
     * mi           = value between indexes for the interpolation of clay
     * mj           = value between indexes for the interpolation of sand
     * fi           = fraction of distance between grid cells for clay
     * fj           = fraction of distance between grid cells for sand
     * mbd          = mineral mulk density
     * mbd_hi_hj    = value for mbdtv[hi][hj] if outside triangular part of table it is reflected 
      from mbdtv[li][lj] otherwise it is just the real point
     */
    double mi = clay * 10.0d;
    int li = (int) mi;
    double fi = mi - li;
    int hi = Math.min(li + 1, 10);

    double mj = sand * 10.0d;
    int lj = (int) mj;
    double fj = mj - lj;
    int hj = Math.min(lj + 1, 10);

    double mbd, mbd_hi_hj;
    if (li + lj == 10) {
      //On table edge, no interpolation necessary
      mbd = mbdtv[li][lj];
    } else {
      if (hi + hj > 10) {
        //Interpolation on the triangular edge of the table
        //mirror li, lj value to make using grid interpolation possible
        mbd_hi_hj = mbdtv[li][hj] + mbdtv[hi][lj] - mbdtv[li][lj];
      } else {
        //interpolation within the table, use grid point
        mbd_hi_hj = mbdtv[hi][hj];
      }
      mbd = (1 - fi) * (1 - fj) * mbdtv[li][lj]
          + (1 - fi) * fj * mbdtv[li][hj]
          + fi * (1 - fj) * mbdtv[hi][lj]
          + fi * fj * mbd_hi_hj;
    }
    return mbd;
  }


  /**
   * The following method estimates settled soil bulk density from intrinsic
   * properties. See Rawls (1983) Soil Science 135, 123-125
   *
   * If the sum of the sand and clay fractions exceed 1.0 then the fractions are
   * averaged.
   *
   * @param clay fraction of soil clay content
   * @param sand fraction of soil sand content
   * @param om organic matter
   * @return settled bulk density
   */
  public static double estimateSettledBulkDensity(double clay, double sand, double om) {
    //Average om bds gm/cm^3
    double mbd = estimateMineralBulkDensity(clay, sand);
    return 1.0d / ((om / averageOrganicMatterBulkDensity) + ((1.0d - om) / mbd));
  }


  /**
   *
   * @param settledBulkDensity
   * @param wetBulkDensity
   * @return
   */
  public static double estimateLinearExtensibility(double settledBulkDensity, double wetBulkDensity) {
    if (wetBulkDensity > settledBulkDensity) {
      throw new IllegalArgumentException("Wet bulk density greater than settled bulk density.");
    }
    return (Math.pow(settledBulkDensity / wetBulkDensity, 1 / 3) - 1) * 100;
  }

}