using System;
using System.Collections;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Text;
using System.Windows.Forms;
using System.IO;
using ZedGraph;


namespace Grad_Project
{
    public partial class Form1 : Form
    {
        private static int N;        // Number of grid
        private static double[] dX;  // Leght of each grid
        private static double[] phi; // Porosity of each grid
        private static double[] k;   // Permeability of each grid
        private static double[] q;   // Production rate
        private static double[] W;   // West Transmisibility
        private static double[] E;   // East Transmisibility
        private static double TMAX;  //end of simulation
        private static string[] outpuText = new string[1];// will hold all dt and Pressure to save to file
        ArrayList Output = new ArrayList();
        private static double tol=0.00000001;
        double qq, NTs;
        private static string[] textt;
        private static double[] P;
        ColorSymbolRotator rotator = new ColorSymbolRotator();
        private static int enoORfinal=0;
        private static double L = 0;
        private static double t, dt, NODES, TAU1, tau, TauNODES, TimeNew;
        private static double alphac = 5.614583;
        private static int Pgraph = 0;//Pressure graph ciziledi
  

        public Form1()
        {
            InitializeComponent();

        }     

        private void Form1_Load(object sender, EventArgs e)
        {
            buttonCalculate.Visible = false;
            textBoxTdi.Visible = false;
            buttonAdd.Visible = false;
            buttonResetGraph.Visible = false;

            GraphPane myPaneP = zg1.GraphPane;//birinci graph whole reservoir
            myPaneP.Title.Text = "Pressure Distrubition Over Reservoir";
            myPaneP.XAxis.Title.Text = "Lenght ft";
            myPaneP.YAxis.Title.Text = "Pressure Psi";


            GraphPane myPaneCS = zg2.GraphPane;//ikinci graph at specific cell
            myPaneCS.Title.Text = "Pressure At Cell";
            myPaneCS.XAxis.Title.Text = "Time day";
            myPaneCS.YAxis.Title.Text = "Pressure Psi";
           
        }

        private void genGraphCell(int eORf)
        {
            int DesiredN;

            GraphPane myPaneCS = zg2.GraphPane;//ikinci graph at specific cell
            myPaneCS.Title.Text = "Pressure At Cell";
            myPaneCS.XAxis.Title.Text = "Time day";
            myPaneCS.YAxis.Title.Text = "Pressure Psi";

            // Make up some data points from the Sine function
            PointPairList list = new PointPairList();
            //int LenghtOfMatrix =x
            list.Clear();
            LineItem Pressure;
            
            if (eORf == 0)
            {
                list.Clear();
                list.Add(0, Convert.ToDouble(textBoxPressureInitial.Text));
                list.Add(Convert.ToDouble(textBoxTMAX.Text), Convert.ToDouble(textBoxPressureInitial.Text));
                Pressure = myPaneCS.AddCurve("Initial Pressure", list, rotator.NextColor, rotator.NextSymbol);
                myPaneCS.CurveList.Clear();
            }
            else
            {

                DesiredN = Convert.ToInt32(textBoxTdi.Text);
                for (int ij = 1; ij < textt.Length; ij++)
                {
                    list.Add(Convert.ToDouble(textt[ij].Split(',')[0]), Convert.ToDouble(textt[ij].Split(',')[DesiredN]));
                }

                Pressure = myPaneCS.AddCurve("Pressure at" + DesiredN.ToString() + "Cell", list, rotator.NextColor, rotator.NextSymbol);
            }
            //myPaneCS.CurveList.Remove(Pressure);

            Pressure.Symbol.Fill = new Fill(Color.White);
            // Calculate the Axis Scale Ranges
            zg2.AxisChange();
            // OPTIONAL: Show tooltips when the mouse hovers over a point
            zg2.IsShowPointValues = true;
            // zg2.PointValueEvent += new ZedGraphControl.PointValueHandler(MyPointValueHandler);
            zg2.Invalidate();
        }
        private void genGraphPressure(int eORf)
        {
            // Get a reference to the GraphPane instance in the ZedGraphControl
            GraphPane myPaneP = zg1.GraphPane;//birinci graph whole reservoir
            
            // Set the titles and axis labels
            myPaneP.Title.Text = "Pressure Distrubition Over Reservoir ";
            myPaneP.XAxis.Title.Text = "Lenght ft";
            myPaneP.YAxis.Title.Text = "Pressure Psi";
            
            double LG = dX[0]/2;//Lenght of grid
            PointPairList list = new PointPairList();
           // PointPairList listInitial = new PointPairList();

            LineItem Pressure;

            if (eORf == 0)
            {
                list.Clear();
                list.Add(0, Convert.ToDouble(textBoxPressureInitial.Text));
                list.Add(Convert.ToDouble(textBoxLength.Text), Convert.ToDouble(textBoxPressureInitial.Text));
                Pressure = myPaneP.AddCurve("Initial Pressure", list, rotator.NextColor, rotator.NextSymbol);
                
            }
            //myPaneP.CurveList.Remove(Pressure);
            else
            {
                list.Clear();
                list.Add(LG, P[0]);
                for (int ij = 1; ij < N-1; ij++)
                {

                    LG = LG + (dX[ij - 1] + dX[ij]) / 2;
                    list.Add(LG, P[ij]);
                }
                list.Add(LG+(dX[N-2]+dX[N-1])/2, P[N-1]);
                 Pressure = myPaneP.AddCurve("Pressure after " + t.ToString() + " Days", list, rotator.NextColor, rotator.NextSymbol);

            }


            Pressure.Symbol.Fill = new Fill(Color.White);
            // Calculate the Axis Scale Ranges
            zg1.AxisChange();
            // OPTIONAL: Show tooltips when the mouse hovers over a point
            zg1.IsShowPointValues = true;
            //zg1.PointValueEvent += new ZedGraphControl.PointValueHandler(MyPointValueHandler);
            zg1.Invalidate();
        }
        private void buttonReadData_Click(object sender, EventArgs e)//read data from input.dat file
        {
            //read data from input.dat file
            
            string[] lines= File.ReadAllLines("input.dat");

            N = lines.Length; //Numbers of Grids
            
            labelNumb.Text = N.ToString();
            
            dX = new double[N];
            phi = new double[N];
            k = new double[N];
            q = new double[N];
            dt = Convert.ToDouble(textBoxDt.Text);
            
            //Getting Reservuar Properties
            for (int i = 0; i < N; i++)
            {
                dX[i] = Convert.ToDouble(lines[i].Split(',')[0]);
                phi[i] = Convert.ToDouble(lines[i].Split(',')[1]);
                k[i] = Convert.ToDouble(lines[i].Split(',')[2])/1000;
                q[i] = 0;// Convert.ToDouble(lines[i].Split(',')[3]);
                L = L + dX[i];
            }
            textBoxLength.Text = L.ToString();
            buttonCalculate.Visible = true;
            genGraphPressure(0);
            if(Pgraph==0)
                zg1.GraphPane.CurveList.Clear();
            zg1.GraphPane.IsAlignGrids = true;
            zg1.AxisChange();
            // OPTIONAL: Show tooltips when the mouse hovers over a point
            zg1.IsShowPointValues = true;
            // zg2.PointValueEvent += new ZedGraphControl.PointValueHandler(MyPointValueHandler);
            zg1.Invalidate();
            genGraphCell(0);

        }
        private void buttonCalculate_Click(object sender, EventArgs e)
        {
            W = new double[N]; // West side
            double[] C = new double[N]; // Centre
            E = new double[N]; // East side
            double[] Q = new double[N]; // Right hand side
            P = new double[N]; // Pressure
            double[] NTS = new double[N];//next time step/stop
            

            double dy = Convert.ToDouble(textBoxWidthW.Text);
            double dz = Convert.ToDouble(textBoxHeightH.Text);
            double A = dy * dz;
            double B = Convert.ToDouble(textBoxFVF.Text);
            double mu = Convert.ToDouble(textBoxViscosity.Text);
            double ct = Convert.ToDouble(textBoxTotalCompressibility.Text);
            double Pl = Convert.ToDouble(textBoxPL.Text); //Left side Pressure
            double Pr = Convert.ToDouble(textBoxPR.Text); //Right Side Pressure
            double ls, rs;//left and right side controller
            bool davam=true;
     
            TMAX = Convert.ToDouble(textBoxTMAX.Text);//tmax
            //NODES = Convert.ToDouble(textBoxNODES.Text);
            //dt = Convert.ToDouble(textBoxDTi.Text);
            
            // t = 0;
            // TauNODES = Math.Log(TMAX);
            // TAU1 = Math.Log(dt);
            // t = Math.Exp(TAU1);
            // int step = 2;
            
            

            // Ilk addim ucun P matrisini initial P ile doldurmaq lazimdi.

            for (int i = 0; i <= N - 1; i = i + 1)
            {
                P[i] = Convert.ToDouble(textBoxPressureInitial.Text);
            }

            #region//Time Grid No sozunu ve grid nolari arrayliste ata
            outpuText[0] = "Time/Grid No";
            for (int ii = 0; ii < N; ii++)
                outpuText[0] = outpuText[0] + "," + (ii+1).ToString();
            Output.Add(outpuText[0]);
            //initial basinci dosyaya yazmaq ucun Output array liste t=0 ve Pi degerini ata
            outpuText[0] = t.ToString();
            for (int ii = 0; ii < N; ii++)
                outpuText[0] = outpuText[0] + "," + P[ii].ToString();
            Output.Add(outpuText[0]);
            #endregion

            W[0] = 0.0;   E[N - 1] = 0.0;

            for (int i = 1; i <= N - 1; i = i + 1)
            {
                W[i] = 1.127 * 2 * A / (dX[i - 1] / k[i - 1] + dX[i] / k[i]) / (mu * B);
            }
            //dt isin icine girdi timestepi goze almaq lazimdi



            do
            {
                //ilk td burda atanmalidi

                #region //time ile elaqeli
             
                for (int count = 0; count < N; count++)// butun quyulari yoxla en yaxin Stop belirle NTs
                {
                    RateRead(count + 1, t + tol, out qq, out NTs);
                    q[count] = qq;
                    NTS[count] = NTs;
                }

                //for (int count = 0; count < N - 1; count++)
                //{
                //    if (NTS[count] <= NTS[count + 1])
                //        NTs = NTS[count];//Enkucuk Nt al
                //    else
                //        NTs = NTS[count + 1];//En kucuk Nt al
                    
                //}
                
                 //t NTs den boyuk olmadigini yoxla
                //if (t > NTs)
                //    {
                //        t = t - dt;
                //        dt = NTs - t;
                //    }

                //    if (t == NTs)
                //        t = t + tol;
                

                #endregion
                    for (int i = 0; i <= N - 2; i = i + 1)
                {
                    E[i] = 1.127 * 2 * A / (dX[i] / k[i] + dX[i + 1] / k[i + 1]) /( mu * B);
                }
                  // sag sol uclrada P olub olmadigini yoxla
            if (Pl > 0) { ls = 1.0; } else { ls = 0.0; }
            if (Pr > 0) { rs = 1.0; } else { rs = 0.0; }

            C[0] = -( E[0] 
                + (dX[0] * A * phi[0] * ct) / (alphac * B * dt) 
                + (ls * 1.127 * 2 * A * k[0] / ((dX[0]) * mu * B)));

            C[N - 1] = 
                -( 
                   W[N - 1] 
                + (dX[N - 1] * A * phi[N - 1] * ct) / (alphac * B * dt)
                + (rs * 1.127 * 2 * A * k[N - 1] / ((dX[N - 1]) * mu * B))
                  );

            for (int i = 1; i <= N - 2; i = i + 1)
                { C[i] = -( W[i] + E[i] + (dX[i] * A * phi[i] * ct) / (alphac * B * dt));  }
            
            for (int i = 0; i < N ; i = i + 1)
                {Q[i] = -((dX[i] * A * phi[i] * ct) / (alphac * B * dt)) * P[i] - q[i]; }
            
       // Left and Rightd side pressures Pl Pr
            
            Q[0] = Q[0] - ((ls * 1.127 * 2 * A * k[0] / ((dX[0]) * mu * B)) * Pl);
            Q[N - 1] = Q[N - 1] - (rs * 1.127 * 2 * A * k[N - 1] / ((dX[N - 1]) * mu * B)) * Pr;
            
            P = TDMASolve(W, C, E, Q);
            t = t + dt;
        
            #region outputTexte at dt ve P lerei
            outpuText[0] = t.ToString();
                for (int ii = 0; ii < N; ii++)
                    outpuText[0] = outpuText[0] + "," + P[ii].ToString();
            Output.Add(outpuText[0]);

            #endregion
             
                //--------------------
            //tau = TAU1 + (TauNODES - TAU1) * ((step - 1) / (NODES - 1));
            //TimeNew = Math.Exp(tau);
            //dt = TimeNew - t;

            //t = TimeNew;
            //step++;
            if (t >= TMAX)
                davam = false;
            } //do while sonu
             while (davam);

            //
            textt = (string[])Output.ToArray(typeof(string)); 
            File.WriteAllLines("output.csv", textt);
            Output.Clear();
            enoORfinal = 1;
            genGraphPressure(enoORfinal);
            
            textBoxTdi.Visible = true;
            int centre = (int)(N / 2+1);
            textBoxTdi.Text = centre.ToString();
            buttonAdd.Visible = true;
        }
        private void RateRead(int wellNo, double time, out double Rate, out double Nt)
        {
            Rate = 0;
            Nt = 10;
            string[] RateList = File.ReadAllLines("rate.dat");
            int Row = RateList.Length; //Numbers of lines in file
            double ts, tf, cr, tfOi, tsOi;//ts start time, tf final time, cr current rate
            int ok = 0;
            for (int j = 0; j < Row; j++)
            {
                if (ok == 1) break;
                int xx = Convert.ToInt32(RateList[j].Split(',')[0]);
                if (wellNo == xx)//kuyunun oldugu satiri bul;
                {
                    int nos = Convert.ToInt32(RateList[j].Split(',')[1]);//number of rate steps
                    int eos;//number of digits sepwrated by ',' up to end of steps
                    eos = nos * 3 + 1;


                    for (int i = 2; i < eos; i = i + 3)//satir icinde irlei get
                    {

                        ts = Convert.ToDouble(RateList[j].Split(',')[i]);
                        tf = Convert.ToDouble(RateList[j].Split(',')[i + 1]);
                        cr = Convert.ToDouble(RateList[j].Split(',')[i + 2]);

                        if (ts <= time && tf >= time)
                        {
                            Rate = cr;
                            Nt = tf;
                            ok = 1;
                            //return Rate;
                        }
                        else //if there is no rate at the time 
                        {
                            Rate = 0;
                            for (int ti = 3; ti < eos - 1; ti = ti + 3)
                            {
                                tfOi = Convert.ToDouble(RateList[j].Split(',')[ti]);//tfinal of i
                                tsOi = Convert.ToDouble(RateList[j].Split(',')[ti + 2]);//t start of i+1
                                if (tfOi <= time && tsOi >= time)
                                {
                                    Nt = tsOi;
                                    ok = 1;
                                }
                            }//for
                        }//else
                        if (ok == 1) break;
                    }//setr icinde gez

                }//kolonlarda axdar


            }//for

            if (ok == 0) // if there is no info about well 
            {
                Rate = 0;
                Nt = TMAX;//TMAX burda geri donecek
                for (int j = 0; j < Row; j++)
                {
                    int kuyuNo = Convert.ToInt32(RateList[j].Split(',')[0]);
                    double tstar = Convert.ToInt32(RateList[j].Split(',')[2]);
                    if (kuyuNo == wellNo & time < tstar)
                        Nt = Convert.ToInt32(RateList[j].Split(',')[2]);

                }
            }


        }
        double[] TDMASolve(double[] a, double[] b, double[] c, double[] f)
        {
            /* Warning: will modify c and d! */
            int n = f.Length; /* All arrays should be the same length */
            double[] x = new double[n];
            double[] d = new double[n];
            d = f;
            double id;

            /* Modify the coefficients. */
            c[0] /= b[0]; /* Division by zero risk. */
            d[0] /= b[0]; /* Division by zero would imply a singular matrix. */
            for (int i = 1; i < n; i++)
            {
                id = b[i] - c[i - 1] * a[i];
                c[i] /= id; /* This calculation during the last iteration is redundant. */
                d[i] = (d[i] - d[i - 1] * a[i]) / id;
            }

            /* Now back substitute. */
            x[n - 1] = d[n - 1];
            for (int i = n - 2; i >= 0; i--)
            {
                x[i] = d[i] - c[i] * x[i + 1];
            }

            return x;
        }
        private void button1_Click(object sender, EventArgs e)
        {
            
            buttonResetGraph.Visible = true;
            genGraphCell(1);
            zg2.GraphPane.IsAlignGrids = true;
            zg2.AxisChange();
            // OPTIONAL: Show tooltips when the mouse hovers over a point
            zg2.IsShowPointValues = true;
            // zg2.PointValueEvent += new ZedGraphControl.PointValueHandler(MyPointValueHandler);
            zg2.Invalidate();
        }
        private void button1_Click_1(object sender, EventArgs e)
        {
            buttonResetGraph.Visible = false;

        zg2.GraphPane.CurveList.Clear();
            zg2.GraphPane.IsAlignGrids = true;
            zg2.AxisChange();
            // OPTIONAL: Show tooltips when the mouse hovers over a point
            zg2.IsShowPointValues = true;
            // zg2.PointValueEvent += new ZedGraphControl.PointValueHandler(MyPointValueHandler);
            zg2.Invalidate();
        }
    }
}