using System;
using System.Collections.Generic;
using System.Text;

namespace GeoFly
{
    public class GreenAmptInfil
    {

        public GreenAmptInfil(GridLayerPara g_GridLayerPara)
        {
            this.g_GridLayerPara = g_GridLayerPara;
        }

        //设置栅格计算参数
        public void SetGridPara(int currow,int curcol,double timelen,double Prerateinf,double Precumr,
                             double Prerintns,double Precuminf,double Preexcum,double pcp,double Cp)
        {
            m_row = currow;
            m_col = curcol;
            m_timeLen = timelen;
            m_Pcp = pcp;
            m_dPreRateinf = Prerateinf;
            m_dPrecumr = Precumr;
            m_dPrerintns = Prerintns;
            m_dPrecuminf = Precuminf;
            m_dPreexcum = Preexcum;

            m_dEhc = EffHydroConductivity(g_GridLayerPara, Cp);
        }

        
        /// <summary>
        /// Green-Ampt法计算地表超渗产流
        /// </summary>
        /// <param name="SP_Sw"></param>
        public void GreenAmptExcessRunoff(double SP_Sw)
        {
            m_dSoilw = 0;
            m_dSurfQ = 0;
            int SoilID = (int)g_GridLayerPara.g_SoilLayer.Values[m_row, m_col];
            double dthet=g_GridLayerPara.SoilTypes[SoilID].SoilWaterDeficitPercent(SP_Sw);            

            double psidt = dthet * g_GridLayerPara.SoilTypes[SoilID].SP_WFCS;

            double drateinf=0;
            double drintns=0;
            double dcumr=0;
            double dcuminf=0;
            double dexinc=0;
            double dexcum=0;

            double dt = m_timeLen / 60;

            dcumr = m_dPrecumr + m_Pcp;
            drintns = m_Pcp / dt;

            if (m_dPreRateinf <= 0) 
                m_dPreRateinf = m_dEhc;
            if (m_dPrecumr <= 0)
                m_dPrecumr = 0;
            if (m_dPrerintns <= 0)
                m_dPrerintns = drintns;
            if (m_dPrecuminf <= 0)
                m_dPrecuminf = 0;
            if (m_dPreexcum <= 0)
                m_dPreexcum = 0;

            if (drintns <= m_dPreRateinf)
            {
                m_dSoilw = m_dPrerintns * dt;
                dcuminf = m_dPrecuminf + m_dSoilw;
                m_dSurfQ = 0;

                if (m_dPreexcum > 0)
                {
                    dexcum = m_dPreexcum;
                    dexinc = 0;
                }
                else
                {
                    dexcum = 0;
                    dexinc = 0;
                }
            }
            else
            {
                double dtmp, df;
                dtmp = m_dEhc * dt;
                bool bExistLoop = false;
                do
                {
                    df = 0;
                    df = m_dPrecuminf + m_dEhc * dt + psidt * Math.Log((dtmp + psidt) / (m_dPrecuminf + psidt));
                    if (Math.Abs(df - dtmp) <= 0.001)
                    {
                        dcuminf = df;
                        dexcum = dcumr - dcuminf;
                        if (dexcum <= 0)
                            dexcum = 0;
                        dexinc = dexcum - m_dPreexcum;
                        if (dexinc < 0)
                            dexinc = 0;
                        m_dSurfQ = m_dSurfQ + dexinc;
                        bExistLoop = true;
                    }
                    else
                        dtmp = df;
                } while (!bExistLoop);

                if (m_dPrecuminf > 0)
                    m_dSoilw = dcuminf - m_dPrecuminf;

            }

            if (dcuminf <= 0.00001)
                drateinf = m_dEhc;
            else
                drateinf = m_dEhc * (psidt / (dcuminf + 0.00001) + 1);

            m_dPreRateinf = drateinf;
            m_dPrecumr = dcumr;
            m_dPrerintns = drintns;
            m_dPrecuminf = dcuminf;
            m_dPreexcum = dexcum;
        }


        /// <summary>
        /// 计算有效水力传导率系数
        /// </summary>
        /// <param name="g_GridLayerPara"></param>
        /// <param name="Cp">默认为-1</param>
        /// <returns></returns>

        private double EffHydroConductivity(GridLayerPara g_GridLayerPara, double Cp)
        {
            double dehc=0;
            int SoilID= (int)g_GridLayerPara.g_SoilLayer.Values[m_row, m_col];
            if (Cp == -1)
                dehc = g_GridLayerPara.SoilTypes[SoilID].SP_Sat_K / 2;
            else
                dehc = g_GridLayerPara.SoilTypes[SoilID].SP_Sat_K / (1 + Math.Exp(0.2 * Cp));
            if (dehc <= 0)
                dehc = 0.001;
            return dehc;
        }


        //计算输出结果
        public double m_dSurfQ;
        public double m_dSoilw;
        /// <summary>
        /// 前一时段下渗率
        /// </summary>
        public double m_dPreRateinf;
        /// <summary>
        /// 前一时段累积下渗量
        /// </summary>
        public double m_dPrecumr;

        public double m_dPrerintns;
        public double m_dPrecuminf;
        public double m_dPreexcum;


        //外部输入
        int m_row;
        int m_col;
        public double m_timeLen;
        public double m_Pcp;
        public bool bDataOK;

        //内部计算
        public double m_dEhc;
        public GridLayerPara g_GridLayerPara;
    }
}