using System;
using System.Collections.Generic;
using System.Text;
using Mapack;
namespace GeoFly
{
    public class KridgingInter:SpatialInterpolate
    {
        /// <summary>
        /// 求半变异函数的值
        /// </summary>
        /// <param name="dist">距离值</param>
        /// <param name="a">变程</param>
        /// <param name="c0"></param>
        /// <param name="c"></param>
        /// <returns></returns>        
        public double Variogram(double dist, double a, double c0, double c)
        {
            if (dist < 0)
                throw new Exception("dist不能为负值");
            double result = 0;
            if (dist <= a)
                result = c0 + c * (3 / 2 * dist / a) - Math.Pow(dist / a, 3) / 2;
            else if (dist > a)
                result = c0 + c;
            if (result < 0)
                throw new Exception("");
            return result;
        }
        /// <summary>
        /// 克立格插值法
        /// </summary>
        /// <param name="rowCount">待插值的栅格行数</param>
        /// <param name="colCount">待插值的栅格列数</param>
        /// <param name="resolution">每个栅格的分辨率</param>
        /// <param name="Xsite">样本点的X坐标(相对于栅格的左上角点的X坐标)</param>
        /// <param name="Ysite">样本点的Y坐标(相对于栅格的左上角点的Y坐标)</param>
        /// <param name="Zsite">样本点的属性值(用于插值)</param>
        /// <param name="a">半变异函数的变程</param>
        /// <param name="c0">块金值</param>
        /// <param name="c">基台值减去块金值的部分</param>
        /// <returns></returns>
        public GridLayer SpatialGridOut(DateTime date, double a, double c0, double c)
        {
            if (this.pStationinfo == null || this.pMeteoData == null)
                return null;

            //栅格数据结构初始化
            m_gridLayer = HydroSimulate.g_GridLayerPara.g_DemLayer.AttributesCopy();

            ProgressBar bar = new ProgressBar();
            bar.Show();
            bar.Text = "正在内插栅格数据";
            if (this.m_CoordType == CoordType.UTM_Coord)
            {
                this.PrepareFromUTM(date);
            }
            else
            {
                this.PrepareFromLongLa(date);
            }

            
            //下面构建方程组：K*lamda=D
            //K矩阵
            Matrix K = new Matrix(X.Length + 1, X.Length + 1);
            //D向量
            Matrix D = new Matrix(X.Length + 1, 1);

            //给K矩阵赋值
            for (int i = 0; i < X.Length + 1; i++)
            {
                for (int j = 0; j < X.Length + 1; j++)
                {
                    if (i < X.Length && j < X.Length)
                    {
                        //计算两个站点间的距离
                        double distX = X[i] - X[j];
                        double distY = Y[i] - Y[j];
                        double dist = Math.Sqrt(distX * distX + distY * distY);
                        K[i, j] = this.Variogram(dist, a, c0, c);
                    }
                    else
                    {
                        K[i, j] = 1.0;
                    }
                }
            }
            K[X.Length, X.Length] = 0;
            //K.save("c:\\LiuMatrix.txt");
            IMatrix inverseK = K.Inverse;
            //给D矩阵赋值并插值
            for (int row = 0; row < m_gridLayer.rowCount; row++)		//开始新一日的空间插值计算
            {
                for (int col = 0; col < m_gridLayer.colCount; col++)
                {
                    //获取当前计算栅格中心坐标
                    double CurrentX = col * 1.0;
                    double CurrentY = row * 1.0;
                    for (int i = 0; i < X.Length; i++)
                    {
                        double distX = X[i] - CurrentX;
                        double distY = Y[i] - CurrentY;
                        double dist = Math.Sqrt(distX * distX + distY * distY);
                        D[i, 0] = this.Variogram(dist, a, c0, c);
                    }
                    D[X.Length, 0] = 1;
                    IMatrix lambda = inverseK.Multiply(D);
                    double sum_Z = 0;
                    for (int i = 0; i < X.Length; i++)
                    {
                        sum_Z += lambda[i, 0] * Z[i];
                    }
                    //if (sum_Z < 0)
                    //    sum_Z = 0;
                    m_gridLayer[row, col] = sum_Z;
                }
            }
            return this.m_gridLayer;
        }
        
    }
}