/*
 * $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 d.util;

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

    /**
     *
     * @param clay
     * @param sand
     * @return
     */
    public static double estimateMineralBulkDensity(double clay, double sand) {
        int li, lj, hi, hj;
        double mi, mj, fi, fj;
        double mbd, mbd_hi_hj;

        /* 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
         */
        //Data table of settled bulk density as a function of sand (across the top) and clay (down the side)
        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}
        };

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

        mi = clay * 10.0d;
        li = (int) mi;
        fi = mi - li;
        hi = Math.min(li + 1, 10);

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

        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) {
        double bds, mbd, mbd_hi_hj;

        //Average om bds gm/cm^3
        double averageOrganicMatterBulkDensity = 0.224d;

        mbd = estimateMineralBulkDensity(clay, sand);
        bds = 1.0d / ((om / averageOrganicMatterBulkDensity) + ((1.0d - om) / mbd));

        return bds;
    }

    /**
     *
     * @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;
    }

}