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

import csip.Config;
import csip.SessionLogger;
import d.soils.wwe02_wepssoilinput.StratifiedTextures;
import java.io.File;
import java.io.IOException;
import soils.Component;
import soils.Fragments;
import soils.Horizon;
import soils.MapUnit;
import soils.TextureGroup;
import soils.utils.EvalResult;

/**
 *
 * @author Brad
 * Only use this class if we are suppose to estimate the IFC values
 */
public class SoilCalc {

    private static String FORCE_CALC = "force_calculations";

    
    public static void convertUnits(MapUnit mapUnit,Component comp){
        
        mapUnit.brockdepmin(mapUnit.brockdepmin() *10);
        comp.calculated_resdeptmin_r(comp.calculated_resdeptmin_r() *10);
        comp.slope_r(comp.slope_r()/100);
        
        for(Horizon h : comp.horizons.values()){
            h.hzdepb_r(Utils.metricConversion(h.hzdepb_r(), Utils.Metric.CM, Utils.Metric.MM));
            h.hzdept_r(Utils.metricConversion(h.hzdept_r(), Utils.Metric.CM, Utils.Metric.MM));
            h.sandtotal_r(Utils.pctToFraction(h.sandtotal_r()));
            h.claytotal_r(Utils.pctToFraction(h.claytotal_r()));
            h.silttotal_r(Utils.pctToFraction(h.silttotal_r()));
            h.sandvc_r(Utils.pctToFraction(h.sandvc_r()));
            h.sandco_r(Utils.pctToFraction(h.sandco_r()));
            h.sandmed_r(Utils.pctToFraction(h.sandmed_r()));
            h.sandfine_r(Utils.pctToFraction(h.sandfine_r()));
            h.sandvf_r(Utils.pctToFraction(h.sandvf_r()));
            h.om_r(Utils.pctToFraction(h.om_r()));
            h.caco3_r(Utils.pctToFraction(h.caco3_r()));
        }      
    }
    
    public static double getSurfaceAlbedo(double surfaceAlbedo,double organicMaterial) {
        if (surfaceAlbedo == 0.0 || EvalResult.testDefaultDouble(surfaceAlbedo) || surfaceAlbedo < 0 || surfaceAlbedo > 1) {
            surfaceAlbedo = 0.6 / (Math.exp(0.4 * organicMaterial * 100.0));
        }
        
        return surfaceAlbedo;
    }

    /**
     * This variable in the WEPS code is only set by the GUI but I created this in case we hear otherwise
     * @return 0
     */
    public static double getSurfaceFragmentCover() {
        return 0;
    }
    
    public static double getSoilpH(Horizon horizon){
        double soilPh;
        
        soilPh = !EvalResult.testDefaultDouble(horizon.ph1to1h2o_r()) ? horizon.ph1to1h2o_r() : horizon.ph01mcacl2_r();
        if (soilPh == 0.0 || EvalResult.testDefaultDouble(soilPh))
            soilPh = 7;
        return soilPh;
    }

    static double getCationExchangeCapacity(Horizon horizon ) {
        double cec7;
        
        cec7 = !EvalResult.testDefaultDouble(horizon.cec7_r()) ? horizon.cec7_r() : horizon.ecec_r();
        if(cec7 == 0 || EvalResult.testDefaultDouble(cec7)){
            cec7 = horizon.claytotal_r() * 100 * 0.5 + horizon.om_r() * 100 * 2.0;
            
            cec7 = Utils.checkBoundaries(cec7, 0.0, 400);
        }
        horizon.cec7_r(cec7);
        return cec7;
    }

    static double getLinearExtensibility(Horizon horizon) {
        double lep;
        double bsd;
        
        lep = horizon.lep_r();
        
        if(lep == 0 || EvalResult.testDefaultDouble(lep)){
            bsd = SoilUtil.estimateSettledBulkDensity(horizon.claytotal_r(), horizon.sandtotal_r(), horizon.om_r());
            
            if(horizon.wthirdbar_r() > bsd)
                bsd = horizon.wthirdbar_r();
            lep = SoilUtil.estimateLinearExtensibility(bsd, horizon.wthirdbar_r());
            
            lep = Utils.checkBoundaries(lep,0,30);
        }
                
        horizon.lep_r(lep);
        return lep;
    }

    static double getAggregateMeanDiameter(Horizon horizon) {
        double mean;
        
        mean = Math.exp(1.343 - 2.235 * horizon.sandtotal_r() - // aggregate geometric mean diameter
                    1.226 * horizon.silttotal_r() - 0.0238 * horizon.sandtotal_r() / horizon.claytotal_r()
                    + 33.6 * horizon.om_r() + 6.85
                    * horizon.caco3_r()) * (1 + 0.006 * horizon.hzdepb_r());
        
        mean = Utils.checkBoundaries(mean, 0.03, 30.0);     
        
        return mean;
    }
    
    static double getAggregateStdDeviation(Horizon horizon, double mean){
        double deviation;
        
        if(EvalResult.testDefaultDouble(mean))
            mean = getAggregateMeanDiameter(horizon);
        
        deviation = 1 / (0.012448 + 0.002463 * mean
                + (0.093467 / (Math.pow(mean, 0.5))));
        
        deviation = Utils.checkBoundaries(deviation, 1.0, 20.0);        
        
        return deviation;
    }

    static double getAggregateStdDeviation(Horizon horizon) {  
        return getAggregateStdDeviation(horizon,Double.NaN);
    }

    static double getMaxAggregateSize(Horizon horizon) {
        double maxsize;
        double mean;
        double dev;
        
        mean = getAggregateMeanDiameter(horizon);
        dev = getAggregateStdDeviation(horizon,mean);
        
        maxsize = Math.pow(dev,(1.52 * (Math.pow(mean, -0.449)))) * mean;
        
        maxsize = Utils.checkBoundaries(maxsize, 1.0, 1000.0);
        
        return maxsize;
    }

    /**
     * The WEPS code set this value to 0.01.
     * @return 
     */
    static double getMinAggregateSize() {
        return 0.01;
    }

    /**
     * The WEPS code sets this to 1.8
     * It then checks to make sure it is between 0.6 and 2.5?
     * @return 
     */
    static double getAggregateDensity(Horizon horizon) {
        return 1.8;
    }

    static double getAggregateStability(Horizon horizon) {
        double stability;
        
        stability = (horizon.claytotal_r() > 0.5) ? 2.73 : 0.83 + 15.7 * horizon.claytotal_r()
                    - 23.8 * Math.pow(horizon.claytotal_r(), 2);	
        
        stability = Utils.checkBoundaries(stability, 0.1, 7.0);
        
        return stability;
    }

    /**
     * The WEPS code sets this to 0.01
     * @return 
     */
    static double getCrustThickness() {
        return 0.01;
    }

    static double getCrustDensity(Horizon horizon) {
        return getAggregateDensity(horizon);
    }

    static double getCrustStability(Horizon horizon) {
        return getAggregateStability(horizon);
    }

    /**
     * Only set by WEPS GUI
     * @return 
     */
    static double getCrustFraction() {
        return 0.0;
    }

    /**
     * Only set by WEPS GUI
     * @return 
     */
    static double getCrustLooseMaterialMass() {
        return 0.0;
    }

    /**
     * Only set by WEPS GUI
     * @return 
     */
    static double getCrustLooseMaterialFraction() {
        return 0.0;
    }

    /**
     * WEPS code sets this to 4
     * @return 
     */
    static double getRandomRoughness() {
        return 4;
    }

    /**
     * Only set by WEPS GUI
     * @return 
     */
    static double getRoughnessOrientation() {
        return 0.0;
    }

    /**
     * Only set by WEPS GUI
     * @return 
     */
    static double getRoughnessHeight() {
        return 0.0;
    }

    /**
     * Only set by WEPS GUI
     * @return 
     */
    static double getRoughnessSpacing() {
        return 10.0;
    }

    /**
     * Only set by WEPS GUI
     * @return 
     */
    static double getRoughnessWidth(){
        return 10.0;
    }

    static double getInitialSWC(Horizon horizon) {
        return (getFieldCapacitySWC(horizon) + getWiltingPointSWC(horizon)) / 2;
    }

    static double getSaturatedSWC(Horizon horizon) {
        double saturated;
        double sax_x = 0.332, sax_y = -7.251e-4, sax_z = 0.1276;			// coef used to calc THETAS
        
        saturated = sax_x + sax_y * horizon.sandtotal_r() * 100. + // saturated soil water content
                    sax_z * (Math.log(horizon.claytotal_r() * 100.) / Math.log(10.));
        
        saturated = Utils.checkBoundaries(saturated, 0.208,0.550);
        
        return saturated;
    }
    
    static double getSoilCB(Horizon horizon) {
        double cb;
        double sax_e = -3.140, sax_f = -2.22e-3, sax_g = -3.484e-5;			// coef used to calc CB
        
        cb = -(sax_e + // cb(k) == ifc3_layer[][23]
            (sax_f * Math.pow(horizon.claytotal_r() * 100., 2.0))
            + (sax_g * Math.pow(horizon.sandtotal_r() * 100., 2.0) * horizon.claytotal_r() * 100.));
        
        cb = Utils.checkBoundaries(cb, 0.917, 27.027);
        
        return cb;
    }

    static double getFieldCapacitySWC(Horizon horizon) {
        double capacity;
        boolean force = Config.getBoolean(FORCE_CALC, false);

        
        capacity = Utils.pctToFraction(horizon.wthirdbar_r());
        
        if (force || capacity == 0.0 || EvalResult.testDefaultDouble(capacity)) {
            capacity = Math.pow((33.0 / get_sax_a(horizon)), (1 / -getSoilCB(horizon)));
            
            capacity = Utils.checkBoundaries(capacity, 0.012, .0335);
        }
        
        return capacity;
    }

    static double getWiltingPointSWC(Horizon horizon) {
        double wilting;
        boolean force = Config.getBoolean(FORCE_CALC, false);
        
        wilting = Utils.pctToFraction(horizon.wfifteenbar_r());
        
        if (force || wilting == 0.0 || EvalResult.testDefaultDouble(wilting)) {
            wilting = Math.pow((1500.0 / get_sax_a(horizon)), (1 / -getSoilCB(horizon)));
            
            wilting = Utils.checkBoundaries(wilting, 0.005, 0.242);
        }
        
        return wilting;
    }

    static double getAirEntryPotential(Horizon horizon) {
        double aep;
        
        aep = -get_sax_a(horizon) * Math.pow(getSaturatedSWC(horizon), -getSoilCB(horizon)); 
        
        aep = Utils.checkBoundaries(aep, -17.91, 0.0);
        
        return aep;
    }

    static double getSaturatedHydraulicConductivity(Horizon horizon) {
        double shc;
        boolean force = Config.getBoolean(FORCE_CALC, false);
        double sk;
        
        shc = horizon.ksat_r() / 1000000; // cvt from ppm to p
        
        if (force || shc == 0.0 || EvalResult.testDefaultDouble(shc)) {
            sk = 9.6 - 0.81 * Math.log(horizon.silttotal_r() * 100.) / Math.log(10.)
                        - 1.09 * Math.log(horizon.claytotal_r() * 100.) / Math.log(10.) - 4.64 * horizon.dbthirdbar_r();
            
            shc = Math.pow(10, sk) * 0.24 * (1.0 / 86400);
            
            shc = Utils.checkBoundaries(shc, 0.0, 0.001);
        }
        
        return shc;
    }
    
    private static double get_sax_a(Horizon horizon){
        double sax_a = 100 * Math.exp(-4.396 - // calc A factor for pote
            0.0715 * (horizon.claytotal_r() * 100.) - // using two steps because that is how
            (4.880e-4) * Math.pow(horizon.sandtotal_r() * 100., 2) - // saxton95 does it
            (4.285e-5) * Math.pow(horizon.sandtotal_r() * 100., 2) * (horizon.claytotal_r() * 100.));
        
        return sax_a;
    }

    public static void averageStratifiedLayers(Component comp,StratifiedTextures st) throws IOException {
        boolean first = true;
        double sand = 0;
        double clay = 0;
        double vsf = 0;
        int count = 0;
        String[] textures;
        double silt;

        Layers:
        for (Horizon h : comp.horizons.values()) {
            for(TextureGroup t : h.textureGroups.values()){
                if(t.stratextsflag() && t.rvindicator() && t.texture().trim().toUpperCase().startsWith("SR")){ // stratified key
                    textures = t.texture().split(" ", -1);
                    
                first = true;
                sand = 0;
                clay = 0;
                vsf = 0;
                count = 0;
                
                    for (String subTexture : textures) {
                        subTexture = subTexture.trim();
                        if (first) {
                            //skip the sr marker
                            first = false;
                            continue;
                        }
                        if (st.isTextureValid(subTexture)) {
                            count++;
                            sand += st.getSand(subTexture);
                            sand /= count;

                            clay += st.getClay(subTexture);
                            clay /= count;

                            vsf += st.getVeryFineSand(subTexture);
                            vsf /= count;
                        } else {
                            //can't fix up this layer
                            continue Layers;
                        }
                    }
                    if (sand + clay > 1) {
                        //TODO: log error
                        continue Layers;
                    }
                    silt = 1 - sand - clay;
                    //if(EvalResult.testDefaultDouble(h.sandtotal_r()))
                        h.sandtotal_r(sand * 100);
                    //if(EvalResult.testDefaultDouble(h.silttotal_r()))
                        h.silttotal_r(silt * 100);
                    //if(EvalResult.testDefaultDouble(h.claytotal_r()))
                        h.claytotal_r(clay * 100);
                    //if(EvalResult.testDefaultDouble(h.sandvf_r()))
                        h.sandvf_r(vsf * 100);
                }
            }       
        }
    }
    
    public static String getSurfaceTexture(Component comp){
        
        return getTexture(Utils.getSurfaceHorizon(comp));
    }
    
    public static String getTexture(Horizon horizon){
        for(TextureGroup t : horizon.textureGroups.values())
            if(t.rvindicator())
                return t.texture();
        return "";
    }
    
    public static double getRockFrags(Horizon horizon){
        double rockFrag = 0;
        
        for(Fragments f : horizon.fragments.values()){
            rockFrag += Utils.averageMissingValue(f.fragvol_h(), f.fragvol_l(), f.fragvol_r());
        }
        
        return rockFrag / 100;
    }
    
    public static double testNumericSoilElement(double high,double low,double r, double defaultValue) {

        if(EvalResult.testDefaultDouble(r))
            r = defaultValue;
        
        if(r < low)
            r = defaultValue;
        if(r > high)
            r= defaultValue;
        
        return r;
    }
    
    public static int testNumericSoilElement(int high,int low,int r, int defaultValue) {

        if(EvalResult.testDefaultInteger(r))
            r = defaultValue;
        
        if(r < low)
            r = defaultValue;
        if(r > high)
            r = defaultValue;
        
        return r;
    }
}