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

namespace vcsFAPP
{
    public partial class Form1 : Form
    {

        struct DRAWPARA
        {          //作図用データ構造体
            public int kpt;       //拡大倍率（画面表示：1，ファイル保存：2）
            public string sxjiku; //ｘ軸名
            public string syjiku; //ｙ軸名
            public double xmin;   //ｘ軸最小値
            public double xmax;   //ｘ軸最大値
            public double dx;     //ｘ軸増分
            public double ymin;   //ｙ軸最小値
            public double ymax;   //ｙ軸最大値
            public double dy;     //ｙ軸増分
            public int kxi0;      //ｘ軸最小値座標（px）
            public int kxf0;      //ｘ軸最大値座標（px）
            public int kyi0;      //ｙ軸最小値座標（px）
            public int kyf0;      //ｙ軸最大値座標（px）
            public int HSIZE;     //画像横寸法
            public int DHL;       //画像左余白（縦軸描画スペース：px）
            public int DHR;       //画像右余白（px）
            public int VSIZE;     //画像縦サイズ（px）
            public int DVL;       //画像下余白（横軸描画スペース：px）
            public int DVU;       //画像上余白（px）
        }
        //--------------------------------------------------------------------------------
        
        public Form1()
        {
            InitializeComponent();
        }
        //--------------------------------------------------------------------------------

        private void toolStripButton1_Click(object sender, EventArgs e)
        {
            main_part();
        }
        //--------------------------------------------------------------------------------

        private void main_part()
        {
            int i, k;
            string fname = "";
            string fnameR = ""; //入力ファイル名
            string fnameW = ""; //出力ファイル名
            string fnameP = ""; //png画像出力ファイル名
            System.IO.StreamReader sr;
            System.IO.StreamWriter sw;
            string dat;
            string[] sbuf;
            char delim = ',';
            int ndata, mm, nn;
            string scom;
            DRAWPARA pltdata;
            double[] datax;
            double[] datay;
            double[] func;
            double[] xr;
            double[] xi;
            double[] ak;
            double[] bk;

            double dt;

            INI_PLTDATA(out pltdata);

            openFileDialog1.InitialDirectory = System.IO.Directory.GetCurrentDirectory();
            if (openFileDialog1.ShowDialog() == DialogResult.OK) { fnameR = openFileDialog1.FileName; }
            if (saveFileDialog1.ShowDialog() == DialogResult.OK) { fnameW = saveFileDialog1.FileName; }
            fname = System.IO.Path.GetFileNameWithoutExtension(fnameW);
            fnameP = System.IO.Path.GetDirectoryName(fnameW) + "\\fig_" + fname + ".png";

            //データ入力
            sr = new System.IO.StreamReader(fnameR, System.Text.Encoding.Default);
            dat = sr.ReadLine(); sbuf = dat.Split(delim); scom = sbuf[0];
            dat = sr.ReadLine(); sbuf = dat.Split(delim); pltdata.sxjiku = sbuf[0]; pltdata.syjiku = sbuf[1];
            dat = sr.ReadLine(); sbuf = dat.Split(delim); pltdata.HSIZE = int.Parse(sbuf[0]); pltdata.DHL = int.Parse(sbuf[1]); pltdata.DHR = int.Parse(sbuf[2]);
            dat = sr.ReadLine(); sbuf = dat.Split(delim); pltdata.VSIZE = int.Parse(sbuf[0]); pltdata.DVL = int.Parse(sbuf[1]); pltdata.DVU = int.Parse(sbuf[2]);
            dat = sr.ReadLine(); sbuf = dat.Split(delim); pltdata.xmin = double.Parse(sbuf[0]); pltdata.xmax = double.Parse(sbuf[1]); pltdata.dx = double.Parse(sbuf[2]);
            dat = sr.ReadLine(); sbuf = dat.Split(delim); pltdata.ymin = double.Parse(sbuf[0]); pltdata.ymax = double.Parse(sbuf[1]); pltdata.dy = double.Parse(sbuf[2]);
            dat = sr.ReadLine(); sbuf = dat.Split(delim); ndata = int.Parse(sbuf[0]); mm = int.Parse(sbuf[1]);
            datax = new double[ndata];
            datay = new double[ndata];
            func = new double[ndata];
            for (i = 0; i <= ndata - 1; i++)
            {
                dat = sr.ReadLine(); sbuf = dat.Split(delim);
                datax[i] = double.Parse(sbuf[0]);
                datay[i] = double.Parse(sbuf[1]);
            }
            sr.Close();

            //フーリエ変換の準備と実行
            nn = 2;
            do
            {
                nn = nn * 2;
            } while (nn < ndata * 1); //必要に応じて後続の0を追加
            xr = new double[nn];
            xi = new double[nn];
            for (i = 1; i <= nn; i++)
            {
                xr[i - 1] = 0.0; xi[i - 1] = 0.0;
            }
            for (i = 1; i <= ndata; i++)
            {
                xr[i - 1] = datay[i - 1] - datay[ndata - 1];
            }
            FFT(nn, xr, xi);//フーリエ変換
            for (i = 1; i <= nn; i++)
            { //戻り値をデータ数で除す
                xr[i - 1] = xr[i - 1] / (double)nn;
                xi[i - 1] = xi[i - 1] / (double)nn;
            }

            dt = datax[1] - datax[0];
            ak = new double[(int)(nn / 2) + 1];
            bk = new double[(int)(nn / 2) + 1];
            for (i = 0; i <= (int)(nn / 2); i++)
            {
                ak[i] = 2.0 * xr[i];
                bk[i] = -2.0 * xi[i];
            }
            for (i = 0; i <= ndata - 1; i++)
            {
                func[i] = 0.5 * ak[0];
                for (k = 1; k <= mm; k++)
                {
                    func[i] = func[i] + ak[k] * Math.Cos(2.0 * Math.PI * k * i / nn / dt);
                    func[i] = func[i] + bk[k] * Math.Sin(2.0 * Math.PI * k * i / nn / dt);
                }
                func[i] = func[i] + datay[ndata - 1];
            }

            //結果出力
            sw = new System.IO.StreamWriter(fnameW, false, System.Text.Encoding.Default);
            sw.WriteLine(scom);
            sw.WriteLine("ndata," + ndata.ToString("0"));
            sw.WriteLine("x,y,func");
            for (i = 0; i <= ndata - 1; i++)
            {
                dat = datax[i].ToString("E");
                dat = dat + "," + datay[i].ToString("E");
                dat = dat + "," + func[i].ToString("E");
                sw.WriteLine(dat);
            }
            sw.Close();

            //画像出力
            PLTCTL(pltdata, ndata - 1, datax, datay, func, scom);
            pictureBox2.Image.Save(fnameP, System.Drawing.Imaging.ImageFormat.Png);

            MessageBox.Show("計算完了");

        }
        //--------------------------------------------------------------------------------
        private void FFT(int nn, double[] xr, double[] xi)
        {
            //**************************************************
            //高速フーリエ変換・逆変換
            //**************************************************
            //nn  ：データ数（2の累乗）
            //xr[]：入出力データ実数部
            //xi[]：入出力データ虚数部

            int g, h, i, j, k, l, m, n, p, q;
            double a, b, xd;

            n = nn;
            double[] s = new double[(int)(n / 2) + 1];
            double[] c = new double[(int)(n / 2) + 1];

            i = 0; j = 0; k = 0; l = 0; p = 0; h = 0; g = 0; q = 0;
            m = (int)(Math.Log((double)n) / Math.Log(2.0) + 1.0);
            a = 0.0;
            b = Math.PI * 2.0 / (double)n;
            for (i = 0; i <= (int)(n / 2); i++)
            {
                s[i] = Math.Sin(a);
                c[i] = Math.Cos(a);
                a = a + b;
            }
            l = n; h = 1;
            for (g = 1; g <= m; g++)
            {
                l = (int)(l / 2); k = 0;
                for (q = 1; q <= h; q++)
                {
                    p = 0;
                    for (i = k; i <= l + k - 1; i++)
                    {
                        j = i + l;
                        a = xr[i] - xr[j];
                        b = xi[i] - xi[j];
                        xr[i] = xr[i] + xr[j];
                        xi[i] = xi[i] + xi[j];
                        if (p == 0)
                        {
                            xr[j] = a;
                            xi[j] = b;
                        }
                        else
                        {
                            xr[j] = a * c[p] + b * s[p];
                            xi[j] = b * c[p] - a * s[p];
                        }
                        p = p + h;
                    }
                    k = k + l + l;
                }
                h = h + h;
            }
            j = (int)(n / 2);
            for (i = 1; i <= n - 1; i++)
            {
                k = n;
                if (j < i)
                {
                    xd = xr[i]; xr[i] = xr[j]; xr[j] = xd;
                    xd = xi[i]; xi[i] = xi[j]; xi[j] = xd;
                }
                k = (int)(k / 2);
                while (j >= k)
                {
                    j = j - k; k = (int)(k / 2);
                    if (k == 0) break;
                }
                j = j + k;
            }
        }
        //---------------------------------------------------------------------------

        private void PLTCTL(DRAWPARA pltdata, int nd, double[] datax, double[] datay, double[] func,string scom)
        {
            //nd：データ数（入力はデータ数，プログラム上数値はデータ数−１）
            //datax[]：プロットデータｘ座標
            //datay[]：プロットデータｙ座標
            Bitmap bmp;
            Graphics g;
            PictureBox thePicBox = pictureBox1;
            int i;

            //プロット領域定義
            pltdata.kxi0 = pltdata.DHL;
            pltdata.kxf0 = pltdata.HSIZE - pltdata.DHR;
            pltdata.kyi0 = pltdata.DVU;
            pltdata.kyf0 = pltdata.VSIZE - pltdata.DVL;

            for (i = 0; i <= 1; i++)
            {
                switch (i)
                {
                    case 0: //画面描画
                        pltdata.kpt = 1;
                        thePicBox = pictureBox1;
                        thePicBox.Visible = true;
                        break;
                    case 1: //画像書き出し
                        pltdata.kpt = 2;
                        thePicBox = pictureBox2;
                        thePicBox.Visible = false;
                        break;
                }
                thePicBox.Size = new Size(pltdata.kpt * pltdata.HSIZE, pltdata.kpt * pltdata.VSIZE);
                bmp = new Bitmap(thePicBox.Width, thePicBox.Height);
                thePicBox.Image = bmp;
                g = Graphics.FromImage(thePicBox.Image);
                g.FillRectangle(Brushes.White, 0, 0, thePicBox.Width, thePicBox.Height);
                PLOT(g, pltdata, nd, datax, datay, func,scom);
                g.Dispose();
            }
            pictureBox1.Left = 0;
            pictureBox1.Top = toolStrip1.Height;
            ClientSize = new Size(pictureBox1.Width, pictureBox1.Height + toolStrip1.Height);
        }

        //---------------------------------------------------------------------------------
        private void PLOT(Graphics g, DRAWPARA pltdata, int nd, double[] datax, double[] datay, double[] func,string scom)
        {
            int i;
            int kxi, kxf, kyi, kyf;
            double xx, yy;
            int kxx, kyy;
            int kxx1, kyy1, kxx2, kyy2;
            string str;
            System.Drawing.Font f = new Font("ＭＳゴシック", pltdata.kpt * 10);
            System.Drawing.SizeF TextSize1;
            System.Drawing.SizeF TextSize2;

            kxi = pltdata.kpt * pltdata.kxi0;
            kxf = pltdata.kpt * pltdata.kxf0;
            kyi = pltdata.kpt * pltdata.kyi0;
            kyf = pltdata.kpt * pltdata.kyf0;
            //座標軸
            System.Drawing.Pen LPen = new Pen(System.Drawing.Color.Black);
            LPen.DashStyle = System.Drawing.Drawing2D.DashStyle.Dot;
            DRN_XJIKU(g, f, LPen, pltdata, kxi, kxf, kyi, kyf);
            DRN_YJIKU(g, f, LPen, pltdata, kxi, kxf, kyi, kyf);
            //枠線描画
            g.DrawRectangle(new Pen(Color.Black, 2), kxi, kyi, kxf - kxi, kyf - kyi);
            //ｙ軸名描画
            str = pltdata.syjiku;
            TextSize2 = g.MeasureString(str, f);
            TextSize1 = g.MeasureString("-1.0", f);
            kxx = (int)(kxi - pltdata.kpt * 10 - TextSize1.Width - TextSize2.Height);
            kyy = (int)((kyi + kyf) / 2 + TextSize2.Width / 2);
            INC_STR(g, f, str, kxx, kyy, -90);
            //ｘ軸名描画
            str = pltdata.sxjiku;
            TextSize2 = g.MeasureString(str, f);
            kxx = (int)((kxi + kxf) / 2 - TextSize2.Width / 2);
            kyy = (int)(kyf + pltdata.kpt * 5 + TextSize1.Height);
            INC_STR(g, f, str, kxx, kyy, 0);
            //データプロット
            int dc;
            dc = pltdata.kpt * 5;
            for (i = 0; i <= nd; i++)
            {
                xx = datax[i];
                yy = datay[i];
                kxx = kxi + (int)((xx - pltdata.xmin) * (kxf - kxi) / (pltdata.xmax - pltdata.xmin)) - (int)(dc / 2);
                kyy = kyf - (int)((yy - pltdata.ymin) * (kyf - kyi) / (pltdata.ymax - pltdata.ymin)) - (int)(dc / 2);
                g.DrawEllipse(new Pen(Color.Red, 1), kxx, kyy, dc, dc);
            }
            //回帰推定線
            xx = datax[0];
            yy = func[0];
            kxx1 = kxi + (int)((xx - pltdata.xmin) * (kxf - kxi) / (pltdata.xmax - pltdata.xmin));
            kyy1 = kyf - (int)((yy - pltdata.ymin) * (kyf - kyi) / (pltdata.ymax - pltdata.ymin));
            for (i = 1; i <= nd; i++)
            {
                xx = datax[i];
                yy = func[i];
                kxx2 = kxi + (int)((xx - pltdata.xmin) * (kxf - kxi) / (pltdata.xmax - pltdata.xmin));
                kyy2 = kyf - (int)((yy - pltdata.ymin) * (kyf - kyi) / (pltdata.ymax - pltdata.ymin));
                g.DrawLine(new Pen(Color.Blue, 2), kxx1, kyy1, kxx2, kyy2);
                kxx1 = kxx2;
                kyy1 = kyy2;
            }
            //残差
            dc = pltdata.kpt * 5;
            for (i = 0; i <= nd; i++)
            {
                xx = datax[i];
                yy = datay[i]-func[i];
                kxx = kxi + (int)((xx - pltdata.xmin) * (kxf - kxi) / (pltdata.xmax - pltdata.xmin)) - (int)(dc / 2);
                kyy = kyf - (int)((yy - pltdata.ymin) * (kyf - kyi) / (pltdata.ymax - pltdata.ymin)) - (int)(dc / 2);
                g.DrawEllipse(new Pen(Color.Green, 1), kxx, kyy, dc, dc);
            }
            //右上コメント
            str = scom;
            TextSize1 = g.MeasureString(str, f);
            kxx = kxf - 2 * pltdata.kpt - (int)TextSize1.Width;
            kyy = kyi - 2 * pltdata.kpt - (int)TextSize1.Height;
            g.DrawString(str, f, Brushes.Black, kxx, kyy);
            f.Dispose();
        }
        //---------------------------------------------------------------------------------
        private void INC_STR(Graphics g, System.Drawing.Font f, string str, int kxx, int kyy, float ang)
        {
            //軸名描画
            g.ScaleTransform(1, 1);//横・縦の表示比率を設定
            g.TranslateTransform(kxx, kyy);//表示位置の設定（表示位置を原点とする座標移動）
            g.RotateTransform(ang);//表示角度を指定
            g.DrawString(str, f, Brushes.Black, 0, 0);//描画実行
            g.ResetTransform();//単位行列にリセット
        }
        //---------------------------------------------------------------------------------
        private void DRN_XJIKU(Graphics g, System.Drawing.Font f, System.Drawing.Pen LPen,
                               DRAWPARA pltdata, int kxi, int kxf, int kyi, int kyf)
        {
            int i, ix, kxx;
            double xx, wv;
            string str;
            System.Drawing.SizeF TextSize1;
            //普通ｘ軸描画
            ix = (int)((pltdata.xmax - pltdata.xmin) / pltdata.dx);
            for (i = 0; i <= ix; i++)
            {
                xx = pltdata.xmin + (double)(i) * pltdata.dx;
                kxx = kxi + (int)((xx - pltdata.xmin) * (kxf - kxi) / (pltdata.xmax - pltdata.xmin));
                wv = pltdata.xmin + (double)(i) * pltdata.dx;
                str = wv.ToString("0");
                if ((int)(pltdata.dx * 1000.0) % 10 != 0) { str = wv.ToString("0.000"); }
                if ((int)(pltdata.dx * 1000.0) % 10 == 0 && (int)(pltdata.dx * 100.0) % 10 != 0) { str = wv.ToString("0.00"); }
                if ((int)(pltdata.dx * 1000.0) % 10 == 0 && (int)(pltdata.dx * 100.0) % 10 == 0 && (int)(pltdata.dx * 10.0) % 10 != 0) { str = wv.ToString("0.0"); }
                TextSize1 = g.MeasureString(str, f);
                g.DrawLine(LPen, kxx, kyi, kxx, kyf);
                g.DrawString(str, f, Brushes.Black, kxx - TextSize1.Width / 2, kyf + pltdata.kpt * 3);
            }
        }
        //---------------------------------------------------------------------------------
        private void DRN_YJIKU(Graphics g, System.Drawing.Font f, System.Drawing.Pen LPen,
                               DRAWPARA pltdata, int kxi, int kxf, int kyi, int kyf)
        {
            int i, iy, kyy;
            double yy, wv;
            string str;
            System.Drawing.SizeF TextSize1;
            //普通ｙ軸描画
            iy = (int)((pltdata.ymax - pltdata.ymin) / pltdata.dy);
            for (i = 0; i <= iy; i++)
            {
                yy = pltdata.ymin + (double)(i) * pltdata.dy;
                kyy = kyf - (int)((yy - pltdata.ymin) * (kyf - kyi) / (pltdata.ymax - pltdata.ymin));
                wv = pltdata.ymin + (int)(i) * pltdata.dy;
                str = wv.ToString("0");
                if ((int)(pltdata.dy * 1000.0) % 10 != 0) { str = wv.ToString("0.000"); }
                if ((int)(pltdata.dy * 1000.0) % 10 == 0 && (int)(pltdata.dy * 100.0) % 10 != 0) { str = wv.ToString("0.00"); }
                if ((int)(pltdata.dy * 1000.0) % 10 == 0 && (int)(pltdata.dy * 100.0) % 10 == 0 && (int)(pltdata.dy * 10.0) % 10 != 0) { str = wv.ToString("0.0"); }
                TextSize1 = g.MeasureString(str, f);
                g.DrawLine(LPen, kxi, kyy, kxf, kyy);
                g.DrawString(str, f, Brushes.Black, kxi - TextSize1.Width - pltdata.kpt * 2, kyy - TextSize1.Height / 2);
            }
        }
        //---------------------------------------------------------------------------------

        private void INI_PLTDATA(out DRAWPARA pltdata)
        {
            //pltdata初期化
            pltdata.kpt = 0;
            pltdata.HSIZE = 0;
            pltdata.DHL = 0;
            pltdata.DHR = 0;
            pltdata.VSIZE = 0;
            pltdata.DVL = 0;
            pltdata.DVU = 0;
            pltdata.sxjiku = "";
            pltdata.xmin = 0.0;
            pltdata.xmax = 0.0;
            pltdata.dx = 0.0;
            pltdata.syjiku = "";
            pltdata.ymin = 0.0;
            pltdata.ymax = 0.0;
            pltdata.dy = 0.0;
            pltdata.kxi0 = 0;
            pltdata.kxf0 = 0;
            pltdata.kyi0 = 0;
            pltdata.kyf0 = 0;
        }
        //---------------------------------------------------------------------------------
    
    }
}