/*
 * SoilUtil.java
 *
 * Created on February 28, 2007, 2:59 PM
 *
 * This class provides calculations for soil estimates.
 */

package usda.weru.soil;

/**
 *
 * @author joelevin
 */
public class SoilUtil {
    
    public static double estimateMineralBulkDensity(double clay, double sand){
        int li, lj, hi, hj;
        double mi, mj, fi, fj;
        double mbd, mbd_hi_hj;
        
        if (Double.isNaN(clay))
            clay = 0;
        
        if (Double.isNaN(sand))
            sand = 0;
        
        /* 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;
    }
    
    
    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, 1d/3d) - 1) * 100;        
    }
    
}