// <summary>ソースコード：テスト用流量発生モデルクラス</summary>
// <author>CommonMP</author>

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

using CommonMP.HYSSOP.CoreImpl;
using CommonMP.HYSSOP.CoreImpl.HSData;
using CommonMP.HYSSOP.CoreImpl.HSTools;

using CommonMP.HYMCO.Interface;
using CommonMP.HYMCO.Interface.Model;
using CommonMP.HYMCO.Interface.Data;
using CommonMP.HYMCO.CoreImpl.Data;
using CommonMP.HYMCO.CoreImpl.Tool;


namespace CommonMP.HYMCO.OptionImple.McHyHdFormulaSampleH13
{
    /// <summary><para>class outline:</para>
    /// <para>水理公式集例題１－９より　Kinematic wave 法による流出計算 （直接差分法）クラス</para>
    /// </summary>
    /// <remarks><para>history:</para>
    /// <para>[CommonMP][ver 1.2.0][2011/11/11][新規作成]</para>
    /// </remarks>
    public class McKinematicWaveBasin : McForecastModelBase
    {
        /// <summary>流域計算情報</summary>
        McKinematicWaveBasinCalInfo m_KMdlInf = null;

        /// <summary><para>method outline:</para>
        /// <para>入力側の接続情報チェック</para>
        /// </summary>
        /// <example><para>usage:</para>
        /// <para>bool bRtn = ReceiveConnectionCheck(ref csErrorInf)</para>
        /// </example>
        /// <param name="csErrorInf">エラー出力</param>
        /// <param name="lInputDataNum">入力情報数</param>
        /// <param name="csInputCellData">演算に必要な入力情報配列</param>
        /// <returns>=true:正常、=false:異常</returns>
        /// <exception cref="">無し</exception>
        /// <remarks><para>remarks:</para>
        /// <para>受信するデータが自モデルが期待している情報か否かをチェックする</para>
        /// </remarks>
        protected override bool ReceiveConnectionCheck(ref McStructErrorInfo csErrorInf, long lInputDataNum, McReceiveCellDataIF[] csInputCellData)
        {
            bool bRtn = true;
            if (lInputDataNum == 1)
            {
                //McModelCellCorrespondenceInf csCellMng = csInputCellData[0].GetCellIDData() as McModelCellCorrespondenceInf;
                // == 入力パターンに応じて何らかの処理を行いたい場合の例 ===
                HySID csInputPtnID = csInputCellData[0].GetReceivePatternID();
                if (csInputPtnID.Equals(McHyHdFormulaSampleH13Define.IN_PATTERN_RAIN) == true)
                {
                    // 正常
                }
                else
                {
                    // エラー
                    csErrorInf.AddCheckErrorData(this.GetID(), McHyHdFormulaSampleH13Define.KINEMATIC_WAVE_BASIN_KIND,
                                            "Unexpected input data pattern. (From " + csInputCellData[0].GetUpperElementID().ToString() + ")");
                    bRtn = false;
                }
                /*
                // == 伝送データのセル配列及びセル内の変数の数を取得する例 =====
                if (McTranInfoDefine.SINGLE_CELL_SERIAL.Equals(csInputCellData[0].GetTranInfoKind()) == true)
                {
                    // 正常
                }
                else
                {
                    // エラー
                    csErrorInf.AddCheckErrorData(this.GetID(), McHyHdFormulaSampleH13Define.KINEMATIC_WAVE_BASIN_KIND,
                                            "Unexpected input data type. (From " + csInputCellData[0].GetUpperElementID().ToString() + ")");
                    bRtn = false;
                }
                */
            }
            else
            {
                // エラー
                csErrorInf.AddCheckErrorData(this.GetID(), McHyHdFormulaSampleH13Define.KINEMATIC_WAVE_BASIN_KIND, "Receive data number must be one.");
                bRtn = false;
            }
            return bRtn;
        }
        /// <summary><para>method outline:</para>
        /// <para>出力側の接続情報チェック</para>
        /// </summary>
        /// <example><para>usage:</para>
        /// <para>bool bRtn = SendConnectionCheck(ref csErrorInf)</para>
        /// </example>
        /// <param name="csErrorInf">エラー出力</param>
        /// <param name="lOutputDataNum">出力情報数</param>
        /// <param name="csOutputCellData">出力情報配列</param>
        /// <returns>=true:正常、=false:異常</returns>
        /// <exception cref="">無し</exception>
        /// <remarks><para>remarks:</para>
        /// <para>送信端子に設定されている伝送データが自モデルが期待している情報か否かをチェックする</para>
        /// </remarks>
        protected override bool SendConnectionCheck(ref McStructErrorInfo csErrorInf, long lOutputDataNum, McSendCellDataIF[] csOutputCellData)
        {
            bool bRtn = true;

            if (lOutputDataNum == 0)
            {
                // 警告（エラーではない）
                csErrorInf.AddCheckWarningData(this.GetID(), McHyHdFormulaSampleH13Define.KINEMATIC_WAVE_BASIN_KIND, "No send port.");
            }
            else
            {
                m_KMdlInf.m_dNextTQTimeOut = new double[lOutputDataNum];
                m_KMdlInf.m_dNextTHTimeOut = new double[lOutputDataNum];
            }
            for (long lP = 0; lP < lOutputDataNum; lP++)
            {   // 出力側接続線数分チェックを行う
                if (McHyHdFormulaSampleH13Define.OUT_PATTERN_Q.Equals(csOutputCellData[lP].GetSendPatternID()) == true)
                {   // 流量の出力ならば
                    long lD1 = 0; long lD2 = 0; long lD3 = 0; long lCellDataNum = 0;
                    long lCellDim = csOutputCellData[lP].GetDimension(ref lD1, ref lD2, ref lD3, ref lCellDataNum);
                    if (lCellDim == 1 && lD1 == 1)
                    {   // １次元配列で１セルならば
                        if (lCellDataNum < 1)
                        {   // セル内の変数の数
                            csErrorInf.AddCheckErrorData(this.GetID(), McHyHdFormulaSampleH13Define.KINEMATIC_WAVE_BASIN_KIND,
                                "Variable number in the Cells is too short. (Send To " + csOutputCellData[lP].GetLowerElementID().ToString() + ")");
                            bRtn = false;
                        }
                    }
                    else
                    {
                        // エラー
                        csErrorInf.AddCheckErrorData(this.GetID(), McHyHdFormulaSampleH13Define.KINEMATIC_WAVE_BASIN_KIND,
                                                "Unexpected send data type. (Send To " + csOutputCellData[lP].GetLowerElementID().ToString() + ")");
                        bRtn = false;
                    }
                }
                else if (McHyHdFormulaSampleH13Define.OUT_PATTERN_DEPTH.Equals(csOutputCellData[lP].GetSendPatternID()) == true)
                {
                    long lD1 = 0; long lD2 = 0; long lD3 = 0; long lCellDataNum = 0;
                    long lCellDim = csOutputCellData[lP].GetDimension(ref lD1, ref lD2, ref lD3, ref lCellDataNum);
                    if (m_KMdlInf.nx + 1 > lD1)
                    {
                        csErrorInf.AddCheckErrorData(this.GetID(), McHyHdFormulaSampleH13Define.DIFFUSION_WAVE_LANE_KIND,
                                                "Cell dimention(" + csOutputCellData[lP].GetTranInfoKind().ToString() + ") is too short. (send to " + csOutputCellData[lP].GetLowerElementID().ToString() + ")");
                        bRtn = false;
                    }
                }
                else
                {
                    // エラー
                    csErrorInf.AddCheckErrorData(this.GetID(), McHyHdFormulaSampleH13Define.KINEMATIC_WAVE_BASIN_KIND,
                                            "Unexpected send data type. (Send To " + csOutputCellData[lP].GetLowerElementID().ToString() + ")");
                    bRtn = false;
                }
            }
            return bRtn;
        }

        //=========================
        // 演算実行前処理関連メソッド
        //=========================

        /// <summary><para>method outline:</para>
        /// <para>モデルを初期化する</para>
        /// </summary>
        /// <example><para>usage:</para>
        /// <para>Initialize(csInitialData)</para>
        /// </example>
        /// <param name="csInitialData">初期化設定情報</param>
        /// <param name="lInputDataNum">入力情報数</param>
        /// <param name="csInputCellData">演算に必要な入力情報配列</param>
        /// <returns>無し</returns>
        /// <exception cref="">無し</exception>
        /// <remarks><para>remarks:</para>
        /// <para>無し</para>
        /// </remarks>
        protected override bool Initialize(ref McPropertyInfoRoot csInitialData, long lInputDataNum, ref McReceiveCellDataIF[] csInputCellData)
        {
            //if (csInitialData != null)
            //{
            //}
            for (long lLp = 0; lLp < m_KMdlInf.m_dNextTQTimeOut.Length; lLp++)
            {
                m_KMdlInf.m_dNextTQTimeOut[lLp] = 0.0;
                m_KMdlInf.m_dNextTHTimeOut[lLp] = 0.0;
            }
            long lQIdx = 0;  // 降雨量は０番目
            HySCellData csCell = null;  // 取得するセル
            csCell = csInputCellData[0].GetInterpolatedCell(0); // セル配列の０番目のセルを取得
            m_KMdlInf.rr = csInputCellData[0].Data(csCell, lQIdx);   // セル内の０番目の変数値を取得



            m_KMdlInf.Qi = 0.0;  /* すべての斜面上の地点で単位幅流量を 0 とする */

            for (int i = 0; i <= m_KMdlInf.nx; i++)
            {
                m_KMdlInf.q1[i] = m_KMdlInf.g[2 * i] / m_KMdlInf.L * m_KMdlInf.Qi;
                if (m_KMdlInf.q1[i] == 0.0)
                {
                    m_KMdlInf.depth[i] = 0.0;
                }
                else
                {
                    m_KMdlInf.depth[i] = hq_equation(m_KMdlInf.type_qh_equation, m_KMdlInf.Qi / m_KMdlInf.L, m_KMdlInf.alpha_a[2 * i], m_KMdlInf.alpha_s[2 * i], m_KMdlInf.m, m_KMdlInf.d);
                }
                m_KMdlInf.h1[i] = m_KMdlInf.g[2 * i] * m_KMdlInf.depth[i];
            }

            m_KMdlInf.qout = m_KMdlInf.q1[m_KMdlInf.nx];
            
            return true;
        }
        

        //=======================
        // 演算実行処理関連メソッド
        //=======================

        /// <summary><para>method outline:</para>
        /// <para>モデル演算</para>
        /// </summary>
        /// <example><para>usage:</para>
        /// <para>long lRtn = Calculate(ref csInputCellData)</para>
        /// </example>
        /// <param name="lInputDataNum">入力情報数</param>
        /// <param name="csInputCellData">演算に必要な入力情報配列</param>
        /// <returns>=0:正常 -1:異常</returns>
        /// <exception cref="">無し</exception>
        /// <remarks><para>remarks:</para>
        /// <para>無し</para>
        /// </remarks>
        protected override long Calculate(long lInputDataNum, ref McReceiveCellDataIF[] csInputCellData)
        {
            long lQIdx = 0;  // 降雨量は０番目
            HySCellData csCell = null;  // 取得するセル
            csCell = csInputCellData[0].GetInterpolatedCell(0); // セル配列の０番目のセルを取得
            m_KMdlInf.rr = csInputCellData[0].Data(csCell, lQIdx);   // セル内の０番目の変数値を取得

            /* 降雨強度の設定 */
            for (int j = 0; j <= 2 * m_KMdlInf.nx; j++)
            {
                m_KMdlInf.rr1[j] = m_KMdlInf.rr * m_KMdlInf.xcos[j]; /* 差分開始時刻の降雨強度 */
                m_KMdlInf.rr2[j] = m_KMdlInf.rr * m_KMdlInf.xcos[j]; /*  1 差分時間間隔後の降雨強度 */
            }

            /* 境界条件の設定 */
            m_KMdlInf.q2[0] = m_KMdlInf.g[0] / m_KMdlInf.L * m_KMdlInf.Qb;
            if (m_KMdlInf.q2[0] == 0.0)
            {
                m_KMdlInf.Hb = 0.0;
            }
            else
            {
                m_KMdlInf.Hb = hq_equation(m_KMdlInf.type_qh_equation, m_KMdlInf.Qb / m_KMdlInf.L, m_KMdlInf.alpha_a[0], m_KMdlInf.alpha_s[0], m_KMdlInf.m, m_KMdlInf.d);
            }
            m_KMdlInf.h2[0] = m_KMdlInf.g[0] * m_KMdlInf.Hb;


            /* 差分計算 */
            if (m_KMdlInf.type_scheme == 1)
            { /* Lax Wendroff scheme */
                LaxWendroff(m_KMdlInf.nx, m_KMdlInf.delta_t, m_KMdlInf.h1, m_KMdlInf.q1, m_KMdlInf.h2, m_KMdlInf.q2, m_KMdlInf.rr1, m_KMdlInf.rr2,
                    m_KMdlInf.type_qh_equation, m_KMdlInf.g, m_KMdlInf.alpha_a, m_KMdlInf.alpha_s, m_KMdlInf.m, m_KMdlInf.d);
            }
            else if (m_KMdlInf.type_scheme == 2)
            { /* Nonlinear scheme */
                Nfdm(m_KMdlInf.nx, m_KMdlInf.delta_t, m_KMdlInf.h1, m_KMdlInf.q1, m_KMdlInf.h2, m_KMdlInf.q2, m_KMdlInf.rr1, m_KMdlInf.rr2,
                 m_KMdlInf.type_qh_equation, m_KMdlInf.g, m_KMdlInf.alpha_a, m_KMdlInf.alpha_s, m_KMdlInf.m, m_KMdlInf.d);
            }

            /* 流出量の計算結果のメモリへの格納 */
            m_KMdlInf.qout = m_KMdlInf.q2[m_KMdlInf.nx];

            /* 水深の計算結果のメモリへの格納 */
            m_KMdlInf.depth[0] = m_KMdlInf.Hb; /* 斜面上端での水深 */
            for (int j = 1; j <= m_KMdlInf.nx; j++)
            {
                m_KMdlInf.depth[j] = m_KMdlInf.h2[j] / m_KMdlInf.g[2 * j];
            }

            return 0;
        }


        /// <summary><para>method outline:</para>
        /// <para>モデル演算結果を外部のエレメントに対して公開する</para>
        /// </summary>
        /// <example><para>usage:</para>
        /// <para>long lRtn = DataFusion( lOutputDataNum, ref csOutputCellData)</para>
        /// </example>
        /// <param name="lOutputDataNum">出力情報数</param>
        /// <param name="csOutputCellData">出力情報配列</param>
        /// <returns>=0:正常 -1:異常</returns>
        /// <exception cref="">無し</exception>
        /// <remarks><para>remarks:</para>
        /// <para>複数接続があった場合、全て同じ情報が出ることに注意（流量を分配して出しているわけではない）</para>
        /// </remarks>
        protected override long DataFusion(long lOutputDataNum, ref McSendCellDataIF[] csOutputCellData)
        {
            /*  Do Nothing
            double dCalTtlTime = this.m_csTotalPassingTime.GetTime();
            HySCellData[] csSndCellData = null;
            for (long lP = 0; lP < m_lOutputDataNum; lP++)
            {
                if (McHyHdFormulaSampleH13Define.OUT_PATTERN_Q.Equals(csOutputCellData[lP].GetSendPatternID()) == true)
                {   // 流量出力ならば
                        csSndCellData = csOutputCellData[lP].PrepareSendCellD1();
                        csSndCellData[0].m_dData[0] = m_KMdlInf.qout;
                        m_KMdlInf.m_dNextTQTimeOut += m_KMdlInf.delta_tq;
                }
                else if (McHyHdFormulaSampleH13Define.OUT_PATTERN_DEPTH.Equals(csOutputCellData[lP].GetSendPatternID()) == true)
                {   // 水深出力ならば
                        for (long lLp = 0; lLp <= m_KMdlInf.nx; lLp++)
                        {
                            csSndCellData = csOutputCellData[lP].PrepareSendCellD1();
                            csSndCellData[lLp].m_dData[0] = m_KMdlInf.depth[lLp];
                        }
                        m_KMdlInf.m_dNextTHTimeOut += m_KMdlInf.delta_th;
                }
            }
            */
            return 0;
        }

        /// <summary><para>method outline:</para>
        /// <para>モデル演算結果を外部のエレメントに対して公開する</para>
        /// </summary>
        /// <example><para>usage:</para>
        /// <para>long lRtn = DataFusion(ref csOutputDataList)</para>
        /// </example>
        /// <param name="csOutputDataList">演算結果を公開する出力情報リスト</param>
        /// <returns>=0:正常 -1:異常</returns>
        /// <exception cref="">無し</exception>
        /// <remarks><para>remarks:</para>
        /// <para>無し</para>
        /// </remarks>
        public override long DataFusion(ref HySDataLinkedList csOutputDataList)
        {
            long lRtn = 0;

            double dCalTtlTime = this.m_csTotalPassingTime.GetTime();
            HySCellData[] csSndCellData = null;
            for (long lP = 0; lP < m_lOutputDataNum; lP++)
            {
                if (McHyHdFormulaSampleH13Define.OUT_PATTERN_Q.Equals(m_OutputCellData[lP].GetSendPatternID()) == true)
                {   // 流量出力ならば
                    if (dCalTtlTime >= m_KMdlInf.m_dNextTQTimeOut[lP])
                    {
                        // 現在時刻を指定して伝送するレコード情報を作成する
                        m_OutputCellData[lP].SetCurrentTime(this.m_csSimTime);

                        // レコードに計算結果を設定する
                        csSndCellData = m_OutputCellData[lP].PrepareSendCellD1();
                        csSndCellData[0].m_dData[0] = m_KMdlInf.qout;
                        m_KMdlInf.m_dNextTQTimeOut[lP] += m_KMdlInf.delta_tq;

                        // モデルがセットした伝送情報を　実際に伝送データとして流す
                        m_OutputCellData[lP].SendData();

                    }
                }
                else if (McHyHdFormulaSampleH13Define.OUT_PATTERN_DEPTH.Equals(m_OutputCellData[lP].GetSendPatternID()) == true)
                {   // 水深出力ならば
                    if (dCalTtlTime >= m_KMdlInf.m_dNextTHTimeOut[lP])
                    {
                        // 現在時刻を指定して伝送するレコード情報を作成する
                        m_OutputCellData[lP].SetCurrentTime(this.m_csSimTime);

                        // レコードに計算結果を設定する
                        for (long lLp = 0; lLp <= m_KMdlInf.nx; lLp++)
                        {
                            csSndCellData = m_OutputCellData[lP].PrepareSendCellD1();
                            csSndCellData[lLp].m_dData[0] = m_KMdlInf.depth[lLp];
                        }

                        // モデルがセットした伝送情報を　実際に伝送データとして流す
                        m_OutputCellData[lP].SendData();

                        m_KMdlInf.m_dNextTHTimeOut[lP] += m_KMdlInf.delta_th;
                    }
                }
            }

            return lRtn;
        }
        //====================
        // その他必要なメソッド
        //====================

        /// <summary><para>method outline:</para>
        /// <para>プロパティ情報を設定する</para>
        /// </summary>
        /// <example><para>usage:</para>
        /// <para>SetProperty(csCellMdlPropertyInfo)</para>
        /// </example>
        /// <param name="csCellMdlPropertyInfo">セル型プロパティ情報</param>
        /// <returns>無し</returns>
        /// <exception cref="">無し</exception>
        /// <remarks><para>remarks:</para>
        /// <para>無し</para>
        /// </remarks>
        public override bool SetProperty(McCellModelPropertyIF csCellMdlPropertyInfo)
        {
            // 使用しやすいようにキャストしておく
            m_KMdlInf = (McKinematicWaveBasinCalInfo)m_csCalInfo;


            McCellModelPropertyInfo csPrptyInfo = csCellMdlPropertyInfo as McCellModelPropertyInfo;
            if (csPrptyInfo != null)
            {
                // 演算ステップ時刻設定
                this.m_csDltTime = new HySTime(csPrptyInfo.GetStepTime());

                // 例題：S1-9:"calcon.dat"
                csPrptyInfo.GetInfo("type_scheme", ref m_KMdlInf.type_scheme);
                csPrptyInfo.GetInfo("_nx", ref m_KMdlInf.nx);
                //csPrptyInfo.GetInfo("delta_t", ref m_KMdlInf.delta_t);
                m_KMdlInf.delta_t = csPrptyInfo.GetStepTime() / 3600.0;  // Hr換算

                csPrptyInfo.GetInfo("delta_tq", ref m_KMdlInf.delta_tq);
                csPrptyInfo.GetInfo("delta_th", ref m_KMdlInf.delta_th);

                /*==================*/
                /*  斜面形状の設定  */
                /*==================*/
                /*---------------------------------*/
                /*  斜面形状のパラメータの読み込み */
                /*---------------------------------*/
                // 例題：S1-9:"slope.dat"
                csPrptyInfo.GetInfo("_L", ref m_KMdlInf.L);
                csPrptyInfo.GetInfo("_theta", ref m_KMdlInf.theta);
                m_KMdlInf.theta = m_KMdlInf.pi * m_KMdlInf.theta / 180.0;           /* 度を rad に変換する。 */
                csPrptyInfo.GetInfo("_bu", ref m_KMdlInf.bu);
                csPrptyInfo.GetInfo("_bl", ref m_KMdlInf.bl);
                m_KMdlInf.sall = m_KMdlInf.L * (m_KMdlInf.bu + m_KMdlInf.bl) / 2.0 * Math.Cos(m_KMdlInf.theta);  /* 斜面の面積 (m*m) */

                /*----------------------------------*/
                /* 地形パターン関数と斜面勾配の設定 */
                /*----------------------------------*/
                m_KMdlInf.g = new double[2 * m_KMdlInf.nx + 1];
                m_KMdlInf.xcos = new double[2 * m_KMdlInf.nx + 1];
                for (int i = 0; i <= 2 * m_KMdlInf.nx; i++)
                {
                    m_KMdlInf.g[i] = (m_KMdlInf.bu + (m_KMdlInf.bl - m_KMdlInf.bu) / (2 * m_KMdlInf.nx) * i) / m_KMdlInf.sall * m_KMdlInf.L; /* 地形パターン関数 g の設定 */
                    m_KMdlInf.xcos[i] = Math.Cos(m_KMdlInf.theta);                   /* 斜面勾配の設定 */
                }


                /*====================================*/
                /*  流量流積関係式の設定と単位の変換  */
                /*====================================*/
                /*------------------------*/
                /*  パラメータの読み込み  */
                /*------------------------*/
                // 例題：S1-9:"qh_eqn.dat"
                csPrptyInfo.GetInfo("type_qh_equation", ref m_KMdlInf.type_qh_equation);
                csPrptyInfo.GetInfo("_k", ref m_KMdlInf.k);
                m_KMdlInf.k = m_KMdlInf.k * 36000.0; /* 36000 倍しているのは，k の単位 cm/sec */
                                                     /* を mm/hr に変換するため */
                csPrptyInfo.GetInfo("_n", ref m_KMdlInf.n);
                m_KMdlInf.n = m_KMdlInf.n / 36000.0; /* 36000 で割るのは，n の単位 m^{-1/3} sec */
                                                     /* を mm^{-1/3} hr に変換するため */

                csPrptyInfo.GetInfo("_m", ref m_KMdlInf.m);
                csPrptyInfo.GetInfo("_D", ref m_KMdlInf.D);
                csPrptyInfo.GetInfo("_gamma", ref m_KMdlInf.gamma);
                m_KMdlInf.d = m_KMdlInf.gamma * m_KMdlInf.D;     /* 実質の A 層厚の設定 */
                /*--------------------*/
                /*  パラメータの設定  */
                /*--------------------*/
                m_KMdlInf.alpha_a = new double[2 * m_KMdlInf.nx + 1];
                m_KMdlInf.alpha_s = new double[2 * m_KMdlInf.nx + 1];
                if (m_KMdlInf.type_qh_equation == 1)
                { /* ダルシー式の場合 */
                    for (int i = 0; i <= 2 * m_KMdlInf.nx; i++)
                    {
                        m_KMdlInf.alpha_a[i] = m_KMdlInf.k * Math.Sin(m_KMdlInf.theta) / m_KMdlInf.gamma / (m_KMdlInf.L * 1000.0);
                        /* L を 1000 倍しているのは単位を m から mm に変換するため */
                    }
                }
                else if (m_KMdlInf.type_qh_equation == 2)
                { /* マニング式の場合 */
                    for (int i = 0; i <= 2 * m_KMdlInf.nx; i++)
                    {
                        m_KMdlInf.alpha_s[i] = Math.Sqrt(Math.Sin(m_KMdlInf.theta)) / m_KMdlInf.n / (m_KMdlInf.L * 1000.0);
                    }
                }
                else if (m_KMdlInf.type_qh_equation == 3)
                { /* 椎葉式の場合 */
                    for (int i = 0; i <= 2 * m_KMdlInf.nx; i++)
                    {
                        m_KMdlInf.alpha_a[i] = m_KMdlInf.k * Math.Sin(m_KMdlInf.theta) / m_KMdlInf.gamma / (m_KMdlInf.L * 1000.0);
                        m_KMdlInf.alpha_s[i] = Math.Sqrt(Math.Sin(m_KMdlInf.theta)) / m_KMdlInf.n / (m_KMdlInf.L * 1000.0);
                    }
                }


                /*============*/
                /*  流出計算  */
                /*============*/
                /*--------------------------------------*/
                /*  差分計算の結果を入れるメモリの取得  */
                /*--------------------------------------*/
                m_KMdlInf.h1 = new double[m_KMdlInf.nx + 1];
                m_KMdlInf.q1 = new double[m_KMdlInf.nx + 1];
                m_KMdlInf.h2 = new double[m_KMdlInf.nx + 1];
                m_KMdlInf.q2 = new double[m_KMdlInf.nx + 1];
                m_KMdlInf.depth = new double[m_KMdlInf.nx + 1];

                /*----------------------------------------*/
                /*  降雨の空間分布を格納するメモリの取得  */
                /*----------------------------------------*/
                m_KMdlInf.rr1 = new double[2 * m_KMdlInf.nx + 1];
                m_KMdlInf.rr2 = new double[2 * m_KMdlInf.nx + 1];



                /*------------------*/
                /*  境界条件の設定  */
                /*------------------*/
                m_KMdlInf.Qb = 0.0; /* すべての時刻において斜面上端での単位幅流量を 0 とする */
            }
            return true;
        }

        /// <summary>
        /// 流量流積関係式を表す関数．水深を入力としてそれに対する流量を返す．
        /// </summary>
        /// <param name="type"></param>
        /// <param name="z"></param>
        /// <param name="alpha_a"></param>
        /// <param name="alpha_s"></param>
        /// <param name="m"></param>
        /// <param name="d"></param>
        /// <returns></returns>
        static double qh_equation
                    (
                     long type,       /* 流量流積関係式のタイプ */
                                /* 1 のとき，ダルシー式(中間流) */
                                /* 2 のとき，マニング式(地表面流) */
                                /* 3 のとき，高棹・椎葉式(中間流＋地表面流) */
                     double z,       /* みかけの水深 */
                     double alpha_a, /* 中間流に対する流量流積関係式の係数 */
                     double alpha_s, /* 地表面流に対する流量流積関係式の係数 */
                     double m,       /* 流量流積関係式のべき数 */
                     double d        /* 実質の A 層厚 (mm) */
                    )
        {
            if (type == 1)
            { /* ダルシー式の場合 */
                if (z <= 0.0)
                {
                    return 0.0;
                }
                else
                {
                    return alpha_a * z;
                }
            }
            else if (type == 2)
            { /* マニング式の場合 */
                if (z <= 0.0)
                {
                    return 0.0;
                }
                else
                {
                    return alpha_s * Math.Pow(z, m);
                }
            }
            else if (type == 3)
            { /* 高棹・椎葉式の場合 */
                if (z <= 0.0)
                {
                    return 0.0;
                }
                else if (0.0 < z && z <= d)
                { /* 中間流のみの場合 */
                    return alpha_a * z;
                }
                else
                {        /* 中間流+地表面流の場合 */
                    return alpha_a * z + alpha_s * Math.Pow(z - d, m);
                }
            }
            return 0.0;
        }

        /// <summary>
        /// 流量流積関係式の逆関数．流量を入力としてそれに対する水深を返す．
        /// </summary>
        /// <param name="type"></param>
        /// <param name="q"></param>
        /// <param name="alpha_a"></param>
        /// <param name="alpha_s"></param>
        /// <param name="m"></param>
        /// <param name="d"></param>
        /// <returns></returns>
        static double hq_equation
                    (
                     long type,       /* 流量流積関係式のタイプ */
                                /* 1 のとき，ダルシー式(中間流) */
                                /* 2 のとき，マニング式(地表面流) */
                                /* 3 のとき，高棹・椎葉式(中間流＋地表面流) */
                     double q,       /* 流量 */
                     double alpha_a, /* 中間流に対する流量流積関係式の係数 */
                     double alpha_s, /* 地表面流に対する流量流積関係式の係数 */
                     double m,       /* 流量流積関係式のべき数 */
                     double d        /* 実質の A 層厚 (mm) */
                    )
        {
            double h0, q0;  /* 作業用 */

            if (type == 1)
            { /* ダルシー式の場合 */
                if (q <= 0.0)
                {
                    return 0.0;
                }
                else
                {
                    return q / alpha_a;
                }
            }
            else if (type == 2)
            { /* マニング式の場合 */
                if (q <= 0.0)
                {
                    return 0.0;
                }
                else
                {
                    return Math.Pow(q / alpha_s, 1.0 / m);
                }
            }
            else if (type == 3)
            { /* 高棹・椎葉式の場合 */
                if (q <= 0.0)
                {
                    return 0.0;
                }
                else if (0 < q && q <= alpha_a * d)
                {     /* 中間流のみの場合 */
                    return q / alpha_a;
                }
                else
                {        /* 中間流+地表面流の場合 */
                    h0 = Math.Pow((q - alpha_a * d) / alpha_s, 1.0 / m) + d;
                    q0 = alpha_s * Math.Pow(h0 - d, m) + alpha_a * h0;
                    for(;;)
                    {
                        h0 = h0 - (q0 - q) / (m * alpha_s * Math.Pow(h0 - d, m - 1.0) + alpha_a);
                        q0 = alpha_s * Math.Pow(h0 - d, m) + alpha_a * h0;
                        if (Math.Abs(q - q0) < 0.001) break;
                    }
                    return h0;
                }
            }
            return 0.0;
        }

        /// <summary>
        /// 流量流積関係式を水深で微分した値(伝播速度) dq/dh を返す．
        /// </summary>
        /// <param name="type"></param>
        /// <param name="z"></param>
        /// <param name="alpha_a"></param>
        /// <param name="alpha_s"></param>
        /// <param name="m"></param>
        /// <param name="d"></param>
        /// <returns></returns>
        static double d_qh_equation
                    (
                     long type,       /* 流量流積関係式のタイプ */
                                /* 1 のとき，ダルシー式(中間流) */
                                /* 2 のとき，マニング式(地表面流) */
                                /* 3 のとき，高棹・椎葉式(中間流＋地表面流) */
                     double z,       /* 実質の水深 */
                     double alpha_a, /* 中間流に対する流量流積関係式の係数 */
                     double alpha_s, /* 地表面流に対する流量流積関係式の係数 */
                     double m,       /* 流量流積関係式のべき数 */
                     double d        /* 実質の A 層厚 (mm) */
                    )
        {
            if (type == 1)
            { /* ダルシー式の場合 */
                return alpha_a;
            }
            else if (type == 2)
            { /* マニング式の場合 */
                if (z <= 0.0)
                {
                    return 0.0;
                }
                else
                {
                    return m * alpha_s * Math.Pow(z, m - 1.0);
                }
            }
            else if (type == 3)
            { /* 高棹・椎葉式の場合 */
                if (z <= d)
                {
                    return alpha_a;
                }
                else
                {
                    return alpha_a + m * alpha_s * Math.Pow(z - d, m - 1.0);
                }
            }
            return 0.0;
        }

        /// <summary>
        /// 4 点陰形式スキームを用いて 1 差分時間間隔先の時刻の水位・流量を求める。
        /// </summary>
        /// <param name="nx"></param>
        /// <param name="dt"></param>
        /// <param name="h1"></param>
        /// <param name="q1"></param>
        /// <param name="h2"></param>
        /// <param name="q2"></param>
        /// <param name="r1"></param>
        /// <param name="r2"></param>
        /// <param name="type"></param>
        /// <param name="g"></param>
        /// <param name="alpha_a"></param>
        /// <param name="alpha_s"></param>
        /// <param name="m"></param>
        /// <param name="d"></param>
        static void Nfdm
                    (
                     long nx,            /* 空間差分間隔の個数 (差分の節点数は nx+1 である) */
                     double dt,         /* 差分時間間隔 (hr) */
                     double[] h1,        /* 既知(現在時刻)の水深 */
                     double[] q1,        /* 既知(現在時刻)の流量 */
                     double[] h2,        /*  1 差分時間間隔先の水深 */
                     double[] q2,        /*  1 差分時間間隔先の流量 */
                     double[] r1,        /* 現在時刻の入力強度 */
                     double[] r2,        /*  1 差分時間間隔先の時刻での入力強度 */
                     long type,          /* 流量集積関係式のタイプ */
                        /*  1 のとき，ダルシー式(中間流) */
                        /*  2 のとき，マニング式(地表面流) */
                        /*  3 のとき，椎葉式(中間流＋地表面流) */
                     double[] g,         /* 地形パターン関数 */
                     double[] alpha_a,   /* 中間流に対する流量流積関係式の係数 */
                     double[] alpha_s,   /* 地表面流に対する流量流積関係式の係数 */
                     double m,          /* 流量流積関係式のべき数 */
                     double d           /* 実質の A 層厚 (mm) */
                    )
        {
            long
              i;              /* 作業用 */
            double
              dx,             /* 変数変換後の空間差分間隔 */
              c,              /* 伝搬速度 dq/dh */
              cmin = 0.0000001, /* 伝搬速度の最小値 */
              qest,           /* 計算途中の流量の値 */
              hh,             /* 水深の内挿値 */
              gg,             /* 地形パターン関数の内挿値 */
              aa,             /* 中間流に対する流量流積関係式の係数の内挿値 */
              ss,             /* 地表面流に対する流量流積関係式の係数の内挿値 */
              gr,             /* 差分対象時間・区間での平均降雨強度 */
              theta = 0.6,    /* 時間差分近似の重み係数 */
              lambda = 0.5;   /* 空間差分近似の重み係数 */

            /*============*/
            /*  初期設定  */
            /*============*/

            dx = 1.0 / (double)nx;      /* 空間差分間隔の設定 */

            /*==================*/
            /*  差分計算の開始  */
            /*==================*/

            for (i = 1; i <= nx; i++)
            {
                gg = (1.0 - lambda) * g[2 * (i - 1)] + lambda * g[2 * i];
                aa = (1.0 - lambda) * alpha_a[2 * (i - 1)] + lambda * alpha_a[2 * i];
                ss = (1.0 - lambda) * alpha_s[2 * (i - 1)] + lambda * alpha_s[2 * i];
                gr = 0.25 * (g[2 * (i - 1)] * r1[2 * (i - 1)] + g[2 * i] * r1[2 * i]
                       + g[2 * (i - 1)] * r2[2 * (i - 1)] + g[2 * i] * r2[2 * i]);

                /* 初期値の設定 */
                qest = 0.5 * (q1[i] + q2[i - 1]);

                hh = g[2 * i] * hq_equation(type, qest / g[2 * i], alpha_a[2 * i], alpha_s[2 * i], m, d);
                hh = (1.0 - theta) * (lambda * h1[i] + (1.0 - lambda) * h1[i - 1])
                       + theta * (lambda * hh + (1.0 - lambda) * h2[i - 1]);

                /* 伝播速度の計算 */
                c = d_qh_equation(type, hh / gg, aa, ss, m, d);
                if (c < cmin) c = cmin;

                for (; ; )
                {
                    q2[i] = (lambda / dt * q1[i] - (1.0 - lambda) / dt * (q2[i - 1] - q1[i - 1])
                         + c * (theta / dx * q2[i - 1] - (1.0 - theta) / dx * (q1[i] - q1[i - 1]) + gr))
                         / (lambda / dt + c * theta / dx);
                    if (q2[i] < 0.0) q2[i] = 0.0;

                    /* 収束判定 */
                    if (Math.Abs(q2[i] - qest) <= 0.0001 * qest)
                    {
                        h2[i] = g[2 * i] * hq_equation(type, q2[i] / g[2 * i], alpha_a[2 * i], alpha_s[2 * i], m, d);
                        break;
                    }

                    qest = q2[i];
                    hh = (1.0 - theta) * (lambda * h1[i] + (1.0 - lambda) * h1[i - 1])
                         + theta * (lambda * h2[i] + (1.0 - lambda) * h2[i - 1]);

                    /* 伝播速度の計算 */
                    c = d_qh_equation(type, hh / gg, aa, ss, m, d);
                    if (c < cmin) c = cmin;
                }
            }

            for (i = 0; i <= nx; i++)
            { /* 次の時刻の初期条件を設定する */
                q1[i] = q2[i];
                h1[i] = h2[i];
            }
        }

        /// <summary>
        ///  Lax Wendroff スキームを用いて 1 差分時間間隔先の時刻の水位・流量を求める．
        /// </summary>
        /// <param name="nx"></param>
        /// <param name="delta_t"></param>
        /// <param name="h1"></param>
        /// <param name="q1"></param>
        /// <param name="h2"></param>
        /// <param name="q2"></param>
        /// <param name="r1"></param>
        /// <param name="r2"></param>
        /// <param name="type"></param>
        /// <param name="g"></param>
        /// <param name="alpha_a"></param>
        /// <param name="alpha_s"></param>
        /// <param name="m"></param>
        /// <param name="d"></param>
        static void LaxWendroff
                    (
                     long nx,          /* 空間差分間隔の個数 (差分の節点数は nx+1 である) */
                     double delta_t,  /* 差分時間間隔 */
                     double[] h1,      /* 既知(現在時刻)の水深 */
                     double[] q1,      /* 既知(現在時刻)の流量 */
                     double[] h2,      /*  1 時間差分間隔先の水深 (計算結果) */
                     double[] q2,      /*  1 差分時間間隔先の流量 (計算結果) */
                     double[] r1,      /* 現在時刻の入力強度 */
                     double[] r2,      /* 1 差分時間間隔先の時刻での入力強度 */
                     long type,        /* 流量集積関係式のタイプ */
                        /*  1 のとき，ダルシー式(中間流) */
                        /*  2 のとき，マニング式(地表面流) */
                        /*  3 のとき，椎葉式(中間流＋地表面流) */
                     double[] g,       /* 地形パターン関数 */
                     double[] alpha_a, /* 中間流に対する流量流積関係式の係数 */
                     double[] alpha_s, /* 地表面流に対する流量流積関係式の係数 */
                     double m,        /* 流量流積関係式のべき数 */
                     double d         /* 実質の A 層厚 (mm) */
                    )
        {
            long
              j=0;            /* 作業用 */
            double
              dx,           /* 変数変換後の空間差分間隔 */
              dt,           /* 差分時間間隔 */
              ds,           /* 差分計算終了までの時間間隔 */
              rt,           /* 差分開始時点での降水強度 */
              qb1,          /* 与えられた q1[0] の値 */
              qb2,          /* 与えられた q2[0] の値 */
              drdt,         /* dr/dt を差分で近似した値 */
              dqdx,         /* dq/dx を差分で近似した値 */
              dqdh,         /* dq/dh を差分で近似した値 */
              H1,           /* dr/dt を差分で近似した値 */
              H2,           /* d2h/dt2 を差分で近似した値 */
              G1, G2, z;    /* 作業用 */

            /*============*/
            /*  初期設定  */
            /*============*/

            dx = 1.0 / (double)nx;      /* 空間差分間隔の設定 */
            ds = delta_t;            /* 差分計算終了までの時間間隔 */
            qb1 = q1[0];             /* 境界条件設定のために値を保存しておく */
            qb2 = q2[0];             /* 境界条件設定のために値を保存しておく */

            /*==================*/
            /*  差分計算の開始  */
            /*==================*/

            do
            {
                dt = ds;

                if (nx == 1){ goto LASTNODE;}

            START:

                z = 0.5 * (h1[0] + h1[1]) / g[1];
                dqdh = d_qh_equation(type, z, alpha_a[1], alpha_s[1], m, d);

                if (dx < dt * dqdh)
                {   /* Courant 条件を満たしているかチェック */
                    dt = dx / dqdh;       /* 満たしていなければ，Courant 条件を満たす */
                    /* ように差分時間間隔を設定しなおす．       */
                }

                dqdx = (q1[1] - q1[0]) / dx;
                rt = (r1[1] * ds + r2[1] * (delta_t - ds)) / delta_t;
                G2 = (g[1] * rt - dqdx) * dqdh;

                for (j = 1; j <= nx - 1; j++)
                {
                    G1 = G2;

                    z = 0.5 * (h1[j] + h1[j + 1]) / g[2 * j + 1];
                    dqdh = d_qh_equation(type, z, alpha_a[2 * j + 1], alpha_s[2 * j + 1], m, d);

                    if (dx < dt * dqdh)
                    {  /* Courant 条件を満たしているかチェック */
                        dt = dx / dqdh;      /* 満たしていなければ，Courant 条件を満たすように */
                        goto START;        /* 差分時間間隔を設定しなおして最初から計算しなおす．*/
                    }

                    rt = (r1[2 * j] * ds + r2[2 * j] * (delta_t - ds)) / delta_t;
                    H1 = g[2 * j] * rt - 0.5 * (q1[j + 1] - q1[j - 1]) / dx;

                    dqdx = (q1[j + 1] - q1[j]) / dx;
                    rt = (r1[2 * j + 1] * ds + r2[2 * j + 1] * (delta_t - ds)) / delta_t;
                    G2 = (g[2 * j + 1] * rt - dqdx) * dqdh;
                    drdt = (r2[2 * j] - r1[2 * j]) / delta_t;
                    H2 = g[2 * j] * drdt - (G2 - G1) / dx;

                    h2[j] = h1[j] + dt * H1 + 0.5 * dt * dt * H2;

                    if (h2[j] < 0.0) h2[j] = 0.0;
                    z = h2[j] / g[2 * j];
                    q2[j] = g[2 * j] * qh_equation(type, z, alpha_a[2 * j], alpha_s[2 * j], m, d);
                }

            LASTNODE:                          /* 斜面最下端の水深の後退差分近似 */
                j = nx - 1;
                dqdx = (q1[nx] - q1[nx - 1]) / dx;
                rt = (r1[2 * nx] * ds + r2[2 * nx] * (delta_t - ds)) / delta_t;
                H1 = g[2 * nx] * rt - dqdx;
                h2[nx] = h1[nx] + dt * H1;

                if (h2[j] < 0.0) { h2[j] = 0.0; }
                z = h2[nx] / g[2 * nx];
                q2[nx] = g[2 * nx] * qh_equation(type, z, alpha_a[2 * nx], alpha_s[2 * nx], m, d);

                ds = ds - dt;                  /* 1 差分時間間隔終了までの残り時間の計算 */
                if (ds < 0.001 * delta_t)
                {
                    ds = 0.0;
                }

                /* 次の時刻の初期条件を設定する */
                q1[0] = qb2 + (qb1 - qb2) * ds / delta_t;
                if (g[0] == 0.0)
                {
                    h1[0] = 0.0;
                }
                else
                {
                    h1[0] = g[0] * hq_equation(type, q1[0] / g[0], alpha_a[0], alpha_s[0], m, d);
                }
                for (j = 1; j <= nx; j++)
                {
                    h1[j] = h2[j];
                    q1[j] = q2[j];
                }
            } while (ds > 0.0);
        }


    }
}
